電離層三維層析成像的自適應(yīng)聯(lián)合迭代重構(gòu)算法

姚宜斌，湯俊，張良，何暢勇，張順

1武漢大學(xué)測繪學(xué)院，武漢 430079

2武漢大學(xué)地球空間環(huán)境與大地測量教育部重點實驗室，武漢 430079

1 引言

電離層是日地空間環(huán)境的一個重要組成部分，與人類的高新技術(shù)和日常生活息息相關(guān).為了深入研究電離層電子密度時空分布的精細結(jié)構(gòu)，Austen等（1986，1988）首次提出了電離層層析成像（Computerized Ionospheric Tomography，CIT）技術(shù)的設(shè)想.隨后，國際上利用該技術(shù)相繼開展了許多實驗和理論研究（Pryse，2003；Raymund et al.，1990；Raymund，1994；Rius et al.，1997；Spencer et al.，1998）.顧及電離層層析成像過程中投影視角的有限性和測站布設(shè)的稀疏性，電離層層析成像中需要解決的核心問題是由于觀測數(shù)據(jù)的缺失而引起的不適定問題，近年來，許多學(xué)者通過引入電離層電子密度的先驗信息，來克服反演過程中的不適定性，由此出現(xiàn)了各種不同的反演算法（Arikan et al.，2007；Chartier et al.，2012；Hobiger et al.，2008；Lee and Kamalabadi，2009；Li et al.，2012；Liu et al.，2010；Ma et al.，2005；Nesterov and Kunitsyn，2011；Wen et al.，2007，2008）.其中，聯(lián)合迭代重構(gòu)算法（Simultaneous Iteration Reconstruction Technique，SIRT）由于計算簡便且易于實現(xiàn)，在電離層層析成像中得到廣泛應(yīng)用.Pryse和Kersley（1992）最早將聯(lián)合迭代重構(gòu)算法應(yīng)用于電離層層析成像；Hobiger等（2008）利用先驗約束的聯(lián)合迭代重構(gòu)算法反演電離層電子密度，并與常用的聯(lián)合迭代重構(gòu)算法進行對比，發(fā)現(xiàn)其迭代收斂速度和反演結(jié)果精度均有所提高；Liu等（2010）針對邊緣缺少約束的問題，提出了改進的先驗約束聯(lián)合迭代重構(gòu)算法，并對其優(yōu)越性進行了驗證；Nesterov和Kunitsyn（2011）提出了基于正則化的聯(lián)合迭代重構(gòu)算法，分別利用全球和區(qū)域的數(shù)據(jù)驗證了該方法的有效性.這些算法都有效的反演了電離層電子密度，并能夠反映電離層電子密度的時空特性變化規(guī)律，然而，由于在每輪迭代過程中，松弛因子和加權(quán)參數(shù)固定不變，從而使得電離層電子密度反演過程中迭代收斂較慢，反演結(jié)果精度不高.

針對上述問題，本文發(fā)展了一種自適應(yīng)的聯(lián)合迭代算法（Adaptive Simultaneous Iteration Reconstruction Technique，ASIRT），并通過模擬數(shù)據(jù)和實測數(shù)據(jù)實驗來驗證該方法的有效性和優(yōu)越性.

2 電離層層析成像原理

電離層層析成像就是根據(jù)反演區(qū)域內(nèi)的一系列衛(wèi)星信號傳播路徑上的TEC信息來反演所選定區(qū)域內(nèi)電離層電子密度的時空分布.如圖1所示，在利用衛(wèi)星信號觀測進行的電離層層析成像過程中，所獲得的電離層電子總含量（Total Electron Content，TEC）是電離層電子密度沿衛(wèi)星至接收機射線路徑上的積分（Liu，2004；鄒玉華，2003），可表示為

式中，Ne為電離層電子密度，l為衛(wèi)星信號傳播路徑，r為t時刻經(jīng)度、緯度和高度所組成的位置向量.電離層電子密度與電離層TEC之間是非線性的.在實際反演過程中，為了反演方便，通常采用離散的反演方法將待反演的電離層空間離散化.對于離散化的像素基層析反演模型，選取像素指標函數(shù)bj作為基函數(shù)，如果射線穿過某像素，則bj為1，否則為0；并將電離層按經(jīng)度、緯度以及高度方向上離散化為三維的格網(wǎng)，其公式表示為

