使用EOF垂直基函數的像素基電離層層析重構

范玉榮 符杰林* 安 濤 王俊義 林基明

1(桂林電子科技大學認知無線電與信息處理教育部重點實驗室 廣西 桂林 541004)2(中國科學院上海天文臺 上海 200030)3(廣西高校衛星導航與位置感知重點實驗室 廣西 桂林 541004)

0 引 言

電離層是距離地面高度60～2 000 km的大氣電離區域，其電子密度分布會影響無線電信號在電離層中的傳播條件，因而電離層電子密度分布信息不僅對電離層物理具有重要意義，而且對星地通信和導航也具有重要意義[1]。使用GNSS觀測數據和電離層層析技術對電離層進行三維建模，不僅克服了二維電離層模型的局限性，而且具有全天候觀測、全球覆蓋、低成本和高效率等優勢[2]。

電離層層析成像是根據反演區域內的大量信號傳播路徑上的總電子含量來反演特定區域的電子密度分布[3]。其通常分為兩類：函數基電離層層析和像素基電離層層析。函數基電離層模型采用一組函數來表示電離層電子密度分布，其優點在于可以用很少的參數即可表述區域內的電離層分布情況，但由于地面GNSS站的分布和衛星星座的幾何特性，直接對觀測資料進行求解往往十分困難[4]。對于像素基電離層層析，通常將反演區域劃分為一系列格網并假定每個格網的電子密度為常數且均勻分布。然而，由于GNSS觀測視角和地面測站數量及分布的限制，基于GNSS的像素基電離層層析通常存在數據不足而引起的不適定問題。為了解決不適定問題帶來的影響，國內外眾多學者先后提出了各種解決方法，文獻[5]使用同時迭代算法(Simultaneous Iteration Reconstruction Technique，SIRT)進行電離層電子密度的重構，減小了噪聲對解算結果的影響；Wen等[6]對代數重構(Algebraic reconstruction Technique，ART)算法松弛因子進行了改進，提出了IART算法；姚宜斌等[7]針對聯合迭代重構算法迭代收斂慢且易受噪聲影響的問題，利用上一輪迭代的電子密度反演結果，自適應地調整松弛因子和加權參數，提出了ASIRT算法；文獻[8]提出的TV-MART算法，使用總變差(Total Variation，TV)最小化結合MART算法對電子密度進行重構，減小了由噪聲引起的不穩定性；文獻[9]對IART算法使用自搜索松弛因子改進，減小了由于松弛因子不合適而引起的誤差；文獻[10]將局部加權線性回歸算法應用在電離層層析中，在一定程度上減弱了沒有射線通過的格網對初始值的依賴。以上方法在一定程度上有效克服了重構過程中的不適定問題，提高了反演精度。但是，以上算法大多計算比較復雜，尤其在格網數目較大時，以上迭代算法會由于反演未知數過多和計算復雜的原因導致反演效率較低[11]，不便于實現電離層的實時預報預警研究。

本文通過從IRI2016中提取EOF垂直基函數，顯著減少了未知數的個數，然后結合MART算法重構電離層電子密度，提高了反演效率，同時在一定程度上減小了不適定問題的影響。通過數值模擬仿真和實測數據實驗結果表明該算法相對于傳統層析算法在計算效率上有顯著提升，同時在反演精度上也有一定的改進。

1 算法原理

本文算法在傳統電離層層析成像的基礎上，從IRI模型中提取垂直電子密度剖面獲取EOF作為垂直方向的基函數。使用幾個EOF的線性組合來表示每一個垂直電子密度剖面，從而將反演的未知數由每個格網的電子密度轉變為垂直剖面的EOF系數，顯著減少了需要反演的未知數。之后利用MART將重新構造的反演矩陣進行迭代計算得到EOF系數值，再結合EOF獲得所有格網的電子密度。

1.1 電離層層析成像原理

電離層層析成像使用GPS射線的電離層斜向總電子含量(Slant Total Electron Content，STEC)重構電離層電子密度分布，其中STEC可以從雙頻GPS(Global Positioning System)接收機的偽距和載波相位觀測值中提取出來，如式(1)所示，具體方法的詳細描述見文獻[12]。

(1)

式中：P4,sm是載波相位平滑偽距之后的觀測值；c表示光速；DCBi和DCBj分別表示衛星和接收機的差分碼偏差；f1和f2分別對應GNSS兩個頻點的頻率。獲得的STEC可以表示成電子密度沿信號路徑上的積分，即：

