CN105573963A - 一种电离层水平不均匀结构重构方法 - Google Patents

Abstract

本发明公开一种电离层水平不均匀结构重构方法，包括如下步骤：（一）建立区域电离层电子浓度网格，在返回散射探测方位上，按照地面距离和高度分别划分网格；（二）准备反演算法的输入，（三）利用频率范围的实测前沿数据进行反演，采用求解非线性问题的Newton-Kontorovich方法进行迭代求解，（四）接着利用频率进行反演，（五）经过n次反演，最终得到二维电子浓度剖面，即我们的最优反演结果。本发明所公开的电离层水平不均匀结构重构方法，克服了现有技术中的缺点，基于解空间约束的思想，提出了频段递增逐步逼近反演算法，对解空间做出了合理限制。

Description

一种电离层水平不均匀结构重构方法

技术领域

本发明涉及电离层研究及应用领域，尤其涉及一种电离层水平不均匀结构重构方法。

背景技术

高频返回散射探测作为重要的电离层探测手段，能够实现遥远区域电离层空间上的连续监测，探测获取的高频返回散射扫频电离图显示了探测频率-群路径-回波能量三者之间的关系。电离图包含了探测路径上的电离层状态信息，通过对其反演可以实时获取大面积范围的电离层参数。由于电离层的时间聚焦和球形聚焦等效应，电离图上有着较清晰、陡峭的前沿(也称最小群时延)，一般能够准确判读。返回散射前沿除了受电离层电子浓度分布影响之外，几乎不受任何其他因素(如天线波束、地面特性等)的影响。因此，返回散射前沿被广泛用来进行电离层的反演研究。

不少学者已经致力于利用返回散射前沿反演电离层参数的研究，目前反演方法大致可以归纳为如下几类：(1)基于迭代的拟合算法。首先是Rao(1974)利用返回散射电离图前沿任意三组(p′，f)层数据反演准抛物模型的三个参数，随后，Rao(1975)提出了一种利用返回散射电离图离散散射源回波描迹反演电离层参数和散射源地面距离的方法。这两种反演算法都假设电离层球形对称，不考虑地磁场的影响，并且都是反演单层准抛物模型参数。1979年DuBroff将这种反演算法推广到电离层非球形对称的情况，假设了一个简单的梯度电离层模型，反演包括准抛物模型的3个参数及其梯度共6个参数。(2)引入求解不适定问题的理论与方法以解决返回散射电离图反演不稳定性问题。Chuang等(1977)利用地球物理反演中发展起来的BG理论建立了斜向探测和返回散射探测反演模型，并用简单模型的垂测反演结果验证了该方法的有效性。Fridman等(1994)以返回散射发射站处垂直探测得到的电离层电子浓度剖面和返回散射前沿作为输入，利用求解非线性问题的Newton-Kontorovich方法和线性不适定问题的Tikhonov正则化方法，反演电离层电子浓度的二维分布，获得电离层的水平不均匀结构。遗传算法是近些年来发展起来的一种高效的非线性全局优化算法，也有学者利用遗传算法来进行返回散射电离图的反演(谢树果，2005)。宋君等(2011)使用模拟退火方法对返回散射电离图前沿进行了反演，但是该方法只适用于电离层较平稳时期。现有的反演方法一般采用全频段反演的方式，在电离层存在较复杂水平不均匀特性时不能很好地反映电离层电子浓度的变化趋势。

发明内容

本发明所要解决的技术问题就是提供一种可以如实反映电离层水平方向的不均匀结构的电离层水平不均匀结构重构方法。

本发明采用如下技术方案：

一种电离层水平不均匀结构重构方法，包括如下步骤：

(一)建立区域电离层电子浓度网格，在返回散射探测方位上，按照地面距离和高度分别划分网格；

(二)准备反演算法的输入，包括：

①实测的返回散射前沿G_实测(f)；

②利用发射站的垂测电离图反演得到返回散射发射站上空的电子浓度剖面N₀(h)；

③将实测前沿G_实测(f)对应的探测频率范围[f_a，f_b]划分为n个频段，每一频段的长度记为Δf_i(i＝1，2，…，n)，获取频率范围[f_a，f_a+Δf₁]，[f_a，f_a+Δf₁+Δf₂]，…，[f_a，f_a+Δf₁+Δf₂+…+Δf_n]＝[f_a，f_b]对应的实测前沿的测量误差先验值δ₁，δ₂，…，δ_n，一般有δ₁≤δ₂≤…≤δ_n；

(三)利用频率范围[f_a，f_a+Δf₁]的实测前沿数据G_实测([f_a，f_a+Δf₁])进行反演，采用求解非线性问题的Newton-Kontorovich方法进行迭代求解，具体步骤如下：

