CN114741839A - 一种分析甚低频电磁波在地-电离层中传播的fdtd方法 - Google Patents

Abstract

本发明公开了一种分析甚低频电磁波在地‑电离层中传播的FDTD方法。本发明首先模拟出地‑电离层电磁环境，再通过两次坐标系转换简化运算，接着用SO‑FDTD算法计算具体传播情况，最后控制各地‑电离层参数分析总结传播特性。本发明从单纯的辐射场计算扩展到传播特性，更进一步地研究了甚低频电磁波在地‑电离层波导中的传播，对于超远程导航、潜艇通信、天气预测等领域的研究有着重要参考意义。本发明可以解决和分析有关甚低频电磁波在地‑电离层波导中的传播特性问题。

Description

一种分析甚低频电磁波在地-电离层中传播的FDTD方法

技术领域

本发明属于电磁场数值计算技术领域，具体涉及一种电磁波在地-电离层中传播的FDTD方法。

背景技术

甚低频电磁波是频率在3kHz～30kHz的电磁波，其具有传播距离远、传播损耗小、幅度和相位稳定、渗透性强等优点，因此广泛应用于超远程导航、潜艇通信、天气预测等领域。地-电离层对甚低频电磁波具有良好的反射特性，再加上甚低频电磁波波长与地-电离层间的距离相近，所以甚低频电磁波在地-电离层之间的传播可类似于在波导中的传播，因此也称为地-电离层波导传播模式。但实际上地面，电离层以及地面和电离层之间的电磁时空变化异常复杂，所以甚低频电磁波在地-电离层中传播时会表现出非常复杂的特性，因此实现更准确的预测对上述应用领域有着十分重要的意义。

在计算甚低频电磁波在地-电离层中传播时，考虑到需要进一步提高预测的准确性，所以解析法不再适用。而在数值方法中，通常使用FDTD(Finite Difference TimeDomain，时域有限差分算法)来模拟甚低频电磁波在地-电离层波导中的传播。而现有的成果中，大多只通过FDTD计算甚低频电磁波在地-电离层波导中传播时的场分布情况，对于通过FDTD分析甚低频电磁波在地-电离层波导中传播特性，尤其是各地-电离层参数对传播带来的影响，目前的研究还较少，而这些却对于提高预测速度和精度十分重要。因此，非常有必要设计一种能够分析甚低频电磁波在地-电离层波导中传播特性的算法。

发明内容

为了克服现有技术的不足，本发明提供了一种分析甚低频电磁波在地-电离层中传播的FDTD方法。本发明首先模拟出地-电离层电磁环境，再通过两次坐标系转换简化运算，接着用SO-FDTD算法计算具体传播情况，最后控制各地-电离层参数分析总结传播特性。本发明从单纯的辐射场计算扩展到传播特性，更进一步地研究了甚低频电磁波在地-电离层波导中的传播，对于超远程导航、潜艇通信、天气预测等领域的研究有着重要参考意义。本发明可以解决和分析有关甚低频电磁波在地-电离层波导中的传播特性问题。

本发明解决其技术问题所采用的技术方案包括如下步骤：

步骤1：以地心为极点，地心与电磁波发射点连线延长所得到的射线为极轴，建立极坐标系，将三维球坐标系下的问题转换到二维极坐标系；

从电离层IRI2016模型和NRLMSISE-00大气模型中爬取数据，建立甚低频电磁波在地-电离层波导中传播路径上粒子的密度和温度的实时数据集；

步骤2：由步骤1得到的实时数据集，计算参数——电子密度N_e和碰撞频率v，实时模拟传播路径上的地-电离层环境；

步骤3：将极坐标系下电子密度N_e和碰撞频率v经过双线性插值算法处理后转换到直角坐标系；

步骤4：根据步骤3转换到直角坐标系的电子密度N_e和碰撞频率v，利用SO-FDTD算法计算甚低频电磁波在传播路径上的场分布情况，将电磁波接收点一天内不同时间的场强变换情况与甚低频台站上实测数据进行比对和验证；

