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

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

電離層內(nèi)大量存在的自由電子會(huì)影響無(wú)線電信號(hào)在電離層中的傳播，因而電離層電子密度分布信息不僅對(duì)電離層物理具有重要意義，而且對(duì)星地通信和衛(wèi)星導(dǎo)航也具有重要意義[1]。傳統(tǒng)的二維電離層建模假定電離層區(qū)域所有的自由電子都集中在一定高度的球殼薄層上，雖然也能達(dá)到較好的精度，但無(wú)法反映電離層的空間分布情況[2]。

GNSS電離層層析成像根據(jù)反演區(qū)域內(nèi)大量信號(hào)傳播路徑上的斜向總電子含量(slant total electronic density, 簡(jiǎn)稱STEC)來(lái)反演特定區(qū)域的電子密度分布，能夠有效反演得到三維電離層信息。相比傳統(tǒng)電離層探測(cè)方法，GNSS電離層層析具有探測(cè)持續(xù)時(shí)間長(zhǎng)、探測(cè)范圍廣及成本低等優(yōu)勢(shì)，近年來(lái)在電離層探測(cè)中被廣泛應(yīng)用[2-3]。然而，由于GNSS衛(wèi)星觀測(cè)視角和地面觀測(cè)站數(shù)量及分布的限制，基于GNSS的電離層層析通常存在因數(shù)據(jù)不足而引起的不適定問題[2-10]。為解決不適定問題帶來(lái)的影響，國(guó)內(nèi)外學(xué)者提出了各種解決方法。文獻(xiàn)[4]采用同時(shí)迭代重建技術(shù)(SIRT)進(jìn)行電離層電子密度的反演，減少了噪聲對(duì)解算結(jié)果的影響；聞德保等[5]通過(guò)代數(shù)重構(gòu)算法(ART)對(duì)松弛因子進(jìn)行了改進(jìn)，提出了IART算法；姚宜斌等[6]針對(duì)聯(lián)合迭代重構(gòu)算法迭代收斂慢且易受噪聲影響的問題，利用上一輪迭代的電子密度反演結(jié)果，自適應(yīng)地調(diào)整松弛因子和加權(quán)參數(shù)，提出了ASIRT算法；文獻(xiàn)[7]提出了一種TV-MART算法，用總變差(TV)最小化結(jié)合MART算法對(duì)電子密度進(jìn)行反演，減少了由噪聲引起的不穩(wěn)定性；文獻(xiàn)[8]用自搜索松弛因子對(duì)IART算法進(jìn)行改進(jìn)，減少了由于松弛因子不合適而引起的誤差。以上方法在一定程度上解決了反演過(guò)程中的不適定問題，提高了反演精度，但仍存在由于觀測(cè)站分布不均和觀測(cè)角度有限帶來(lái)的很多格網(wǎng)依舊缺乏射線穿過(guò)而得不到有效修正的問題，尤其是靠近邊緣部分的格網(wǎng)，電離層層析成像的病態(tài)性依然突出[9]。

文獻(xiàn)[10]提出一種約束自適應(yīng)聯(lián)合迭代反演方法，通過(guò)對(duì)觀測(cè)矩陣附加平滑約束，對(duì)無(wú)射線穿過(guò)的格網(wǎng)電子密度進(jìn)行了一定的修正，但當(dāng)無(wú)射線穿過(guò)的格網(wǎng)較多時(shí)，其精度提升并不明顯。文獻(xiàn)[11]將邊緣射線用于電離層層析中，提高了射線格網(wǎng)覆蓋率，顯著提高了反演精度，但該反演精度取決于處于層析范圍內(nèi)的邊緣射線的總電子含量觀測(cè)值的精度，而這個(gè)精度較難控制。

鑒于此，通過(guò)在地面缺乏觀測(cè)站的區(qū)域設(shè)立虛擬觀測(cè)站，建立一系列虛擬觀測(cè)射線，再建立二維薄殼電離層模型，并根據(jù)投影函數(shù)計(jì)算出這些虛擬觀測(cè)射線的STEC值。這些虛擬觀測(cè)射線能夠有效擴(kuò)大射線的覆蓋范圍，使得原本無(wú)射線穿過(guò)的格網(wǎng)的電子密度也能得到有效修正，從而得到更加準(zhǔn)確的電離層電子密度分布。

1 系統(tǒng)模型及問題分析

電離層層析成像使用GPS射線的電離層STEC反演電離層電子密度分布。STEC可從雙頻GPS接收機(jī)的偽距和載波相位觀測(cè)值中提取[12]：

