电离层层析图像重建的改进算法＊

肖宏波，门守强

（西安工业大学 理学院，西安710021）

电离层层析图像重建（Computerized Ionospheric Tomography，CT）是电离层反演的一种新技术，可以有效检测电离层大尺度结构（如赤道异常、中纬槽等），且探测范围广、易于实现的电离层探测技术.1986年University of Illinois的Austen et al［1］.首次提出了将计算机辅助层析成像技术与无线电信标结合应用于反演电离层电子密度的设想，这一新设想立即引发了国际上相继开展电离层层析图像成像的实验和理论研究，为电离层探测提供了一种重要的遥感方法，并在近来的地震预测应用中显示了广阔的应用前景［2］.

在实际电离层CT问题中，待测区域、GPS卫星轨道和地面接收台站之间构成特定的几何关系，这种地基电离层CT实验缺乏水平或接近水平方向的扫描射线以及有限的台站密度，使得采集的投影数据严重不完整，属于不完全数据重建问题，这是影响图像反演重建质量的主要因素［3］.乘法代数重建算法（Multiplicative Algebraic Reconstruction Techniques，MART）［4］是最常用的电离层CT图像重建算法.对电离层反演区域进行离散，采用滤波反投影和经过归一化处理的理论值或经验电离层模式（如IRI－2001模式）相结合的办法来确定迭代初始值［5－6］，MART 反演重建结果较好，但对于迭代初始值非常敏感，初始值的选取对反演结果的影响很大［7－8］.所以重建算法的合理选择是图像重建必须考虑的问题.文中利用单目标无约束优化方法（Steepest Descent Method，SD）［9－10］对重建数据进行迭代预处理，然后再利用MART算法进行反演重建，既可以减少电离层CT图像重建算法MART对初始值的敏感，又没有太多增加计算的复杂度与运算量，数值模拟及实测数据结果表明，减少了电离层电子密度反演误差，提高了重建结果的可靠性.

1 反演算法

根据电波传播理论［2］，电离层CT图像重建技术测量的是沿信号射线路径上的总电子含量（Total Electron Contents，TEC）［11－12］

其中：Ne为沿传播路径Path上任意一点的电子密度；r，θ为对应的坐标.该方程表明只要测出一定区域内不同射线路径上的TEC值，即可重建出该区域的电子密度分布.在基于GPS卫星信标测量的电离层CT图像重建中，TEC可利用GPS双频信标的测码伪距观测值求得.

由台站接收到的GPS卫星信标可以确定TEC的值，为了在计算机上用数值方法求解式（1）的积分方程，可以采用分块的方法将反演区域离散化成N个网格，这样上式就改写成

其中：yi某一表示路径上接收到的TEC；Xn为象素点（电子密度值）构成的向量；Ym为TEC构成的向量；m为射线路径数；Amn即投影矩阵.在实际的电离层CT图像重建中，由于投影数据严重不完整，因此利用迭代优化算法求解式（1），重建出电离层电子密度分布，MART算法是最常用的一种.

在利用MART算法进行电离层电子密度反演重建时，由于该算法对初始值的敏感，使得迭代初始值的选取尤为重要，初始值的稍微偏差或不合适都会造成反演误差过大甚至不符合实际［5－8］.为了减少迭代初始值选取对电离层对演重建质量的影响，提高电离层CT图像重建的质量，本文将无约束优化图像重建算法SD与MART算法结合起来，以减少迭代初始值对反演结果的影响，提高电离层CT重建质量.

在电离层CT反演重建问题中，反演区域沿垂直和纬向被均匀划分为n个象素点，在每一个象素点上假设电子密度是均匀分布的，则可由卫星位置及台站位置计算出投影矩阵Amn.其迭代过程为

在目标优化为基础的图像重建问题中，首先考虑的是目标函数应当使图像的再投影尽可能的接近于实际的投影数据，这类问题归结为范数极小问题，迭代算法将反演重建问题转换为极小值问题.为了减少迭代初始值选取不合适造成的反演误差，在获得投影数据后，利用无约束图像优化重建算法中的最速下降法SD求解式（2），然后将其结果作为式（3）的迭代初始值进行反演重建.SD算法的迭代过程为

其中：k个迭代点为Xk，λk为最优步长，令k＝0，1，2，3，…，就可以得到一个点列 X0，X1，X2，…，其中X0是初始点.

由于最速下降SD的无约束特征，所以首先选择适当的初始电子密度分布作为迭代初始值，从式（4）出发进行迭代求解.采用的是IRI－2001［13－14］国际电离层参考模型给出的电子密度作为初始分布，由这些分布出发利用SD进行迭代可以得到一簇重建结果，记为XSD；然后将利用SD算法反演得到的结果作为MART的初始迭代值，利用式（3）进行二次迭代求解，得到最终的电离层电子密度反演重建结果.即