步骤5：控制各地-电离层参数变化，对计算得到的场分布情况进行分析，得到甚低频电磁波在地-电离层波导中传播特性。

进一步地，所述步骤1具体包括以下过程：

步骤1-1：利用爬虫爬取电离层IRI2016模型传播路径上的电子密度N_e和电子温度T_e，爬取NRLMSISE-00大气模型传播路径上的氧气分子密度

氧原子密度N_o和氮气分子密度

；

步骤1-2：建立甚低频电磁波在地-电离层波导中传播路径上粒子密度和温度的实时数据集；

以地心为极点，地心与电磁波发射点连线延长所得到的射线为极轴，建立极坐标系，将三维球坐标系下的问题转换到二维极坐标系；对于二维极坐标系内任一点P(ρ,θ)，极径ρ对应该点到地心的距离，极角θ代表到该点与地心连线和电磁波发射点与地心连线之间的夹角；从电磁波发射点开始直到电磁波接收点结束，每间隔一个角度，就通过爬虫爬取一次对应位置从地面到低电离层中每隔一定高度的粒子密度和温度的实时数据，建立实时数据集；

进一步地，所述步骤2具体包括以下过程：

步骤2-1：模拟地面环境；

将地面近似为均匀的电磁媒质，具体如表1所示：

表1地面近似电磁媒质表

	介电常数ε	磁导率μ
			平均陆地状态	10	0.003
平均海洋状态	80	4

步骤2-2：模拟电离层环境；

电离层的两个参数——电子密度N_e和碰撞频率v，其中电子密度N_e直接从电离层IRI2016模型中获取，碰撞频率v则通过经验公式(1)计算，具体如下：

其中，

式中，V_e,i代表电子和离子的碰撞频率，

代表电子和氧气分子的碰撞频率，V_e,O代表电子和氧原子的碰撞频率，

代表电子和氮气分子的碰撞频率，N_e代表电子密度，T_e代表电子温度，

代表氧气分子密度，N_O代表氧原子密度，

代表氮气分子密度；

步骤2-3：模拟地-电离层之间环境；

对于地面和电离层之间来说，采用经验指数模型进行表示，具体如下：

v(z)＝1.82×10¹¹e^-0.15z (3)

N(z)＝1.43×10⁷e^-0.15He^{[(β-0.15)(z-H)]} (4)

式中，v(z)代表高度为z时的碰撞频率，N(z)代表高度为z时的电子密度，β代表梯度系数(km^-1)，H代表电离层参考高度(km)；

其中经验指数模型内中低纬度地区β和H的推荐值见表2，其中，f为工作频率；

表2中低纬度地区β和H取值

	夏季	冬季
			白天	β＝0.3,H＝70	β＝0.3,H＝72
晚上	β＝0.0077f+0.31,H＝87	β＝0.0077f+0.31,H＝87

进一步地，所述步骤3具体包括以下过程：

对于FDTD来说，极坐标系和直角坐标系中的Yee元胞不是一一对应的，所以不能直接转换，需要通过插值进行处理；采用双线性插值算法将极坐标系下电子密度N_e和碰撞频率v转换到直角坐标系；双线性插值算法的核心思想是在x方向和y方向分别进行一次插值，即直角坐标系中每一点的值是由邻近四个极坐标系中的点插值得到的；对于点P(x,y)的值，已知Q₁₁(x₁,y₁)，Q₁₂(x₂,y₂)，Q₁₃(x₃,y₃)和Q₁₄(x₄,y₄)四个点的值，插值公式为

进一步地，所述步骤4具体包括以下过程：

步骤4-1：利用SO-FDTD计算传播路径上的场分布；

线性各向同性介质麦克斯韦旋度方程为：

式中，G代表磁场强度，D代表电通量密度，σ代表电导率，E代表电场强度，B代表磁通量密度，σ_m代表等效磁导率；

电离层的介电常数随频率变化，所以电离层也属于色散介质，对于此类介质有：

B＝μG (7)

D(ω)＝ε(ω)E(ω) (8)