(1)

其中：P4,sm為載波相位平滑偽距之后的觀測(cè)值；c為光速；DCBi、DCBj分別為衛(wèi)星和接收機(jī)的差分碼偏差(differential code biases, 簡(jiǎn)稱DCB)。

獲得的STEC觀測(cè)值可表示為電子密度沿信號(hào)路徑上的積分，

(2)

其中：Ne為電子密度；l為信號(hào)路徑；r為t時(shí)刻經(jīng)度、緯度和高度所組成的位置向量。將反演區(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è)噪聲與其他觀測(cè)誤差的和。電離層層析模型的簡(jiǎn)化示意圖如圖1所示。

圖1 電離層層析模型簡(jiǎn)化示意圖

式(3)的矩陣形式可表示為

y=Ax+e，

(4)

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

由于觀測(cè)衛(wèi)星數(shù)有限，式(4)中的投影矩陣A通常是一個(gè)稀疏矩陣，運(yùn)用迭代重構(gòu)算法對(duì)式(4)進(jìn)行解算。由于MART迭代速度較快，且能解決反演值為負(fù)的問題，在電離層層析中被廣泛使用，其反演迭代公式為

(5)

在實(shí)際觀測(cè)中，受地面觀測(cè)站分布不均和衛(wèi)星觀測(cè)角度的限制，部分格網(wǎng)無(wú)任何射線穿過(guò)，因而這些格網(wǎng)的電子密度依賴于初始經(jīng)驗(yàn)值，得不到有效修正。通過(guò)增加虛擬射線的方式增加射線穿過(guò)格網(wǎng)的覆蓋率，通過(guò)事先建立的薄殼電離層模型結(jié)合投影函數(shù)獲得虛擬射線的STEC觀測(cè)值。

2 利用虛擬射線改進(jìn)的反演方法

薄殼電離層模型由于其結(jié)構(gòu)簡(jiǎn)單且擁有不錯(cuò)的精度而在電離層建模方法中被廣泛應(yīng)用[13]。本研究采用低階球諧函數(shù)模型進(jìn)行電離層建模。電離層球諧函數(shù)模型[14]為

(6)

根據(jù)式(6)用足量的STEC觀測(cè)值建立方程組，利用最小二乘法可解算出球諧函數(shù)未知球諧系數(shù)。通過(guò)二維薄殼模型和投影函數(shù)可得到區(qū)域內(nèi)任意位置任意仰角的STEC值，且觀測(cè)仰角越高，其精度越高[15]。

在地面觀測(cè)站分布較少的區(qū)域設(shè)置虛擬觀測(cè)站，將其周圍仰角較高的一系列射線作為虛擬射線，利用二維薄殼模型和投影函數(shù)計(jì)算虛擬射線的STEC值，與實(shí)際觀測(cè)射線一起依據(jù)式(4)、(5)進(jìn)行反演重構(gòu)。設(shè)置虛擬射線的具體操作流程如下：

1)將反演區(qū)域底部按經(jīng)緯度劃分為2°×1.5°的格網(wǎng)，將上下左右的格網(wǎng)內(nèi)均無(wú)觀測(cè)站的格網(wǎng)的中心位置設(shè)為虛擬觀測(cè)站位置。

2)將反演區(qū)域頂部按經(jīng)緯度劃分為0.5°×0.5°的格網(wǎng)，每個(gè)格網(wǎng)中心點(diǎn)與地面虛擬觀測(cè)站連線，從中選擇仰角不低于60°的射線作為虛擬射線。為使虛擬射線通過(guò)盡可能多的格網(wǎng)，將虛擬觀測(cè)站位置高度設(shè)置為-50 km。

3)計(jì)算虛擬射線在450 km高度面的穿刺點(diǎn)位置及天頂角，利用之前得到的薄殼電離層模型和投影函數(shù)計(jì)算得到虛擬射線的STEC值。

4)將虛擬射線與實(shí)測(cè)射線按式(5)構(gòu)建矩陣，用MART進(jìn)行解算。

3 數(shù)據(jù)處理及分析

采用歐洲地區(qū)2018年9月20日的IGS觀測(cè)站的GPS觀測(cè)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。反演區(qū)域設(shè)為緯度35°N～55°N，經(jīng)度5°W～25°E，高度范圍為100～1 000 km。在緯度、經(jīng)度和高度上的空間分辨率分別取1°、2°、20 km。IGS觀測(cè)站及虛擬觀測(cè)站的分布如圖2所示。其中，深色·表示實(shí)際觀測(cè)站，淺色· 表示虛擬觀測(cè)站，深色×表示測(cè)高儀觀測(cè)站。

