CN105184039A - 一种电离层垂直剖面建模及参数反演方法 - Google Patents

Abstract

本发明公开了一种电离层垂直剖面建模及参数反演方法，所述方法包括以下步骤，步骤1：将电离层建模为包含E层、谷层、层、层的四层模型；步骤2：基于建立的电离层模型，推导各层回波虚高的计算公式；步骤3：利用实测电离图数据，结合各层回波虚高计算结果，进行各层电离层参数的反演。本发明所公开的电离层垂直剖面建模及参数反演方法，基于移位切比雪夫多项式模型的约束优化层参数、层参数的垂测电离图反演方法，即得到谷参数后，选取层较高区域回波描迹数据，在保证剖面连续光滑的约束条件下，计算层剖面多项式系数，该方法适用于层未充分发育的情况，可以有效提高该种情况下反演精度和稳定性。

Description

一种电离层垂直剖面建模及参数反演方法

技术领域

本发明涉及电离层研究及应用领域，尤其涉及一种电离层垂直剖面建模及参数反演方法。

背景技术

电离层垂直探测(简称垂测)能够获得反映电离层虚高与频率关系的垂测电离图，所谓虚高，其实质是电磁波传播时间的度量，电磁波从垂测发射机发出，经电离层反射后到达垂测接收机所经历的时间乘以光速就称为虚高。垂测获得的虚高并不是电磁波在电离层中的真实反射高度，真实反射高度获得需要对垂测电离图进行反演，即利用垂测电离图频率-虚高描迹反演电离层剖面(电离层高度与等离子体频率或电子浓度的对应关系)。垂测电离图的反演对研究电离层结构和电离层波传播问题具有重要意义，一直以来受到十分广泛的重视，当然，反演也具有相当大的难度。

目前，应用较为普遍的垂测电离图反演方法是基于模式法思想发展的电离层参数反演方法，其中，黄雪钦等公开了一种基于移位切比雪夫多项式模型反演电离层剖面的方法，该方法中，将F层建模为移位切比雪夫多项式，求解满足计算虚高与实测虚高在最小二乘意义上最佳吻合的多项式系数，从而确定电离层剖面，后来，针对回波描迹中明显有F₁层描迹的情况，做了模型上的一些改进，即将F₁层和F₂层分别用不同的移位切比雪夫多项式表示，但该方法更适用于F₁层充分发育的情况，对于F₁层未充分发育的情况，基于反演结果合成的垂测电离图与实测电离图在F₁层临频附近差别较大，另外，未考虑层与层连接处的光滑性问题也是该方法的不足之处。

发明内容

本发明所要解决的技术问题就是提供一种电离层垂直剖面建模及参数反演方法。

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

一种电离层垂直剖面建模及参数反演方法，其改进之处在于，所述方法包括以下步骤：

步骤1：将电离层建模为包含E层、谷层、F₁层、F₂层的四层模型，其中，E层和谷层剖面用抛物模型表示，F₁层和F₂层剖面用移位切比雪夫多项式模型表示；

步骤2：基于建立的电离层模型，推导各层回波虚高的计算公式；

步骤3：利用实测电离图数据，结合各层回波虚高计算结果，进行各层电离层参数的反演。

进一步的，所述步骤1具体包括：

步骤11：电离层电子浓度剖面具体形式如下式所示：

$\left\{\begin{array}{c}{f}_{\mathrm{NE}}^2=f_{\mathrm{CE}}^2[1-{\left(\frac{h-h_{\mathrm{mE}}}{y_{\mathrm{mE}}}\right)}^2]h_{\mathrm{bE}}\≤h\≤h_{\mathrm{mE}},h_{\mathrm{mE}}\≤h\≤h_1\\ f_{\mathrm{NV}}^2=f_{\mathrm{CV}}^2[1+{\left(\frac{h-h_{\mathrm{mV}}}{y_{\mathrm{mV}}}\right)}^2]h_1\≤h\≤h_2\\ h=A_{I+1}+g^{1/2}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{T}_i\left(g\right)h_2\≤h\≤h_{\mathrm{mF}1}\\ h=C_{N+1}+l^{1/2}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i{T}_i\left(l\right)h_{\mathrm{mF}1}\≤h\≤h_{\mathrm{mF}2}\end{array}\right.$

上式中各符号的具体含义如下：

E层：f_NE表示E层等离子体频率；f_CE表示E层临频；h_mE表示E层峰高；y_mE表示E层半厚；h_bE＝h_mE-y_mE表示E层底高；

谷层：f_NV表示谷层等离子体频率；f_CV表示谷层最小等离子体频率；h_mV表示谷层等离子体频率为f_CV时对应的电离层高度；y_mV表示谷层半厚；h₂＝h_mE+W，W定义为谷层宽度；

F₁层：T_i(g)为移位切比雪夫多项式，具有下式所示形式：

$\left\{\begin{array}{c}{T}_i\left(g\right)=2(2g-1)T_{i-1}\left(g\right)-T_{i-2}\left(g\right)\\ T_0\left(g\right)=1,T_1\left(g\right)=2g-1\end{array}\right.,g=\frac{\ln (f_{\mathrm{NF}1}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))\ln (f_{\mathrm{CF}1}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))}{\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))}\≤g\≤1$

f_NF1表示F₁层等离子体频率；f_CF1表示F₁层临频；A_i(i＝0～I+1)为移位切比雪夫多项式系数，

且： $A_{I+1}=h_2-\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i,$

h_mF1表示F₁层峰高，且： $h_{\mathrm{mF}1}={[A_{I+1}+g^{1/2}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{T}_i\left(g\right)]|}_{g=\frac{\ln (f_{\mathrm{CF}1}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))}{\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))}},$

