一种基于共轭梯度法的电离层层析成像混合反演方法

摘要

一种基于共轭梯度法的电离层层析成像混合反演方法。该发明基于考虑散射效应的绝对和相对总电子含量计算，采用乘法代数重构法和共轭梯度法结合的混合反演方法，基于电离层层析成像技术对电子密度空间分布进行高精度反演。首先，采用三频信标数据，考虑电离层散射效应，在绝对和相对总电子密度反演中均加入散射效应时延项，通过三频差分获得电离层总电子含量；其次，选定电离层待反演区，结合空间几何计算，得到层析成像的投影矩阵；最后，采用乘法代数重构算法和共轭梯度法相结合的混合反演方法得到电离层电子密度分布。本发明具有精度高，可靠性好等优点，可用于针对电离层电子密度的大气空间监测与地震灾害预测等。

说明书

技术领域

本发明涉及一种基于共轭梯度法的电离层层析成像混合反演方法，可用于投影矩阵 秩亏情况下的电离层电子密度高精度反演，如大气探测与电离层异常监测系统等。

背景技术

电离层是地球表面以上约60km至1000km的大气层，由分子、原子、电子及正、 负离子组成，属于日地空间环境的一部分。它的状态受到太阳及地球活动的影响，与人 类的生存、通信等活动息息相关。电子密度作为电离层的主要参数之一，其监测对了解 电离层状态具有重要意义，而基于卫星信标技术的电离层层析成像技术是获取电离层数 据的关键手段。但由于观测环境和观测视角有限等条件的限制，在待重构的区域内，信 标的地面接收机布设常存在数量稀少、分布不均匀等问题，导致观测信息的不足，投影 矩阵秩亏，因此对于弥补反演数据稀疏、解决层析成像的不适定性一直是电离层层析成 像的一个研究热点。当前广泛使用的基于卫星信标的电离层层析成像反演方法主要包 括：迭代重构算法(代数迭代重构算法、乘法代数重构算法、同时迭代重构算法等)， 模式参数拟合算法，傅里叶变换算法、卡尔曼滤波算法等。其中，傅里叶变换法由于投 影不完整而难以取得理想结果，卡尔曼滤波法则由于观测中存在的数值不稳定和观测数 据粗差等容易导致滤波发散而失效。对于观测信息不足导致的投影矩阵秩亏问题，模式 参数拟合法和迭代重构法理论上取得了比较完善的解决方法，但前者依赖于高精度的电 离层模型，后者在观测信息不足较大，数据稀疏严重的情况下，有限的迭代修正处理难 以实现高精度的电离层电子密度反演。

发明内容

本发明的技术解决问题是：针对观测数据不足、投影矩阵秩亏造成的电离层层析成 像的不适定性和可靠性等问题，基于考虑散射效应的绝对和相对总电子含量反演及乘法 代数重构法和共轭梯度法相结合的混合反演方法，反演计算高精度的电离层电子密度空 间分布。

本发明的技术解决方案为：一种基于共轭梯度法的电离层层析成像混合反演方法， 包括以下步骤：

首先，采用三频信标数据，考虑电离层散射效应，在绝对和相对总电子含量反演中 均加入散射效应时延项，通过三频差分获得电离层总电子含量；其次，选定电离层待反 演区，结合空间几何计算，得到层析成像的投影矩阵；最后，采用乘法代数重构算法和 共轭梯度法相结合的混合反演方法得到电离层电子密度空间分布；具体步骤如下：

第一步，基于三频信标观测的电离层总电子含量计算：

三频信标包括特高频、甚高频和低频三个异频载波，其频率分别为 n₁＝3,n₂＝8,f₀＝16.67MHz；首先，将三个频率载波两两差 分得时延和相位延迟方程，联立方程进行数学解算得到考虑电离层散射效应的绝对总电 子含量TEC_p和相对总电子含量TEC_φ；接着，根据观测信号的伪距和相位延迟数据求得 TEC_p与TEC_φ之间的差值：

$\Delta TEC=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}({\mathrm{TEC}}_{pi}-{\mathrm{TEC}}_{\phi i})$

其中，N为信号观测次数；

计算相位积分常数：

$k_0=\frac{40.3{\left(n_2^2-n_1^2\right)}^2}{{\mathrm{cf}}_0{n}_2^2\left(n_2^4+n_1^4\right)}\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\left({\mathrm{TEC}}_{pi}-{\mathrm{TEC}}_{\phi i}\right)$

