一种短波射线追踪技术中的电离层混合建模方法

摘要

一种短波射线追踪技术中的电离层混合建模方法，它由输入模块、IRI模块、QPS模块和射线追踪模块组成八个步骤：一：输入待测时间、工作频率、发射机坐标、接收机坐标、天线仰角范围等初始参数；二：由发射机、接收机坐标计算出反射中点位置、地心角、大圆距离、方位角；三：由IRI模型软件计算待测时间电波反射中点处的E层、F2层的有关参数；四：将由IRI模型算出的各层参数传递给QPS模型作为输入参量建立反射空间的多层准抛物等离子体频率分布；五：利用反射空间多层准抛物等离子体频率分布以群路径为自变量从发射机位置开始逐步求解射线追踪方程组，得到由发射机到接收机的电波群路径仿真值；六：在IRI模型等离子体频率分布下进行射线追踪计算，得到IRI模型下的电波群路径仿真值；七：开展试验，得到实测电波群路径；八：比较仿真值与实测值间的差距，验证混合建模方法的有效性。

说明书

(一)技术领域

本发明涉及一种基于国际参考电离层模型和多层准抛物模型的电离层混合建模方法，尤其涉及一种短波射线追踪技术中的电离层混合建模方法，属于短波通信技术领域。

(二)背景技术

在短波波段，电离层在一个波长内的变化很小，可以把电波近似为射线，利用射线追踪技术来研究电波在电离层中的传播和反射。电离层中的电子浓度分布称为电离层模型，射线追踪技术必须结合电离层模型进行，目前，在我国电离层射线追踪技术中使用的典型电离层模型主要有两种：一是国际参考电离层(International Reference Ionosphere)模型，简称IRI模型；二是多层准抛物(Quasi-Parabolic Segments)模型，也称QPS模型。电离层模型直接影响着短波射线追踪的效果，因此电离层模型的选择在短波通信、超视距雷达、卫星导航系统的设计与应用中具有重要意义。

电离层是地球大气层受太阳辐射电离产生的一种电中性媒质，由于电离层受太阳辐射、宇宙射线、地磁场等多种因素共同影响，其作用机理复杂，至今还没有一种模型能准确预报电离层变化。目前，在国内理论和工程领域应用的电离层模型主要是国际参考电离层和多层准抛物模型。

国际参考电离层(IRI)由国际空间委员会和国际无线电科学联合会共同推荐，是一种基于大量观测数据的电离层统计预报模式，反映了宁静电离层的平均状态。其电离层参量的高度变化由简单的数学函数描述，经度、纬度变化由正交多项式描述，采用CCIR和URSI系数来计算F2层临界频率和F2层最大电子浓度高度。该模型以一定精度普适于全球，是一种通用的电离层经验模型，IRI模型每5年更新一次，目前最新的是IRI-2007。同时国际无线电联合会还编制了IRI模型的计算软件，只要输入日期、具体时间和观测点经纬度，就可算出对应时间观测点处电离层各层的最大电子浓度、临界频率、最大电子浓度高度、半厚等多种参数，通过该模型可以预测未来几年内地球上任意位置处的电离层状态，在与电离层状态密切相关的短波通信、超视距雷达、卫星导航等领域具有重要指导意义。IRI电离层模型和相关计算软件可以从有关网站下载得到。虽然在模型编制过程中，因中国地区电离层观测数据较少，对中国地区电离层状态的预报精度不如对欧美地区电离层状态预报精度高，但该模型输入简单，使用方便，并且描述了未来几年全球的电离层状态，在我国仍然得到了广泛应用。

另一种常用的电离层模型是准抛物模型，因其采用准抛物线方程来描述空间中的电子浓度分布而得名。因为电子浓度N与等离子体频率f_N间存在以下关系

$f_N^2=80.6N$

所以常用等离子体频率分布表示电子浓度分布，下式为单层准抛物模型：

$f_N^2=\left(\begin{array}{cc}a-b{(1-\frac{r_m}{r})}^2& r_b<r<r_m\left(\frac{r_b}{r_b-y_m}\right)\\ 0& \mathrm{others}\end{array}\right)$

式中，r为空间位置距地心的距离(单位km)，r_b为电离层底高(距地心的距离，单位km)，y_m为层的半厚度(单位km)，r_m＝r_b+y_m为最大电子浓度对应的高度(单位km)，b＝a(r_b/y_m)²，f_N为等离子体频率(单位MHz)，大小只与电子浓度有关，f_o为电离层临界频率(单位MHz)，代表电离层能够反射的最高频率。

1988年，Dyson和Bennett引入了多层准抛物模型(QPS模型)。设电离层由E层、F1层、F2层组成，且各层均为宁静电离层，不考虑地磁、扰动、梯度等影响，则各层的等离子体频率剖面可表示如下：

$f_{\mathrm{NE}}^2=a_E-b_E{(1-\frac{r_{\mathrm{mE}}}{r})}^2$E层

$f_{\mathrm{Nj}1}^2=a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r})}^2$第一连接层，连接E层和F1层

$f_{\mathrm{NF}1}^2=a_{F1}-b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r})}^2$F1层

$f_{\mathrm{Nj}2}^2=a_{j2}-b_{j2}{(1-\frac{r_{j2}}{r})}^2$第二连接层，连接F1层和F2层

$f_{\mathrm{NF}2}^2=a_{F2}-b_{F2}{(1-\frac{r_{\mathrm{mF}2}}{r})}^2$F2层

上式中，r为空间位置距地心的距离，f_NE为E层等离子体频率，f_oE为E层临界频率，r_bE为E层底高，y_mE为E层的半厚度，r_mE＝r_bE+y_mE为E层最大电子浓度对应高度，b_E＝a_E(r_bE/y_mE)²；f_Nj1为第一连接层等离子体频率，f_oj1为第一连接层临界频率，r_bj1为第一连接层底高，y_mj1为第一连接层的半厚度，r_j1＝r_bj1+y_mj1为第一连接层最大电子浓度对应高度，b_j1＝a_j1(r_bj1/y_mj1)²；f_NF1为F1层等离子体频率，f_oF1为F1层临界频率，r_bF1为F1层底高，y_mF1为F1层的半厚度，r_mF1＝r_bF1+y_mF1为F1层最大电子浓度对应高度，b_F1＝a_F1(r_bF1/y_mF1)²；f_Nj2为第二连接层等离子体频率，f_oj2为第二连接层临界频率，r_bj2为第二连接层底高，y_mj2为第二连接层的半厚度，r_j2＝r_bj2+y_mj2为第二连接层最大电子浓度对应高度，b_j2＝a_j2(r_bj2/y_mj2)²；f_NF2为F2层等离子体频率，f_oF2为F2层临界频率，r_bF2为F2层底高，y_mF2为F2层的半厚度，r_mF2＝r_bF2+y_mF2为F2层最大电子浓度对应高度，b_F2＝a_F2(r_bF2/y_mF2)²。其中，两个连接层仅起平滑连接E层、F1层、F2层的作用，厚度很薄。

上面介绍的两种电离层模型是使用较多的两种通用模型，两种模型采用不同的方程来描述空间中的电子浓度分布，在实际使用中各有优缺点。

IRI模型输入参数简单，但在我国使用该模型进行射线追踪结果精度较低，特别是在F2层估算结果与实测值差异较大，不适用于如超视距雷达等对计算精度要求较高的场合。虽然采用解析式的方法来描述电离层，只需要知道各层的特征参数(底高、半厚、临界频率)就可建立相应的电离层模型，结构简单，使用方便，同时经过中国电波传播研究所的验证，在我国中纬度地区使用QPS模型进行射线追踪能够得到较为准确的结果，但在使用该模型进行射线追踪时必须先输入电波反射区上空各层的电离层基本参数(如：底高、半厚、临界频率等特征参数)用来建立模型，而这些重要的参数往往事先并不知道，需要对电波反射区进行探测获得或输入经验值。在无法开展电离层探测和无法取得准确经验值的情况下，使用该模型进行射线追踪的误差反而会增大，不便进行实时追踪和预测，这样就限制了QPS模型的使用。而IRI模型能对过去几十年以及未来五年(自发布之日算起)内全球各地的电离层状态进行预测，通过IRI模型能够得到给定时间地点下电离层各层底高、半厚、临界频率等参数，恰好弥补了QPS模型无法预测电离层参数的不足，因此，IRI模型和QPS模型混合建模的数字射线追踪技术就具有了重要意义。为了完善IRI模型和QPS模型在工程应用中的不足，在此基础上，提出了混合使用IRI模型和QPS模型进行电离层建模的方法。