式中，μ代表磁导率，ε(ω)代表介电常数，ω代表角频率；

由式(6)进行FDTD离散得到：

式中，n代表迭代步数；

将式(8)由频域转换到时域，对于x分量，得：

式中，ε₀代表真空介电常数，ε_r代表相对介电常数，E_x(t)代表E(t)的x分量，D_x(t)代表D(t)的x分量。

电离层是Drude介质，所以有：

式中，p_l和q_l是多项式系数，M和N代表多项式总项数。

即有：

设函数：

式(13)左端平均值近似，右端中心差分近似得：

引进离散移位算子z_l，定义为：

z_lfⁿ＝fⁿ⁺¹ (15)

联立式(14)和式(15)，得：

联立式(13)和式(16)，得：

式(17)代入式(12)，得：

两边同乘(z_l+1)^N，得：

因为移位算子的作用使得

所以得到步进公式为：

式中，a_l和b_l由有理分式式(11)中的系数p₀,p₁,...,p_N和q₀,q₁,...,q_M表示；

当M＝N＝2时，整理得：

其中，

电离层是等离子体，所以有：

其中ω_p是等离子体频率，v_c是电子碰撞频率，ε_∞＝1；

将SO-FDTD的步进计算步骤归结如下：(1)由E→G，用式(9)第一式计算；(2)由G→D，用式(9)第二式计算；(3)由D→E，用式(20)计算；由此能够在时间空间轴上逐步推进求解甚低频电磁波在地-电离层中传播产生的辐射场分布情况；

步骤4-2：比对计算数据与实测数据，验证方法可行性。

进一步地，所述步骤5具体包括以下过程：

步骤5-1：设置直角坐标值为变量，求得电磁场分布的矩阵数值；绘图得到传播路径上各个电磁场分量的振幅值和相位值分布，以此为基础分析甚低频电磁波在地-电离层波导中的传播特性；

步骤5-2：考虑各地-电离层参数对传播过程产生的影响，设置各地-电离层参数为变量，通过对比各地-电离层参数不同时的电磁场分量的振幅和相位值分布情况，得出甚低频电磁波在地-电离层波导中的传播特性结论。

本发明的有益效果如下：

本发明与现有的甚低频在地-电离层波导中传播计算方法相比，将三维球坐标系中的问题转换到二维极坐标系，再转换到二维直角坐标系，大大减小了计算量，计算效率更高。现有方法多对于甚低频电磁波在地-电离层波导中传播产生的辐射场分布进行仿真计算，本发明将辐射场的计算拓展为甚低频在地-电离层波导中的传播特性研究，尤其通过控制各地-电离层参数，探究了地-电离层参数对传播的影响。本发明思路新颖，具有创新性，对于甚低频电磁波的应用方面具有重要意义。

附图说明

图1是本发明方法流程图。

图2是本发明方法的问题模型示意图。

图3是本发明方法爬取数据的爬虫示意图。

图4是本发明方法实时模拟地-电离层电磁环境方法示意图。

图5是本发明方法坐标系转换示意图。

图6是本发明方法三维球坐标系到二维极坐标系转换方法示意图。

图7是本发明方法二维极坐标系到二维直角坐标系转换方法示意图。

图8是本发明方法的SO-FDTD运行原理示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

本发明解决的技术问题是：为提高现有技术中计算甚低频电磁波在地-电离层波导中传播产生的辐射场时空分布的速度和精度，本发明设计一种分析甚低频电磁波在地-电离层中传播的FDTD方法。

本发明基于Visual studio编译平台，使用Fortran语言，实时模拟地-电离层电磁环境，首先从电离层IRI2016模型和NRLMSISE-00大气模型中获取各种粒子的密度和温度，进而通过计算可以获取对传播起到直接作用的两个电离层参数——电子密度N_e和碰撞频率v，这样就实时模拟出了地-电离层环境。接着利用SO-FDTD算法，计算出甚低频电磁波在地-电离层波导中产生的辐射场分布情况，以其中不同时间内接收点的场强变化情况与甚低频台站获得的实际数据进行比对，验证算法的可行性。最后通过控制各地-电离层参数变化来具体研究其对甚低频电磁波在地-电离层波导中传播的影响，总结甚低频电磁波在地-电离层波导中传播特性。