①初始二维电子浓度剖面认为是水平均匀的，即步骤(一)中的电子浓度网格每一地面距离上空的电子浓度剖面都等同于返回散射发射站处的电子浓度剖面N₀(h)，则有N(h，x)＝N₀(h)，在此电离层模型下，利用电离层短波数字三维射线追踪技术合成频率范围[f_a，f_a+Δf₁]的理论前沿G_理论([f_a，f_a+Δf₁])；

②使用符号||||表示均方根，定义理论前沿值与实测前沿值之间的均方根误差为：

其中M表示参与反演的前沿点个数；

如果下列条件满足

||G_理论([f_a，f_a+Δf₁])-G_实测([f_a，f_a+Δf₁])||≤δ₁(9)

则认为步骤(三)的①中的初始电子浓度剖面就是反演问题的最优解，否则进行下面的迭代过程；

③求解频率范围[f_a，f_a+Δf₁]的前沿数据对应的反演问题，上式中K(f，x)称为核函数，f表示工作频率，x表示地面距离，表示均匀电离层模型下返回散射前沿对应的地面散射点到发射站的地面距离，G₁(f)＝G(f)-G₀(f)，其中G(f)表示实测的返回散射前沿，G₀(f)表示均匀电离层模型下的返回散射前沿，u(x)是待求函数；根据Tikhonov正则化方法，求解反演问题等价于求解 ${\mathrm{\Φ}}_{\mathrm{\α}}\[u,{\tilde{G}}_1\]={\∫}_{f_1}^{f_2}{\[Au-{\tilde{G}}_1\left(f\right)\]}^{2}df+\mathrm{\α}\mathrm{\Ω}\[u\]$ 的极小值问题，后者又等价于求解下列由(10)式积分微分方程和(11)式边界条件组成的问题：

${\∫}_0^{x_0^m}\stackrel{\‾}{K}(\mathrm{\ξ},t)u\left(t\right)dt-\mathrm{\α}\{\frac{d}{d\mathrm{\ξ}}\[q\×\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}\]-u\left(\mathrm{\ξ}\right)\}=b\left(\mathrm{\ξ}\right)---\left(10\right)$

满足边界条件

${\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}|}_{\mathrm{\ξ}=0}={\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}|}_{\mathrm{\ξ}=x_0^m}=0---\left(11\right)$

其中

$\stackrel{\‾}{K}(\mathrm{\ξ},t)={\∫}_{f_1}^{f_2}K(f,\mathrm{\ξ})K(f,t)df$

$b\left(\mathrm{\ξ}\right)={\∫}_{f_1}^{f_2}K(f,\mathrm{\ξ})G_1\left(f\right)df$

G₁(f)＝G_实测(f)-G_理论(f)

这里，f₁＝f_a，f₂＝f_a+Δf₁，G_实测(f)＝G_实测([f_a，f_a+Δf₁])，G_理论(f)＝G_理论([f_a，f_a+Δf₁])；