(三)发明内容

(1)发明目的：

IRI模型和QPS模型是在短波射线追踪技术中使用的两种典型模型，IRI模型使用方便，只需输入时间、发射/接收机位置即可预测地球任意区域上空的电离层变化情况，在无法获知电波反射区域电离层参数时，可以使用该模型进行射线追踪，预测电波群路径、地面距离、电波反射高度、发射仰角等信息，但该模型在我国中纬度地区的预测精度较低，特别是在F2层与实测结果相差较大，单一使用该模型进行射线追踪，估算结果无法满足诸如导航、超视距雷达等对电离层参数精度要求较高的应用领域。QPS模型简单，易于构造，只需知道各层的底高、半厚和临频就可建立电离层电子浓度分布，相对IRI模型，其仿真值与我国中纬地区实测值更为接近，但在使用该模型进行射线追踪时，必须获得反射区域电离层临频、底高、最大电子浓度高度等基本参数，在无法获知电离层信息时使用不便。本发明提供一种短波射线追踪技术中的电离层混合建模方法，它将IRI模型和QPS模型混合应用，由IRI模型计算出反射区电离层参数后代入QPS模型中进行射线追踪，既确保了仿真输入条件的简洁方便，解决了QPS模型输入参数获取不便的难题，同时又提高了仿真估值的精确程度，避免了单纯使用IRI模型仿真结果精度低的问题。经过与实测结果的比对，使用混合模型进行射线追踪对白天F2层的计算精度要高于由IRI模型进行射线追踪的计算精度，验证了该方法的可行性与有效性。

(2)技术方案：

该方法是由四个模块组成方法步骤：

本发明一种短波射线追踪技术中的电离层混合建模方法，它是基于国际参考电离层和多层准抛物电离层模型的混合应用来实现一种新的电离层建模方法，该方法兼具两种电离层模型的优点，图1是电离层混合建模方法的流程图。该方法由四个模块组成方法步骤：它们是输入模块、IRI模块，QPS模块，射线追踪模块。四个模块间的位置、连接关系及信号走向如下述：用户在输入模块中输入待测时间、发射机与接收机地理坐标、天线仰角范围及工作频率等初始参数，由该模块自动算出电波反射中点地理坐标、地心角、大圆距离及方位角，该模块将待测时间和算出的电波反射点坐标传递给IRI模块供生成相应的电离层参数使用，同时将除待测时间和反射点位置外的其他参数传递给射线追踪模块供射线追踪计算使用；IRI模块根据待测时间和电波反射中点位置计算出电波发射区域的E层临界频率、E层最大电子浓度高度、E层半厚、F2层临界频率、F2层最大电子浓度高度、F2层半厚等参数；将由IRI模块算出的电离层参数输入到QPS模块中，建立多层准抛物线形式的等离子体频率分布；射线追踪模块使用工作频率、方位角、天线仰角范围、多层准抛物等离子体频率分布作为初始参数以固定步长逐步计算由发射机到接收机的电波群路径。

所述输入模块是用户与算法的主要交互接口，以界面形式内嵌在电离层混合建模方法中，该模块流程图参见图2所示。该模块与用户的主要接口包括：时间、发射机坐标、接收机坐标、工作频率、天线仰角范围。它供用户输入相关使用参数，用户在此模块中可设置待测时间、工作频率、发射机坐标、接收机坐标、天线仰角范围等初始参数，该模块可以根据发射机、接收机坐标算出电波反射中点位置、地心角、大圆距离及从发射机到接收机的方位角，以上参数是射线追踪方程组的初始值，该模块一方面将电波反射点位置和待测时间传递给IRI模块供生成相应的电离层参数使用，同时将除去待测时间和反射点位置外的其他参数传递给射线追踪模块供射线追踪计算使用。

所述IRI模块是由国际参考电离层模型来构建电子浓度分布的计算模块，该模块所用到的国际参考电离层(IRI)模型和计算程序都可从美国宇航局网站下载获得。该模块流程图见图3所示，它以软件封装的形式存在，模型所用到的公式和CCIR系数都已封装在软件中，使用由上一模块中传来的反射点位置和时间即可生成反射区域各层电离层的关键参数并传递给QPS模块用以建立反射区域等离子体频率分布，同时IRI模块根据IRI模型自动生成反射区域等离子体频率分布并直接传递给射线追踪模块。

所述QPS模块是由多层准抛物模型来构建等离子体频率分布的计算模块，是自行开发的程序，流程图见图4所示，该模块以封装的形式内嵌在电离层混合建模方法中，作为混合建模方法的一个模块，QPS模型所用到的公式都封装在该程序模块中。该模块使用上一步中由IRI模块估计的各层电离层关键参数来构建反射区域上空等离子体频率分布，并将各点等离子体频率传递给射线追踪模块。

该模块中设电离层由E层、F1层、F2层组成，且各层均为宁静电离层，不考虑地磁、扰动、梯度等影响，则各层的等离子体频率剖面可表示如下：

$f_{\mathrm{NE}}^2=a_E-b_E{(1-\frac{r_{\mathrm{mE}}}{r})}^2$E层

$f_{\mathrm{Nj}1}^2=a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r})}^2$第一连接层，连接E层和F1层

$f_{\mathrm{NF}1}^2=a_{F1}-b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r})}^2$F1层

$f_{\mathrm{Nj}2}^2=a_{j2}-b_{j2}{(1-\frac{r_{j2}}{r})}^2$第二连接层，连接F1层和F2层

$f_{\mathrm{NF}2}^2=a_{F2}-b_{F2}{(1-\frac{r_{\mathrm{mF}2}}{r})}^2$F2层

上式中，r为空间位置距地心的距离，f_NE为E层等离子体频率，f_oE为E层临界频率，r_bE为E层底高，y_mE为E层的半厚度，r_mE＝r_bE+y_mE为E层最大电子浓度对应高度，b_E＝a_E(r_bE/y_mE)²；f_Nj1为第一连接层等离子体频率，f_oj1为第一连接层临界频率，r_bj1为第一连接层底高，y_mj1为第一连接层的半厚度，r_j1＝r_bj1+y_mj1为第一连接层最大电子浓度对应高度，b_j1＝a_j1(r_bj1/y_mj1)²；f_NF1为F1层等离子体频率，f_oF1为F1层临界频率，r_bF1为F1层底高，y_mF1为F1层的半厚度，r_mF1＝r_bF1+y_mF1为F1层最大电子浓度对应高度，b_F＝a_F1(r_bF1/y_mF1)²；f_Nj2为第二连接层等离子体频率，f_oj2为第二连接层临界频率，r_bj2为第二连接层底高，y_mj2为第二连接层的半厚度，r_j2＝r_bj2+y_mj2为第二连接层最大电子浓度对应高度，b_j2＝a_j2(r_bj2/y_mj2)²；f_NF2为F2层等离子体频率，f_oF2为F2层临界频率，r_bF2为F2层底高，y_mF2为F2层的半厚度，r_mF2＝r_bF2+y_mF2为F2层最大电子浓度对应高度，b_F2＝a_F2(r_bF2/y_mF2)²。其中，两个连接层仅起平滑连接E层、F1层、F2层的作用，厚度很薄。

对各层等离子体频率分布求导，可以得到各层等离子体频率随高度的变化率

$\frac{{\mathrm{df}}_{\mathrm{NE}}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mE}}{b}_E}{r^2}(1-\frac{r_{\mathrm{mE}}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}1}^2}{\mathrm{dr}}=-\frac{{2r}_{j1}{b}_{j1}}{r^2}(1-\frac{r_{j1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}1}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{r^2}(1-\frac{r_{\mathrm{mF}1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}2}^2}{\mathrm{dr}}=-\frac{{2r}_{j2}{b}_{j2}}{r^2}(1-\frac{r_{j2}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}2}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}2}{b}_{F2}}{r^2}(1-\frac{r_{\mathrm{mF}2}}{r})$