式中，n為離散化的格網(wǎng)數(shù)（總的像素數(shù)），xj（j＝1，…，n）為模型參數(shù)（離散化后的電離層格網(wǎng)電子密度）.那么對每條射線路徑上的TEC測量值可以表示為

式中，m為電離層TEC觀測值總數(shù)，aij為投影矩陣元素，即第i條射線在第j個格網(wǎng)內(nèi)的截距.考慮到測量中觀測噪聲和離散誤差的影響，且假定在一定時間段格網(wǎng)內(nèi)電子密度是不變的，則每條射線傳播路徑上的電離層TEC測量數(shù)據(jù)可表示為

將（5）式用矩陣形式表示如下：

式中，y為電離層TEC觀測值組成的m維列向量，A為投影矩陣（射線在對應(yīng)像素內(nèi)的截距構(gòu)成的m個n維的行向量，由于可視衛(wèi)星數(shù)有限，其通常是一個稀疏矩陣），x為未知參數(shù)組成的n維列向量，e為觀測噪聲和離散誤差組成的m維列向量.

3 自適應(yīng)聯(lián)合迭代重構(gòu)算法

3.1 算法原理

聯(lián)合迭代重構(gòu)算法是基于平方優(yōu)化技術(shù)的一種并行迭代方法.每一次對所有觀測路徑上的電離層電子密度進行迭代，并根據(jù)每一次迭代的修正量再對電子密度分布做整體修正（Lytle，1979）.其迭代式表示如下：

從式（8）中可以看出，加權(quán)參數(shù)只與投影矩陣A相關(guān)，而投影矩陣A本身又附加了很多的前提假設(shè)，如假定一定時間段內(nèi)格網(wǎng)電子密度不變，由此，這一方法并不完全符合實際的電離層層析成像過程.因此，本文將自適應(yīng)的思想引入聯(lián)合迭代重構(gòu)算法，即考慮將上一輪的迭代電子密度值引入到權(quán)值中，以便在某種程度上改正投影矩陣的不精確性.則改正后的加權(quán)參數(shù)為

式中，M和R為對稱正定矩陣.該算法中的松弛因子λk對于迭代過程中結(jié)果的精度及收斂速度有很重要的影響.本文中松弛因子（Elfving et al.，2012）的選擇如下：

由此，本文在聯(lián)合迭代重構(gòu)算法的基礎(chǔ)上，通過利用上一輪的電離層電子密度反演結(jié)果，自適應(yīng)地調(diào)整松弛因子和加權(quán)參數(shù)，發(fā)展了一種新的自適應(yīng)的聯(lián)合迭代重構(gòu)算法.

3.2 收斂性分析

考慮到如下線性方程的形式：式中，A具有非負性.以下推導(dǎo)過程中，向量都是列向量，ai為系數(shù)矩陣A的第i行的轉(zhuǎn)置，內(nèi)積的標準形式為〈x，y〉＝xTy，ρ（·）為譜半徑（最大正特征值）.對于每一個線性方程（14），定義超平面Hi＝｛x∈Rn｜〈ai，x〉＝y(tǒng)i｝.

如果

則迭代方程（12）收斂于解x＊（minx‖Ax－y‖M）.其中，τ是一個任意小的固定常數(shù).

由收斂條件可知：y＝Ax＊；假定rk＝y(tǒng)－Axk，令dk＝Mkrk，ek＝xk－x＊，則

4 計算結(jié)果與分析

4.1 模擬數(shù)據(jù)實驗

為驗證算法ASIRT的可行性以及相對于SIRT的優(yōu)越性，本文首先利用模擬的TEC數(shù)據(jù)來進行實驗.利用模擬數(shù)據(jù)進行反演時，所選取的經(jīng)度范圍為0°—20°E，緯度范圍為40°N—60°N，高度范圍為100～1000km，格網(wǎng)間隔在經(jīng)度和緯度方向上為1°，在高度方向上為50km.為了更加接近實際情況，測站數(shù)據(jù)和衛(wèi)星數(shù)據(jù)均取真實值，選取歐洲地區(qū)的IGS跟蹤站并計算衛(wèi)星在反演時段內(nèi)的坐標，獲得相應(yīng)反演時段內(nèi)射線穿過格網(wǎng)內(nèi)的截距，并構(gòu)成觀測方程（6）中的投影矩陣A.利用IRI2007模型得到2012年4月8日UT10∶00時待反演區(qū)域內(nèi)各網(wǎng)格中心點處的電離層電子密度xIRI，并且將各條射線傳播路徑上的電離層TEC用ysimu表示.考慮到實際觀測中觀測誤差和離散誤差的存在，進行模擬時，加入一定量的隨機誤差e，即