(2)

式中：Ne為電子密度；l為信號路徑；r為t時刻經度、緯度和高度所組成的位置向量；s為信號傳播路徑。將反演區域按照經度、緯度、高度方向上劃分為三維格網，那么每條路徑上的STEC測量值可以表示為：

(3)

式中：m為穿過電離層的射線總數；n為總格網數目；aij為第i條射線在第j個格網內的截距；xj為第j個格網的電子密度；ei為觀測噪聲和誤差。如圖1所示為電離層層析模型的簡化示意圖。

圖1 電離層層析模型示意圖

將式(3)寫成矩陣形式可以表示為：

ym×1=Am×n·xn×1+em×1

(4)

式中：y表示GNSS信號射線傳播路徑上電離層STEC構成的m維列向量；A為GNSS信號射線穿過格網時的截距構成的m×n維投影矩陣，由于觀測衛星數有限，其通常是一個稀疏矩陣；x為所有格網像素中心電子密度構成的n維列向量；e為噪聲和觀測誤差構成的列向量。

1.2 EOF基函數

經驗正交函數可以提取一組簡化的變量來減少數據的維數，這些變量可以解釋數據中的大部分特性[13]。因而原始數據可以表示為提取出的少量經驗正交函數的線性組合，被廣泛應用在氣象和電離層數據分析和建模中。如施闖等[14]使用EOF和球諧函數分別作為垂直方向和水平方向基函數建立了3D電離層模型；Ansari等[15]使用EOF建立了韓國區域電離層總電子含量(Total Electron Content，TEC)模型。

為了獲取電離層垂直剖面的EOF，通過IRI2016模型提取電子密度剖面矩陣，參數如表1所示。IRI2016模型是由空間研究協會(CORPAR)和國際無線電科學聯盟(URSI)發起的經驗電離層模型，其綜合了全球范圍內的電離層測高儀、散射雷達等電離層探測儀器的觀測數據，因而在電離層建模中多使用其作為電離層背景模式。

表1 從IRI2016獲取電子密度剖面的參數范圍

產生的電子密度剖面矩陣用Np×q表示，其中:p表示垂直像素個數；q表示從IRI2016中獲取的電子密度垂直剖面數。通過對垂直剖面數據矩陣Np×q進行奇異值分解(SVD)，構造一組EOF基函數，具體如下：

(5)

式中：奇異值包含在對角矩陣S中；EOF基函數在矩陣U中。我們可以通過從奇異值中計算每個EOF的貢獻百分比來選擇占主導地位的EOF[15]。由于大部分情況下前三個EOF即可以代表96%以上的信息量，因而本文選取前三個占比較大的EOF。圖2為UT02:00時刻的EOF示意圖。

圖2 UT02:00時刻的前三個EOF

這三個EOF可以構成矩陣Ep×3，那么對于反演區域內任一垂直剖面的電子密度都可以用EOF來線性表示。

1.3 EOF垂直基函數與MART結合

MART由于其迭代速度較快且可以克服反演值為負的問題，因而在電離層層析中被廣泛使用。MART反演迭代公式如下：

(6)

式(3)可以用EOF基函數表示為：

(7)

式中：Fn×c是由Ep×3按照對角組合得到的，其擴展方法如下：

(8)

式中：n為格網總數；c為每個剖面對應的EOF系數個數。式(8)可以進一步寫為：

Hm×c·Yc×1=ym×1

(9)

值得注意的是，由于EOF有負數值，為防止迭代過程中解發散的問題，對Hij做了絕對值處理，因而最終的反演迭代公式為：

(10)

2 實驗與結果分析

2.1 數值模擬實驗

由于對全部反演區域真實的電離層狀態未知，直接使用實測數據進行電離層電子密度反演無法全面地評判算法的有效性和穩定性，因而本文設計了數值模擬實驗以驗證算法的有效性和穩定性。選取的反演區域為35°N至55°N和5°W至25°E，高度范圍為100至1 000 km。在緯度、經度和高度上的空間分辨率分別取2°、2°、20 km。選取2018年3月20日歐洲地區45個IGS觀測站的GPS雙頻觀測數據進行建模。具體觀測站的選取和分布如圖3所示，其中三角形標志為PQ052電離層測高儀觀測站的位置。在實驗過程中，使用10 min GPS觀測數據進行實驗，測站和衛星位置均取真實值。