一种分析甚低频电磁波在地-电离层中传播的FDTD方法，包括如下步骤：

步骤1：以地心为极点，地心与电磁波发射点连线延长所得到的射线为极轴，建立极坐标系，将三维球坐标系下的问题转换到二维极坐标系；

利用Python编程语言中的Requests库实现从电离层IRI2016模型(https://ccmc.gsfc.nasa.gov/modelweb/models/iri2016_vitmo.php)和NRLMSISE-00大气模型(https://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php)大规模爬取数据，建立甚低频电磁波在地-电离层波导中传播路径上的各种粒子的密度和温度的实时数据集；

步骤2：由步骤1得到的实时数据集，利用经验公式计算对电离层起到直接影响作用的两个参数——电子密度N_e和碰撞频率v，实时模拟传播路径上的地-电离层环境；

步骤3：将极坐标系下电子密度N_e和碰撞频率v经过双线性插值算法处理后转换到直角坐标系；

步骤4：根据步骤3转换到直角坐标系的电子密度N_e和碰撞频率v，利用SO-FDTD算法计算甚低频电磁波在传播路径上的场分布情况，将电磁波接收点一天内不同时间的场强变换情况与甚低频台站上实测数据进行比对和验证；

步骤5：控制各地-电离层参数变化，对计算得到的场分布情况进行分析总结，得到甚低频电磁波在地-电离层波导中传播特性。

进一步地，所述步骤1具体包括以下过程：

步骤1-1：使用Python编程语言中的Requests库实现爬虫来爬取所需数据。对于电离层IRI2016模型(https://ccmc.gsfc.nasa.gov/modelweb/models/iri2016_vitmo.php)，需要爬取传播路径上的电子密度N_e和电子温度T_e。对于NRLMSISE-00大气模型(https://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php)，需要爬取传播路径上的氧气分子密度

氧原子密度N_o和氮气分子密度

。

步骤1-2：建立甚低频电磁波在地-电离层波导中传播路径上粒子密度和温度的实时数据集；

以地心为极点，地心与发射点连线延长所得到的射线为极轴，建立极坐标系，从而将三维球坐标系下的问题转换到二维极坐标系。对于坐标系内任一点P(ρ,θ)来说，极径ρ对应该点到地心的距离，极角θ代表到该点与地心连线和发射点与地心连线之间的夹角。如图6，从发射点开始直到接收点结束(实际上考虑到FDTD中的边界，所以是从发射点稍外开始到接收点稍外结束)，每隔一个很小的角度，就通过爬虫爬取一次对应位置从地面到低电离层中每隔一个很小的高度，各种粒子密度和温度的变化情况，这样仅用一个二维数组就存储了传播路径上各种粒子密度和温度的数据集。

进一步地，所述步骤2具体包括以下过程：

步骤2-1：模拟地面环境；

将地面近似为均匀的电磁媒质，具体如表1所示：

表1地面近似电磁媒质表

	介电常数ε	磁导率μ
			平均陆地状态	10	0.003
平均海洋状态	80	4

步骤2-2：模拟电离层环境；

对于电离层来说，起到直接影响作用的有两个参数——电子密度N_e和碰撞频率v，其中电子密度N_e直接从电离层IRI2016模型中获取，碰撞频率v则通过经验公式(1)计算，具体如下：

其中，

式中，V_e,i代表电子和离子的碰撞频率，

代表电子和氧气分子的碰撞频率，V_e,O代表电子和氧原子的碰撞频率，

代表电子和氮气分子的碰撞频率，N_e代表电子密度，T_e代表电子温度，

代表氧气分子密度，N_O代表氧原子密度，

代表氮气分子密度；

步骤2-3：模拟地-电离层之间环境；

如图4，对于地面和电离层之间来说，由于白天电离层下边界高度在65km左右，晚上电离层下边界高度在80km左右，从地面到电离层下边界这一部分是空缺的，所以采用经验指数模型进行表示，具体如下：

v(z)＝1.82×10¹¹e^-0.15z (3)

N(z)＝1.43×10⁷e^-0.15He^{[(β-0.15)(z-H)]} (4)

式中，v(z)代表高度为z时的碰撞频率，N(z)代表高度为z时的电子密度，β代表梯度系数(km^-1)，H代表电离层参考高度(km)；

其中经验指数模型内中低纬度地区β和H的推荐值见表2，其中，f为工作频率；

表2中低纬度地区β和H取值

	夏季	冬季
			白天	β＝0.3,H＝70	β＝0.3,H＝72
晚上	β＝0.0077f+0.31,H＝87	β＝0.0077f+0.31,H＝87

进一步地，所述步骤3具体包括以下过程：

对于FDTD来说，极坐标系和直角坐标系中的Yee元胞不是一一对应的，所以不能直接转换，需要通过插值进行处理；采用双线性插值算法将极坐标系下电子密度N_e和碰撞频率v转换到直角坐标系；双线性插值算法的核心思想是在x方向和y方向分别进行一次插值，即直角坐标系中每一点的值是由邻近四个极坐标系中的点插值得到的；对于点P(x,y)的值，已知Q₁₁(x₁,y₁)，Q₁₂(x₂,y₂)，Q₁₃(x₃,y₃)和Q₁₄(x₄,y₄)四个点的值，插值公式为

进一步地，所述步骤4具体包括以下过程：

步骤4-1：利用SO-FDTD计算传播路径上的场分布；

线性各向同性介质麦克斯韦旋度方程为：

式中，G代表磁场强度，D代表电通量密度，σ代表电导率，E代表电场强度，B代表磁通量密度，σ_m代表等效磁导率；

电离层的介电常数随频率变化，所以电离层也属于色散介质，对于此类介质有：

B＝μG (7)

D(ω)＝ε(ω)E(ω) (8)

式中，μ代表磁导率，ε(ω)代表介电常数，ω代表角频率；

由式(6)进行FDTD离散得到：

将式(8)由频域转换到时域，对于x分量，得：

电离层是非常典型的Drude介质，所以有：

即有：

设函数：

式(13)左端平均值近似，右端中心差分近似得：

引进离散移位算子z_l，定义为：

z_lfⁿ＝fⁿ⁺¹ (15)

联立式(14)和式(15)，得：

联立式(13)和式(16)，得：

式(17)代入式(12)，得：

两边同乘(z_l+1)^N，得：

因为移位算子的作用使得

所以得到步进公式为：

当M＝N＝2时，整理得：

其中，

电离层是等离子体，所以有：

其中ω_p是等离子体频率，v_c是电子碰撞频率，ε_∞＝1；

将SO-FDTD的步进计算步骤归结如下：(1)由E→G，用式(9)第一式计算；(2)由G→D，用式(9)第二式计算；(3)由D→E，用式(20)计算；由此能够在时间空间轴上逐步推进求解甚低频电磁波在地-电离层中传播产生的辐射场分布情况；

步骤4-2：比对计算数据与实测数据，验证方法可行性。

进一步地，所述步骤5具体包括以下过程：

步骤5-1：为了方便分析甚低频电磁波在地-电离层波导中的传播特性，在Visualstudio平台利用Fortran语言进行编程，设置统一坐标系下的直角坐标值为变量，求得电磁场分布的矩阵数值。利用Origin软件进行绘图，得到传播路径上各个电磁场分量的振幅值和相位值分布，以此为基础分析甚低频电磁波在地-电离层波导中的传播特性。

步骤5-2：考虑各地-电离层参数对传播过程产生的影响，可以设置各地-电离层参数为变量，通过对比各地-电离层参数不同时的电磁场分量的振幅和相位值分布情况，得出甚低频电磁波在地-电离层波导中的传播特性结论。

具体实施例：

本发明方法在于通过Visual studio平台通过Fortran语言编程，对爬取到的初始数据经处理后构建出地-电离层环境，通过两次坐标系变换简化问题，利用SO-FDTD算法计算传播路径上的场分布情况，控制地-电离层参数分析总结传播特性。

如图1所示为计算流程图。首先模拟地-电离层的真实电磁环境，对于电离层来说，从电离层IRI2016模型和NRLMSISE-00大气模型爬取各种粒子的密度和温度，利用经验公式求得电子密度和碰撞频率，对于地面来说，近似成某种均匀的电磁媒质，对于地面和电离层之间来说，利用指数经验模型近似得到电子密度和碰撞频率；在此基础上已经将三维球坐标系下的问题转换到二维极坐标系，再利用双线性插值算法将二维极坐标系下的问题转换到二维直角坐标系；SO-FDTD算法可求出传播路径上的场分布情况，比对计算结果和实测数据，验证方法可行性；控制各地-电离层参数，观察其对传播过程产生的影响，分析总结传播特性。

如图2所示为问题示意图。从内到外依次为地面、地面和电离层之间以及电离层，发射点和接收点均位于地面上。实际传播过程中，甚低频电磁波在地-电离层中的传播可近似于在波导中的传播。

如图3所示为爬虫示意图。爬虫向电离层IRI2016模型发出一个携带时间，经纬度，高度等参数的GET请求，返回电子密度N_e和电子温度T_e的信息。同理，爬虫向NRLMSISE-00大气模型发出一个携带时间，经纬度，高度等参数的GET请求，返回氧气分子密度

氧原子密度N_o和氮气分子密度

信息。

如图4所示为实时模拟地-电离层电磁环境方法示意图。从内到外依次为地面、地面和电离层之间以及电离层，对于电离层来说，从电离层IRI2016模型和NRLMSISE-00大气模型爬取各种粒子的密度和温度，利用经验公式求得电子密度和碰撞频率，对于地面来说，近似成某种均匀的电磁媒质，对于地面和电离层之间来说，利用指数经验模型近似得到电子密度和碰撞频率，这样就能够实时模拟出地-电离层电磁环境。

如图5所示为坐标系转换示意图。从左到右依次为三维球坐标系示意图，二维极坐标系示意图和二维直角坐标系示意图。甚低频电磁波在地-电离层波导中的传播是一个三维球坐标系下的问题。当以地心为极点，地心与发射点连线延长所得到的射线为极轴，建立极坐标系，三维球坐标系下的问题可以转换到二维极坐标系中。在极坐标系中取一个能够同时包含发射点和接收点的“矩形”，二维极坐标系下的问题可以转换到二维直角坐标系中。

如图6所示为三维球坐标系到二维极坐标系转换方法示意图(或初始数据集构建示意图)。左为二维极坐标系示意图，右为初始数据集示意图。从发射点开始直到接收点结束(实际上考虑到FDTD中的边界，所以是从发射点稍外开始到接收点稍外结束)，每隔一个很小的角度，就通过爬虫爬取一次对应位置从地面到低电离层中每隔一个很小的高度，各种粒子密度和温度的变化情况，这样仅用一个二维数组就存储了传播路径上各种粒子密度和温度的数据集。

如图7所示为二维极坐标系到二维直角坐标系转换方法示意图。点P(x,y)对应FDTD直角坐标系上的某个Yee元胞，点Q₁₁(x₁,y₁)，Q₁₂(x₂,y₂)，Q₁₃(x₃,y₃)和Q₁₄(x₄,y₄)对应距点P(x,y)最近的FDTD极坐标上四个的Yee元胞，根据双线性插值算法，点P(x,y)的值由Q₁₁(x₁,y₁)，Q₁₂(x₂,y₂)，Q₁₃(x₃,y₃)和Q₁₄(x₄,y₄)四个点的值来决定，将极坐标系下的电子密度和碰撞频率转换到直角坐标系。

如图8所示为SO-FDTD运行原理示意图。图中，每一个电场分量都有四个磁场分量环绕，每一个磁场分量都有四个电场分量环绕；电场和磁场在时间顺序上交替抽样，抽样时间间隔相差半个时间步；E→H→D→E每一步均有明确的公式，所以可以在时间空间轴上通过迭代逐步推进求解甚低频电磁波在地-电离层中传播产生的辐射场分布情况。

Claims (6)

1.一种分析甚低频电磁波在地-电离层中传播的FDTD方法，其特征在于，包括如下步骤：

步骤1：以地心为极点，地心与电磁波发射点连线延长所得到的射线为极轴，建立极坐标系，将三维球坐标系下的问题转换到二维极坐标系；

从电离层IRI2016模型和NRLMSISE-00大气模型中爬取数据，建立甚低频电磁波在地-电离层波导中传播路径上粒子的密度和温度的实时数据集；

步骤2：由步骤1得到的实时数据集，计算参数——电子密度N_e和碰撞频率v，实时模拟传播路径上的地-电离层环境；

步骤3：将极坐标系下电子密度N_e和碰撞频率v经过双线性插值算法处理后转换到直角坐标系；

步骤4：根据步骤3转换到直角坐标系的电子密度N_e和碰撞频率v，利用SO-FDTD算法计算甚低频电磁波在传播路径上的场分布情况，将电磁波接收点一天内不同时间的场强变换情况与甚低频台站上实测数据进行比对和验证；

步骤5：控制各地-电离层参数变化，对计算得到的场分布情况进行分析，得到甚低频电磁波在地-电离层波导中传播特性。

2.根据权利要求1所述的一种分析甚低频电磁波在地-电离层中传播的FDTD方法，其特征在于，所述步骤1具体包括以下过程：

步骤1-1：利用爬虫爬取电离层IRI2016模型传播路径上的电子密度N_e和电子温度T_e，爬取NRLMSISE-00大气模型传播路径上的氧气分子密度

氧原子密度N_o和氮气分子密度

步骤1-2：建立甚低频电磁波在地-电离层波导中传播路径上粒子密度和温度的实时数据集；

以地心为极点，地心与电磁波发射点连线延长所得到的射线为极轴，建立极坐标系，将三维球坐标系下的问题转换到二维极坐标系；对于二维极坐标系内任一点P(ρ,θ)，极径ρ对应该点到地心的距离，极角θ代表到该点与地心连线和电磁波发射点与地心连线之间的夹角；从电磁波发射点开始直到电磁波接收点结束，每间隔一个角度，就通过爬虫爬取一次对应位置从地面到低电离层中每隔一定高度的粒子密度和温度的实时数据，建立实时数据集。

3.根据权利要求2所述的一种分析甚低频电磁波在地-电离层中传播的FDTD方法，其特征在于，所述步骤2具体包括以下过程：

步骤2-1：模拟地面环境；

将地面近似为均匀的电磁媒质，具体如表1所示：

表1 地面近似电磁媒质表

介电常数ε 磁导率μ 平均陆地状态 10 0.003 平均海洋状态 80 4

步骤2-2：模拟电离层环境；

电离层的两个参数——电子密度N_e和碰撞频率v，其中电子密度N_e直接从电离层IRI2016模型中获取，碰撞频率v则通过经验公式(1)计算，具体如下：

其中，

式中，V_e,i代表电子和离子的碰撞频率，

代表电子和氧气分子的碰撞频率，V_e,O代表电子和氧原子的碰撞频率，

代表电子和氮气分子的碰撞频率，N_e代表电子密度，T_e代表电子温度，

代表氧气分子密度，N_O代表氧原子密度，

代表氮气分子密度；

步骤2-3：模拟地-电离层之间环境；

对于地面和电离层之间来说，采用经验指数模型进行表示，具体如下：

v(z)＝1.82×10¹¹e^-0.15z (3)

N(z)＝1.43×10⁷e^-0.15He^{[(β-0.15)(z-H)]} (4)

式中，v(z)代表高度为z时的碰撞频率，N(z)代表高度为z时的电子密度，β代表梯度系数(km^-1)，H代表电离层参考高度(km)；

其中经验指数模型内中低纬度地区β和H的推荐值见表2，其中，f为工作频率；

表2 中低纬度地区β和H取值

夏季 冬季 白天 β＝0.3,H＝70 β＝0.3,H＝72 晚上 β＝0.0077f+0.31,H＝87 β＝0.0077f+0.31,H＝87

4.根据权利要求3所述的一种分析甚低频电磁波在地-电离层中传播的FDTD方法，其特征在于，所述步骤3具体包括以下过程：

对于FDTD来说，极坐标系和直角坐标系中的Yee元胞不是一一对应的，所以不能直接转换，需要通过插值进行处理；采用双线性插值算法将极坐标系下电子密度N_e和碰撞频率v转换到直角坐标系；双线性插值算法的核心思想是在x方向和y方向分别进行一次插值，即直角坐标系中每一点的值是由邻近四个极坐标系中的点插值得到的；对于点P(x,y)的值，已知Q₁₁(x₁,y₁)，Q₁₂(x₂,y₂)，Q₁₃(x₃,y₃)和Q₁₄(x₄,y₄)四个点的值，插值公式为

5.根据权利要求4所述的一种分析甚低频电磁波在地-电离层中传播的FDTD方法，其特征在于，所述步骤4具体包括以下过程：

步骤4-1：利用SO-FDTD计算传播路径上的场分布；

线性各向同性介质麦克斯韦旋度方程为：

式中，G代表磁场强度，D代表电通量密度，σ代表电导率，E代表电场强度，B代表磁通量密度，σ_m代表等效磁导率；

电离层的介电常数随频率变化，所以电离层也属于色散介质，对于此类介质有：

B＝μG (7)

D(ω)＝ε(ω)E(ω) (8)

式中，μ代表磁导率，ε(ω)代表介电常数，ω代表角频率；

由式(6)进行FDTD离散得到：

式中，n代表迭代步数；

将式(8)由频域转换到时域，对于x分量，得：

式中，ε₀代表真空介电常数，ε_r代表相对介电常数，E_x(t)代表E(t)的x分量，D_x(t)代表D(t)的x分量；

电离层是Drude介质，所以有：

式中，p_l和q_l是多项式系数，M和N代表多项式总项数；

即有：

设函数：

式(13)左端平均值近似，右端中心差分近似得：

引进离散移位算子z_l，定义为：

z_lfⁿ＝fⁿ⁺¹ (15)

联立式(14)和式(15)，得：

联立式(13)和式(16)，得：

式(17)代入式(12)，得：

两边同乘(z_l+1)^N，得：

因为移位算子的作用使得

所以得到步进公式为：

式中，a_l和b_l由有理分式式(11)中的系数p₀,p₁,...,p_N和q₀,q₁,...,q_M表示；

当M＝N＝2时，整理得：

其中，

电离层是等离子体，所以有：

其中ω_p是等离子体频率，v_c是电子碰撞频率，ε_∞＝1；

将SO-FDTD的步进计算步骤归结如下：(1)由E→G，用式(9)第一式计算；(2)由G→D，用式(9)第二式计算；(3)由D→E，用式(20)计算；由此能够在时间空间轴上逐步推进求解甚低频电磁波在地-电离层中传播产生的辐射场分布情况；

步骤4-2：比对计算数据与实测数据，验证方法可行性。

6.根据权利要求5所述的一种分析甚低频电磁波在地-电离层中传播的FDTD方法，其特征在于，所述步骤5具体包括以下过程：

步骤5-1：设置直角坐标值为变量，求得电磁场分布的矩阵数值；绘图得到传播路径上各个电磁场分量的振幅值和相位值分布，以此为基础分析甚低频电磁波在地-电离层波导中的传播特性；

步骤5-2：考虑各地-电离层参数对传播过程产生的影响，设置各地-电离层参数为变量，通过对比各地-电离层参数不同时的电磁场分量的振幅和相位值分布情况，得出甚低频电磁波在地-电离层波导中的传播特性结论。

Images