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

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

2武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430079

1 引言

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

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

2 電離層層析成像原理

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

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

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

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

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

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

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

3.1 算法原理

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

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

式中，M和R為對(duì)稱正定矩陣.該算法中的松弛因子λk對(duì)于迭代過程中結(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具有非負(fù)性.以下推導(dǎo)過程中，向量都是列向量，ai為系數(shù)矩陣A的第i行的轉(zhuǎn)置，內(nèi)積的標(biāo)準(zhǔn)形式為〈x，y〉＝xTy，ρ（·）為譜半徑（最大正特征值）.對(duì)于每一個(gè)線性方程（14），定義超平面Hi＝｛x∈Rn｜〈ai，x〉＝y(tǒng)i｝.

如果

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 重構(gòu)電離層電子密度剖面與電離層測(cè)高儀測(cè)量剖面的比較（（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層電子密度峰值以及峰值高度與測(cè)高儀結(jié)果的比較Fig.6 Comparison of NmF2and hmF2values from CIT reconstruction by two algorithms and ionosonde data

5 結(jié)論

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

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.

附中文參考文獻(xiàn)

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