圖3 IGS站點及測高儀觀測站位置分布圖

根據衛星位置和觀測站位置計算射線在每個格網中的截距并構成式(4)中的投影矩陣A。考慮到EOF是由IRI模型中獲取得到的，因而本文設置的真實電子密度為IRI2016的基礎上添加一定的偏移量，如式(11)所示；同時在模擬STEC數據上增加了最大為1.5TECu(total electron content unit)的噪聲；此外，初始電子密度設為xIRI。

(11)

(12)

(13)

圖4中(a)、(b)和(c)分別是UT07:00、UT12:00和UT19:00三個時刻兩種算法解算結果在(51°N,11°E)處的電子密度剖面對比圖?？梢钥闯?，在高度小于200 km的區域，MART反演結果比使用EOF基函數的反演結果要更加接近真實值，而在高度高于200 km的區域則反之。此外，表2給出了以上三個時刻兩種算法的均方根誤差和平均電子密度誤差，可以看出，兩種算法反演得到的電離層電子密度都比較接近真實值，在電子密度較大的時候反演誤差也隨之增大，但都遠小于此時的電子密度峰值?？傮w來看，使用EOF基函數的算法相較于MART在精度上有一定提升。

(a) UT07:00

圖5 兩種算法的迭代次數與精度對比

(a) UT05:00

表2 兩種算法反演精度對比

在實驗中，使用MART算法，需要反演的未知數個數為6 750個，而通過使用EOF作為垂直基函數，其需要反演的未知數個數只有450個，顯著減少了反演未知數個數，大大加快了反演速度。圖5是UT02:00時刻兩種算法的迭代次數與精度對比?？梢钥闯觯捎肊OF作為垂直基函數，迭代7次左右就已經達到收斂，而傳統MART算法則需要25次左右才達到收斂，且使用EOF作為垂直基函數的算法精度更高。這是因為采用EOF作為垂直基函數，當垂直方向的一系列格網中格網電子密度有修正時，可以帶動整條垂直方向的格網都進行修正，從而收斂速度更快。這也證實了使用EOF作為垂直基函數的算法在計算效率上的優越性。

2.2 實測數據實驗

為進一步評估使用EOF作為基函數反演電離層電子密度的精度，使用2018年3月20日歐洲地區45個IGS觀測站的實際GPS雙頻觀測數據獲取STEC觀測值進行計算。反演區域和格網劃分與上述模擬實驗一致。利用反演區域內電離層測高儀PQ052(50°N，14.6°E)的觀測數據對反演結果進行核驗。

圖6給出了UT05:00、UT11:00和UT19:00三個時刻兩種算法反演得到的電離層電子密度剖面與測高儀觀測數據的對比。從比較結果來看，使用EOF作為垂直基函數的算法的反演結果與測高儀觀測數據更為接近。說明使用EOF作為垂直基函數的算法的精度相較于傳統MART反演算法確實有一定提升。但是，兩種算法在電子密度峰值高度上均與測高儀觀測值存在一定差異，一方面是由于采用的GPS射線仰角較高而導致的垂直方向分辨率不高，另一方面是由于采用的EOF來源于IRI2016，因而反演得到的電子密度峰值高度與IRI2016接近。

圖7是一天內UT 01:00至UT23:00各個時刻兩種算法反演得到的電子密度的F2層峰值電子密度(NmF2)與測高儀觀測值的對比。可以看出，兩種算法反演獲得的NmF2與測高儀觀測數據總體上符合都較好，體現了電離層一天內的峰值電子密度隨著時間推移先增大后逐漸減少的特性。大部分情況下，使用EOF垂直基函數的算法更接近測高儀觀測結果。

圖7 兩種算法反演得到的NmF2與測高儀結果的對比

3 結 語

本文采用從IRI中提取出的EOF作為垂直基函數，結合MART算法實現對電離層電子密度進行快速精確層析反演。不同于傳統電離層層析方法中將格網電子密度作為未知數反演，本文算法采用EOF作為垂直方向基函數，將EOF系數作為反演未知數，顯著提高了反演效率，解決了傳統方法中由于未知數數目較多而導致的反演效率低的問題。使用歐洲地區45個觀測站10 min內GPS數據對電離層進行了建模實驗。通過模擬數據實驗和實測數據實驗表明，相對于傳統MART反演算法，本文算法在反演效率和精度上都有所提升。