其中，c＝3×10⁸m/s为电波的传播速度；

最后，将相位积分常数带入计算得到电离层总电子含量；

第二步，层析成像投影矩阵的计算：

首先根据待反演的电离层区域的经、纬度及高度位置信息，将其分别沿经、纬度和 高度进行等间隔的网格划分，使连续的电离层离散化为n个类方形小块即像素，并假设 每一个像素范围内的电子密度为同一值；然后，信号穿过电离层射向地面时，信号的传 播射线被各像素分割，根据卫星和地面接收机的精密定位，通过空间几何推算出m条射 线与经、纬面及高度球层的交点坐标；接着，根据每一像素的位置坐标和所求交点坐标 的关系，对各射线是否穿过该像素进行判断，并求出m条射线在n个像素中的截距d； 最后，分别以同一射线在n个像素中的截距作为行向量，以同一个像素截取m条射线的 截距作列向量，得到投影矩阵A＝[d_ij]_m×n；

第三步，基于共轭梯度法的电离层层析成像混合反演方法：

电离层的电子密度与电离层的总电子含量关系式为：

y＝Ax

其中，x为电子密度，A为投影矩阵，y为电离层总电子含量的测量值；由于投影矩阵 受观测数据有限的影响，常常为秩亏矩阵，无法直接对上式进行逆运算得到电子密度； 首先采用乘法代数迭代法进行电子密度的初步反演，利用国际参考电离层模型生成各个 像素的电子密度迭代初始值X⁽¹⁾；然后，根据乘法代数迭代公式，依次对每个像素的电 子密度进行m次的迭代修正，得到电离层电子密度的乘法代数重构结果；在此基础上， 又将电离层的电子密度反演不适定问题正则化为无约束最优化问题，以乘法代数重构结 果作为初值，采用共轭梯度法对电离层电子密度进行优化反演，考虑投影后等式左右两 端的误差最小化，得到正定二次目标函数：

f(x)＝x^TA^TAx-2y^TAx+y^Ty

最终，通过极小化目标函数，求电子密度反演的无约束最优化问题的最优解，得到 电子密度的最终反演结果。

所述第一步中电离层总电子含量的计算过程为：

在一阶电离层折射影响的基础上，考虑电离层散射效应，并将其影响项同时加入到 信号时延和相位延迟的计算中，通过提高数据信息来得到更高精度的反演结果，考虑电 离层散射效应的电离层时延如下：

$t\left(f\right)=t_0+\frac{40.3}{{\mathrm{cf}}^2}\underset{s}{\int }{N}_e\left(s\right)ds+\frac{K}{f^4}$

其中，t₀为自由空间与中性大气产生的时延，f为载波频率，c为电波的传播速度，N_e为 电子密度，K为电离层散射效应参数；通过差分处理得到时延和相位延迟表达式如下：

${\mathrm{\Delta t}}_{ij}=\frac{40.3}{c}\left(\frac{1}{f_i^2}-\frac{1}{f_j^2}\right)\underset{s}{\int }{N}_{e}ds+\left(\frac{1}{f_i^4}-\frac{1}{f_j^4}\right)K$

${\mathrm{\Delta \phi }}_{ij}=\frac{40.3}{c}\frac{\left(f_j^2-f_i^2\right)}{f_i^2{f}_j^2}\underset{s}{\int }{N}_{e}ds+\frac{f_j^4-f_i^4}{f_i^3{f}_j^4}K$

其中，△t_ij＝t_i-t_j，与三频信号对应；

对三频信号联立方程进行数学解算，得到考虑散射效应影响后，采用时延和相位延 迟求电离层绝对总电子含量的表达式：

${\mathrm{TEC}}_p=\frac{{\mathrm{cf}}_0^2{n}_1^2{n}_2^2}{40.3{\left(n_2^2-n_1^2\right)}^3\left(n_2^2+n_1^2\right)}\left[\left(n_2^8-n_1^4{n}_2^4\right){\mathrm{\Delta t}}_{13}-\left(n_2^8-n_1^8\right){\mathrm{\Delta t}}_{12}\right]$

$TEC=\frac{{\mathrm{cf}}_0{n}_2^2\left(n_2^4+n_1^4\right)}{40.3{\left(n_2^2-n_1^2\right)}^2}\left[\left(\frac{n_2^4}{n_2^4+n_1^4}{\mathrm{\Delta \phi }}_{13}^{\prime }-{\mathrm{\Delta \phi }}_{12}^{\prime }\right)+k_0\right]={\mathrm{TEC}}_{\phi }+\Delta TEC$