将下述步骤三中得到的E层临界频率f_oE、E层最大电子浓度高度r_mE、E层半厚y_mE、F2层临界频率f_oF2、F2层最大电子浓度高度r_mF2、F2层半厚y_mF2及连接层参数代入上式中，得到电波反射区域等离子体频率的梯度分布。

第一连接层和第二连接层的参数可由下式求得，假设第一连接层与E层的交点在E层的最大电子浓度处，则根据两层的等离子体频率及其梯度分别相等，有

a_j1＝a_E

r_j1＝r_mE

设第一连接层与F1层的交点在r＝r_c1处，则在r_c1处两层的等离子体频率及其梯度分别相等，有

$-\frac{{2r}_{j1}{b}_{j1}}{{r_{c1}}^2}(1-\frac{r_{j1}}{r_{c1}})=\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{{r_{c1}}^2}(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})$

$a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r_{c1}})}^2=a_{F1}+b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})}^2$

联立上两式解得

$r_{c1}=\frac{r_{\mathrm{mF}1}{b}_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}{a_{F1}-a_E+b_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}$

$b_{j1}=\frac{r_{\mathrm{mF}1}{r}_{c1}{b}_{F1}-r_{\mathrm{mF}1}^2{b}_{F1}}{r_{\mathrm{mE}}^2-r_{\mathrm{mE}}{r}_{c1}}$

同理得到

a_j2＝a_F1，r_j2＝r_mF1，$r_{c2}=\frac{r_{\mathrm{mF}2}{b}_{F2}(1-r_{\mathrm{mF}2}/r_{j2})}{a_{j2}-a_{F2}+b_{F2}(1-r_{\mathrm{mF}2}/r_{j2})},$$b_{j2}=\frac{r_{\mathrm{mF}2}^2{b}_{F2}-r_{\mathrm{mF}2}{r}_{c2}{b}_{F2}}{r_{j2}^2-r_{j2}{r}_{c2}}$

其中，r_c1为第一连接层与E层交点处的高度，r_c2为第二连接层与F2层交点处的高度。

所述射线追踪模块是电波群路径的计算模块，该模块系自行开发的程序，流程图可参考图5所示，所用到的射线追踪方程组都封装在该模块中。它以输入模块中的工作频率、天线仰角范围、方位角等参数作为射线追踪方程组的初值，根据反射区域上空等离子体频率分布由发射机开始迭代计算电波传播位置及方向，最后取距接收机最近的一点作为从发射机到接收机的传播轨迹，将传播轨迹描迹显示并得到电波群路径。

该模块从球坐标系下以电波传播群路径P′为自变量的射线追踪方程组出发

方程组用于计算射线路径上点的坐标及该点处的波矢量。式中，P′为群路径；r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；ω为圆频率，c为光速，n为折射指数，与电子浓度分布有关。

$n^2=1-\frac{2X(1-\mathrm{iZ}-X)}{2(1-\mathrm{iZ})(1-\mathrm{iZ}-X)-Y_T^2\pm \sqrt{T_T^4+4Y_L^2{(1-\mathrm{iZ}-X)}^2}}$

式中，Y＝f_H/f，Y_T＝Ysinψ，Y_L＝Ycosψ，Z＝υ/2πf，f_H＝eB/2πm。其中，f_N为等离子体频率，是一个与电子浓度有关的物理量，由上一步的QPS模型构建；f为电波频率；f_H为磁旋频率；ψ为波矢量与地磁场间的夹角；υ为电子碰撞频率；N为电子浓度；e为电子电量；m为电子质量；ε₀为真空中介电常数；B为地磁场强度。“+”对应寻常波(0波)，“-”对应非常波(X波)。

因碰撞只造成能量的吸收，不影响电波传播方向，故在射线追踪中无须考虑碰撞效应，取Z＝0，在不考虑地磁场影响时，取Y＝0，则折射指数简化为

$n^2=1-x=1-\frac{f_N^2}{f^2}$

则有

$\frac{\mathrm{dX}}{\mathrm{dr}}=\frac{1}{f^2}\frac{{\mathrm{df}}_N^2}{\mathrm{dr}}$

若在射线追踪过程中不考虑电离层梯度变化的影响，则为0。

设电波发射仰角为β，方位角为α，工作圆频率为ω，发射机纬度、经度分别为λ₀和则射线的初始位置为：

r＝r₀

θ＝π/2-λ₀

$k_r=\left|\overrightarrow{k}\right|\cos (\pi /2-\beta )=\omega \sin \beta /c$

$k_{\theta }=-\left|\overrightarrow{k}\right|\cos \beta \cos \alpha =-\omega \cos \beta \cos \alpha /c$

式中，r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；r₀为地球半径，根据r、θ、k_r、k_θ、六个参数的初始值，即可开始射线轨迹的计算。

在射线路径上，以群路径为自变量，射线的坐标矢量和波矢量随着群路径的变化而变化，利用常微分方程数值解法Runge-Kutta法以固定步长逐点求解该射线追踪方程组，得到射线路径上各点的r、θ、k_r、k_θ、其中每换一点，都要计算一次该点的等离子体参量变化率将这些点连接起来，即得到一条球坐标系中的射线轨迹，在终点处，可得到电波群路径、地面距离、接收点经纬度等信息。

本发明一种短波射线追踪技术中的电离层混合建模方法，该方法具体步骤如下：

步骤一：在输入模块中首先输入待测时间、工作频率、发射机坐标、接收机坐标、天线仰角范围等初始参数。

步骤二：由发射机、接收机坐标计算出反射中点位置M、地心角d、大圆距离D和方位角为α，计算公式参见下式

$M_{\mathrm{Lat}}=\frac{T_{\mathrm{Lat}}+R_{\mathrm{Lat}}}{2},$$M_{\mathrm{Lon}}=\frac{T_{\mathrm{Lon}}+R_{\mathrm{Lon}}}{2},$$\mathrm{\alpha =arccos}\frac{\sin R_{\mathrm{Lat}}-{\sin T}_{\mathrm{Lat}}\cos d}{\cos T_{\mathrm{Lat}}\sin d}$

d＝arccos(sinT_Lat sinR_Lat+cosT_Lat cosR_Lat cos(T_Lon-R_Lon))，D＝r₀·d

上式中，M_Lat为反射中点纬度，M_Lon为反射中点经度，T_Lat为发射机纬度，T_Lon为发射机经度，R_Lat为接收机纬度，R_Lon为接收机经度，α为从发射机到接收机的方位角，以正北为基准东偏为正，d为地心角，是收发两点与地心连线间的夹角，单位为弧度，r₀为地球半径，D为大圆距离，是收发两点间的地面距离，单位为km。将算出的反射点地理坐标和待测时间传递给IRI模块，将其他参数送入射线追踪模块。

步骤三：在IRI模块中，待测时间地点上空的电离层参数由国际电信联盟ITU推荐使用的CCIR系数得到，其具体系数和算法已封装到IRI模型计算程序中，由国际无线电联合会推荐使用，该模型和计算程序可从美国宇航局网站下载获得。在模型中输入待测时间和电波反射中点地理坐标可以由软件自动生成对应位置上空的E层临界频率f_oE、E层最大电子浓度高度r_mE、E层半厚y_mE、F2层临界频率f_oF2、F2层最大电子浓度高度r_mF2、F2层半厚y_mF2等参数，得到待测时间电波反射点区域的等离子体频率分布。将以上算出的电离层各层关键参数传递给QPS模块用以在步骤四中建立反射区域等离子体频率分布，同时IRI模块根据IRI模型算出反射区域等离子体频率并直接传递给射线追踪模块供步骤六使用。

步骤四：将由IRI模型算出的各时刻电离层参数传递给QPS模块作为输入参量建立反射空间的多层准抛物等离子体频率分布和等离子体频率梯度分布。

