肖宏波，門守強

（西安工業(yè)大學(xué) 理學(xué)院，西安710021）

電離層層析圖像重建（Computerized Ionospheric Tomography，CT）是電離層反演的一種新技術(shù)，可以有效檢測電離層大尺度結(jié)構(gòu)（如赤道異常、中緯槽等），且探測范圍廣、易于實現(xiàn)的電離層探測技術(shù).1986年University of Illinois的Austen et al［1］.首次提出了將計算機輔助層析成像技術(shù)與無線電信標(biāo)結(jié)合應(yīng)用于反演電離層電子密度的設(shè)想，這一新設(shè)想立即引發(fā)了國際上相繼開展電離層層析圖像成像的實驗和理論研究，為電離層探測提供了一種重要的遙感方法，并在近來的地震預(yù)測應(yīng)用中顯示了廣闊的應(yīng)用前景［2］.

在實際電離層CT問題中，待測區(qū)域、GPS衛(wèi)星軌道和地面接收臺站之間構(gòu)成特定的幾何關(guān)系，這種地基電離層CT實驗缺乏水平或接近水平方向的掃描射線以及有限的臺站密度，使得采集的投影數(shù)據(jù)嚴重不完整，屬于不完全數(shù)據(jù)重建問題，這是影響圖像反演重建質(zhì)量的主要因素［3］.乘法代數(shù)重建算法（Multiplicative Algebraic Reconstruction Techniques，MART）［4］是最常用的電離層CT圖像重建算法.對電離層反演區(qū)域進行離散，采用濾波反投影和經(jīng)過歸一化處理的理論值或經(jīng)驗電離層模式（如IRI－2001模式）相結(jié)合的辦法來確定迭代初始值［5－6］，MART 反演重建結(jié)果較好，但對于迭代初始值非常敏感，初始值的選取對反演結(jié)果的影響很大［7－8］.所以重建算法的合理選擇是圖像重建必須考慮的問題.文中利用單目標(biāo)無約束優(yōu)化方法（Steepest Descent Method，SD）［9－10］對重建數(shù)據(jù)進行迭代預(yù)處理，然后再利用MART算法進行反演重建，既可以減少電離層CT圖像重建算法MART對初始值的敏感，又沒有太多增加計算的復(fù)雜度與運算量，數(shù)值模擬及實測數(shù)據(jù)結(jié)果表明，減少了電離層電子密度反演誤差，提高了重建結(jié)果的可靠性.

1 反演算法

根據(jù)電波傳播理論［2］，電離層CT圖像重建技術(shù)測量的是沿信號射線路徑上的總電子含量（Total Electron Contents，TEC）［11－12］

其中：Ne為沿傳播路徑Path上任意一點的電子密度；r，θ為對應(yīng)的坐標(biāo).該方程表明只要測出一定區(qū)域內(nèi)不同射線路徑上的TEC值，即可重建出該區(qū)域的電子密度分布.在基于GPS衛(wèi)星信標(biāo)測量的電離層CT圖像重建中，TEC可利用GPS雙頻信標(biāo)的測碼偽距觀測值求得.

由臺站接收到的GPS衛(wèi)星信標(biāo)可以確定TEC的值，為了在計算機上用數(shù)值方法求解式（1）的積分方程，可以采用分塊的方法將反演區(qū)域離散化成N個網(wǎng)格，這樣上式就改寫成

其中：yi某一表示路徑上接收到的TEC；Xn為象素點（電子密度值）構(gòu)成的向量；Ym為TEC構(gòu)成的向量；m為射線路徑數(shù)；Amn即投影矩陣.在實際的電離層CT圖像重建中，由于投影數(shù)據(jù)嚴重不完整，因此利用迭代優(yōu)化算法求解式（1），重建出電離層電子密度分布，MART算法是最常用的一種.

