0 引言

随着数值模式能力以及计算机水平的不断提高，数值天气预报对人类生活的指导作用也日益加强.在这种大趋势下，GPS掩星探测技术应运而生，这是一种新兴的探测地球环境参数的方法，具有全天候、低费用、无系统长期漂移等优点，其不受天气状况影响的特点在天气预报及气候的应用中具有独特的优势.其基本原理是，由于大气的不均匀性，电磁波传播时会发生信号延迟(相对于真空环境)，而这种延迟与大气状态参数有关，因此，通过对信号延迟的观测，就可以反演出大气环境参数.与其他卫星资料相比，GPS掩星观测资料能提供高精度、高分辨率、全球准均匀覆盖的地球电离层和中性层大气参数剖面(Melbourne et al., 1994；严豪健等，2007；Jin et al., 2011).

在掩星探测数据处理流程中，Abel反演是一个重要的环节，它把弯角转化成折射指数，由给定的折射指数数据n=n(r)(r为半径)通过下面的公式(称为Abel积分，其中a为影响参数，x=n(r))来计算电磁波传播的弯角序列：

(1)

其逆变换称为Abel变换，即通过给定的弯角序列α(a)来计算折射指数n：

(2)

Abel积分方程属于第一类Volterra型积分方程，最早是由Abel在19世纪末20世纪初在研究质点力学问题时提出的，后来在力学、光谱学、地震学、等离子物理及散射理论等许多领域中都得到了广泛的应用(刘克玲，1982；Gorenflo and Vessella, 1991；Groetsch，1993；Bukhgeim，1999；林剑等，2009；周义炎，2012)).众所周知，第一类Volterra型积分方程是不适定的(黄思训和伍荣生，2011)，尤其其不稳定性为实际应用带来了很大的困难，在实际问题中常常利用Tikhonov正则化进行处理(Vauhkonen et al., 1998；Zhao et al., 2011；Yang and Fu, 2013；崔岩和王彦飞，2015)；同时，无论是Abel积分方程还是Abel变换，在实际应用中都面临着如何进行高精度离散的问题.正因为如此，Abel积分方程(包括Abel变换)引发了许多学者的兴趣(Minerbo and Levy, 1969；Wing，1992；普小云等，1996；Steiner et al., 1999；Avazzadeh et al., 2011；吕涛和黄晋，2013).

从数学物理反问题的角度来看，在弯角α(a)已知的情况下，反演折射指数n(r)有两条途径：①利用Abel变换(2)直接计算；②把Abel积分方程(1)作为反问题求解.由于弯角α(a)是通过观测计算得到的，不可避免地含有误差，这部分误差在积分计算过程中很可能会被放大.因此，如果直接利用Abel变换(2)进行计算，所得到的结果是有问题的(肖庭延和宋金来，2000).如果利用Abel积分方程，首先需要将积分方程离散化，变成线性方程组再进行求解.这里有两个问题需要解决：一是被积函数具有弱奇性，同时，观测数据是离散的，导致了ln(n)可能不够光滑，因此，不能进行简单离散，需要找到合适的方法来定义导数dln(n)/dx以获得积分的高精度数值离散结果; 二是由于Abel积分方程(1)的不适定性，导致离散化后得到的线性方程组即线性反问题仍然具有不适定性，为解决该问题，常用的是正则化技术，其中以Tikhonov正则化技术(Tikhonov and Arsenin, 1977)的应用最为广泛.该技术通过添加正则化项η‖x‖²来实现，但当解具有不连续点时，该方法反演得到的近似解将会与真解产生较大偏差.在实际观测中，由于电离层及水汽作用的影响，弯角数据可能存在高频成分，从而导致剧烈的梯度变化，另外，电离层的影响可能导致解间断，此时，采用Tikhonov正则化方法会产生较大误差.作为处理这类问题有效的方法，Bayes算法被提了出来.

本文利用Bayes算法，基于Abel积分方程来反演折射指数.首先利用EGOPS软件仿真出150 km以下的弯角、碰撞参数及折射率等数据信息，并将60 km以下的折射率数据作为真实值; 然后采用Tikhonov正则化(通过L-曲线来确定最优的正则化参数η)来反演折射指数，并将该技术与Bayes算法的反演结果进行对比，依此来探究当原始数据具有误差时，Bayes算法是否可以得到优于Tikhonov正则化的反演结果.

1 Bayes算法的实施

Bayes估计是将Bayes理论应用于参数估计的结果，是Bayes统计的主要部分.其基本思想是：通过“归纳推理”的统计方法，利用研究对象的先验信息形成先验分布，并结合观测得到的样本信息利用Bayes公式得到未知量的后验概率，再根据后验概率大小进行决策分类.

在求解第一类积分方程的过程中，相比于Tikhonov正则化以及全变差方法，Bayes算法具有很明显的优势，即无论待定的函数是连续的还是具有间断点，均能得到比较令人满意的反演结果(陆璐，2008).

我们将Abel积分方程(1)中的α(a_i)记为α_i，i=1, …, m，积分上限记为a₁，根据大气折射指数与折射率的关系式

(3)

利用等价无穷小代换，将lnn近似为10^－6·N，并记为，将积分区间分为n层，每一层的折射指数对碰撞参数的梯度变化设为常数，因此，方程(1)可以离散表示为

(4)

经过一系列处理，我们将方程(1)离散化为

(5)

其中，矩阵A的分量为A_ik，

(6)

向量e表示观测误差矩阵，满足Gauss型分布：

(7)