Δf_C为垂测电离图智能判读软件自动给出的f_CF1相对于F₁层模型设定临频的偏差；

F₂层：移位切比雪夫多项式T_i(l)中的l具有下式所示形式：

$l=\frac{\ln (f_{\mathrm{NF}2}/f_{\mathrm{CF}2})}{\ln (f_{\mathrm{CF}1}/f_{\mathrm{CF}2})}0\≤l\≤1$

f_NF2表示F₂层等离子体频率；f_CF2表示F₂层临频；C_i(i＝0～N+1)为移位切比雪夫多项式系数，且： $C_{N+1}=h_{\mathrm{mF}1}-\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i$

h_mF2表示F₂层峰高，且：h_mF2＝C_N+1，

E层和谷层的连接点位于E层峰高h_mE，谷层与F₁层的连接点位于高度h₂处，并且在高度h₂处的等离子体频率等于E层临频f_CE，谷层包括两个部分：与E层的连接部分和与F₁层的连接部分，这两部分的连接点位于高度h₁处，F₁层与F₂层连接点位于F₁层峰高h_mF1处；

步骤12：为了使建立的电子浓度剖面满足连续光滑特性，在层与层的连接点处，基于连接点以上及以下电离层模型分别计算的等离子体频率的平方值以及剖面梯度应该相等，根据这一条件，限定相关参数之间的内在关系。

进一步的，所述步骤2具体包括：

步骤21：对于频率小于等于f_CE的电波将在E层反射，进行E层回波虚高计算公式的推导；

步骤22：对于频率大于f_CE、小于等于f_CF1的电波将在F₁层反射，进行F₁层回波虚高计算公式的推导，主要包括在E层传播的群距离Δh_E(f)、在谷层中与E层连接部分传播的群距离Δh′_J(f)、在谷层中与F₁层连接部分传播的群距离Δh′_V(f)和在F₁层传播的群距离Δh′_F1(f)的计算公式；

步骤23、对于频率大于f_CF1、小于等于f_CF2的电波将在F₂层反射，进行F₂层回波虚高计算公式的推导，主要包括Δh′_E(f)、Δh′_J(f)、Δh′_V(f)、Δh′_F1(f)和在F₂层传播的群距离Δh′_F2(f)的计算公式。

进一步的，所述步骤3具体包括：

步骤31：选取实测E层描迹数据，采用区域搜索的方法实现E层f_CE、h_mE(或h_bE)、y_mE三个参数的反演；

步骤32：选取F₁层描迹中大于E层临频和F₁层描迹最小虚高对应的频率之间的数据，基于搜索、迭代方法实现谷层参数的反演；

步骤33：选取F₁层描迹中，F₁层描迹最小虚高对应的频率到f_CF1之间的数据，利用约束寻优法实现F₁层参数的反演；

步骤34：选取F₂层描迹数据，利用约束寻优法实现F₂层参数的反演；

步骤35、对根据不同Δf_C得到的谷层、F₁层和F₂层参数，寻找使基于这些参数计算的所有数据点的计算虚高和实测虚高的误差和最小的一组参数，该组参数确定为最终的谷层、F₁层和F₂层参数。

本发明的有益效果在于：

本发明所公开的电离层垂直剖面建模及参数反演方法，基于移位切比雪夫多项式模型的约束优化F₁层参数、F₂层参数的垂测电离图反演方法，即得到谷参数后，选取F₁层较高区域回波描迹数据，在保证剖面连续光滑的约束条件下，计算F₁层剖面多项式系数，同样，在保证剖面连续光滑的约束条件下，选取F₂层回波描迹数据，计算F₂层剖面多项式系数，最后，基于所有数据点计算虚高和实测虚高误差和最小准则，选取对应初始设置下得到的剖面参数最终确定电离层剖面。该方法适用于F₁层未充分发育的情况，可以有效提高该种情况下反演精度和稳定性。

附图说明

图1为本发明所公开的电离层垂直剖面建模及参数反演方法的流程图；

图2为三层电离层反演实例。

具体实施方式

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

实施例1，如图1所示，本实施例公开了一种电离层垂直剖面建模及参数反演方法，包括以下步骤：

(1)建立电离层剖面数学模型：

本发明基于模式法的思想，将电离层建模为包含E层、谷层、F₁层、F₂层的四层模型，E层和谷层剖面用抛物模型表示，F₁层和F₂层剖面用移位切比雪夫多项式模型表示，电离层电子浓度剖面具有式(1)所示形式：