设电离层由E层、F1层、F2层组成，且各层均为宁静电离层，不考虑地磁、扰动、梯度等影响，则各层的等离子体频率剖面可表示如下：

$f_{\mathrm{NE}}^2=a_E-b_E{(1-\frac{r_{\mathrm{mE}}}{r})}^2$E层

$f_{\mathrm{Nj}1}^2=a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r})}^2$第一连接层，连接E层和F1层

$f_{\mathrm{NF}1}^2=a_{F1}-b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r})}^2$F1层

$f_{\mathrm{Nj}2}^2=a_{j2}-b_{j2}{(1-\frac{r_{j2}}{r})}^2$第二连接层，连接F1层和F2层

$f_{\mathrm{NF}2}^2=a_{F2}-b_{F2}{(1-\frac{r_{\mathrm{mF}2}}{r})}^2$F2层

上式中各参数物理意义可参见背景技术中的说明。

对各层等离子体频率分布求导，可以得到各层等离子体频率随高度的变化率

$\frac{{\mathrm{df}}_{\mathrm{NE}}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mE}}{b}_E}{r^2}(1-\frac{r_{\mathrm{mE}}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}1}^2}{\mathrm{dr}}=-\frac{{2r}_{j1}{b}_{j1}}{r^2}(1-\frac{r_{j1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}1}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{r^2}(1-\frac{r_{\mathrm{mF}1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}2}^2}{\mathrm{dr}}=-\frac{{2r}_{j2}{b}_{j2}}{r^2}(1-\frac{r_{j2}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}2}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}2}{b}_{F2}}{r^2}(1-\frac{r_{\mathrm{mF}2}}{r})$

将步骤三中得到的E层临界频率f_oE、E层最大电子浓度高度r_mE、E层半厚y_mE、F2层临界频率f_oF2、F2层最大电子浓度高度r_mF2、F2层半厚y_mF2及连接层参数代入上式中，得到电波反射区域等离子体频率的梯度分布。

第一连接层和第二连接层的参数可由下式求得，假设第一连接层与E层的交点在E层的最大电子浓度处，则根据两层的等离子体频率及其梯度分别相等，有

a_j1＝a_E

r_j1＝r_mE

设第一连接层与F1层的交点在r＝r_c1处，则在r_c1处两层的等离子体频率及其梯度分别相等，有

$-\frac{{2r}_{j1}{b}_{j1}}{{r_{c1}}^2}(1-\frac{r_{j1}}{r_{c1}})=\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{{r_{c1}}^2}(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})$

$a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r_{c1}})}^2=a_{F1}+b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})}^2$

联立上两式解得

$r_{c1}=\frac{r_{\mathrm{mF}1}{b}_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}{a_{F1}-a_E+b_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}$

$b_{j1}=\frac{r_{\mathrm{mF}1}{r}_{c1}{b}_{F1}-r_{\mathrm{mF}1}^2{b}_{F1}}{r_{\mathrm{mE}}^2-r_{\mathrm{mE}}{r}_{c1}}$

同理得到

a_j2＝a_F1，r_j2＝r_mF1，$r_{c2}=\frac{r_{\mathrm{mF}2}{b}_{F2}(1-r_{\mathrm{mF}2}/r_{j2})}{a_{j2}-a_{F2}+b_{F2}(1-r_{\mathrm{mF}2}/r_{j2})},$$b_{j2}=\frac{r_{\mathrm{mF}2}^2{b}_{F2}-r_{\mathrm{mF}2}{r}_{c2}{b}_{F2}}{r_{j2}^2-r_{j2}{r}_{c2}}$

其中，r_c1为第一连接层与E层交点处的高度，r_c2为第二连接层与F2层交点处的高度。

步骤五：利用反射空间多层准抛物等离子体频率分布以群路径为自变量从发射机位置开始逐点求解射线追踪方程组，得到由发射机到接收机的电波群路径。

从球坐标系下以电波传播群路径P′为自变量的射线追踪方程组出发

方程组用于计算射线路径上点的坐标及该点处的波矢量。式中，P′为群路径；r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；ω为圆频率，c为光速，n为折射指数，与等离子体频率分布有关，在不考虑碰撞和地磁场时，折射指数可表示为

n²＝1-X

式中，f_N为等离子体频率，f为电波频率。

$\frac{\mathrm{dX}}{\mathrm{dr}}=\frac{1}{f^2}\frac{{\mathrm{df}}_N^2}{\mathrm{dr}}$

若在射线追踪过程中不考虑电离层梯度变化的影响，则为0。

由步骤一和步骤二得到的参量，电波发射仰角β、方位角α、工作圆频率ω、发射机纬度λ₀和经度可以确定射线的初始位置为：

r＝r₀

θ＝π/2-λ₀

$k_r=\left|\overrightarrow{k}\right|\cos (\pi /2-\beta )=\omega \sin \beta /c$

$k_{\theta }=-\left|\overrightarrow{k}\right|\cos \beta \cos \alpha =-\omega \cos \beta \cos \alpha /c$

式中，r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；r₀为地球半径，根据r、θ、k_r、k_θ、六个参数的初始值，即可开始射线轨迹的计算。

采用Runge-Kutta法求解射线追踪方程组，群路径P′取固定步长，每次求解方程组得到当前电波位置r、θ、和波矢量k_r、k_θ、计算该位置处的等离子体浓度fN和折射指数n，再将解得的r、θ、k_r、k_θ、f_N、n代入方程组中求得下一步电波位置和波矢量，如此循环，直到电波到达地面。电波仰角以固定步长在一定范围内扫描，得到射线轨迹后，改变电波初始仰角，每个仰角对应不同的电波轨迹，分别选取由E层和F层反射的落点位置距接收机最接近的一条射线，得到从发射机到接收机的E层反射群路径和F层反射群路径。

步骤六：用步骤三中IRI模块直接计算出的(所需CCIR系数和算法已内嵌到模型中)待测时间反射区域IRI模型下的等离子体频率代替步骤四中由QPS模块算出的等离子体频率f_N，然后重复步骤五的射线追踪过程，得到IRI模型下的电波群路径。

步骤七：按照步骤一中的仿真条件，利用中国电波传播研究所现有实验设备和观测站点，进行电波斜向探测试验，得到实测电波群路径；

步骤八：比较由步骤五得到的混合模型下仿真群路径和由步骤六得到的IRI模型下仿真群路径与实测值间均方距离的大小，如果混合模型下仿真群路径与实测值间的均方距离小于IRI模型仿真群路径与实测值间均方距离，说明采用混合模型进行射线追踪的精确程度更高。

为了描述模型仿真值与实测值间的差异程度，定义仿真值与实测值间的均方距离S和均方距离差异率，均方距离S为

式中，x_i是i时刻对应的仿真群路径，y_i是i时刻对应的实测群路径，i取从1到N的整数，对应N个不同的试验时间。S越小，说明仿真值与实测值距离越小，总体上越接近实测值。

$S\%=S/\overline{P^{\prime }}\times 100\%$

上式中，S％为均方距离差异率，S为均方距离，为实测群路径均值。S％越小，仿真值与实测值越接近。

(3)优点及功效：

本发明一种短波射线追踪技术中的电离层混合建模方法，它由IRI模型计算出反射区电离层参数后代入QPS模型中进行射线追踪，既确保了仿真输入条件的简洁方便，解决了QPS模型输入参数获取不便的难题，同时又提高了仿真估值的精确程度，避免了单纯使用IRI模型仿真结果精度低的问题。经过与实测结果的对比，使用混合模型进行射线追踪对白天F2层的计算精度要高于由IRI模型进行射线追踪的计算精度，验证了该方法的可行性与有效性。同时所需设备简单，仅需要一台可以顺利运行Visual C++开发软件的普通个人计算机(如：P4 1.6GHz主频CPU、1G内存、Windows XP操作系统)就可进行电离层混合建模仿真工作，便于广泛开展。

(四)附图说明

图1混合建模流程示意图

图2输入模块流程示意图

图3IRI模块流程示意图

图4QPS模块流程示意图

图5射线追踪模块流程示意图

图6电离层斜向探测示意图

