肖宏波，門守強

（西安工業(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ù)嚴(yán)重不完整，屬于不完全數(shù)據(jù)重建問題，這是影響圖像反演重建質(zhì)量的主要因素［3］.乘法代數(shù)重建算法（Multiplicative Algebraic Reconstruction Techniques，MART）［4］是最常用的電離層CT圖像重建算法.對電離層反演區(qū)域進(jìn)行離散，采用濾波反投影和經(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ù)進(jìn)行迭代預(yù)處理，然后再利用MART算法進(jìn)行反演重建，既可以減少電離層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án)重不完整，因此利用迭代優(yōu)化算法求解式（1），重建出電離層電子密度分布，MART算法是最常用的一種.

在利用MART算法進(jìn)行電離層電子密度反演重建時，由于該算法對初始值的敏感，使得迭代初始值的選取尤為重要，初始值的稍微偏差或不合適都會造成反演誤差過大甚至不符合實際［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）的迭代初始值進(jìn)行反演重建.SD算法的迭代過程為

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

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

與MART算法相比，SD算法對電子密度初始分布不敏感，這有利于迭代初始值的選取.改進(jìn)算法首先考慮使用SD算法進(jìn)行迭代求解.考慮到實際的電離層層析成像問題中投影數(shù)據(jù)嚴(yán)重缺少的實際情況，仍然使用MART算法進(jìn)行反演重建，將迭代初值對反演的影響降到最低，使重建結(jié)果更符合實際的電子密度分布，可以提高電離層層析成像的質(zhì)量.由于改進(jìn)算法在重建的過程中，進(jìn)行了兩次迭代求解，增加了計算量，不過并沒有增加計算的復(fù)雜度.所以，與傳統(tǒng)的單獨使用MART算法進(jìn)行反演重建相比，改進(jìn)算法既利用了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算法及文中改進(jìn)算法的反演重建結(jié)果如圖2～3所示，單獨使用MART算法的誤差及使用本文算法的誤差見表1.比較圖1IRI－2001模型給出的電離層原始電子密度分布與圖2單獨使用MART算法重建結(jié)果可以看出，重建圖像在輪廓上反映了原始的電子密度分布，但對細(xì)節(jié)的重建質(zhì)量較差，不能反映電子密度的微小起伏變化.將圖3文中算法的重建結(jié)果與圖2比較可以看出，新算法重建結(jié)果不但反映了原始電子密度的大范圍變化，而且對電子密度分布的細(xì)節(jié)重建結(jié)果也較為理想.

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

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

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

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

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

利用文獻(xiàn)實測TEC數(shù)據(jù)進(jìn)行了計算，并與實際觀測結(jié)果進(jìn)行了比較.TEC數(shù)據(jù)見文獻(xiàn)［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)，通過將反演的電離層二維剖面與同時的非相干散射雷達(dá)探測得到的剖面［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算法，提出了一種改進(jìn)的電離層層析圖像重建算法.該算法利用SD算法的無約束特征對投影數(shù)據(jù)進(jìn)行預(yù)處理，然后將處理結(jié)果進(jìn)行MART迭代以減少初始值不完全對重建結(jié)果的影響.數(shù)值模擬結(jié)果對比分析可以看出，利用相同的迭代初始值，文中算法的反演誤差從約14%減少到11%.通過實測數(shù)據(jù)的兩種反演結(jié)果與同期雷達(dá)探測數(shù)據(jù)比較，相對于單獨使用MART算法，迭代初始值對改進(jìn)算法的影響減小，電子密度反演重建的精度及與原始分布的一致性都高于MART算法，特別對中緯槽位置的分布與實測結(jié)果吻合較好.不過反演誤差的減少還不夠理想，還需要對理論模型做進(jìn)一步分析，并結(jié)合更多的實測結(jié)果進(jìn)行研究.

［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］／／濮祖蔭.空間物理前沿進(jìn)展.北京：氣象出版社，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)濟(jì)，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.