$\left\{\begin{array}{c}{f}_{\mathrm{NE}}^2=f_{\mathrm{CE}}^2\[1-{\left(\frac{h-h_{\mathrm{mE}}}{y_{\mathrm{mE}}}\right)}^2\]h_{\mathrm{bE}}\≤h\≤h_{\mathrm{mE}},h_{\mathrm{mE}}\≤h\≤h_1\\ f_{\mathrm{NV}}^2=f_{\mathrm{CV}}^2\[1+{\left(\frac{h-h_{\mathrm{mV}}}{y_{\mathrm{mV}}}\right)}^2\]h_1\≤h\≤h_2\\ h=A_{I+1}+g^{1/2}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{T}_i\left(g\right)h_2\≤h\≤h_{\mathrm{mF}1}\\ h=C_{N+1}+l^{1/2}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i{T}_i\left(l\right)h_{\mathrm{mF}1}\≤h\≤h_{\mathrm{mF}2}\end{array}\right.---\left(1\right)$

E层和谷层的连接点位于E层峰高h_mE，谷层与F₁层的连接点位于高度h₂处，并且在高度h₂处的等离子体频率等于E层临频f_CE，谷层包括两个部分：与E层的连接部分和与F₁层的连接部分，这两部分的连接点位于高度h₁处，F₁层与F₂层连接点位于F₁层峰高h_mF1处，式(1)中各符号的具体含义如下：

E层：

f_NE表示E层等离子体频率；f_CE表示E层临频；h_mE表示E层峰高；y_mE表示E层半厚；h_bE＝h_mE-y_mE表示E层底高；

谷层：

f_NV表示谷层等离子体频率；f_CV表示谷层最小等离子体频率；h_mV表示谷层等离子体频率为f_CV时对应的电离层高度；y_mV表示谷层半厚；h₂＝h_mE+W，W定义为谷层宽度；

F₁层：

T_i(g)为移位切比雪夫多项式，具有式(2)所示形式：

$\left\{\begin{array}{c}{T}_i\left(g\right)=2(2g-1)T_{i-1}\left(g\right)-T_{i-2}\left(g\right)\\ T_0\left(g\right)=1,T_1\left(g\right)=2g-1\end{array}\right.---\left(2\right)$

$g=\frac{\ln (f_{\mathrm{NF}1}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}{\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}\frac{\ln (f_{\mathrm{CF}1}/(F_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}{\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}\≤g\≤1---\left(3\right)$

f_NF1表示F₁层等离子体频率；f_CF1表示F₁层临频；A_i(i＝0～I+1)为移位切比雪夫多项式系数，且：

$A_{I+1}=h_2-\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i---\left(4\right)$

h_mF1表示F₁层峰高，且：

$h_{\mathrm{mF}1}=[A_{I+1}+g^{1/2}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{T}_i\left(g\right)]{|}_{g=\frac{\ln (f_{\mathrm{CF}1}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}{\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}}---\left(5\right)$

Δf_C为垂测电离图智能判读软件自动给出的f_CF1相对于F₁层模型设定临频的偏差；

F₂层：

移位切比雪夫多项式T_i(l)中的l具有式(6)所示形式：

$l=\frac{\ln (f_{\mathrm{NF}2}/f_{\mathrm{CF}2})}{\ln (f_{\mathrm{CF}1}/f_{\mathrm{CF}2})}0\≤l\≤1---\left(6\right)$

f_NF2表示F₂层等离子体频率；f_CF2表示F₂层临频；C_i(i＝0～N+1)为移位切比雪夫多项式系数，且：

$C_{N+1}=h_{\mathrm{mF}1}-\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i---\left(7\right)$

h_mF2表示F₂层峰高，且：

h_mF2＝C_N+1(8)

为了使建立的电子浓度剖面满足连续光滑特性，在层与层的连接点处，基于连接点以上及以下电离层模型分别计算的等离子体频率的平方值以及剖面梯度应该相等，根据这一条件，限定相关参数之间的内在关系，即：

1)在h＝h_mF1处，有：

$\left\{\begin{array}{c}{f}_{\mathrm{NF}2}^2{|}_{h=h_{\mathrm{mF}1}}=f_{\mathrm{NF}1}^2{|}_{h=h_{\mathrm{mF}1}}\\ \frac{\∂f_{\mathrm{NF}2}^2}{\∂h}{|}_{h=h_{\mathrm{mF}1}}=\frac{\∂f_{\mathrm{NF}1}^2}{\∂h}{|}_{h=h_{\mathrm{mF}1}}\end{array}\right.---\left(9\right)$

利用式(9)中的上式，可以由A_i(i＝0～I+1)计算得到C_N+1，如式(5)和式(7)所示。令则由式(9)中的下式可得：

$\frac{1}{4f_{\mathrm{CF}1}^2\ln (f_{\mathrm{CF}1}/f_{\mathrm{CF}2})}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i(1+2\frac{{\mathrm{dT}}_i\left(l\right)}{\mathrm{dl}}{|}_{l=1})=a---\left(10\right)$

式(10)为后续反演F₂层系数C_i(i＝0～N)的一个约束条件。

2)在h＝h₂处，有：

$\left\{\begin{array}{c}{f}_{\mathrm{NF}1}^2{|}_{h=h_2}=f_{\mathrm{NV}}^2{|}_{h=h_2}\\ \frac{\∂f_{\mathrm{NF}1}^2}{\∂h}{|}_{h=h_2}=\frac{\∂f_{\mathrm{NV}}^2}{\∂h}{|}_{h=h_2}\end{array}\right.---\left(11\right)$

考虑到 $f_{\mathrm{NF}1}^2{|}_{h=h_2}=f_{\mathrm{CE}}^2,$ 并且令 $D\stackrel{\mathrm{\Δ}}{=}{f}_{\mathrm{CE}}^2-f_{\mathrm{CV}}^2$ 以及 $\frac{\∂f_{\mathrm{NF}1}^2}{\∂h}{|}_{h=h_2}\stackrel{\mathrm{\Δ}}{=}\frac{1}{B},$ 即：

$B\stackrel{\mathrm{\Δ}}{=}\frac{1}{4f_{\mathrm{CE}}^2\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i(1+2\frac{{\mathrm{dT}}_i\left(g\right)}{\mathrm{dg}}{|}_{g=1})---\left(12\right)$

式(12)也是后续反演F₁层系数A_i(i＝0～I)的一个约束条件，则由式(11)得到：

$\left\{\begin{array}{c}{h}_{\mathrm{mV}}=h_2-2\mathrm{BD}\\ y_{\mathrm{mV}}=2B{f}_{\mathrm{CV}}{D}^{1/2}\end{array}\right.---\left(13\right)$

3)在h＝h₁处，有：

$\left\{\begin{array}{c}{f}_{\mathrm{NE}}^2{|}_{h=h_1}=f_{\mathrm{NV}}^2{|}_{h=h_1}\\ \frac{\∂f_{\mathrm{NE}}^2}{\∂h}{|}_{h=h_1}=\frac{\∂f_{\mathrm{NV}}^2}{\∂h}{|}_{h=h_1}\end{array}\right.---\left(14\right)$

由式(14)中的上式可以得到：

h_mV＝h_mE+[(4B²D+Q)D]^1/2(15)

其中，式(15)结合式(13)中的上式以及W可以进一步得到：

$D=\frac{W^2}{Q+4\mathrm{BW}}---\left(16\right)$

由式(14)中的下式可以给出h₁的计算结果为：

$h_1=\frac{{4B}^{2}D{h}_{\mathrm{mE}}+h_{\mathrm{mV}}Q}{Q+{4B}^{2}D}---\left(17\right)$

由以上推导可以看出，得到E层参数后，谷层只要确定B和W两个参数，则谷层剖面就可以确定，因此，在后续谷层参数反演中，我们只需确定B和W即可。

(2)各层反射回波虚高的计算：

通过对各层反射回波虚高的计算，得到计算虚高与实测虚高的误差量，从而可用于实现后续基于最小化误差量的各层参数反演。为了简化计算过程，同时又不会引入很大的误差，在计算电波在E层和谷层传播的群距离时没有考虑地磁场的影响，在计算电波在F₁层和F₂层传播的群距离时假设地磁场为一定值，即任意位置处的地磁场与垂测站上空300km高度处的地磁场一致。

1)E层回波虚高的计算

对于频率小于等于f_CE的电波将在E层反射，回波虚高计算公式为：

$h^{\′}\left(f\right)=h_{\mathrm{bE}}+{\∫}_{h_{\mathrm{bE}}}^{h_r}{\mathrm{\μ}}^{\′}\mathrm{dt}---\left(18\right)$

其中，f为电波频率，h_r为电波反射点处高度，μ′为群折射指数，不考虑地磁场时，具有如下形式：

${\mathrm{\μ}}^{\′}={(1-f_N^2/f^2)}^{-1/2}---\left(19\right)$

式中f_N表示对应位置处等离子体频率。

基于建立的E层电离层模型，则式(18)可以进一步计算得到：

$h^{\′}\left(f\right)=h_{\mathrm{bE}}+\frac{1}{2}{y}_{\mathrm{mE}}\frac{f}{f_{\mathrm{CE}}}\ln \frac{f_{\mathrm{CE}}+f}{f_{\mathrm{CE}}-f}---\left(20\right)$

2)F₁层回波虚高的计算

对于频率大于f_CE、小于等于f_CF1的电波将在F₁层反射，回波虚高计算公式为：

$h^{\′}\left(f\right)=h_{\mathrm{bE}}+{\∫}_{h_{\mathrm{bE}}}^{h_{\mathrm{mE}}}{\mathrm{\μ}}^{\′}\mathrm{dh}+{\∫}_{h_{\mathrm{mE}}}^{h_1}{\mathrm{\μ}}^{\′}\mathrm{dh}+{\∫}_{h_1}^{h_2}{\mathrm{\μ}}^{\′}\mathrm{dh}+{\∫}_{h_2}^{h_r}{\mathrm{\μ}}^{\′}\mathrm{dh}---\left(21\right)$

其中，式(21)中第二项为电波在E层传播的群距离，记为Δh′_E(f)，第三项为电波在谷层中与E层连接部分传播的群距离，记为Δh′_J(f)，第四项为电波在谷层中与F₁层连接部分传播的群距离，记为Δh′_V(f)，第五项为电波在F₁层传播的群距离，记为Δh′_F1(f)。

计算Δh′_E(f)、Δh′_J(f)、Δh′_V(f)时用到的μ′仍然具有式(19)的形式，计算结果分别为：

$\mathrm{\Δ}{h}_E^{\′}\left(f\right)=\frac{1}{2}{y}_{\mathrm{mE}}\frac{f}{f_{\mathrm{CE}}}\ln \frac{f+f_{\mathrm{CE}}}{f-f_{\mathrm{CE}}}---\left(22\right)$

$\mathrm{\Δ}{h}_J^{\′}\left(f\right)=f{Q}^{1/2}\ln \left\{\right[\sqrt{\frac{\mathrm{QD}}{Q+4B^{2}D}}+\sqrt{f^2-f_{\mathrm{CE}}^2+\frac{\mathrm{QD}}{Q+4B^{2}D}}\left]{(f^2-f_{\mathrm{CE}}^2)}^{-1/2}\right\}---\left(23\right)$

$\mathrm{\Δ}{h}_V^{\′}\left(f\right)=2\mathrm{fB}\sqrt{D}[\arcsin \left(\sqrt{\frac{D}{f^2-f_{\mathrm{CE}}^2+D}}\right)+\arcsin \left(\frac{2\mathrm{BD}}{\sqrt{(Q+4B^{2}D)(f^2-f_{\mathrm{CE}}^2+D)}}\right)]---\left(24\right)$

计算Δh′_F1(f)时用到的μ′具有以下形式：

${\mathrm{\μ}}^{\′}=\frac{G_o}{{\mathrm{\μ}}_o}---\left(25\right)$

其中，

${\mathrm{\μ}}_o=\sqrt{1-X_o}---\left(26\right)$

$X_o=f_N^2/f^2---\left(27\right)$

$G_o=\frac{{\mathrm{\μ}}_o}{n_o}\{1+\frac{X_o{\tan }^2\mathrm{\θ}}{M^2}[\frac{1+X_o}{{(1+\mathrm{\γ}{\mathrm{\μ}}_o^4)}^{1/2}}-\frac{2}{1+{(1+\mathrm{\γ}{\mathrm{\μ}}_o^4)}^{1/2}}\left]\right\}---\left(28\right)$

$\mathrm{\γ}=\frac{4{\tan }^2\mathrm{\θ}}{Y_o^2{\cos }^2\mathrm{\θ}}---\left(29\right)$

Y_o＝f_H/f(30)

$M=1+{\mathrm{\μ}}_o^2\frac{2{\tan }^2\mathrm{\θ}}{1+{(1+\mathrm{\γ}{\mathrm{\μ}}_o^4)}^{1/2}}---\left(31\right)$

${\left(\frac{{\mathrm{\μ}}_o}{n_o}\right)}^2=\frac{M}{1+2{\tan }^2\mathrm{\θ}/[1+{(1+\mathrm{\γ}{\mathrm{\μ}}_o^4)}^{1/2}]}---\left(32\right)$

式中，f_H为垂测站上空300km处磁旋频率，θ为垂测站上空300km处磁倾角。此时，无法直接给出Δh′_F1(f)的解析表达式，只能采用数值积分方法计算Δh′_F1(f)，考虑到在反射点处，μ′为无穷大，为了能够进行数值计算，作以下变量替换，即：

$f_{\mathrm{NF}1}^2=f^2-t^2(f^2-f_{\mathrm{CE}}^2)0\≤t\≤1---\left(33\right)$

则Δh′_F1(f)可以写为：

$\mathrm{\Δ}{h}_{F1}^{\′}\left(f\right)={\∫}_{h_2}^{h_r}{\mathrm{\μ}}^{\′}\mathrm{dh}={\∫}_{f_{\mathrm{CE}}^2}^{f^2}{\mathrm{\μ}}^{\′}\frac{\mathrm{dh}}{d{f}_{\mathrm{NF}}^2}d{f}_{\mathrm{NF}}^2={\∫}_0^{1}2{\mathrm{\μ}}^{\′}t(f^2-f_{\mathrm{CE}}^2)\frac{\mathrm{dh}}{d{f}_{\mathrm{NF}}^2}\mathrm{dt}---\left(34\right)$

根据建立的F₁层电离层模型(式(1)中的第三个式子)，可得到：

$\frac{\mathrm{dh}}{d{f}_{\mathrm{NF}1}^2}=\frac{1}{2\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))f_{\mathrm{NF}}^2}\left\{\frac{1}{2}{g}^{-1/2}\right[\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{T}_i\left(g\right)+2g\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i\frac{d{T}_i\left(g\right)}{\mathrm{dg}}\left]\right\}---\left(35\right)$

将式(35)代入到式(34)，则Δh′_F1(f)可进一步写为：

$\mathrm{\Δ}{h}_{F1}^{\′}\left(f\right)=\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{S}_i\left(f\right)---\left(36\right)$

式中，

$S_i\left(f\right)=\frac{f^2-f_{\mathrm{CE}}^2}{2\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}{\∫}_0^1\frac{{\mathrm{\μ}}^{\′}t}{f_{\mathrm{NF}1}^2{g}^{1/2}}[T_i\left(g\right)+2g\frac{d{T}_i\left(g\right)}{\mathrm{dg}}]\mathrm{dt}---\left(37\right)$

此时，在临近反射点(t→0)处，由式(26)～式(32)，可以得到μ_o→0，M→1，从而由式(25)得到所以，只要磁倾角S_i(f)就可以由式(37)计算得到(其中f_CF1可由垂测电离图智能判读软件自动给出)，从而可根据式(36)得到Δh′_F1(f)。

3)F₂层回波虚高的计算

对于频率大于f_CF1、小于等于f_CF2的电波将在F₂层反射，回波虚高计算公式为：

$h^{\′}\left(f\right)=h_{\mathrm{bE}}+\mathrm{\Δ}{h}_E^{\′}\left(f\right)+\mathrm{\Δ}{h}_J^{\′}\left(f\right)+\mathrm{\Δ}{h}_V^{\′}\left(f\right)+{\∫}_{h_2}^{h_{\mathrm{mF}1}}{\mathrm{\μ}}^{\′}\mathrm{dh}+{\∫}_{h_{\mathrm{mF}1}}^{h_r}{\mathrm{\μ}}^{\′}\mathrm{dh}---\left(38\right)$

其中，式(38)中第五项为电波在F₁层传播的群距离，记为Δh′_F1(f)，第六项为电波在F₂层传播的群距离，记为Δh′_F2(f)。

此时，Δh′_E(f)、Δh′_J(f)、Δh′_V(f)仍然可分别用式(22)、式(23)、式(24)计算，Δh′_F1(f)仍然可以采用式(36)计算，但由于此时电波在F₂层反射，其中的S_i(f)的计算公式变为：

$S_i\left(f\right)=\frac{f_{\mathrm{CF}1}^2-f_{\mathrm{CE}}^2}{2\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))}{\∫}_0^1\frac{{\mathrm{\μ}}^{\′}{t}^{\′}}{f_{\mathrm{NF}1}^2{g}^{1/2}}[T_i\left(g\right)+2g\frac{{\mathrm{dT}}_i\left(g\right)}{\mathrm{dg}}]{\mathrm{dt}}^{\′}---\left(39\right)$

其中，

$f_{\mathrm{NF}1}^2=f_{\mathrm{CF}1}^2-t^{\′2}(f_{\mathrm{CF}1}^2-f_{\mathrm{CE}}^2)0\≤t^{\′}\≤1---\left(40\right)$

Δh′_F2(f)采用与F₁层回波虚高中Δh′_F1(f)相同的计算方法得到：

${\mathrm{\Δh}}_{F2}^{\′}\left(f\right)=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i{P}_i\left(f\right)---\left(41\right)$

式中，

$P_i\left(f\right)=\frac{f^2-f_{\mathrm{CF}1}^2}{2\ln (f_{\mathrm{CF}1}/f_{\mathrm{CF}2})}{\∫}_0^1\frac{{\mathrm{\μ}}^{\′}s}{f_{\mathrm{NF}2}^2{l}^{1/2}}[T_i\left(l\right)+2l\frac{{\mathrm{dT}}_i\left(l\right)}{\mathrm{dl}}]\mathrm{ds}---\left(42\right)$

$f_{\mathrm{NF}2}^2=f^2-s^2(f^2-f_{\mathrm{CF}1}^2)0\≤s\≤1---\left(43\right)$

此时，在临近反射点(s→0)处，由式(25)～式(32)，得到f_CF2可由垂测电离图智能判读软件自动给出。

(3)各层参数的反演：

基于计算得到的虚高和实测虚高的差最小化准则，实现各层参数的反演。

1)E层参数的反演