图710月28日E层群路径仿真值与实测值对比示意图

图810月28日F2层群路径仿真值与实测值对比示意图

图910月31日E层群路径仿真值与实测值对比示意图

图1011月1日F2层群路径仿真值与实测值对比示意图

图中符号说明如下：

IRI为国际参考电离层；QPS模型为多层准抛物模型；

r为球坐标下射线距地心的距离；  6370为地球半径；

T表示探测时间；  G表示电波群路径，单位为km；  E表示电离层的E层；

F2表示电离层的F2层；  f表示电波频率；

图7、图8、图9、图10中的图示符号说明：

表示在IRI模型下仿真得到的群路径；

表示在IRI和QPS混合建模下仿真得到的群路径；

表示由实测得到的群路径。

(五)具体实施方式

1.组成

本发明一种短波射线追踪技术中的电离层混合建模方法，它是基于国际参考电离层和多层准抛物电离层模型的混合应用来实现一种新的电离层建模方法，该方法兼具两种电离层模型的优点，图1是电离层混合建模方法的流程图。该方法由四个模块组成：它们是输入模块、IRI模块，QPS模块，射线追踪模块。四个模块间的位置、连接关系及信号走向如下述：用户在输入模块中输入待测时间、发射机与接收机地理坐标、天线仰角范围及工作频率等初始参数，由该模块自动算出电波反射中点地理坐标、地心角、大圆距离及方位角，该模块将待测时间和算出的电波反射点坐标传递给IRI模块供生成相应的电离层参数使用，同时将除待测时间和反射点位置外的其他参数传递给射线追踪模块供射线追踪计算使用；IRI模块根据待测时间和电波反射中点位置计算出电波发射区域的E层临界频率、E层最大电子浓度高度、E层半厚、F2层临界频率、F2层最大电子浓度高度、F2层半厚等参数；将由IRI模块算出的电离层参数输入到QPS模块中，建立多层准抛物线形式的等离子体频率分布；射线追踪模块使用工作频率、方位角、天线仰角范围、多层准抛物等离子体频率分布作为初始参数以固定步长逐步计算由发射机到接收机的电波群路径。

所述输入模块是用户与算法的主要交互接口，以软件封装形式内嵌在电离层混合建模方法中，该模块流程图见图2所示。该模块与用户的主要接口包括：时间、发射机坐标、接收机坐标、工作频率、天线仰角范围。它供用户输入相关使用参数，用户在此模块中可设置待测时间、工作频率、发射机坐标、接收机坐标、天线仰角范围等初始参数，该模块可以根据发射机、接收机坐标算出电波反射中点位置、地心角、大圆距离及从发射机到接收机的方位角，以上参数是射线追踪方程组的初始值，该模块一方面将电波反射点位置和待测时间传递给IRI模块供生成相应的电离层参数使用，同时将除去待测时间和反射点位置外的其他参数传递给射线追踪模块供射线追踪计算使用。

所述IRI模块是由国际参考电离层模型来构建电子浓度分布的计算模块，该模块所用到的国际参考电离层(IRI)模型和计算程序都可从美国宇航局网站公开下载获得。该模块流程图见图3所示，它以软件封装的形式存在，模型所用到的公式和CCIR系数都已封装在软件中，使用由上一模块中传来的反射点位置和时间即可生成反射区域各层电离层的关键参数并传递给QPS模块用以建立反射区域等离子体频率分布，同时IRI模块根据IRI模型自动生成反射区域等离子体频率分布并直接传递给射线追踪模块。

所述QPS模块是由多层准抛物模型来构建等离子体频率分布的计算模块，该模块是自行开发的计算程序，该模块流程图见图4所示，该模块以封装的形式内嵌在电离层混合建模方法中，作为混合建模方法的一个模块，QPS模型所用到的公式都封装在该程序模块中。该模块使用上一步中由IRI模块估计的各层电离层关键参数来构建反射区域上空等离子体频率分布，并将各点等离子体频率传递给射线追踪模块。

该模块中设电离层由E层、F1层、F2层组成，且各层均为宁静电离层，不考虑地磁、扰动、梯度等影响，则各层的等离子体频率剖面可表示如下：

$f_{\mathrm{NE}}^2=a_E-b_E{(1-\frac{r_{\mathrm{mE}}}{r})}^2$E层

$f_{\mathrm{Nj}1}^2=a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r})}^2$第一连接层，连接E层和F1层

$f_{\mathrm{NF}1}^2=a_{F1}-b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r})}^2$F1层

$f_{\mathrm{Nj}2}^2=a_{j2}-b_{j2}{(1-\frac{r_{j2}}{r})}^2$第二连接层，连接F1层和F2层

$f_{\mathrm{NF}2}^2=a_{F2}-b_{F2}{(1-\frac{r_{\mathrm{mF}2}}{r})}^2$F2层

上式中，r为空间位置距地心的距离，f_NE为E层等离子体频率，f_oE为E层临界频率，r_bE为E层底高，y_mE为E层的半厚度，r_mE＝r_bE+y_mE为E层最大电子浓度对应高度，b_E＝a_E(r_bE/y_mE)²；f_Nj1为第一连接层等离子体频率，f_oj1为第一连接层临界频率，r_bj1为第一连接层底高，y_mj1为第一连接层的半厚度，r_j1＝r_bj1+y_mj1为第一连接层最大电子浓度对应高度，b_j1＝a_j1(r_bj1/y_mj1)²；f_NF1为F1层等离子体频率，f_oF1为F1层临界频率，r_bF1为F1层底高，y_mF1为F1层的半厚度，r_mF1＝r_bF1+y_mF1为F1层最大电子浓度对应高度，b_F1＝a_F1(r_bF1/y_mF1)²；f_Nj2为第二连接层等离子体频率，f_oj2为第二连接层临界频率，r_bj2为第二连接层底高，y_mj2为第二连接层的半厚度，r_j2＝r_bj2+y_mj2为第二连接层最大电子浓度对应高度，b_j2＝a_j2(r_bj2/y_mj2)²；f_NF2为F2层等离子体频率，f_oF2为F2层临界频率，r_bF2为F2层底高，y_mF2为F2层的半厚度，r_mF2＝r_bF2+y_mF2为F2层最大电子浓度对应高度，b_F2＝a_F2(r_bF2/y_mF2)²。

对各层等离子体频率分布求导，可以得到各层等离子体频率随高度的变化率

$\frac{{\mathrm{df}}_{\mathrm{NE}}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mE}}{b}_E}{r^2}(1-\frac{r_{\mathrm{mE}}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}1}^2}{\mathrm{dr}}=-\frac{{2r}_{j1}{b}_{j1}}{r^2}(1-\frac{r_{j1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}1}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{r^2}(1-\frac{r_{\mathrm{mF}1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}2}^2}{\mathrm{dr}}=-\frac{{2r}_{j2}{b}_{j2}}{r^2}(1-\frac{r_{j2}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}2}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}2}{b}_{F2}}{r^2}(1-\frac{r_{\mathrm{mF}2}}{r})$

将步骤三中得到的E层临界频率r_oE、E层最大电子浓度高度r_mE、E层半厚y_mE、F2层临界频率f_oF2、F2层最大电子浓度高度r_mF2、F2层半厚y_mF2及连接层参数代入上式中，得到电波反射区域等离子体频率的梯度分布。

第一连接层和第二连接层的参数可由下式求得，假设第一连接层与E层的交点在E层的最大电子浓度处，则根据两层的等离子体频率及其梯度分别相等，有

a_j1＝a_E

r_j1＝r_mE

设第一连接层与F1层的交点在r＝r_c1处，则在r_c1处两层的等离子体频率及其梯度分别相等，有

$-\frac{{2r}_{j1}{b}_{j1}}{{r_{c1}}^2}(1-\frac{r_{j1}}{r_{c1}})=\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{{r_{c1}}^2}(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})$

$a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r_{c1}})}^2=a_{F1}+b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})}^2$

联立上两式解得

$r_{c1}=\frac{r_{\mathrm{mF}1}{b}_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}{a_{F1}-a_E+b_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}$