分別采用ASIRT和SIRT算法進行反演，并將反演的結(jié)果進行對比.圖2給出了10°E子午面上電離層電子密度隨緯度和高度的變化特性，其中圖2a為IRI2007模型計算出的結(jié)果，圖2b和圖2c分別為ASIRT和SIRT算法迭代相同次數(shù)后的反演結(jié)果.從中可以看出，ASIRT算法反演的電離層電子密度分布整體上與IRI2007模型計算的電離層電子密度分布符合得較好，且優(yōu)于SIRT算法，這說明利用ASIRT算法進行電離層電子密度反演是可行、有效的.

表1給出了兩種算法在不同經(jīng)度面上的反演結(jié)果誤差統(tǒng)計，從中可以看出，在不同經(jīng)度面上，利用ASIRT算法反演結(jié)果誤差均小于SIRT算法.統(tǒng)計結(jié)果顯示，利用ASIRT反演的重構(gòu)區(qū)域內(nèi)所有像素的電離層電子密度的平均絕對誤差約為2.1×109el／m3，最大電離層電子密度誤差的絕對值約為5.8×1010el／m3，SIRT反演的重構(gòu)區(qū)域內(nèi)所有像素的電離層電子密度的平均絕對誤差約為5.0×109el／m3，最大電離層電子密度誤差的絕對值約為9.8×1010el／m3，且ASIRT算法經(jīng)過9次迭代后收斂，SIRT算法經(jīng)過16次迭代后收斂，這說明了本文算法在收斂速度和反演結(jié)果精度上均有提高，從而證實了該算法的優(yōu)越性.

圖1 電離層層析幾何分布示意圖Fig.1 Schematic diagram of the computation domain of IED tomography

圖2 10°E子午面上電離層電子密度剖面圖（a）IRI2007模型獲取的；（b）ASIRT算法重構(gòu)的；（c）SIRT算法重構(gòu)的.Fig.2 IED profiles at 10°E by different algorithms（a）IRI2007model；（b）ASIRT algorithm；（c）SIRT algorithm.

為了進一步檢驗該算法，給出了高度面上的反演結(jié)果.圖3給出了兩種算法在300km高度面上電子密度分布及反演誤差分布，圖3b和3d分別為利用ASIRT算法反演的電子密度分布和誤差分布，圖3c和3e分別為利用SIRT算法反演的電子密度分布和誤差分布.從中可以看出，相對于SIRT算法，ASIRT算法反演的結(jié)果總體來說更加接近于真值.表2給出了兩種算法在不同高度面上的反演結(jié)果誤差統(tǒng)計，從中可以看出，在不同高度面上，利用ASIRT算法反演結(jié)果誤差均小于SIRT算法，這進一步說明了該ASIRT的優(yōu)越性.

表1 不同經(jīng)度面上兩種算法誤差統(tǒng)計分析表（單位：1010el／m3）Table 1 Error analysis of IED reconstructed by two algorithms at different longitudes（unit：1010el／m3）

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

利用實測數(shù)據(jù)對ASIRT算法進行了檢驗.重構(gòu)時間為2012年4月8日，重構(gòu)區(qū)域及格網(wǎng)劃分與上述模擬實驗一致.本文采用的觀測數(shù)據(jù)來自歐洲地區(qū)IGS觀測網(wǎng)絡(luò)，選取重構(gòu)區(qū)域內(nèi)的臺站觀測信息進行反演，并利用該區(qū)域內(nèi)的一個垂直探測站PRUHONICE（50.00°N，14.60°E）的觀測數(shù)據(jù)進行獨立檢核.該電離層垂直站安裝的是數(shù)字測高儀DPS－4D，發(fā)射天線為兩只交叉的雙三角天線，高36m；接收天線場為四交叉天線循環(huán)；標準測定設(shè)置為每15分鐘產(chǎn)生一幅0.05MHz分辨率的電離層圖.