由式(1)中上式可知，决定E层剖面的三个参数主要是f_CE、h_mE(或h_bE)、y_mE，其中f_CE可由垂测电离图智能判读软件自动给出，误差小于0.2MHz，本发明中采用一种区域搜索的方法实现E层参数的反演，具体为：

假设垂测电离图智能判读得到的E层描迹有K个点，其对应的工作频率和虚高分别为f_k和h″(f_k)，读出的E层临频和最小虚高分别记为为和h″_minE，则对参数f_CE、h_bE、y_mE分别在[h″_minE-δ₁，h″_minE+δ₂]、[0，δ₃](其中δ₁、δ₂和δ₃是搜索范围控制量)以一定步进取值得到不同组参数，每一组参数根据式(20)计算得到K个点的h′(f_k)，然后计算实测虚高和计算虚高的误差平方和：

$\mathrm{\ϵ}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(h^{\′\′}\left(f_k\right)-h^{\′}\left(f_k\right))}^2---\left(44\right)$

使ε达到最小的那组参数即确定为E层参数。

2)谷层参数的反演

本发明中基于搜索、迭代方法实现谷层B和W两个参数的反演。

在F₁层描迹中，大于E层临频和F₁层描迹最小虚高(记为h″_minF1)对应的频率之间的数据对谷层参数比较敏感，因此，在谷层参数的反演过程中，选择这部分描迹点用于确定谷层参数，假设共有K个点，其对应的工作频率和虚高分别为f_k和h″(f_k)。