$b_{j1}=\frac{r_{\mathrm{mF}1}{r}_{c1}{b}_{F1}-r_{\mathrm{mF}1}^2{b}_{F1}}{r_{\mathrm{mE}}^2-r_{\mathrm{mE}}{r}_{c1}}$

同理得到

a_j2＝a_F1，r_j2＝r_mF1，$r_{c2}=\frac{r_{\mathrm{mF}2}{b}_{F2}(1-r_{\mathrm{mF}2}/r_{j2})}{a_{j2}-a_{F2}+b_{F2}(1-r_{\mathrm{mF}2}/r_{j2})},$$b_{j2}=\frac{r_{\mathrm{mF}2}^2{b}_{F2}-r_{\mathrm{mF}2}{r}_{c2}{b}_{F2}}{r_{j2}^2-r_{j2}{r}_{c2}}$

其中，r_c1为第一连接层与E层交点处的高度，r_c2为第二连接层与F2层交点处的高度。

所述射线追踪模块是电波群路径的计算模块，该模块是自行开发的计算程序，所用到的射线追踪方程组都封装在该模块中。该模块流程图参见图5所示，它以输入模块中的工作频率、天线仰角范围、方位角等参数作为射线追踪方程组的初值，根据反射区域上空等离子体频率分布由发射机开始迭代计算电波传播位置及方向，最后取距接收机最近的一点作为从发射机到接收机的传播轨迹，将传播轨迹描迹显示并得到电波群路径。

该模块从球坐标系下以电波传播群路径P′为自变量的射线追踪方程组出发

方程组用于计算射线路径上点的坐标及该点处的波矢量。式中，P′为群路径；r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；ω为圆频率，c为光速，n为折射指数，与电子浓度分布有关。

$n^2=1-\frac{2X(1-\mathrm{iZ}-X)}{2(1-\mathrm{iZ})(1-\mathrm{iZ}-X)-Y_T^2\pm \sqrt{T_T^4+4Y_L^2{(1-\mathrm{iZ}-X)}^2}}$

式中，Y＝f_H/f，Y_T＝Ysinψ，Y_L＝Ycosψ，Z＝υ/2πf，f_H＝eB/2πm。其中，f_N为等离子体频率，是一个与电子浓度有关的物理量，由上一步的QPS模型构建；f为电波频率；f_H为磁旋频率；ψ为波矢量与地磁场间的夹角；υ为电子碰撞频率；N为电子浓度；e为电子电量；m为电子质量；ε₀为真空中介电常数；B为地磁场强度。“+”对应寻常波(0波)，“-”对应非常波(X波)。

因碰撞只造成能量的吸收，不影响电波传播方向，故在射线追踪中无须考虑碰撞效应，取Z＝0，在不考虑地磁场影响时，取Y＝0，则折射指数简化为

n²＝1-X

有

$\frac{\mathrm{dX}}{\mathrm{dr}}=\frac{1}{f^2}\frac{{\mathrm{df}}_N^2}{\mathrm{dr}}$

若在射线追踪过程中不考虑电离层梯度变化的影响，则为0。

设电波发射仰角为β，方位角为α，工作圆频率为ω，发射机纬度、经度分别为λ₀和则射线的初始位置为：

r＝r₀

θ＝π/2-λ₀

$k_r=\left|\overrightarrow{k}\right|\cos (\pi /2-\beta )=\omega \sin \beta /c$

$k_{\theta }=-\left|\overrightarrow{k}\right|\cos \beta \cos \alpha =-\omega \cos \beta \cos \alpha /c$

式中，r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；r₀为地球半径，根据r、θ、k_r、k_θ、六个参数的初始值，即可开始射线轨迹的计算。

在射线路径上，以群路径为自变量，射线的坐标矢量和波矢量随着群路径的变化而变化，利用常微分方程数值解法Runge-Kutta法以固定步长逐点求解该射线追踪方程组，得到射线路径上各点的r、θ、k_r、k_θ、其中每前进一点，都要计算一次该点的等离子体参量变化率将这些点连接起来，即得到一条球坐标系中的射线轨迹，在终点处，可得到电波群路径、地面距离、接收点经纬度等信息。

本发明具有较高实用性，根据两种常用电离层模型，开展混合建模工作，弥补了各自的不足，使用简单方便，只需要一台可正常运行Visual C++开发软件的普通个人计算机即可实现(如：P4 1.6GHz主频CPU、1G内存、Windows XP操作系统)，不需要特定设备和环境条件。

2.步骤

本发明一种短波射线追踪技术中的电离层混合建模方法，该方法具体步骤如下：

步骤一：在输入模块中首先输入发射机、接收机的地理坐标、待测时间、工作频率、天线仰角范围等初始参数。发射机位于青岛，坐标为东经120.30°，北纬36.10°，接收机位于长春，坐标为东经125.28°，北纬43.84°，待测时间设为2008年10月28日、10月31日和11月1日的白天，三天选用不同的工作频率，分别为9.322MHz、8MHz和10MHz，天线仰角范围20°-60°；

步骤二：由发射机、接收机地理坐标按以下公式计算出反射中点M地理坐标(东经122.79°，北纬39.97°)、地心角d(0.15rad)、大圆距离D(959km)和方位角为α(24.67°)；

$M_{\mathrm{Lat}}=\frac{T_{\mathrm{Lat}}+R_{\mathrm{Lat}}}{2},$$M_{\mathrm{Lon}}=\frac{T_{\mathrm{Lon}}+R_{\mathrm{Lon}}}{2},$$\mathrm{\alpha =arccos}\frac{\sin R_{\mathrm{Lat}}-{\sin T}_{\mathrm{Lat}}\cos d}{\cos T_{\mathrm{Lat}}\sin d}$

d＝arccos(sinT_Lat sinR_Lat+cosT_Lat cosR_Lat cos(T_Lon-R_Lon))，D＝r₀·d

上式中，M_Lat为反射中点纬度，M_Lon为反射中点经度，T_Lat为发射机纬度，T_Lon为发射机经度，R_Lat为接收机纬度，R_Lon为接收机经度，α为从发射机到接收机的方位角，以正北为基准东偏为正，d为地心角，是收发两点与地心连线间的夹角，以弧度为单位，r₀为地球半径，为6370km，D为大圆距离，是收发两点间的地面距离，单位为km。

步骤三：在IRI模块中输入待测时间和反射点坐标，由IRI模型自动生成待测时间电波反射点处的电离层参数值和反射区域的等离子体频率分布，计算所需的CCIR系数和算法已封装在IRI模型中，直接使用即可，该模型可从美国宇航局网站下载获得。生成的电离层参数包括E层临界频率f_oE、E层最大电子浓度高度r_mE、E层半厚y_mE、F2层临界频率f_oF2、F2层最大电子浓度高度r_mF2、F2层半厚y_mF2。表1-表3是由IRI模型推算出的青岛至长春短波电路于2008年10月28日、10月31日和11月1日白天不同时刻所对应的电离层参数，因未出现F1层，所以表中没有F1层的相关参数。

表1.10月28日电离层参数

  探测时间  2008-10-28  f_oE  (MHz)  h_mE  (km)  y_mE  (km)  f_oF2  (MHz)  h_mF2  (km)  y_mF2  (km)  8:15  2.4  110  30  7.7  282  173  8:45  2.6  110  30  8.2  290  90  9:15  2.7  110  30  8.5  299  184  9:45  2.8  110  30  8.8  308  187  10:15  2.9  110  30  9.1  316  188  10:45  2.9  110  30  9.3  323  189  11:15  2.9  110  30  9.4  326  190  11:45  3.0  110  30  9.6  326  189  12:15  2.9  110  30  9.8  324  190  12:45  2.9  110  30  10.1  319  189  13:15  2.9  110  30  10.3  314  188  13:45  2.8  110  30  10.5  308  186  14:15  2.7  110  30  10.6  303  184  14:45  2.6  110  30  10.6  299  179  15:15  2.4  110  30  10.5  296  173  15:45  2.2  110  30  10.4  294  164