表2 不同高度面上兩種算法誤差統(tǒng)計分析表（單位：1010el／m3）Table 2 Error analysis of IED reconstructed by the two algorithms at different heights（unit：1010el／m3）

圖4 TEC分布及重構(gòu)電離層電子密度分布（a）和（c）是9∶00UT；（b）和（d）是23∶00UT.（a）和（b）為TEC圖；（c）和（d）為15°E的垂直剖面圖.Fig.4 Images of TEC distribution and reconstruction IED（a）and（c）9∶00UT；（b）and（d）23∶00UT.（a）and（b）TEC；（c）and（d）Vertical cross－sections along 15°E.

如圖4所示，（a）和（c）給出了UT9∶00時的分布，（b）和（d）給出了UT23∶00時的分布，分別代表了白天和夜晚電離層電子密度的空間分布.（a）和（b）是實際空間的TEC分布，（c）和（d）是15°E的垂直剖面.從中可看出，利用ASIRT反演的電離層電子密度隨著緯度的增大而逐漸降低，大體上與實際空間的TEC分布走向相符合，說明該算法是可行的.為了驗證算法的可靠性，將反演結(jié)果與測高儀數(shù)據(jù)進行比較.圖5展示了不同時刻的電離層測高儀所得剖面與ASIRT和SIRT兩種算法反演結(jié)果的比較，（a）和（b）分別為UT9∶00時和UT23∶00時的電子密度剖面的比較結(jié)果.從比較結(jié)果來看，ASIRT反演結(jié)果整體上與測高儀數(shù)據(jù)更為接近.說明ASIRT反演的精度確實優(yōu)于SIRT.然而，在垂直方向上，這兩種算法反演的結(jié)果與測高儀數(shù)據(jù)均存在一定差異，原因是觀測信息不足導(dǎo)致垂直分辨率不高.

圖6展示了一天12個反演時段內(nèi)兩種算法層析反演獲得的F2層電子密度峰值（NmF2）和F2層電子密度峰值高度（hmF2）與電離層測高儀站觀測數(shù)據(jù)的比較結(jié)果.從中可以看出，兩種算法反演獲得的F2層電子密度峰值與測高儀觀測結(jié)果總體上符合較好，但ASIRT算法總體上更接近于測高儀觀測結(jié)果.然而，F(xiàn)2層的峰值高度與測高儀存在一定的差異.由于觀測噪聲、電離層空間離散誤差以及測站幾何結(jié)構(gòu)限制等因素，使得反演結(jié)果的垂直分辨率較差.這說明在電離層層析成像過程中僅僅通過改進反演算法來改善電子密度空間結(jié)構(gòu)（特別是垂直分辨率）是不夠的，利用其他技術(shù)增加用于反演的觀測信息是解決這一問題的最根本的手段.

圖5 重構(gòu)電離層電子密度剖面與電離層測高儀測量剖面的比較（（a）UT9∶00；（b）UT23∶00）Fig.5 Comparison of estimated IED profiles by two algorithms and IED profiles measured by ionosonde（（a）UT9∶00；（b）UT23∶00）

圖6 兩種算法反演的F2層電子密度峰值以及峰值高度與測高儀結(jié)果的比較Fig.6 Comparison of NmF2and hmF2values from CIT reconstruction by two algorithms and ionosonde data

5 結(jié)論

本文利用一種新的自適應(yīng)的聯(lián)合迭代重構(gòu)算法研究了電離層層析成像的問題，實驗結(jié)果表明，該算法在一定程度上能夠改善常用的聯(lián)合迭代重構(gòu)算法的反演效果.相對于常用的聯(lián)合迭代重構(gòu)算法，新算法的收斂速度和反演精度均有所提高.盡管該算法改善了電離層層析成像的效果，但仍然不能很好地提高層析成像的垂直分辨率，解決這一問題，需要進一步增加觀測信息.

Arikan O，Arikan F，Erol C B.2007.Computerized ionospheric tomography with the IRI model.Adv.SpaceRes.，39（5）：859－866.