谷层参数反演的基本步骤为：

①设定W＝0；

②设定B＝0，I＝7；

③基于最小二乘法计算F₁层剖面系数A_i(i＝0～I+1)；

④检查计算的系数A_i(i＝0～I)是否满足F₁层剖面单调递增特性，(a)如果不满足，则I＝I-1，如果I＜0，则执行⑤，否则，执行③；(b)如果满足，对A_I+1和前一次迭代记录的值进行比较，如果二者之差小于某一较小值(例如0.5km)，则计算K个点实测虚高和计算虚高的误差平方和，记录当前B、W以及计算的误差平方和的值，否则，在没有超过限定的最大迭代次数情况下，根据A_i(i＝0～I)自动调整I，按照式(12)更新B的值，执行③，若超过限定的最大迭代次数，则执行⑤；

⑤W＝W+1(单位为km)，如果W小于设定的搜索范围(例如0.7h″_minF1-h_mE)，则执行②，否则，执行⑥；

⑥找出记录的误差平方和的最小值，该最小值对应的B、W即确定为谷层参数，如果没有记录下有效的B和W，则W＝0，B＝0。

其中，上述步骤③的具体方法为：

根据式(21)给出K个点实测虚高h″(f_k)和计算虚高h′(f_k)的误差平方和：