圖2 觀測(cè)站分布圖

3.1 DCB及薄殼電離層模型精度評(píng)估

DCB是提取STEC觀測(cè)值中的主要誤差來(lái)源，精確估計(jì)DCB可提高STEC精度，同時(shí)也可得到更加精確的薄殼電離層模型。單獨(dú)用24 h的GPS觀測(cè)數(shù)據(jù)解算得到當(dāng)天的衛(wèi)星和接收機(jī)DCB信息，在之后的薄殼電離層模型解算和電離層層析中都利用這些DCB信息以得到更加精確的STEC觀測(cè)值。為驗(yàn)證解算得到的DCB的精度，將解算得到的DCB與中科院(CAS)發(fā)布的DCB數(shù)據(jù)進(jìn)行對(duì)比，得到如圖3所示的差值折線圖。

圖3 解算得到的DCB與CAS的DCB差值

由于CAS給出的接收機(jī)DCB數(shù)據(jù)有限，圖3(b)只包含CAS給出了數(shù)據(jù)的部分觀測(cè)站的DCB差值。從圖3可看出，解算得到的GPS衛(wèi)星DCB和接收機(jī)DCB與CAS發(fā)布的DCB差值均在1 ns內(nèi)，其中衛(wèi)星DCB的最大誤差為0.537 ns，觀測(cè)站接收機(jī)DCB的最大誤差為0.945 ns。這表明，解算得到的DCB具有較高的精度，從而保證了通過(guò)提取得到的STEC的可靠性。

用30 min內(nèi)GPS觀測(cè)數(shù)據(jù)進(jìn)行電離層薄殼模型解算，將解算得到的薄殼模型與CODE發(fā)布的全球電離層圖GIM進(jìn)行對(duì)比。解算得到的每個(gè)時(shí)刻的電離層薄殼模型與GIM的差值均值和均方根如表1所示。從表1可看出，解算得到的薄殼模型與GIM的誤差基本在1.3 TECU內(nèi)。這表明，得到的薄殼模型精度較高，同時(shí)也保證了添加的虛擬射線的STEC精度。

表1 薄殼模型與GIM的差值均值與均方根

3.2 模擬數(shù)據(jù)實(shí)驗(yàn)

由于對(duì)反演區(qū)域整體真實(shí)的電離層狀態(tài)未知，直接用實(shí)測(cè)數(shù)據(jù)進(jìn)行電離層電子密度反演無(wú)法全面評(píng)判算法的有效性和穩(wěn)定性，因而用數(shù)值模擬實(shí)驗(yàn)驗(yàn)證算法的有效性和穩(wěn)定性。根據(jù)觀測(cè)站位置和衛(wèi)星位置計(jì)算射線在每個(gè)格網(wǎng)中的截距并構(gòu)建投影矩陣A。其中，觀測(cè)站坐標(biāo)和衛(wèi)星坐標(biāo)均為2019年9月20日UT11:45～12:15觀測(cè)間隔內(nèi)的真實(shí)值。利用IRI2016模型獲得2019年9月20日UT12:00時(shí)刻待反演區(qū)域內(nèi)各格網(wǎng)中心點(diǎn)處的電離層電子密度xIRI。將A與xIRI相乘，得到各條射線上的模擬STEC。為了更接近真實(shí)情況，在模擬STEC數(shù)據(jù)上增加了最大為模擬STEC均值10%的隨機(jī)噪聲。同時(shí)，將初始電離層電子密度設(shè)置為0.8xIRI。用xIRI作為真實(shí)值進(jìn)行評(píng)估。

圖4為添加了虛擬射線后的射線分布圖，其中淺色線為實(shí)際觀測(cè)射線，深色線為添加的虛擬觀測(cè)射線。從圖4可看出，虛擬射線填補(bǔ)了由于實(shí)際觀測(cè)站和衛(wèi)星位置分布限制導(dǎo)致的格網(wǎng)覆蓋空白，使得射線格網(wǎng)覆蓋率大大提升。經(jīng)統(tǒng)計(jì)，添加虛擬射線后格網(wǎng)覆蓋率達(dá)到了90%以上，比添加前的格網(wǎng)覆蓋率提高了約30%。