在利用MART算法進行電離層電子密度反演重建時，由于該算法對初始值的敏感，使得迭代初始值的選取尤為重要，初始值的稍微偏差或不合適都會造成反演誤差過大甚至不符合實際［5－8］.為了減少迭代初始值選取對電離層對演重建質(zhì)量的影響，提高電離層CT圖像重建的質(zhì)量，本文將無約束優(yōu)化圖像重建算法SD與MART算法結(jié)合起來，以減少迭代初始值對反演結(jié)果的影響，提高電離層CT重建質(zhì)量.

在電離層CT反演重建問題中，反演區(qū)域沿垂直和緯向被均勻劃分為n個象素點，在每一個象素點上假設(shè)電子密度是均勻分布的，則可由衛(wèi)星位置及臺站位置計算出投影矩陣Amn.其迭代過程為

在目標(biāo)優(yōu)化為基礎(chǔ)的圖像重建問題中，首先考慮的是目標(biāo)函數(shù)應(yīng)當(dāng)使圖像的再投影盡可能的接近于實際的投影數(shù)據(jù)，這類問題歸結(jié)為范數(shù)極小問題，迭代算法將反演重建問題轉(zhuǎn)換為極小值問題.為了減少迭代初始值選取不合適造成的反演誤差，在獲得投影數(shù)據(jù)后，利用無約束圖像優(yōu)化重建算法中的最速下降法SD求解式（2），然后將其結(jié)果作為式（3）的迭代初始值進行反演重建.SD算法的迭代過程為

其中：k個迭代點為Xk，λk為最優(yōu)步長，令k＝0，1，2，3，…，就可以得到一個點列 X0，X1，X2，…，其中X0是初始點.

由于最速下降SD的無約束特征，所以首先選擇適當(dāng)?shù)某跏茧娮用芏确植甲鳛榈跏贾?，從式?）出發(fā)進行迭代求解.采用的是IRI－2001［13－14］國際電離層參考模型給出的電子密度作為初始分布，由這些分布出發(fā)利用SD進行迭代可以得到一簇重建結(jié)果，記為XSD；然后將利用SD算法反演得到的結(jié)果作為MART的初始迭代值，利用式（3）進行二次迭代求解，得到最終的電離層電子密度反演重建結(jié)果.即

與MART算法相比，SD算法對電子密度初始分布不敏感，這有利于迭代初始值的選取.改進算法首先考慮使用SD算法進行迭代求解.考慮到實際的電離層層析成像問題中投影數(shù)據(jù)嚴重缺少的實際情況，仍然使用MART算法進行反演重建，將迭代初值對反演的影響降到最低，使重建結(jié)果更符合實際的電子密度分布，可以提高電離層層析成像的質(zhì)量.由于改進算法在重建的過程中，進行了兩次迭代求解，增加了計算量，不過并沒有增加計算的復(fù)雜度.所以，與傳統(tǒng)的單獨使用MART算法進行反演重建相比，改進算法既利用了MART算法的優(yōu)點，又可以把初始迭代值對反演結(jié)果的影響降低，提高CIT反演重建質(zhì)量.

2 數(shù)值模擬及實測數(shù)據(jù)結(jié)果分析

2.1 數(shù)值模擬分析

在反演算法中誤差函數(shù)是一個重要的參數(shù)，反演誤差［15］可以定義為）

經(jīng)過多次數(shù)值模擬重建檢驗表明，文中提出的混合算法可以提高圖像反演重建的質(zhì)量.這里給出了其中一次的重建結(jié)果.IRI－2001作為參考電離層模式產(chǎn)生的電子密度分布及所用參數(shù)如圖1所示.反演區(qū)域水平方向被分成46個像素，垂直方向分為19個像素，垂直高度范圍為100～1000km，緯度范圍為北緯15°～60°.假設(shè)衛(wèi)星軌跡在給定的區(qū)域是水平的，并且每隔0.125°計算一個位置，另外假設(shè)對于每一個衛(wèi)星位置電離層是時不變的，電子密度在每一個像素內(nèi)是均勻的.