$\begin{array}{c}\mathrm{\ϵ}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(h^{\′\′}\left(f_k\right)-h^{\′}\left(f_k\right))}^2=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(h^{\′\′}\left(f_k\right)-h_{\mathrm{bE}}-\mathrm{\Δ}{h}_E^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_J^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_V^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_{F1}^{\′}\left(f_k\right))}^2\\ =\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(\mathrm{\Δ}{h}^{\′}\left(f_k\right)-\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{S}_i\left(f_k\right))}^2\end{array}---\left(45\right)$

其中， $\mathrm{\Δ}{h}^{\′}\left(f_k\right)\stackrel{\mathrm{\Δ}}{=}{h}^{\′\′}\left(f_k\right)-h_{\mathrm{bE}}-\mathrm{\Δ}{h}_E^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_J^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_V^{\′}\left(f_k\right).$ 要使ε达到最小，即求解满足式(46)的系数A_i(i＝0～I)：

$\frac{\∂\mathrm{\ϵ}}{\∂A_j}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}2(\mathrm{\Δ}{h}^{\′}\left(f_k\right)S_j\left(f_k\right)-\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{S}_i\left(f_k\right)S_j\left(f_k\right))=0(j=0~I)---\left(46\right)$

式(46)进一步化简为：

$\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i\left(\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{S}_i\left(f_k\right)S_j\left(f_k\right)\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\mathrm{\Δ}{h}^{\′}\left(f_k\right)S_j\left(f_k\right)(j=0~I)---\left(47\right)$

求解上述方程组即可得到系数A_i(i＝0～I)，然后按照式(4)可计算A_I+1。

3)F₁层参数的反演