表2.10月31日电离层参数

  探测时间  2008-10-31  f_oE  (MHz)  h_mE  (km)  y_mE  (km)  f_oF2  (MHz)  h_mF2  (km)  y_mF2  (km)  8:15  2.4  110  30  7.8  283  173  8:45  2.6  110  30  8.2  290  179

  探测时间  2008-10-31  f_oE  (MHz)  h_mE  (km)  y_mE  (km)  f_oF2  (MHz)  h_mF2  (km)  y_mF2  (km)  9:15  2.7  110  30  8.6  299  184  9:45  2.8  110  30  8.9  308  186  10:15  2.8  110  30  9.1  316  188  10:45  2.9  110  30  9.3  323  189  11:15  2.9  110  30  9.5  326  190  11:45  2.9  110  30  9.6  327  190  12:15  2.9  110  30  9.8  324  190  12:45  2.9  110  30  10.0  320  180  13:15  2.8  110  30  10.3  314  188  13:45  2.8  110  30  10.5  308  186  14:15  2.7  110  30  10.3  303  183  14:45  2.5  110  30  10.6  299  179  15:15  2.4  110  30  10.6  296  172  15:45  2.1  110  30  10.4  295  163

表3.11月1日电离层参数

  探测时间  2008-11-1  f_oE  (MHz)  h_mE  (km)  y_mE  (km)  f_oF2  (MHz)  h_mF2  (km)  y_mF2  (km)  9:15  2.7  110  30  8.2  292  184  9:45  2.8  110  30  8.4  299  186

  探测时间  2008-11-1  f_oE  (MHz)  h_mE  (km)  y_mE  (km)  f_oF2  (MHz)  h_mF2  (km)  y_mF2  (km)  10:15  2.8  110  30  8.6  305  188  10:45  2.9  110  30  8.7  310  189  11:15  2.9  110  30  8.8  312  190  11:45  2.9  110  30  9.0  313  190  12:15  2.9  110  30  9.2  312  190  12:45  2.9  110  30  9.4  308  189  13:15  2.8  110  30  9.6  302  188  13:45  2.8  110  30  9.8  296  186  14:15  2.7  110  30  9.8  290  183  14:45  2.5  110  30  9.8  285  179  15:15  2.4  110  30  9.6  281  172  15:45  2.1  110  30  9.4  278  162  16:15  1.8  110  30  9.2  277  150  16:45  1.4  110  30  8.9  278  135

步骤四：将步骤三中由IRI模型算出的各时刻电离层参数传递给QPS模块作为输入参量代入QPS公式建立反射空间的多层准抛物等离子体频率分布和等离子体频率梯度分布；

设电离层由E层、F1层、F2层组成，且各层均为宁静电离层，不考虑地磁、扰动、梯度等影响，则各层的等离子体频率剖面可表示如下：

$f_{\mathrm{NE}}^2=a_E-b_E{(1-\frac{r_{\mathrm{mE}}}{r})}^2$E层

$f_{\mathrm{Nj}1}^2=a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r})}^2$第一连接层，连接E层和F1层

$f_{\mathrm{NF}1}^2=a_{F1}-b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r})}^2$F1层

$f_{\mathrm{Nj}2}^2=a_{j2}-b_{j2}{(1-\frac{r_{j2}}{r})}^2$第二连接层，连接F1层和F2层

$f_{\mathrm{NF}2}^2=a_{F2}-b_{F2}{(1-\frac{r_{\mathrm{mF}2}}{r})}^2$F2层

上式中各参数物理意义可参见背景技术中的说明。

对各层等离子体频率分布求导，可以得到各层等离子体频率随高度的变化率

$\frac{{\mathrm{df}}_{\mathrm{NE}}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mE}}{b}_E}{r^2}(1-\frac{r_{\mathrm{mE}}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}1}^2}{\mathrm{dr}}=-\frac{{2r}_{j1}{b}_{j1}}{r^2}(1-\frac{r_{j1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}1}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{r^2}(1-\frac{r_{\mathrm{mF}1}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{Nj}2}^2}{\mathrm{dr}}=-\frac{{2r}_{j2}{b}_{j2}}{r^2}(1-\frac{r_{j2}}{r})$

$\frac{{\mathrm{df}}_{\mathrm{NF}2}^2}{\mathrm{dr}}=-\frac{{2r}_{\mathrm{mF}2}{b}_{F2}}{r^2}(1-\frac{r_{\mathrm{mF}2}}{r})$

将步骤三中得到的E层临界频率f_oE、E层最大电子浓度高度r_mE、E层半厚y_mE、F2层临界频率f_oF2、F2层最大电子浓度高度r_mF2、F2层半厚y_mF2及连接层参数代入上式中，得到电波反射区域等离子体频率的梯度分布。

第一连接层和第二连接层的参数可由下式求得，假设第一连接层与E层的交点在E层的最大电子浓度处，则根据两层的等离子体频率及其梯度分别相等，有

a_j1＝a_E

r_j1＝r_mE

设第一连接层与F1层的交点在r＝r_c1处，则在r_c1处两层的等离子体频率及其梯度分别相等，有

$-\frac{{2r}_{j1}{b}_{j1}}{{r_{c1}}^2}(1-\frac{r_{j1}}{r_{c1}})=\frac{{2r}_{\mathrm{mF}1}{b}_{F1}}{{r_{c1}}^2}(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})$

$a_{j1}-b_{j1}{(1-\frac{r_{j1}}{r_{c1}})}^2=a_{F1}+b_{F1}{(1-\frac{r_{\mathrm{mF}1}}{r_{c1}})}^2$

由上两式解得

$r_{c1}=\frac{r_{\mathrm{mF}1}{b}_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}{a_{F1}-a_E+b_{F1}(1-r_{\mathrm{mF}1}/r_{\mathrm{mE}})}$

$b_{j1}=\frac{r_{\mathrm{mF}1}{r}_{c1}{b}_{F1}-r_{\mathrm{mF}1}^2{b}_{F1}}{r_{\mathrm{mE}}^2-r_{\mathrm{mE}}{r}_{c1}}$

同理

a_j2＝a_F1，r_j2＝r_mF1，$r_{c2}=\frac{r_{\mathrm{mF}2}{b}_{F2}(1-r_{\mathrm{mF}2}/r_{j2})}{a_{j2}-a_{F2}+b_{F2}(1-r_{\mathrm{mF}2}/r_{j2})},$$b_{j2}=\frac{r_{\mathrm{mF}2}^2{b}_{F2}-r_{\mathrm{mF}2}{r}_{c2}{b}_{F2}}{r_{j2}^2-r_{j2}{r}_{c2}}$

其中，r_c1为第一连接层与E层交点处的高度，r_c2为第二连接层与F2层交点处的高度。

步骤五：利用步骤四中建立的反射空间多层准抛物等离子体频率分布以群路径为自变量从发射机位置开始使用Runge-Kutta法以固定步长迭代求解射线追踪方程组，每次求解方程组得到当前电波轨迹位置和传播方向，将以上参数代入方程组中再次计算，得到下一步电波轨迹位置和传播方向，如此迭代，至电波到达接收机为止，最终得到待测时间由发射机到接收机的电波群路径；

从球坐标系下以电波传播群路径P′为自变量的射线追踪方程组出发

方程组用于计算射线路径上点的坐标及该点处的波矢量。式中，P′为群路径；r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；ω为圆频率，c为光速，n为折射指数，与等离子体频率分布有关，在不考虑碰撞和地磁场时，折射指数可表示为

n²＝1-X

式中，f_N为等离子体频率，f为电波频率。

$\frac{\mathrm{dX}}{\mathrm{dr}}=\frac{1}{f^2}\frac{{\mathrm{df}}_N^2}{\mathrm{dr}}$

若在射线追踪过程中不考虑电离层梯度变化的影响，则为0。

由步骤一和步骤二得到的参量，电波发射仰角β＝20°～60°、方位角α＝24.67°、三天中工作圆频率分别对应ω＝2π×9.322/8/10MHz、发射机纬度λ₀＝36.10°和经度可以确定射线的初始位置为：

r＝r₀

θ＝π/2-λ₀

$k_r=\left|\overrightarrow{k}\right|\cos (\pi /2-\beta )=\omega \sin \beta /c$

