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

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

0 引 言

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

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

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

1 算法原理

本文算法在傳統(tǒng)電離層層析成像的基礎(chǔ)上，從IRI模型中提取垂直電子密度剖面獲取EOF作為垂直方向的基函數(shù)。使用幾個EOF的線性組合來表示每一個垂直電子密度剖面，從而將反演的未知數(shù)由每個格網(wǎng)的電子密度轉(zhuǎn)變?yōu)榇怪逼拭娴腅OF系數(shù)，顯著減少了需要反演的未知數(shù)。之后利用MART將重新構(gòu)造的反演矩陣進(jìn)行迭代計算得到EOF系數(shù)值，再結(jié)合EOF獲得所有格網(wǎng)的電子密度。

1.1 電離層層析成像原理

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

(1)

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

(2)

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

(3)

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

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

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

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

(4)

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

1.2 EOF基函數(shù)

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

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

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

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

(5)

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

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

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

1.3 EOF垂直基函數(shù)與MART結(jié)合

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

(6)

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

(7)

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

(8)

式中：n為格網(wǎng)總數(shù)；c為每個剖面對應(yīng)的EOF系數(shù)個數(shù)。式(8)可以進(jìn)一步寫為：

Hm×c·Yc×1=ym×1

(9)

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

(10)

2 實驗與結(jié)果分析

2.1 數(shù)值模擬實驗

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

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

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

(11)

(12)

(13)

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

(a) UT07:00

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

(a) UT05:00

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

在實驗中，使用MART算法，需要反演的未知數(shù)個數(shù)為6 750個，而通過使用EOF作為垂直基函數(shù)，其需要反演的未知數(shù)個數(shù)只有450個，顯著減少了反演未知數(shù)個數(shù)，大大加快了反演速度。圖5是UT02:00時刻兩種算法的迭代次數(shù)與精度對比。可以看出，采用EOF作為垂直基函數(shù)，迭代7次左右就已經(jīng)達(dá)到收斂，而傳統(tǒng)MART算法則需要25次左右才達(dá)到收斂，且使用EOF作為垂直基函數(shù)的算法精度更高。這是因為采用EOF作為垂直基函數(shù)，當(dāng)垂直方向的一系列格網(wǎng)中格網(wǎng)電子密度有修正時，可以帶動整條垂直方向的格網(wǎng)都進(jìn)行修正，從而收斂速度更快。這也證實了使用EOF作為垂直基函數(shù)的算法在計算效率上的優(yōu)越性。

2.2 實測數(shù)據(jù)實驗

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

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

圖7是一天內(nèi)UT 01:00至UT23:00各個時刻兩種算法反演得到的電子密度的F2層峰值電子密度(NmF2)與測高儀觀測值的對比?？梢钥闯觯瑑煞N算法反演獲得的NmF2與測高儀觀測數(shù)據(jù)總體上符合都較好，體現(xiàn)了電離層一天內(nèi)的峰值電子密度隨著時間推移先增大后逐漸減少的特性。大部分情況下，使用EOF垂直基函數(shù)的算法更接近測高儀觀測結(jié)果。

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

3 結(jié) 語

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