与MART算法相比，SD算法对电子密度初始分布不敏感，这有利于迭代初始值的选取.改进算法首先考虑使用SD算法进行迭代求解.考虑到实际的电离层层析成像问题中投影数据严重缺少的实际情况，仍然使用MART算法进行反演重建，将迭代初值对反演的影响降到最低，使重建结果更符合实际的电子密度分布，可以提高电离层层析成像的质量.由于改进算法在重建的过程中，进行了两次迭代求解，增加了计算量，不过并没有增加计算的复杂度.所以，与传统的单独使用MART算法进行反演重建相比，改进算法既利用了MART算法的优点，又可以把初始迭代值对反演结果的影响降低，提高CIT反演重建质量.

2 数值模拟及实测数据结果分析

2.1 数值模拟分析

在反演算法中误差函数是一个重要的参数，反演误差［15］可以定义为）

经过多次数值模拟重建检验表明，文中提出的混合算法可以提高图像反演重建的质量.这里给出了其中一次的重建结果.IRI－2001作为参考电离层模式产生的电子密度分布及所用参数如图1所示.反演区域水平方向被分成46个像素，垂直方向分为19个像素，垂直高度范围为100～1000km，纬度范围为北纬15°～60°.假设卫星轨迹在给定的区域是水平的，并且每隔0.125°计算一个位置，另外假设对于每一个卫星位置电离层是时不变的，电子密度在每一个像素内是均匀的.

图1 IRI－2001国际参考电离层模式产生的电离层分布Fig.1 Electron density created from IRI－2001model

单独使用MART算法及文中改进算法的反演重建结果如图2～3所示，单独使用MART算法的误差及使用本文算法的误差见表1.比较图1IRI－2001模型给出的电离层原始电子密度分布与图2单独使用MART算法重建结果可以看出，重建图像在轮廓上反映了原始的电子密度分布，但对细节的重建质量较差，不能反映电子密度的微小起伏变化.将图3文中算法的重建结果与图2比较可以看出，新算法重建结果不但反映了原始电子密度的大范围变化，而且对电子密度分布的细节重建结果也较为理想.

图2 MART算法反演结果Fig.2 Reconstructed Image from MART

图3 文中改进算法反演结果Fig.3 Reconstructed Image from the improved algorithm

表1 单独使用MART算法和本文改进算法的反演误差比较Tab.1 Comparison of the error norms of MART and of the improved algorithm

在数值检验中，使用的是相同的迭代初始值，比较图2和图3可以看出，相对于单独使用MART算法，改进算法不但对电离层电子密度大尺度分布反演结果很好，而且电子密度分布的细节重建结果也较为理想.文中提出的改进算法得到的重建结果更能精确地重现原始模型中的电子密度分布特征，提高了重建图像的可靠性.

2.2 实测数据分析

利用文献实测TEC数据进行了计算，并与实际观测结果进行了比较.TEC数据见文献［12］，5个接收点由南至北分布，跨越电离层中纬槽.图4为利用MART算法重建的电子密度的二维分布，图5为利用本文提出的混合算法反演得结果.

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

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

从图4～5中可以看出，电离层中纬槽的位置通过两种CT反演算法能够较好的体现，通过将反演的电离层二维剖面与同时的非相干散射雷达探测得到的剖面［12］比较，符合较好.但本文提出的混合算法优于MART算法的反演结果.比较图4和5可见，利用MART反演算法重构的电离层剖面，清楚地显示了中纬槽的位置，但是可以看出，在纬度约67°，500km高度上反演结果中出现了一个电子密度的畸形区域，而且随高度梯度有一定的畸变；使用混合算法后，由图5可以看出，电子密度的畸形趋于明显减少，且中纬槽的幅度也得到了很好的重构.

3 结论

文中基于图像优化重建算法的思想，联合使用SD算法和MART算法，提出了一种改进的电离层层析图像重建算法.该算法利用SD算法的无约束特征对投影数据进行预处理，然后将处理结果进行MART迭代以减少初始值不完全对重建结果的影响.数值模拟结果对比分析可以看出，利用相同的迭代初始值，文中算法的反演误差从约14%减少到11%.通过实测数据的两种反演结果与同期雷达探测数据比较，相对于单独使用MART算法，迭代初始值对改进算法的影响减小，电子密度反演重建的精度及与原始分布的一致性都高于MART算法，特别对中纬槽位置的分布与实测结果吻合较好.不过反演误差的减少还不够理想，还需要对理论模型做进一步分析，并结合更多的实测结果进行研究.

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

［2］杨剑，吴云，周义炎.基于电离层层析成像技术探测汶川地震前电离层异常［J］.大地测量与地球动力学，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］李化欣，潘金孝.最速下降法在图像重建中的应用［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］徐继生，马淑英，阳其罕，等.东亚赤道异常区电离层CT诊断－实验及初步结果［J］.地球物理学报，1995，38（5）：553.XU Ji－sheng，MA Shu－ying，YANG Qi－han，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，FRANKE 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.