选取F₁层描迹中，F₁层描迹最小虚高对应的频率到f_CF1之间的数据用于确定F₁层参数，假设共有K个数据点，其对应的工作频率和虚高分别为f_k和h″(f_k)。在读取了f_CF1的情况下，F₁层剖面由系数A_i(i＝0～I+1)完全确定，可采用类似反演谷层参数时用到的方法来计算这些系数，这里需要注意的是，在谷层参数确定后，谷层与F₁层交点处的剖面梯度也已经确定，那么，根据当前数据计算出的系数A_i(i＝0～I)必须满足式(12)，因此，F₁层参数的反演实际上是一个约束优化问题，即：

$\underset{A_i}{\min }\mathrm{\ϵ}=\underset{A_i}{\min }\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(\mathrm{\Δ}{h}^{\′}\left(f_k\right)-\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{S}_i\left(f_k\right))}^2$

$s.t.=\frac{1}{4f_{\mathrm{CE}}^2\ln (f_{\mathrm{CE}}/(f_{\mathrm{CE}1}+\mathrm{\Δ}{f}_C))}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i(1+2\frac{{\mathrm{dT}}_i\left(g\right)}{\mathrm{dg}}{|}_{g=1})=B---\left(48\right)$

可采用拉格朗日法求解上述问题进行F₁层参数的反演，具体的为：

①根据式(48)建立一个新的函数：

②对式(49)中各个自变量求偏导数，建立方程组：

③求解上述方程组，得到系数A_i(i＝0～I)，检查计算的系数A_i(i＝0～I)是否满足F₁层剖面单调递增特性，(a)如果不满足，则I＝I-1，如果I＜0，执行④，否则，执行①；(b)如果满足，执行④；

④按照式(4)计算A_I+1。

4)F₂层参数的反演

F₂层描迹中的数据用于确定F₂层参数，假设共有K个数据点，其对应的工作频率和虚高分别为f_k和h″(f_k)。在读取了f_CF2的情况下，F₂层剖面由系数C_i(i＝0～N+1)完全确定，在F₁层参数确定后，F₁层与F₂层交点处的剖面梯度也已经确定，那么，根据当前数据计算出的系数C_i(i＝0～N)必须满足式(10)，因此，F₂层参数的反演同样是一个约束优化问题，可采用类似反演F₁层参数时用到的方法来计算这些系数：

$\underset{C_i}{\min }\mathrm{\ϵ}=\underset{C_i}{\min }\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{(\mathrm{\Δ}{\tilde{h}}^{\′}\left(f_k\right)-\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i{P}_i\left(f_k\right))}^2$

$s.t.\frac{1}{4f_{\mathrm{CF}1}^2\ln (f_{\mathrm{CF}1}/f_{\mathrm{CF}2})}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i(1+2\frac{d{T}_i\left(l\right)}{\mathrm{dl}}{|}_{l=1})=a---\left(51\right)$

式中， $\mathrm{\Δ}{\tilde{h}}^{\′}\left(f_k\right)\stackrel{\mathrm{\Δ}}{=}{h}^{\′\′}\left(f_k\right)-h_{\mathrm{bE}}-\mathrm{\Δ}{h}_E^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_J^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_V^{\′}\left(f_k\right)-\mathrm{\Δ}{h}_{F1}^{\′}\left(f_k\right).$ 同样，采用拉格朗日法求解上述问题进行F₂层参数的反演，具体的为：

①根据式(48)建立一个新的函数：

②对式(49)中各个自变量求偏导数，建立方程组：

③求解上述方程组，得到系数C_i(i＝0～N)，检查计算的系数C_i(i＝0～N)是否满足F₂层剖面单调递增特性，(a)如果不满足，则N＝N-1，如果N＜0，执行④，否则，执行①；(b)如果满足，执行④；

④按照式(7)计算C_N+1。

5)谷层、F₁层、F₂层参数的最终确定

对于未充分发育F₁层，垂测电离图智能判读软件自动给出的f_CF1相对于F₁层模型设定临频的偏差Δf_C是一个未知量，理论上，Δf_C的取值介于0～f_CF2-f_CF1，因此，在反演谷层、F₁层、F₂层参数时，将Δf_C在0～f_CF2-f_CF1内遍历，选取使所有数据点计算虚高与实测虚高误差和最小的Δf_C对应的谷层、F₁层、F₂层参数作为最终的谷层、F₁层、F₂层参数。

按照上述方法对三层电离层反演后的结果如图2所示。

Claims (4)

1.一种电离层垂直剖面建模及参数反演方法，其特征在于，所述方法包括以下步骤：

步骤1：将电离层建模为包含E层、谷层、F₁层、F₂层的四层模型，其中，E层和谷层剖面用抛物模型表示，F₁层和F₂层剖面用移位切比雪夫多项式模型表示；

步骤2：基于建立的电离层模型，推导各层回波虚高的计算公式；

步骤3：利用实测电离图数据，结合各层回波虚高计算结果，进行各层电离层参数的反演。

2.根据权利要求1所述的电离层垂直剖面建模及参数反演方法，其特征在于，所述步骤1具体包括：

步骤11：电离层电子浓度剖面具体形式如下式所示：