其中，△φ'₁₃,△φ'₁₂为△φ₁₃,△φ₁₂的小数部分，k₀为相位积分常数，TEC_p表示用时延数据求 得的绝对总电子含量，TEC_φ表示相对总电子含量，TEC表示用相位延迟数据求得的绝对 总电子含量。由于数据观测精度的影响，用相位延迟求信号路径上的电离层绝对总电子 含量精度较高，因此用TEC作为最终求得的绝对总电子含量结果。

所述第三步中基于共轭梯度法的电离层层析成像混合反演方法为：首先采用保证反 演结果非负性的乘法代数重构法，基于最大熵原理进行电子密度重构，接着采用存储量 小、收敛速度快、稳定性高的共轭梯度法进行再反演，两种方法相结合的混合反演方法 利用有限的观测信息，最大化使投影矩阵接近实际，针对观测数据有限情况下的电离层 电子密度反演问题，有效完成电子密度的高精度空间分布反演，具体如下：

首先，利用国际参考电离层模型生成各个像素的电子密度迭代初始值X⁽¹⁾；然后， 根据乘法代数迭代公式：

$X_j^{(k+1)}=X_j^{(k+1)}{\left(\frac{y_i}{<D_i^T,X^{\left(k\right)}>}\right)}^{\frac{{\lambda }_k{d}_{ij}}{\left|\right|D_i\left|\right|}}$

其中，为第k次迭代得到第j个电离层像素的电子密度，y_i为第i条信号射线路径上 的电离层绝对总电子含量，D_i为投影矩阵的第i行向量，X^(k)为第k次迭代得到的n个电 离层像素的电子密度向量，称为松弛因子；依次对每个像素的电子密度进行m 次的迭代修正，得到电离层电子密度的乘法代数重构结果；最后，以电子密度的乘法代 数重构结果作为初值，采用共轭梯度法依次按以下六步对电离层电子密度进行优化反 演：

第一步，取乘法迭代算法的结果作为初值x⁽¹⁾，设置精度要求E；

第二步，计算目标函数负梯度-f'(x^(k))＝2A^TAx^(k)-2A^Ty，其中，x^(k)为第k次迭代的 电子密度，A为投影矩阵，y为电离层总电子含量的测量值；

第三步，以minf(x^(k)+λ^(k)p^(k))为目标，采用一维线性搜索求λ^(k)∈(0,+∞)的最优解， 其中，λ^(k)为第k次迭代的搜索步长，p^(k)为第k次迭代的搜索方向，p⁽¹⁾＝-f'(x⁽¹⁾)；

第四步，计算电子密度x^(k+1)＝x^(k)+λ^(k)p^(k)；

第五步，判断||f'(x^(k+1))||<E是否满足，若满足，令fin＝x^(k+1)，计算结束，得到优化 结果，否则继续计算p^(k+1)＝-f'(x^(k+1))+a^(k+1)p^(k)；

第六步，判断f'(x^(k+1))^Tp^(k+1)≥0，令x⁽¹⁾＝x^(k+1)，k＝1，转第二步，否则令k＝k+1， 转第三步；

得到最优解作为电离层电子密度的最终反演结果。

本发明与现有技术相比的优点在于：

(1)相较于传统基于双频信标和三频信标的时延方程，前者仅考虑电离层一阶折 射影响，后者仅在绝对总电子密度计算中考虑散射效应，本发明利用含散射效应时延的 电离层时延方程，同时在绝对和相对总电子密度计算中考虑散射效应，通过增加数据信 息提高总电子含量的探测精度，以此为基础的电离层电子密度反演，精度相应提高。

(2)对电离层电子密度进行的层析成像反演，采用乘法代数重构法和共轭梯度法 相结合的混合迭代重构算法，乘法代数重构法基于最大熵原理进行电子密度重构，同时 具有保证反演结果非负性的优点，而共轭梯度法则克服了最速下降法收敛慢及牛顿法存 储和计算量大的问题，在解决高维的无约束优化问题中具有存储量小，收敛速度快、稳 定性高等优点。两者结合的混合反演充分利用有限的观测信息，最大化使投影矩阵接近 实际，相较于单一重构算法(代数迭代法、乘法代数迭代法等)，提高了反演结果的精 度和可靠性。