圖1 IRI－2001國際參考電離層模式產(chǎn)生的電離層分布Fig.1 Electron density created from IRI－2001model

單獨使用MART算法及文中改進算法的反演重建結(jié)果如圖2～3所示，單獨使用MART算法的誤差及使用本文算法的誤差見表1.比較圖1IRI－2001模型給出的電離層原始電子密度分布與圖2單獨使用MART算法重建結(jié)果可以看出，重建圖像在輪廓上反映了原始的電子密度分布，但對細節(jié)的重建質(zhì)量較差，不能反映電子密度的微小起伏變化.將圖3文中算法的重建結(jié)果與圖2比較可以看出，新算法重建結(jié)果不但反映了原始電子密度的大范圍變化，而且對電子密度分布的細節(jié)重建結(jié)果也較為理想.

圖2 MART算法反演結(jié)果Fig.2 Reconstructed Image from MART

圖3 文中改進算法反演結(jié)果Fig.3 Reconstructed Image from the improved algorithm

表1 單獨使用MART算法和本文改進算法的反演誤差比較Tab.1 Comparison of the error norms of MART and of the improved algorithm

在數(shù)值檢驗中，使用的是相同的迭代初始值，比較圖2和圖3可以看出，相對于單獨使用MART算法，改進算法不但對電離層電子密度大尺度分布反演結(jié)果很好，而且電子密度分布的細節(jié)重建結(jié)果也較為理想.文中提出的改進算法得到的重建結(jié)果更能精確地重現(xiàn)原始模型中的電子密度分布特征，提高了重建圖像的可靠性.

2.2 實測數(shù)據(jù)分析

利用文獻實測TEC數(shù)據(jù)進行了計算，并與實際觀測結(jié)果進行了比較.TEC數(shù)據(jù)見文獻［12］，5個接收點由南至北分布，跨越電離層中緯槽.圖4為利用MART算法重建的電子密度的二維分布，圖5為利用本文提出的混合算法反演得結(jié)果.

圖4 MART反演重建的電子密度剖面Fig.4 Reconstructed Image from MART

圖5 混合算法反演得到的電子密度剖面Fig.5 Reconstructed Image from the Improved algorithm

從圖4～5中可以看出，電離層中緯槽的位置通過兩種CT反演算法能夠較好的體現(xiàn)，通過將反演的電離層二維剖面與同時的非相干散射雷達探測得到的剖面［12］比較，符合較好.但本文提出的混合算法優(yōu)于MART算法的反演結(jié)果.比較圖4和5可見，利用MART反演算法重構(gòu)的電離層剖面，清楚地顯示了中緯槽的位置，但是可以看出，在緯度約67°，500km高度上反演結(jié)果中出現(xiàn)了一個電子密度的畸形區(qū)域，而且隨高度梯度有一定的畸變；使用混合算法后，由圖5可以看出，電子密度的畸形趨于明顯減少，且中緯槽的幅度也得到了很好的重構(gòu).

3 結(jié)論

文中基于圖像優(yōu)化重建算法的思想，聯(lián)合使用SD算法和MART算法，提出了一種改進的電離層層析圖像重建算法.該算法利用SD算法的無約束特征對投影數(shù)據(jù)進行預(yù)處理，然后將處理結(jié)果進行MART迭代以減少初始值不完全對重建結(jié)果的影響.數(shù)值模擬結(jié)果對比分析可以看出，利用相同的迭代初始值，文中算法的反演誤差從約14%減少到11%.通過實測數(shù)據(jù)的兩種反演結(jié)果與同期雷達探測數(shù)據(jù)比較，相對于單獨使用MART算法，迭代初始值對改進算法的影響減小，電子密度反演重建的精度及與原始分布的一致性都高于MART算法，特別對中緯槽位置的分布與實測結(jié)果吻合較好.不過反演誤差的減少還不夠理想，還需要對理論模型做進一步分析，并結(jié)合更多的實測結(jié)果進行研究.