Austen J R，F(xiàn)ranke S J，Liu C H，et al.1986.Application of computerized tomography techniques to ionospheric research.／／Tauriainen A ed.Radio Beacon Contributions to the Study of Ionization and Dynamics of the Ionosphere and to Corrections to Geodesy and Technical Workshop.Oulu，F(xiàn)inland，Part 1，25－35.

Austen J R，F(xiàn)ranke S J，Liu C H.1988.Ionospheric imaging using computerized tomography.RadioSci.，23（3）：299－307.

Chartier A T，Smith N D，Mitchell C N，et al.2012.The use of ionosondes in GPS ionospheric tomography at low latitudes.J.Geophys.Res.，117：A10326，doi：10.1029／2012JA018054.

Demmel J W.1997.Applied numerical linear algebra.SIAM：Philadelphia，PA.

Elfving T，Hansen P，Relaxation T.2012.Semiconvergence and relaxation parameters for projected SIRT algorithms.SIAMJ.Sci.Comput.，34（4）：A2000－A2017.

Hobiger T，Kondo T，Koyama Y.2008.Constrained simultaneous algebraic reconstruction technique（C－SART）—a new and simple algorithm applied to ionospheric tomography.Earth PlanetsSpace，60：727－735.

Lee J K，Kamalabadi F.2009.GPS－based radio tomography with edge－preserving regularization.IEEETrans.Geosci.RemoteSens.，47（1）：312－324.

Li H，Yuan Y，Li Z，et al.2012.Ionospheric electron concentration imaging using combination of LEO satellite data with groundbased GPS observations over China.IEEETrans.Geosci.RemoteSens.，50（5）：1728－1735.

Liu S，Wang J，Gao J.2010.Inversion of ionospheric electron density based on a constrained simultaneous iteration reconstruction technique.IEEETrans.Geosci.RemoteSens.，48（6）：2455－2459.

Liu Z Z.2004.Ionosphere tomographic modeling and applications using global positioning system（GPS）measurements［Ph.D.thesis］.Calgary：The University of Calgary.

Lytle R J.1979.Computerized geophysical tomography.Proc.IEEE，67（7）：1065－1073.

Ma X F，Maruyama T，Ma G，et al.2005.Three－dimensional ionospheric tomography using observation data of GPS ground receivers and ionosonde by neural network.J.Geophys.Res.，110：A05308，doi：10.1029／2004JA010797.

Nesterov I A，Kunitsyn V E.2011.GNSS radio tomography of the ionosphere：The problem with essentially incomplete data.Adv.SpaceRes.，47（10）：1789－1803.

Pryse S E，Kersley L.1992.A preliminary experimental test of ionospheric tomography.J.Atmos.Terr.Phys.，54（7－8）：1007－1012.

Pryse S E.2003.Radio tomography：a new experimental technique.Surv.Geophys.，24（1）：1－38.

Raymund T D，Austen J R，F(xiàn)ranke S J，et al.1990.Application of computerized tomography to the investigation of ionospheric structures.RadioSci.，25（5）：771－789.

Raymund T D.1994.Ionospheric tomography algorithms.Int.J.Imag.Syst.Tech.，5（2）：75－85.

Rius A，Ruffini G，Cucurull L.1997.Improving the vertical resolution of ionospheric tomography with GPS occultations.Geophys.Res.Lett.，24（18）：2291－2294.

Spencer P，Kersley L，Pryse S E.1998.A new solution to the problem of ionospheric tomography using quadratic programing.RadioSci.，33（3）：607－616.

Wen D B，Yuan Y B，Ou J，et al.2007.Three－dimensional ionospheric tomography by an improved algebraic reconstruction technique.GPSSolut.，11（4）：251－258.

Wen D B，Yuan Y B，Ou J，et al.2008.A hybrid reconstruction algorithm for 3－D ionospheric tomography.IEEETrans.Geosci.RemoteSens.，46（6）：1733－1739.

Zou Y H.2003.A study of time－dependent 3－D ionospheric tomography with ground－based GPS network and occultation observations［Ph.D.thesis］（in Chinese）.Wuhan：Wuhan University.

附中文參考文獻

鄒玉華.2003.GPS地面臺網(wǎng)和掩星觀測結(jié)合的時變?nèi)S電離層層析［博士論文］.武漢：武漢大學(xué).