附图说明

图1为本发明一种基于共轭梯度法的电子密度层析成像混合反演方法的处理流程 图。

具体实施方式

如图1所示，本发明具体实现步骤如下：

首先，采用三频信标数据，考虑电离层散射效应，在绝对和相对总电子密度反演中 均加入散射效应时延项，通过三频差分获得电离层总电子含量；其次，选定电离层待反 演区，结合空间几何计算，得到层析成像的投影矩阵；最后，采用乘法代数重构算法和 共轭梯度法相结合的混合反演方法得到电离层电子密度分布。具体步骤如下：

第一步，基于三频信标观测的电离层总电子含量计算：

三频信标包括特高频、甚高频和低频三个异频载波，其频率分别为 n₁＝3,n₂＝8,f₀＝16.67MHz。在一阶电离层折射影响的基础 上，考虑电离层散射效应，得到电离层时延模型：

$t\left(f\right)=t_0+\frac{40.3}{{\mathrm{cf}}^2}\underset{s}{\int }{N}_e\left(s\right)ds+\frac{K}{f^4}$

其中，t₀为自由空间与中性大气产生的时延，f为载波频率，c为电波的传播速度，N_e为 电子密度，K为电离层散射效应参数。差分处理得到时延和相位延迟方程如下：

${\mathrm{\Delta t}}_{ij}=\frac{40.3}{c}\left(\frac{1}{f_i^2}-\frac{1}{f_j^2}\right)\underset{s}{\int }{N}_{e}ds+\left(\frac{1}{f_i^4}-\frac{1}{f_j^4}\right)K$

${\mathrm{\Delta \phi }}_{ij}=\frac{40.3}{c}\frac{(f_j^2-f_i^2)}{f_i^2{f}_j^2}\underset{s}{\int }{N}_{e}ds+\frac{f_j^4-f_i^4}{f_i^3{f}_j^4}K$

其中，Δt_ij＝t_i-t_j，联立方程进行数学解算得：

${\mathrm{TEC}}_p=\frac{{\mathrm{cf}}_0^2{n}_1^2{n}_2^2}{40.3{\left(n_2^2-n_1^2\right)}^3\left(n_2^2+n_1^2\right)}\left[\left(n_2^8-n_1^4{n}_2^4\right){\mathrm{\Delta t}}_{13}-\left(n_2^8-n_1^8\right){\mathrm{\Delta t}}_{12}\right]$

$TEC=\frac{{\mathrm{cf}}_0{n}_2^2\left(n_2^4+n_1^4\right)}{40.3{\left(n_2^2-n_1^2\right)}^2}\left[\left(\frac{n_2^4}{n_2^4+n_1^4}{\mathrm{\Delta \phi }}_{13}^{\prime }-{\mathrm{\Delta \phi }}_{12}^{\prime }\right)+k_0\right]={\mathrm{TEC}}_{\phi }+\Delta TEC$

其中，△φ'₁₃,△φ'₁₂为△φ₁₃,△φ₁₂的小数部分，k₀为相位积分常数，TEC_p表示用时延数据求 得的绝对总电子含量，TEC_φ表示相对总电子含量，TEC表示用相位延迟数据求得的绝对 总电子含量。由于数据观测精度的影响，用相位延迟求信号路径上的电离层绝对总电子 含量精度较高，因此用TEC作为最终求得的绝对总电子含量结果。根据观测信号的伪距 和相位延迟数据求得：

$\Delta TEC=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}({\mathrm{TEC}}_{pi}-{\mathrm{TEC}}_{\phi i})$

其中，N为信号观测次数；

由此可得相位积分常数：

$k_0=\frac{40.3{\left(n_2^2-n_1^2\right)}^2}{{\mathrm{cf}}_0{n}_2^2\left(n_2^4+n_1^4\right)}\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\left({\mathrm{TEC}}_{pi}-{\mathrm{TEC}}_{\phi i}\right)$

将相位积分常数带入TEC求解式中得到电离层总电子含量。

第二步，层析成像投影矩阵的计算：

首先根据待反演的电离层区域的经、纬度及高度位置信息，以地心为原点，y轴指 向本初子午线与赤道的交点，建立空间直角坐标系，得到电离层待反演区域的经、纬度 和高度面表达式。按照1°×1°×20km的规格对电离层进行等间隔的网格划分，使得连续的 电离层离散化为n个类方形小块(像素)，并假设每一个像素范围内的电子密度为同一 值Ne_j(j＝1,...,n)；然后，由于信号穿过电离层射向地面时，信号的传播射线被各像素分 割，根据卫星和地面接收机的精密定位，得到射线的表达式为：