［1］JEFFREY R A，STEVEN J F，LIU C H.Ionospheric Imaging Using Computerized Tomography［J］.Radio Science，1988，23（3）：299.

［2］楊劍，吳云，周義炎.基于電離層層析成像技術(shù)探測汶川地震前電離層異常［J］.大地測量與地球動力學(xué)，2011，31（1）：9.YANG Jian，WU Yun，ZHOU Yi－yan.Probe into Seismo－Ionospheric Anomaly of Wenchuan Ms8.0 Earthquake Based on Computerized Ionospheric Tomography［J］.Journal of Geodesy and Geodynamics，2011，31（1）：9.（in Chinese）

［3］徐繼生，馬淑英.電離層無線電層析［C］／／濮祖蔭.空間物理前沿進展.北京：氣象出版社，1998.XU Ji－sheng，MA Shu－ying.Radio Tomography of the Ionosphere［C］／／Pu Zu－yin.Developments in Space Physics Front.Beijing：Meteorological Publication，1998.（in Chinese）

［4］CENOR Y.Finite Series－Expansion Reconstruction Methods［J］.Proc IEEE，1983，71（3）：409.

［5］YIZENGAW E，MOLDWIN MB，DYSON P L，et al.First Tomography Image of Ionospheric Outflows［J］.Geophysical Research Letters，2006，33（20）：1.

［6］PRYSE S E.Radio Tomography：A New Experimental Technique［J］.Sury Geophys，2003，24（1）：1.

［7］YAEVZ Y，ARIKAN F，ARKAN O.A Hybrid Reconstruction Algorithm for Computerized Ionospheric Tomography［J］.Recent Advances in Space Technologies，IEEE，2005，9（11）：782.

［8］TSAI L C.Tomographic Imaging of the Ionosphere Using the GPS／MET and NNSS Data［J］.J.Atmos Sol Terr Phys，2002，64（18）：2003.

［9］EGGERMONT P P B，HERMAN G T，LENT A.Iterative Algorithms for Large Partitioned Linear Systems with Applications to Image Reconstruction［J］.Linear Algebra and Its Applications，1981，40：37.

［10］李化欣，潘金孝.最速下降法在圖像重建中的應(yīng)用［J］.科技情報開發(fā)與經(jīng)濟，2006，16（3）：155.LI Hua－xin，PAN Jin－xiao.Application of the Steepest Descent Method in Image Reconstruction［J］.Sci－Tech Information Development ＆ Economy，2006，16（3）：155.（in Chinese）

［11］徐繼生，馬淑英，陽其罕，等.東亞赤道異常區(qū)電離層CT診斷－實驗及初步結(jié)果［J］.地球物理學(xué)報，1995，38（5）：553.XU Ji－sheng，MA Shu－ying，YANG Qi－h(huán)an，et al.Tomographic Diagnosis of Ionosphere in the Region of East－Asian Equatorial Anomaly：An Experimental Campaign and Some Early Results［J］.Chinese Journal of Geophysics，1995，38（5）：553.（in Chinese）

［12］RAYMUND T D，JEFFREY R A，F(xiàn)RANKE S J，et al.Application of Computerized Tomography to the Investigation of Ionospheric Structures［J］.Radio Science，1990，25（5）：771.

［13］ANTONIO R，GIULIO R，AUGUST R.Analysis of Ionospheric Electron Density Distribution from GPS／MET Occupations ［J］.IEEE，1998，36（2），383.

［14］Helen R Na.Resolution Degradation Parameters of Ionospheric Tomography［J］.Radio Science，1994，29（1）：115.

［15］GIULIO R，ALEJANDRO F，ANTONIO R.GPS Tomography of the Ionospheric Electron Content with a Correlation Functiona［J］.IEEE，1998，36（1）：143.