其中，σ表示该分布的标准差，可以是一个数，也可以是一个矩阵.

定义一个向量ξ，其分量为

(8)

ξ与具有如下关系：

(9)

此时，方程(5)化为

(10)

我们把对的估计转化为对ξ的估计，考虑Gauss型的先验密度，针对向量ξ构建先验模型.引入超参数λ_j(hyperparameter)，用来描述不连续点的位置与大小.根据Bayes公式，各参数之间的依赖关系可以用概率表示为(Calvetti and Somersalo, 2007)

(11)

其中，π_prior表示先验密度.

在实际计算中，令k从1开始，给λ赋初值λ¹，然后按如下步骤进行迭代，寻找ξ和λ的最大后验估计值.

第一步：从上式解得ξ^k，使其最大化π(ξ|λ^k, y)，

(12)

即

(13)

其中，Λ=diag(λ^k).

第二步：解λ^k+1，使其最大化π(λ|ξ^k, y)，

(14)

即

(15)

第三步：计算ξ^k+1以及它与上次迭代结果的差|ξ^k+1－ξ^k|.如果小于事先设定的误差水平δ，则输出ξ^k+1和λ^k+1；否则，回到第二步.

整个Bayes算法过程如图 1所示.

2 数值实验 2.1 数据来源

为了探讨不同算法在利用弯角反演折射指数时的效果，我们利用EGOPS软件进行数据仿真.将Abel积分方程(1)中的积分上限截断到150 km处，利用该软件仿真出该高度以下的弯角、碰撞参数以及折射率等数据信息，由于电离层、水汽等的影响, 以及原始弯角数据的数量级变化，我们在弯角数据中加入误差水平为8%的相对误差，以及10^-6与10^-4的绝对误差，然后分别利用Bayes算法以及Tikhonov正则化技术进行反演，将反演结果分别与EGOPS给出的60 km以下的折射率廓线进行比较，并通过50次数值试验，给出了两种方法反演结果的平均误差随高度变化的廓线图.

2.2 反演结果

首先，我们直接利用EGOPS仿真出的弯角数据进行反演，即认为此时的观测是没有误差的，结果如图 2所示.

图 2表示直接利用EGOPS仿真的弯角数据进行反演时，Bayes算法与Tikhonov正则化反演的折射指数的对比图，其中，蓝色实线代表的是EGOPS软件给出的根据其仿真出的弯角数据计算得到的60 km以下的折射指数廓线(我们在这里认为是真实值)，红色的虚线代表Bayes算法的反演结果，绿色的虚线代表Tikhonov正则化的反演结果.上图表明，Bayes算法与Tikhonov正则化方法的反演结果一致，二者的均方根误差均为2.5470×10^-8.

下面我们给出EGOPS仿真出的用于反演计算的弯角数据以及分别加入了8%、10^-6以及10^-4的随机误差的弯角数据示意图，如图 3、图 4、图 5所示.是否能利用该弯角数据较好地反演出折射指数，是我们对Bayes算法以及Tikhonov正则化技术的一个考察方面.

将作为Tikhonov正则化中的罚项，利用图 6所示的L-曲线，确定正则化参数η=0.0001.

利用EGOPS软件模拟的数据，针对所加入的8%的相对误差以及两组绝对误差(10^-6与10^-4)，我们分别给出了Bayes算法及Tikhonov正则化反演的折射指数廓线，如图 7、8、9所示.

图 7、8、9表示在EGOPS仿真的弯角数据中加入随机误差后，Bayes算法与Tikhonov正则化反演出的折射指数廓线对比图，其中图 7对应的是8%的相对误差，图 8、9分别对应的是10^-6与10^-4的绝对误差，可以看出，Bayes算法反演得到的折射指数廓线更接近真实廓线.对于反演结果的评估，我们进行了50次数值试验，给出了在EGOPS仿真的弯角数据中加入8%的相对误差后，Bayes算法与Tikhonov正则化方法反演的折射指数的平均误差随高度变化的廓线，如图 10、11、12所示.

图 10、11分别代表利用Bayes算法以及Tikhonov正则化方法进行反演时，所得结果的平均误差随高度的变化廓线，为了对比的方便，我们绘制了图 12，其中实线代表Bayes算法对应的平均误差，虚线代表Tikhonov正则化技术对应的平均误差，通过这几张廓线图，我们有理由相信，在处理梯度变化很大、甚至有间断情况的数据时，Bayes算法比Tikhonov正则化技术更具优势.

3 总结

对于折射指数的计算，无论是将其作为正问题，直接利用Abel变换(2)直接求解，还是作为反问题，利用Abel积分(1)进行求解，都面临积分方程的不适定问题.文本将Bayes算法应用到Abel积分方程(1)中，利用EGOPS软件仿真出的150 km以下的弯角数据进行反演计算，并将60 km以下的折射指数廓线作为真实廓线.当弯角数据中不含误差时，Bayes算法与Tikhonov正则化技术的反演结果具有很好的一致性，两种方法的均方根误差均为2.5470×10^-8；当在弯角数据中分别加入8%的相对误差以及10^-6与10^-4的绝对误差后，相比于Tikhonov正则化技术，Bayes算法可以反演出更加接近真实廓线的折射指数，反演结果精度也更高.我们知道，在实际情况下，由于各方面客观因素的影响，由观测得到的数据信息往往不可避免地含有误差，因此，可以考虑利用Bayes算法，得到更加令人满意的反演结果.

致谢  感谢两位审稿人提出的宝贵意见.