圖4 射線分布圖

圖5為2種方法迭代30次的均方根誤差曲線。其中：3DCIT表示傳統(tǒng)電離層層析方法；3DCIT-VR表示添加了虛擬射線的電離層層析方法。從圖5可看出，添加虛擬射線后，由于格網(wǎng)覆蓋率的優(yōu)勢(shì)，其一步迭代后的誤差更小，且誤差下降速度更快。最終收斂時(shí)，3DCIT-VR、3DCIT的均方根誤差分別為0.063 9×1011、0.120 6×1011el/m3。

圖5 2種方法迭代誤差

圖6為2種方法反演得到的電離層電子密度誤差經(jīng)度切片。從圖6(a)可看出，未添加虛擬射線時(shí)，由于大量格網(wǎng)無(wú)射線穿過(guò)而未得到任何修正，因而其誤差較大，尤其在地面測(cè)站分布較少的區(qū)域呈現(xiàn)出明顯的不連續(xù)性。從圖6(b)可看出，3DCIT-VR的誤差明顯更小，且更加均勻，這表明添加虛擬射線能夠提高電離層層析精度。

圖6 2種方法反演得到的電離層電子密度誤差經(jīng)度切片

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

為進(jìn)一步驗(yàn)證本方法的優(yōu)越性，用30 min內(nèi)實(shí)際觀測(cè)數(shù)據(jù)進(jìn)行電離層層析反演。先用30 min內(nèi)觀測(cè)數(shù)據(jù)建立區(qū)域薄殼電離層模型，再算得虛擬射線的STEC值，最后將虛擬觀測(cè)數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)一同進(jìn)行迭代反演，得到電離層電子密度分布。將反演結(jié)果與未添加虛擬射線時(shí)的反演結(jié)果進(jìn)行比較，并使用PQ052和AT138兩個(gè)電離層測(cè)高儀的實(shí)測(cè)數(shù)據(jù)對(duì)層析結(jié)果進(jìn)行驗(yàn)證。

圖7為用2種方法在UT6:00、UT12:00和UT18:00三個(gè)時(shí)刻反演得到的電子密度與相應(yīng)時(shí)刻PQ052、AT138兩個(gè)電離層觀測(cè)站的觀測(cè)數(shù)據(jù)對(duì)比曲線。由圖2可知，PQ052位于觀測(cè)站較為豐富的區(qū)域，AT138周圍則缺乏實(shí)際觀測(cè)站。從圖7(a)可看出，2種方法得到的反演結(jié)果較接近。這表明添加虛擬射線并未降低實(shí)際射線豐富區(qū)域的反演精度，且有時(shí)對(duì)反演精度還有一些改善。從圖7(b)可看出，3DCIT-VR反演得到的電離層電子密度剖面更接近電離層實(shí)際觀測(cè)值。這也表明了添加虛擬射線可在一定程度上提高電離層反演精度，從而進(jìn)一步證明了本方法的優(yōu)越性。從圖7還可看出，反演得到的電離層電子密度剖面峰值高度與實(shí)測(cè)值峰值高度存在一些差異，這是因?yàn)橛^測(cè)射線仰角偏高導(dǎo)致垂直分辨率不高，未來(lái)可進(jìn)一步引入掩星觀測(cè)數(shù)據(jù)進(jìn)行反演以改善這個(gè)問題。

4 結(jié)束語(yǔ)

通過(guò)在缺乏地面觀測(cè)站的位置設(shè)置虛擬觀測(cè)站，以添加虛擬射線的方式增加可用射線數(shù)量，進(jìn)而提高格網(wǎng)覆蓋率。本方法在保證射線豐富區(qū)域反演精度的同時(shí)，使無(wú)射線穿過(guò)的格網(wǎng)電離層電子密度也能得到有效修正。通過(guò)模擬數(shù)據(jù)實(shí)驗(yàn)及實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)驗(yàn)證了本方法的有效性，并利用2個(gè)測(cè)高儀觀測(cè)站(PQ052、AT138)分別對(duì)射線豐富區(qū)域和射線欠豐富區(qū)域的電離層反演效果做了驗(yàn)證。實(shí)驗(yàn)結(jié)果表明，本方法效果顯著，解決了缺乏射線穿過(guò)的格網(wǎng)電離層電子密度得不到有效反演的問題。

圖7 2種方法反演得到的電子密度與PQ052和AT138兩個(gè)電離層觀測(cè)站觀測(cè)數(shù)據(jù)對(duì)比