上述问题中，u是我们待求的未知函数；K是已知核函数，按照 $K(f,x)=\frac{{\left(S_0^m\right)}^3}{2\sqrt{{\mathrm{\ρ}}^2-{\left(S_0^m\right)}^2}}B(x,x_0^m)-\frac{{\mathrm{\ρ}}^2}{2S_0^m}\frac{({\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0)}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{{\left(S_0^m\right)}^2}{2}{\∫}_x^{x_0^m}B(x,x^{\′})\frac{{\mathrm{\ρ}}^2}{{\widetilde{\mathrm{\ϵ}}}_0}(\frac{1}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{\∂{\widetilde{\mathrm{\ϵ}}}_0}{\∂z}-$ $\frac{2}{R_E}){\mathrm{dx}}^{\′}$ 描述的方法计算，其中，

$\begin{array}{c}B(x,x_0^m)=\frac{z_{0s}\left(x\right)z_{0x}\left(x_0^m\right)-z_{0x}\left(x\right)z_{0s}\left(x_0^m\right)}{{\left(S_0^m\right)}^2}\×(\frac{\∂({\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0)}{\∂z}-\frac{{\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{\∂{\widetilde{\mathrm{\ϵ}}}_0}{\∂z})\\ +\frac{{\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0}{{\widetilde{\mathrm{\ϵ}}}_0}\\ \×\[z_{0x}\left(x_0^m\right)\left(z_{0xx}\right(x\left)z_{0s}\right(x)-z_{0x}(x\left)z_{0sx}\right(x\left)\right)\\ -2z_{0xx}\left(x\right)z_{0x}\left(x\right)z_{0s}\left(x_0^m\right)\]\end{array}$

$\begin{array}{cccc}{z}_{0x}=\frac{\∂z_0}{\∂x}& z_{0xx}=\frac{{\∂}^2{z}_0}{\∂x^2}& z_{0s}=\frac{\∂z_0}{\∂S}& z_{0sx}=\frac{{\∂}^2{z}_0}{\∂S\∂x}\end{array}$

ρ＝1+h/R_E $\begin{array}{ccc}z={\∫}_0^h{\mathrm{dh}}^{\′}/\mathrm{\ρ}& \widetilde{\mathrm{\ϵ}}={\mathrm{\ρ}}^2\mathrm{\ϵ}& \mathrm{\ϵ}=1-\frac{f_P^2}{f^2}\end{array}$

上式中，f_P表示等离子体频率；S＝sinβ，β表示射线的出射角，即射线传播方向与垂直方向的夹角，有S^m＝sinβ^m，其中β^m表示返回散射前沿对应的射线出射角；R_E表示地球半径；下标“0”表示均匀电离层模型下的参量，参数q一般需要满足条件q≥L²，其中L表示预期的电子浓度水平变化尺度典型值，考虑到中纬地区电离层水平不均匀性的典型尺度通常为上千公里量级或者更大，因此这里取值q＝10⁶km²，由于我们引进了电子浓度网格，可以将上述问题表达式离散化，采用数值方法进行求解，离散化形式如下：

$-\frac{\mathrm{\α}}{d^2}({\mathrm{qu}}_{k-1}+{\mathrm{qu}}_{k+1}-2{\mathrm{qu}}_k)+{\mathrm{\αu}}_k+{\mathrm{\Σ}}_{r=1}^v\left({\stackrel{\‾}{K}}_{k,r}{u}_r{h}_r\right)=b_k,k=1,...,v---\left(12\right)$

满足边界条件：

u₀＝u₁，u_v+1＝u_v(13)

其中

${\stackrel{\‾}{K}}_{k,r}=\underset{j=1}{\overset{w}{\Σ}}\left(K_{j,k}{K}_{j,r}s\right)$

$b_k=\underset{j=1}{\overset{w}{\Σ}}\left(K_{j,k}{G}_{1j}s\right)$

假设网格是均匀划分的，v代表地面距离网格个数，w代表探测频率个数，d表示地面距离网格的宽度，s表示探测频率步进；

固定正则化参数α的值，利用Cholesky分解法求解上述v个方程可以得到v个未知数的解 $u_{{\mathrm{\α}}_k}(k=1,...,v),$ 这是u_α(x)的离散化形式；所有满足 ${\∫}_{f_1}^{f_2}{\[{\tilde{G}}_1\left(f\right)-Au\]}^{2}df\≤{\mathrm{\δ}}^2(f_2-$ $f_1)$ 的α中的最大值称为最优正则化参数α_d，它对应的解称为正则化解，也就是我们要求的解；

其中表示实测前沿与利用射线追踪合成的理论前沿的差值，即f₁和f₂分别表示起始工作频率和终止工作频率，δ称为前沿测量误差的上界，将代入到(1)式，便得到反演后的电子浓度剖面N(h，x)；

④在新的电子浓度剖面N(h，x)下，运用射线追踪技术合成频率范围[f_a，f_a+Δf₁]的理论前沿G_理论([f_a，f_a+Δf₁])，将其代入(9)式，如果满足条件，则N(h，x)就是反演问题的最优解，记为N¹(h，x)，反演过程结束，否则，转入步骤(三)的③；

(四)接着利用频率范围[f_a，f_a+Δf₁+Δf₂]的实测前沿数据G_实测([f_a，f_a+Δf₁+Δf₂])进行反演，其反演方法与步骤(三)相同，只不过需将步骤(三)的①中涉及的初始二维电子浓度剖面取值成N¹(h，x)，即上一频率范围对应的前沿数据反演得到的电子浓度剖面，记该频率范围的反演结果为N²(h，x)；

(五)将频率范围按照步骤(二)的③划分的方法逐步增加，然后按照步骤(四)描述的方法完成反演，经过n次反演，最终得到二维电子浓度剖面Nⁿ(h，x)，即我们的最优反演结果。

本发明的有益效果在于：

本发明所公开的电离层水平不均匀结构重构方法，克服了现有技术中的缺点，基于解空间约束的思想，提出了频段递增逐步逼近反演算法，对解空间做出了合理限制。为了保证解的唯一性和稳定性，又引入了求解非线性问题的Newton-Kontorovich方法和求解线性不适定问题的Tikhonov正则化方法。本发明建立的算法反演精度较高，反演算法稳健，对返回散射前沿判读误差不敏感，能够准确提供返回散射探测方向上地面距离0-2000km的二维电子浓度分布，如实反映电离层水平方向的不均匀结构。

附图说明

图1为本发明所公开的二维电子浓度网格划分示意图；

图2为使用本发明所公开的方法和Fridman等(1994)算法对模拟数据反演结果的比较；

图3为使用本发明所公开的方法和Fridman等(1994)算法对实测数据反演结果的比较。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

反演基本原理：

本发明提出了频段递增逐步逼近反演算法，即将返回散射探测频率划分为若干个频段，每一次反演使用的频段都在上一次反演使用的频段基础上增加一个频段，并且反演的初始电子浓度剖面取值为上一次的反演结果，对解空间做出限制。这样的处理既可以避免确定不同频率区间反演结果的分界点以及随后需要考虑的连续性和光滑性问题，又可以很好地利用先验信息对解做出约束，逐步逼近真实解，提高反演精度。同时针对返回散射反演的非线性问题和不适定问题，又引入了Newton-Kontorovich方法和Tikhonov正则化方法，以保证解的稳定性和唯一性。

电离层模型：

本发明采用垂直剖面形状不变的二维电离层模型，其数学表达式为：

N(h，x)＝N₀(h)[1+u(x)](1)

其中h表示高度，x表示地面距离。该模型对实际电离层做了简化近似，适合描述突变的、不规则的或者水平变化不单一的电离层结构。模型假设电离层电子浓度分布的垂直剖面形状沿水平方向近似不变，且最大电子浓度高度沿水平方向不变，只是最大电子浓度值在水平方向按(1)式变化。式中，N₀(h)表示初始的电子浓度剖面，是已知的；u(x)是待求函数，它表征了临界频率(电子浓度)随地面距离的变化。一旦求解出u(x)，带入(1)式便可求得二维电子浓度分布N(h，x)。

反演方程：

利用返回散射前沿反演得到二维电子浓度剖面需要求解下列第一类Fredholm积分方程

${\∫}_0^{x_0^m}K(f,x)u\left(x\right)dx=G_1\left(f\right)---\left(2\right)$

上式中，K(f，x)称为核函数，f表示工作频率，x表示地面距离，表示均匀电离层模型下返回散射前沿对应的地面散射点到发射站的地面距离，G₁(f)＝G(f)-G₀(f)，其中G(f)表示实测的返回散射前沿，G₀(f)表示均匀电离层模型下的返回散射前沿，u(x)是待求函数。

本文中的核函数采用下列计算公式

$\begin{array}{c}K(f,x)=\frac{{\left(S_0^m\right)}^3}{2\sqrt{{\mathrm{\ρ}}^2-{\left(S_0^m\right)}^2}}B(x,x_0^m)-\frac{{\mathrm{\ρ}}^2}{2S_0^m}\frac{({\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0)}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{{\left(S_0^m\right)}^2}{2}{\∫}_x^{x_0^m}B(x,x^{\′})\frac{{\mathrm{\ρ}}^2}{{\widetilde{\mathrm{\ϵ}}}_0}(\frac{1}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{\∂{\widetilde{\mathrm{\ϵ}}}_0}{\∂z}-\\ \frac{2}{R_E}){\mathrm{dx}}^{\′}---(3)\end{array}$

其中，

$B(x,x_0^m)=\frac{z_{0s}\left(x\right)z_{0x}\left(x_0^m\right)-z_{0x}\left(x\right)z_{0s}\left(x_0^m\right)}{{\left(S_0^m\right)}^2}\×(\frac{\∂({\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0)}{\∂z}-\frac{{\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{\∂{\widetilde{\mathrm{\ϵ}}}_0}{\∂z})$

$\begin{array}{c}+\frac{{\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0}{{\widetilde{\mathrm{\ϵ}}}_0}\\ \×\[z_{0x}\left(x_0^m\right)\left(z_{0xx}\right(x\left)z_{0s}\right(x)-z_{0x}(x\left)z_{0sx}\right(x\left)\right)\\ -2z_{0xx}\left(x\right)z_{0x}\left(x\right)z_{0s}\left(x_0^m\right)\]\end{array}$

$\begin{array}{cccc}{z}_{0x}=\frac{\∂z_0}{\∂x}& z_{0xx}=\frac{{\∂}^2{z}_0}{\∂x^2}& z_{0s}=\frac{\∂z_0}{\∂S}& z_{0sx}=\frac{{\∂}^2{z}_0}{\∂S\∂x}\end{array}$

ρ＝1+h/R_E $\begin{array}{ccc}z={\∫}_0^h{\mathrm{dh}}^{\′}/\mathrm{\ρ}& \widetilde{\mathrm{\ϵ}}={\mathrm{\ρ}}^2\mathrm{\ϵ}& \mathrm{\ϵ}=1-\frac{f_P^2}{f^2}\end{array}$

上式中，f_P表示等离子体频率；S＝sinβ，β表示射线的出射角，即射线传播方向与垂直方向的夹角，有S^m＝sinβ^m，其中β^m表示返回散射前沿对应的射线出射角；R_E表示地球半径；下标“0”表示均匀电离层模型下的参量。

由于返回散射电离图反演属于不适定问题，这里采用Tikhonov正则化方法进行求解，即利用光滑性对解做出限制，以保证解是最优且唯一。

根据Tikhonov正则化理论，满足(4)式的称为反演问题(2)式的一个近似解，

${\∫}_{f_1}^{f_2}{\[{\tilde{G}}_1\left(f\right)-Au\]}^{2}df\≤{\mathrm{\δ}}^2(f_2-f_1)---\left(4\right)$

其中 $Au={\∫}_0^{x_0^m}K(f,x)u\left(x\right)dx---\left(5\right)$

表示实测前沿与利用射线追踪合成的理论前沿的差值，即 f₁和f₂分别表示起始工作频率和终止工作频率，δ称为前沿测量误差的上界。

认为其中最光滑的是最佳的近似，亦即正则化解。使用下列形式的泛函作为光滑度的度量，又称稳定化泛函。

$\mathrm{\Ω}\[u\]={\∫}_0^{x_f}f\[u^2\left(x\right)+q{\left(\frac{du\left(x\right)}{dx}\right)}^2\]dx---\left(6\right)$

这里，q是个常数，q≥0。于是，在(4)式的约束下，使(6)式最小化的解可以认为是反演问题(2)式的最优解。应用拉格朗日乘子法，则反演问题变成求使下列Tikhonov正则化泛函

${\mathrm{\Φ}}_{\mathrm{\α}}\[u,{\tilde{G}}_1\]={\∫}_{f_1}^{f_2}{\[Au-{\tilde{G}}_1\left(f\right)\]}^{2}df+\mathrm{\α}\mathrm{\Ω}\[u\]---\left(7\right)$

取得极小值的u_α。这里α称为正则化参数，α＞0。α越大，稳定化泛函Ω[u]在(7)式中的权重越大，得到的解u_α也越光滑。但如果α取值太大，则会使解过分光滑而偏离真实解。

实施例1，本实施例公开了一种电离层水平不均匀结构重构方法，包括如下步骤：

(1)首先建立我们感兴趣区域的电离层电子浓度网格。在返回散射探测方位上，按照地面距离和高度分别划分网格，如图1所示。

(2)准备反演算法的输入，包括：

①实测的返回散射前沿G_实测(f)；

②返回散射发射站上空的电子浓度剖面N₀(h)；

③将实测前沿G_实测(f)对应的探测频率范围[f_a，f_b]划分为n个频段，每一频段的长度记为Δf_i(i＝1，2，…，n)，获取频率范围[f_a，f_a+Δf₁]，[f_a，f_a+Δf₁+Δf₂]，…，[f_a，f_a+Δf₁+Δf₂+…+Δf_n]＝[f_a，f_b]对应的实测前沿的测量误差先验值δ₁，δ₂，…，δ_n，一般有δ₁≤δ₂≤…≤δ_n。

其中②利用发射站的垂测电离图反演得到。

(3)首先利用频率范围[f_a，f_a+Δf₁]的实测前沿数据G_实测([f_a，f_a+Δf₁])进行反演。由于反演问题(2)式一般是非线性的，这里采用求解非线性问题的Newton-Kontorovich方法进行迭代求解。具体步骤如下：

(I)初始二维电子浓度剖面认为是水平均匀的，即第(1)步中的电子浓度网格每一地面距离上空的电子浓度剖面都等同于返回散射发射站处的电子浓度剖面N₀(h)，则有N(h，x)＝N₀(h)。在此电离层模型下，利用电离层短波数字三维射线追踪技术合成频率范围[f_a，f_a+Δf₁]的理论前沿G_理论([f_a，f_a+Δf₁])。

(II)使用符号||||表示均方根，定义理论前沿值与实测前沿值之间的均方根误差为

其中M表示参与反演的前沿点个数。

如果下列条件满足

||G_理论([f_a，f_a+Δf₁])-G_实测([f_a，f_a+Δf₁])||≤δ₁(15)

则认为第(3)步的步骤(I)中的初始电子浓度剖面就是反演问题的最优解。否则进行下面的迭代过程。

(III)求解频率范围[f_a，f_a+Δf₁]的前沿数据对应(2)式的反演问题。

根据Tikhonov正则化方法，求解反演问题(2)式等价于求解(7)式的极小值问题，后者又等价于求解下列由(16)式积分微分方程和(17)式边界条件组成的问题：

${\∫}_0^{x_0^m}\stackrel{\‾}{K}(\mathrm{\ξ},t)u\left(t\right)dt-\mathrm{\α}\{\frac{d}{d\mathrm{\ξ}}\[q\×\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}\]-u\left(\mathrm{\ξ}\right)\}=b\left(\mathrm{\ξ}\right)---\left(16\right)$

满足边界条件

${\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}|}_{\mathrm{\ξ}=0}={\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}|}_{\mathrm{\ξ}=x_0^m}=0---\left(17\right)$

其中

$\stackrel{\‾}{K}(\mathrm{\ξ},t)={\∫}_{f_1}^{f_2}K(f,\mathrm{\ξ})K(f,t)df$

$b\left(\mathrm{\ξ}\right)={\∫}_{f_1}^{f_2}K(f,\mathrm{\ξ})G_1\left(f\right)df$

G₁(f)＝G_实测(f)-G_理论(f)

这里，f₁＝f_a，f₂＝f_a+Δf₁，G_实测(f)＝G_实测([f_a，f_a+Δf₁])，G_理论(f)＝G_理论([f_a，f_a+Δf₁])。

上述问题中，u是我们待求的未知函数；K是已知核函数，按照(3)式描述的方法计算；参数q一般需要满足条件q≥L²，其中L表示预期的电子浓度水平变化尺度典型值，考虑到中纬地区电离层水平不均匀性的典型尺度通常为上千公里量级或者更大，因此这里取值q＝10⁶km²。需要说明的是，虽然参数q的选择很重要，但只要它在一定范围内取值，反演结果对它是不敏感的。

由于我们引进了电子浓度网格，可以将上述问题表达式离散化，采用数值方法进行求解，离散化形式如下

$-\frac{\mathrm{\α}}{d^2}({\mathrm{qu}}_{k-1}+{\mathrm{qu}}_{k+1}-2{\mathrm{qu}}_k)+{\mathrm{\αu}}_k+{\mathrm{\Σ}}_{r=1}^v\left({\stackrel{\‾}{K}}_{k,r}{u}_r{h}_r\right)=b_k,k=1,...,v---\left(18\right)$

满足边界条件

u₀＝u₁，u_v+1＝u_v(19)

其中

${\stackrel{\‾}{K}}_{k,r}=\underset{j=1}{\overset{w}{\Σ}}\left(K_{j,k}{K}_{j,r}s\right)$

假设网格是均匀划分的，v代表地面距离网格个数，w代表探测频率个数，d表示地面距离网格的宽度，s表示探测频率步进。

固定正则化参数α的值，利用Cholesky分解法求解上述v个方程可以得到v个未知数的解这是u_α(x)的离散化形式。所有满足(4)式的α中的最大值称为最优正则化参数α_d，它对应的解称为正则化解，也就是我们要求的解。将代入到(1)式，便得到反演后的电子浓度剖面N(h，x)。

(IV)在新的电子浓度剖面N(h，x)下，运用射线追踪技术合成频率范围[f_a，f_a+Δf₁]的理论前沿G_理论([f_a，f_a+Δf₁])，将其代入(15)式，如果满足条件，则N(h，x)就是反演问题的最优解，记为N¹(h，x)。反演过程结束。否则，转入第(3)步的步骤(III)。

(4)接着利用频率范围[f_a，f_a+Δf₁+Δf₂]的实测前沿数据G_实测([f_a，f_a+Δf₁+Δf₂])进行反演，其反演方法与第(3)步相同，只不过需将第(3)步的步骤(I)中涉及的初始二维电子浓度剖面取值成N¹(h，x)，即上一频率范围对应的前沿数据反演得到的电子浓度剖面。记该频率范围的反演结果为N²(h，x)。

5)将频率范围按照第(2)步的步骤③划分的方法逐步增加，然后按照第(4)步描述的方法完成反演，经过n次反演，最终得到二维电子浓度剖面Nⁿ(h，x)，即我们的最优反演结果。至此完成本发明的全部反演工作。

利用模拟数据和实测数据分别对本发明建立的算法进行了验证，并与Fridman等1994年提出的反演方法进行了对比，反演结果分别见图2和图3。可以看出，无论是存在严重电离层扰动的模拟数据还是水平方向上具有较大梯度的实测数据，本发明的反演精度都要高于Fridman等(1994)算法，说明本发明的反演效果是非常有效的。

Claims (1)

1.一种电离层水平不均匀结构重构方法，其特征在于，包括如下步骤：

(一)建立区域电离层电子浓度网格，在返回散射探测方位上，按照地面距离和高度分别划分网格；

(二)准备反演算法的输入，包括：

①实测的返回散射前沿G_实测(f)；

②利用发射站的垂测电离图反演得到返回散射发射站上空的电子浓度剖面N₀(h)；

③将实测前沿G_实测(f)对应的探测频率范围[f_a,f_b]划分为n个频段，每一频段的长度记

为Δf_i(i＝1,2,…,n)，获取频率范围[f_a,f_a+Δf₁]，[f_a,f_a+Δf₁+Δf₂]，…，[f_a,f_a+Δf₁+Δf₂+…+Δf_n]＝[f_a,f_b]对应的实测前沿的测量误差先验值δ₁，δ₂，…，δ_n，一般有δ₁≤δ₂≤…≤δ_n；

(三)利用频率范围[f_a,f_a+Δf₁]的实测前沿数据G_实测([f_a,f_a+Δf₁])进行反演，采用求解非线性问题的Newton-Kontorovich方法进行迭代求解，具体步骤如下：

①初始二维电子浓度剖面认为是水平均匀的，即步骤(一)中的电子浓度网格每一地面距离上空的电子浓度剖面都等同于返回散射发射站处的电子浓度剖面N₀(h)，则有N(h,x)＝N₀(h)，在此电离层模型下，利用电离层短波数字三维射线追踪技术合成频率范围[f_a,f_a+Δf₁]的理论前沿G_理论([f_a,f_a+Δf₁])；

②使用符号||||表示均方根，定义理论前沿值与实测前沿值之间的均方根误差为：

其中M表示参与反演的前沿点个数；

如果下列条件满足

||G_理论([f_a,f_a+Δf₁])-G_实测([f_a,f_a+Δf₁])||≤δ₁(9)

则认为步骤(三)的①中的初始电子浓度剖面就是反演问题的最优解，否则进行下面的迭代过程；

③求解频率范围[f_a,f_a+Δf₁]的前沿数据对应的反演问题，上式中K(f,x)称为核函数，f表示工作频率，x表示地面距离，表示均匀电离层模型下返回散射前沿对应的地面散射点到发射站的地面距离，G₁(f)＝G(f)-G₀(f)，其中G(f)表示实测的返回散射前沿，G₀(f)表示均匀电离层模型下的返回散射前沿，u(x)是待求函数；根据Tikhonov正则化方法，求解反演问题等价于求解 ${\mathrm{\Φ}}_{\mathrm{\α}}\[u,{\tilde{G}}_1\]={\∫}_{f_1}^{f_2}{\[Au-{\tilde{G}}_1\left(f\right)\]}^{2}df+\mathrm{\α}\mathrm{\Ω}\[u\]$ 的极小值问题，后者又等价于求解下列由(10)式积分微分方程和(11)式边界条件组成的问题：

${\∫}_0^{x_0^m}\stackrel{\‾}{K}(\mathrm{\ξ},t)u\left(t\right)dt-\mathrm{\α}\{\frac{d}{d\mathrm{\ξ}}\[q\×\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}\]-u\left(\mathrm{\ξ}\right)\}=b\left(\mathrm{\ξ}\right)---\left(10\right)$

满足边界条件

$\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}{|}_{\mathrm{\ξ}=0}=\frac{du\left(\mathrm{\ξ}\right)}{d\mathrm{\ξ}}{|}_{\mathrm{\ξ}=x_0^m}=0---\left(11\right)$

其中

$\stackrel{\‾}{K}(\mathrm{\ξ},t)={\∫}_{f_1}^{f_2}K(f,\mathrm{\ξ})K(f,t)df$

$b\left(\mathrm{\ξ}\right)={\∫}_{f_1}^{f_2}K(f,\mathrm{\ξ})G_1\left(f\right)df$

G₁(f)＝G_实测(f)-G_理论(f)

这里，f₁＝f_a，f₂＝f_a+Δf₁，G_实测(F)＝G_实测([f_a,f_a+Δf₁])，G_理论(f)＝G_理论([f_a,f_a+Δf₁])；

上述问题中，u是我们待求的未知函数；K是已知核函数，按照

$K\left(f,x\right)=\frac{{\left(S_0^m\right)}^3}{2\sqrt{{\mathrm{\ρ}}^2-{\left(S_0^m\right)}^2}}B\left(x,x_0^m\right)-\frac{{\mathrm{\ρ}}^2}{2S_0^m}\frac{\left({\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0\right)}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{{\left(S_0^m\right)}^2}{2}{\∫}_x^{x_0^m}B\left(x,x^{\′}\right)\frac{{\mathrm{\ρ}}^2}{{\widetilde{\mathrm{\ϵ}}}_0}(\frac{1}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{\∂{\widetilde{\mathrm{\ϵ}}}_0}{\∂z}-$ $\frac{2}{R_E}){\mathrm{dx}}^{\′}$ 描述的方法计算，其中，

$\begin{array}{c}B\left(x,x_0^m\right)=\frac{z_{0s}\left(x\right)z_{0x}\left(x_0^m\right)-z_{0x}\left(x\right)z_{0s}\left(x_0^m\right)}{{\left(S_0^m\right)}^2}\×\left(\frac{\∂\left({\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0\right)}{\∂z}-\frac{{\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0}{{\widetilde{\mathrm{\ϵ}}}_0}\frac{\∂{\widetilde{\mathrm{\ϵ}}}_0}{\∂z}\right)\\ +\frac{{\mathrm{\ρ}}^2-{\widetilde{\mathrm{\ϵ}}}_0}{{\widetilde{\mathrm{\ϵ}}}_0}\\ \×\[z_{0x}\left(x_0^m\right)\left(z_{0xx}\left(x\right)z_{0s}\left(x\right)-z_{0x}\left(x\right)z_{0sx}\left(x\right)\right)\\ -2z_{0xx}\left(x\right)z_{0x}\left(x\right)z_{0s}\left(x_0^m\right)\]\end{array}$

$z_{0x}=\frac{\∂z_0}{\∂x}$ $z_{0xx}=\frac{{\∂}^2{z}_0}{\∂x^2}$ $z_{0s}=\frac{\∂z_0}{\∂S}$ $z_{0sx}=\frac{{\∂}^2{z}_0}{\∂S\∂x}$

ρ＝1+h/R_E $z={\∫}_0^h{\mathrm{dh}}^{\′}/\mathrm{\ρ}$ $\widetilde{\mathrm{\ϵ}}={\mathrm{\ρ}}^2\mathrm{\ϵ}$ $\mathrm{\ϵ}=1-\frac{f_P^2}{f^2}$

上式中，f_P表示等离子体频率；S＝sinβ，β表示射线的出射角，即射线传播方向与垂直方向的夹角，有S^m＝sinβ^m，其中β^m表示返回散射前沿对应的射线出射角；R_E表示地球半径；下标“0”表示均匀电离层模型下的参量，参数q一般需要满足条件q≥L²，其中L表示预期的电子浓度水平变化尺度典型值，考虑到中纬地区电离层水平不均匀性的典型尺度通常为上千公里量级或者更大，因此这里取值q＝10⁶km²，由于我们引进了电子浓度网格，可以将上述问题表达式离散化，采用数值方法进行求解，离散化形式如下：

$-\frac{\mathrm{\α}}{d^2}({\mathrm{qu}}_{k-1}+{\mathrm{qu}}_{k+1}-2{\mathrm{qu}}_k)+{\mathrm{\αu}}_k+{\mathrm{\Σ}}_{r=1}^v\left({\stackrel{\‾}{K}}_{k,r}{u}_r{h}_r\right)=b_k,k=1,...,v---\left(12\right)$

满足边界条件：

u₀＝u₁，u_v+1＝u_v(13)

其中

${\stackrel{\‾}{K}}_{k,r}=\underset{j=1}{\overset{w}{\Σ}}\left(K_{j,k}{K}_{j,r}s\right)$

$b_k=\underset{j=1}{\overset{w}{\Σ}}\left(K_{j,k}{G}_{1_j}s\right)$

假设网格是均匀划分的，v代表地面距离网格个数，w代表探测频率个数，d表示地面距离网格的宽度，S表示探测频率步进；

固定正则化参数α的值，利用Cholesky分解法求解上述v个方程可以得到v个未知数的解 $u_{{\mathrm{\α}}_k},(k=1,...,v),$ 这是u_α(x)的离散化形式；所有满足 ${\∫}_{f_1}^{f_2}{\[{\tilde{G}}_1\left(f\right)-Au\]}^{2}df\≤{\mathrm{\δ}}^2(f_2-$ $f_1)$ 的α中的最大值称为最优正则化参数α_d，它对应的解称为正则化解，也就是我们要求的解；

其中 表示实测前沿与利用射线追踪合成的理论前沿的差值，即f₁和f₂分别表示起始工作频率和终止工作频率，δ称为前沿测量误差的上界，将代入到(1)式，便得到反演后的电子浓度剖面N(h,x)；

④在新的电子浓度剖面N(h,x)下，运用射线追踪技术合成频率范围[f_a,f_a+Δf₁]的理论前沿G_理论([f_a,f_a+Δf₁])，将其代入(9)式，如果满足条件，则N(h,x)就是反演问题的最优解，记为N¹(h,x)，反演过程结束，否则，转入步骤(三)的③；

(四)接着利用频率范围[f_a,f_a+Δf₁+Δf₂]的实测前沿数据G_实测([f_a,f_a+Δf₁+Δf₂])进行反演，其反演方法与步骤(三)相同，只不过需将步骤(三)的①中涉及的初始二维电子浓度剖面取值成N¹(h,x)，即上一频率范围对应的前沿数据反演得到的电子浓度剖面，记该频率范围的反演结果为N²(h,x)；

(五)将频率范围按照步骤(二)的③划分的方法逐步增加，然后按照步骤(四)描述的方法完成反演，经过n次反演，最终得到二维电子浓度剖面Nⁿ(h,x)，即我们的最优反演结果。