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

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

0 引 言

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

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

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

1 算法原理

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

1.1 電離層層析成像原理

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

(1)

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

(2)

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

(3)

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

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

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

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

(4)

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

1.2 EOF基函數(shù)

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

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

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

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

(5)

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

圖2 UT02:00時(shí)刻的前三個(gè)EOF

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

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

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

(6)

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

(7)

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

(8)

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

Hm×c·Yc×1=ym×1

(9)

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

(10)

2 實(shí)驗(yàn)與結(jié)果分析

2.1 數(shù)值模擬實(shí)驗(yàn)

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

圖3 IGS站點(diǎn)及測(cè)高儀觀測(cè)站位置分布圖

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

(11)

(12)

(13)

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

(a) UT07:00

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

(a) UT05:00

表2 兩種算法反演精度對(duì)比

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

2.2 實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

3 結(jié) 語(yǔ)

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