$k_{\theta }=-\left|\overrightarrow{k}\right|\cos \beta \cos \alpha =-\omega \cos \beta \cos \alpha /c$

式中，r、θ、是射线路径上点在球坐标系中的坐标；k_r、k_θ、是波矢量在球坐标系中的三个分量；r₀为地球半径，根据r、θ、k_r、k_θ、六个参数的初始值，即可开始射线轨迹的计算。

采用Runge-Kutta法求解射线追踪方程组，群路径P′取固定步长1km，每次求解方程组得到当前电波位置r、θ、和波矢量k_r、k_θ、由电波位置计算等离子体浓度fN和折射指数n，再将解得的r、θ、k_r、k_θ、f_N、n代入方程组中求得下一步电波位置和波矢量，如此迭代，直到电波到达地面，电波仰角以0.5°为步长从20°扫描到60°，每个仰角对应不同的电波轨迹，分别选取由E层反射和F2层反射的落点位置距接收机最接近的一条射线，得到E层群路径和F2层群路径。仿真结果见表4，表中“-”代表在该时刻对应频点没有仿真值。

表4.混合建模仿真群路径

  探测时间  2008-10-28  E层  2008-10-28  F2层  2008-10-31  E层  2008-11-1  F2层  8:15  -  1115  996  -  8:45  -  1121  988  -  9:15  -  1109  986  1134  9:45  994  1110  983  1136  10:15  990  1119  983  1135  10:45  990  1118  983  1137  11:15  990  1107  983  1136  11:45  996  1107  983  1130  12:15  990  1104  983  1121  12:45  990  1098  983  1110  13:15  995  1074  999  1101  13:45  995  1074  999  1092  14:15  -  1080  985  1086

  探测时间  2008-10-28  E层  2008-10-28  F2层  2008-10-31  E层  2008-11-1  F2层  14:45  -  1067  988  1081  15:15  -  1056  994  1081  15:45  -  1062  -  1086

步骤六：用步骤三中IRI模型软件直接计算出的(所需CCIR系数和算法已内嵌到模型中)待测时间反射区域IRI模型下的等离子体频率替换步骤四中由QPS模块算出的等离子体频率fN，然后重复步骤五的射线追踪过程，得到IRI模型下的电波群路径，见表5。

表5.IRI模型仿真群路径

  探测时间  2008-10-28  E层  2008-10-28  F2层  2008-10-31  E层  2008-11-1  F2层  8:15  992  1140  990  1166  8:45  992  1132  990  1166  9:15  990  1144  989  1154  9:45  990  1144  989  1154  10:15  990  1167  989  1170  10:45  991  1168  989  1170  11:15  990  1177  989  1180  11:45  990  1177  989  1180  12:15  988  1157  989  1161  12:45  990  1162  989  1161  13:15  1012  1126  989  1124

  探测时间  2008-10-28  E层  2008-10-28  F2层  2008-10-31  E层  2008-11-1  F2层  13:45  1012  1126  989  1103  14:15  990  1092  989  1088  14:45  991  1093  989  1088  15:15  994  1085  990  1076  15:45  994  1085  990  1076

步骤七：按照步骤一中的仿真条件，利用中国电波传播研究所的实验设备和观测站点进行电波斜向探测试验，得到实测电波群路径；

为了考察混合模型仿真效果的有效性，以电波群路径作为测量参量进行实测验证，如图6所示，按照仿真条件在青岛-长春间建立短波电离层斜向探测电路，发射机位于青岛，坐标为东经120.30°，北纬36.10°，接收机位于长春，坐标为东经125.28°，北纬43.84°；试验时间为2008年10月28日、10月31日和11月1日，从00:15开始，每隔30分钟探测一次，发射信号为脉冲周期信号，脉冲重复周期1s，以扫频方式工作，载波范围从2M到30M，频点间隔50k，测量参量为电波群路径，接收机群路径最小分辨率5km。收发机间地面距离为958.97km，方位角为24.68°。由实测结果发现，试验的三天中未出现F1层，且在夜间F2层消失，实测电波群路径见表6所示，其中10月28日实测电波群路径对应频率9.322MHz，10月31日实测电波群路径对应频率8MHz，11月1日实测电波群路径对应频率10MHz。表中“-”代表在探测时刻没有得到实测值。

表6.实测群路径

  探测时间  2008-10-28  E层  2008-10-28  F2层  2008-10-31  E层  2008-11-1  F2层  8:15  981  1068  974  -  8:45  987  1068  974  -  9:15  981  1068  974  1074  9:45  -  -  974  1068

  探测时间  2008-10-28  E层  2008-10-28  F2层  2008-10-31  E层  2008-11-1  F2层  10:15  974  1074  974  1106  10:45  981  1068  974  1074  11:15  974  1224  974  1081  11:45  974  1224  974  1074  12:15  974  1087  974  1093  12:45  974  1081  974  1093  13:15  -  -  974  1087  13:45  974  1087  974  1068  14:15  981  1074  974  1068  14:45  981  1068  974  1068  15:15  974  1118  981  -  15:45  974  1118  974  -

步骤八：比较由步骤五得到的混合模型下仿真群路径和由步骤六得到的IRI模型下仿真群路径与实测值间均方距离的大小，如果混合模型下仿真群路径与实测值间的均方距离小于IRI模型仿真群路径与实测值间均方距离，说明采用混合模型进行射线追踪的精确程度更高。

图7是10月28日工作频率取9.322MHz下E层群路径仿真值与实测值随探测时间的变化曲线，图8是相同条件下F2层群路径仿真值与实测值随探测时间的变化曲线，从曲线上可以发现混合模型的仿真值与实测值更为接近。除图中所列探测时间，在其他探测时间未能接收到9.322MHz工作频率的反射信号，故其他时间在图中未列出。

图9是10月31日工作频率为8MHz下的E层群路径仿真值与实测值的对比曲线，图10是11月1日工作频率为10MHz下的F2层群路径仿真值与实测值对比曲线。从图中可以发现，由IRI-QPS混合模型得到的仿真值与实测值更加接近。

为了描述模型仿真值与实测值间的差异程度，定义仿真值与实测值间的均方距离S为

$S=\sqrt{\frac{\underset{i=1}{\overset{N}{\Sigma }}{(x_i-y_i)}^2}{N-1}}$

式中，x_i是i时刻对应的仿真群路径，y_i是i时刻对应的实测群路径，i取从1到N的整数，对应N个不同的试验时间。S越小，说明仿真值与实测值距离越小，总体上越接近实测值。

$S\%=S/\overline{P^{\prime }}\times 100\%$

上式中，S％为均方距离差异率，S为均方距离，为实测群路径均值。

表7.群路径仿真值与实测值间的均方距离

表8.群路径均方距离差异率

表7是仿真值与实测值的均方距离，可以发现，IRI模型和混合模型在E层反射电波的群路径与实测值较接近，在F2层反射的群路径与实测值相差较大，可能与IRI模型中国地区实测数据不足有关，也可能与对电离层进行简化，忽略梯度、扰动、地磁影响等因素有关，说明IRI模型对应的宁静电离层只是实际电离层的一种近似，与实际电离层相比仍有较大差距，特别是在F2层。将均方距离除以实测均值可以得到均方距离差异率，反映了仿真值与实测值的总体差异程度，见表8。10月28日在E层，IRI模型的均方距离与实测均值之比为1.82％，混合模型的均方距离与实测均值相比达到1.88％，在F2层，IRI模型的均方距离与实测均值相比达到5.90％，混合模型的均方距离与实测均值相比为4.83％。由表8，在F2层IRI-QPS混合模型的预测值与实测值更为接近，可以看出在E层IRI-QPS混合模型与IRI模型的估计值比较接近，而在F2层使用IRI-QPS混合模型估计值的准确程度要优于仅使用IRI模型的估计值。

由表8发现，IRI-QPS混合建模在白天对F2层的射线追踪仿真结果比仅使用IRI模型射线追踪的仿真结果更精确，对E层的仿真结果精确程度两者相当，同时混合建模的方法还避免了QPS模型输入参数的盲目性，扩展了QPS模型的使用范围。