$\left\{\begin{array}{c}{f}_{\mathrm{NE}}^2=f_{\mathrm{CE}}^2[1-{\left(\frac{h-h_{\mathrm{mE}}}{y_{\mathrm{mE}}}\right)}^2]{,h}_{\mathrm{bE}}\≤h\≤h_{\mathrm{mE}},h_{\mathrm{mE}}\≤h\≤h_1\\ f_{\mathrm{NV}}^2=f_{\mathrm{CV}}^2[1+{\left(\frac{h-h_{\mathrm{mV}}}{y_{\mathrm{mV}}}\right)}^2]{,h}_1\≤h\≤h_2\\ h=A_{I+1}+g^{1/2}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{T}_i\left(g\right),h_2\≤h\≤h_{\mathrm{mF}1}\\ h=C_{N+1}+l^{1/2}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i{T}_i\left(l\right){,h}_{\mathrm{mF}1}\≤h\≤h_{\mathrm{mF}2}\end{array}\right.$

上式中各符号的具体含义如下：

E层：f_NE表示E层等离子体频率；f_CE表示E层临频；h_mE表示E层峰高；y_mE表示E层半厚；h_bE＝h_mE-y_mE表示E层底高；

谷层：f_NV表示谷层等离子体频率；f_CV表示谷层最小等离子体频率；h_mV表示谷层等离子体频率为f_CV时对应的电离层高度；y_mV表示谷层半厚；h₂＝h_mE+W，W定义为谷层宽度；

F₁层：T_i(g)为移位切比雪夫多项式，具有下式所示形式：

$\left\{\begin{array}{c}{T}_i\left(g\right)=2(2g-1)T_{i-1}\left(g\right)-T_{i-2}\left(g\right)\\ T_0\left(g\right)=1,T_1\left(g\right)=2g-1\end{array}\right.,g=\frac{\ln (f_{\mathrm{NF}1}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))\ln (f_{\mathrm{CF}1}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}{\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+\mathrm{\Δ}{f}_C))}\≤g\≤1$

f_NE1表示F₁层等离子体频率；f_CF1表示F₁层临频；A_i(i＝0～I+1)为移位切比雪夫多项式系数，且： $A_{I+1}=h_2-\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i,$

h_mF1表示F₁层峰高，且： $h_{\mathrm{mF}1}=[A_{I+1}+g^{1/2}\underset{i=0}{\overset{I}{\mathrm{\Σ}}}{A}_i{T}_i\left(g\right)]{|}_{g=\frac{\ln (f_{\mathrm{CF}1}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))}{\ln (f_{\mathrm{CE}}/(f_{\mathrm{CF}1}+{\mathrm{\Δf}}_C))}},$ Δf_C为垂测电离图智能判读软件自动给出的f_CF1相对于F₁层模型设定临频的偏差；

F₂层：移位切比雪夫多项式T_i(l)中的l具有下式所示形式：

$l=\frac{\ln (f_{\mathrm{NF}2}/f_{\mathrm{CF}2})}{\ln (f_{\mathrm{CF}1}/f_{\mathrm{CF}2})},0\≤l\≤1$

f_NF2表示F₂层等离子体频率；f_CF2表示F₂层临频；C_i(i＝0～N+1)为移位切比雪夫多项式系数，且： $C_{N+1}=h_{\mathrm{mF}1}-\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{C}_i$

h_mF2表示F₂层峰高，且：h_mF2＝C_N+1，

E层和谷层的连接点位于E层峰高h_mE，谷层与F₁层的连接点位于高度h₂处，并且在高度h₂处的等离子体频率等于E层临频f_CE，谷层包括两个部分：与E层的连接部分和与F₁层的连接部分，这两部分的连接点位于高度h₁处，F₁层与F₂层连接点位于F₁层峰高h_mF1处；

步骤12：为了使建立的电子浓度剖面满足连续光滑特性，在层与层的连接点处，基于连接点以上及以下电离层模型分别计算的等离子体频率的平方值以及剖面梯度应该相等，根据这一条件，限定相关参数之间的内在关系。

3.根据权利要求1所述的电离层垂直剖面建模及参数反演方法，其特征在于，所述步骤2具体包括：

步骤21：对于频率小于等于f_CE的电波将在E层反射，进行E层回波虚高计算公式的推导；

步骤22：对于频率大于f_CE、小于等于f_CF1的电波将在F₁层反射，进行F₁层回波虚高计算公式的推导，主要包括在E层传播的群距离Δh′_E(f)、在谷层中与E层连接部分传播的群距离Δh′_J(f)、在谷层中与F₁层连接部分传播的群距离Δh′_V(f)和在F₁层传播的群距离Δh′_F1(f)的计算公式；

步骤23、对于频率大于f_CF1、小于等于f_CF2的电波将在F₂层反射，进行F₂层回波虚高计算公式的推导，主要包括Δh′_E(f)、Δh′_J(f)、Δh′_V(f)、Δh′_F1(f)和在F₂层传播的群距离Δh′_F2(f)的计算公式。

4.根据权利要求1所述的电离层垂直剖面建模及参数反演方法，其特征在于，所述步骤3具体包括：

步骤31：选取实测E层描迹数据，采用区域搜索的方法实现E层f_CE、h_mE(或h_bE)、y_mE三个参数的反演；

步骤32：选取F₁层描迹中大于E层临频和F₁层描迹最小虚高对应的频率之间的数据，基于搜索、迭代方法实现谷层参数的反演；

步骤33：选取F₁层描迹中，F₁层描迹最小虚高对应的频率到f_CF1之间的数据，利用约束寻优法实现F₁层参数的反演；

步骤34：选取F₂层描迹数据，利用约束寻优法实现F₂层参数的反演；

步骤35、对根据不同Δf_C得到的谷层、F₁层和F₂层参数，寻找使基于这些参数计算的所有数据点的计算虚高和实测虚高的误差和最小的一组参数，该组参数确定为最终的谷层、F₁层和F₂层参数。