$\frac{x-x_1}{x_2-x_1}=\frac{y-y_1}{y_2-y_1}=\frac{z-z_1}{z_2-z_1}$

其中，(x₁,y₁,z₁),(x₂,y₂,z₂)分别表示卫星和地面站在上述坐标系下的直角坐标。

联立方程可得m条射线在穿过电离层时与经面、纬面、高度球面的交点坐标；接着， 根据每一像素的位置坐标和所得交点坐标的关系，判断各射线是否穿过该像素，并求出 m条射线在n个像素中的截距d_ij(i＝1,...,m；j＝1,...,n)：

$d_{ij}=\sqrt{{({\mathrm{xl}}_2-{\mathrm{xl}}_1)}^2+{({\mathrm{yl}}_2-{\mathrm{yl}}_1)}^2+{({\mathrm{zl}}_2-{\mathrm{zl}}_1)}^2}$

其中，(xl₁,yl₁,zl₁),(xl₂,yl₂,zl₂)分别表示第i条射线在第j个像素中的两相邻交点；最后， 以同一射线在n个像素中的截距作为行向量，同一个像素截取m条射线的截距作列向量， 得到投影矩阵A＝[d_ij]_m×n。

第三步，基于共轭梯度法的电离层层析成像混合反演方法：

电离层的电子密度与总电子含量关系式为：

y＝Ax

其中，x为电子密度，A为投影矩阵，y为总电子含量测量值。

投影矩阵A和总电子含量测量值y已知，对上式进行逆运算可得到电子密度x。然 而投影矩阵受观测数据有限影响，常常为秩亏矩阵，无法求逆。首先采用乘法代数迭代 法进行电子密度的初步反演，利用国际参考电离层模型生成各个像素的电子密度迭代初 始值X⁽¹⁾；然后，根据乘法代数迭代公式，依次对每个像素的电子密度进行m次的迭代 修正，得到电离层电子密度的乘法代数重构结果。在此基础上，又将电离层的电子密度 反演问题处理为无约束最优化问题，以乘法代数重构结果作为初值，采用共轭梯度法对 电离层电子密度进行优化反演，考虑投影后等式左右两端的误差最小化，得到正定二次 目标函数：

f(x)＝x^TA^TAx-2y^TAx+y^Ty

通过极小化目标函数，求该无约束问题的最优解，具体步骤为：

第一步，取乘法迭代算法的结果作为初值x⁽¹⁾，设置精度要求E；

第二步，计算目标函数负梯度-f'(x^(k))＝2A^TAx^(k)-2A^Ty，其中，x^(k)为第k次迭代的 电子密度，A为投影矩阵，y为电离层总电子含量的测量值；

第三步，以minf(x^(k)+λ^(k)p^(k))为目标，采用一维线性搜索求λ^(k)∈(0,+∞)的最优解， 其中，λ^(k)为第k次迭代的搜索步长，p^(k)为第k次迭代的搜索方向，p⁽¹⁾＝-f'(x⁽¹⁾)；

第四步，计算电子密度x^(k+1)＝x^(k)+λ^(k)p^(k)；

第五步，判断||f'(x^(k+1))||<E是否满足，若满足，令fin＝x^(k+1)，计算结束，得到优化 结果，否则继续计算p^(k+1)＝-f'(x^(k+1))+a^(k+1)p^(k)；

第六步，判断f'(x^(k+1))^Tp^(k+1)≥0，令x⁽¹⁾＝x^(k+1)，k＝1，转第二步，否则令k＝k+1， 转第三步。

得到最优解作为电离层电子密度的最终反演结果。

其中，一维线性搜索求λ^(k)∈(0,+∞)的最优解采用标准的Armijo搜索，具体步骤如下：

·在ρ∈(0,1)中选取参数ρ的值；

·求搜索步长λ^(k)＝max{ρ^j|j＝1,2,...}，使其满足f(x^(k+1))≤f(x^(k))+δλ^(k)f'(x^(k))^Tp^(k)，其 中，参数δ∈(0,0.5)，f'(x)为目标函数梯度，x^(k)为第k次迭代的电子密度，p^(k)为第k次迭 代的搜索方向；。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。