一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法

摘要

本发明公开了一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法，该方法通过对整个计算区域进行分层处理，将粗网格M(2，4)FDTD算法(用于上层区域)与细网格传统FDTD算法(用于下层区域)相结合，利用修正因子以及亚网格技术，进行低频地波传播时延的高精度快速预测，本发明解决了传统全区域FDTD算法在减小数值色散误差与降低计算机消耗之间的矛盾，提高了低频地波传播时延预测精度的同时，减少了计算机内存资源占用，提高了计算速度。

说明书

技术领域

本发明属于电波传播理论和电磁场数值计算技术领域，具体涉及一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法。

背景技术

低频地波传播时延预测精度是提高陆基无线电定位/导航/授时系统精度的关键。现有的经典理论预测方法主要有：基于均匀光滑地面模型的平地面公式和Fock地波绕射公式等；基于分段不均匀光滑地面模型的Millington经验公式、Wait积分公式、波模转换法等；基于不均匀不光滑地面模型的积分方程方法、抛物方程方法等。然而，一方面，经典的理论预测方法均是在一定的模型近似及理论假设条件下得到的，使用时无法考虑实际传播中的各种因素的复杂变化(如复杂路径的后向反射波的影响、地形横向剧烈变化的影响、大气折射率空时变化的影响等)；另一方面，经典理论预测方法多是针对单一频率的频域预测结果，而对经过调制的脉冲信号(如罗兰C信号)，计算距离越长、结果误差越大。

近些年，随着计算机技术的不断提高以及数值计算理论和方法的进一步发展，时域有限差分方法(Finite-Difference Time-Domain Method，FDTD)被用于该领域来进行低频地波传播时延的高精度预测。它从Maxwell旋度方程出发，利用二阶精度的中心差分近似，将旋度方程中的时间及空间微分算符直接转换为差分形式，从而在有限的空间和时间上对连续的电磁场数据进 行抽样。与经典的地波传播时延理论预测方法相比较，FDTD方法具有较高的预测精度，但是，随着传播距离的增长以及路径复杂度的增大，FDTD方法的预测精度受到显著影响，究其原因主要是FDTD方法存在数值色散误差。理论上，通过进一步增加网格剖分密度或采用更高阶的FDTD算法等措施可以减小其数值色散误差，但这就意味着需要占用更大的计算机内存资源、消耗更长的计算时间，难以用于长距离、大区域甚至三维FDTD在低频地波传播时延预测中的工程推广。

为了解决传统FDTD方法数值色散问题，国内外学者分别从时间离散方式和空间差分格式出发对传统FDTD方法进行了一系列的改进，其中M(2，4)FDTD是Hadi和Piket-May从Maxwell方程的积分形式出发，结合Fang的FDTD(2，4)算法提出的一种高阶FDTD方法，通过引入两个修正系数k₁和k₂，可降低特定频率所有方向的数值色散。但是，该方法无法直接用于非均匀媒质。

目前，尚未见M(2，4)FDTD与FDTD相结合方法在低频地波传播时延预测方面的研究报道和专利。

发明内容

本发明的目的是提供一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法，解决了现有技术中存在的FDTD数值方法在进行低频地波传播时延预测时，随传播距离增长、路径复杂度提高而存在的“数值计算误差增大、计算速度缓慢、内存资源占用较大”的问题。

本发明所采用的技术方案是，一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法，通过对低色散低频地波计算区域进行分层处理，将粗网格M(2，4)FDTD与细网格传统FDTD方法相结合，进行低频地波传 播时延的预测，从而达到在保证预测精度的同时，提高计算速度、减少计算机内存占用，具体步骤如下：

步骤1、设置并生成模型文件、输入模型文件；

步骤2、参数设置及初始化；

步骤3、添加场源；

步骤4、更新整个计算区域的ρ向电场分量E_ρ；

步骤5、更新整个计算区域的z向电场分量E_z；

步骤6、更新整个计算区域的φ向磁场分量H_φ；

步骤7、更新激励源；

步骤8、判断结束条件，循环；

步骤9、观测点观测量计算并输出。

本发明的特点还在于，

步骤1具体为：

步骤(1.1)、设模型文件的上层区域大小N_ρ1×N_z1，下层区域大小N_ρ2×N_z2，CFS-PML层数为N_PML，其中N_ρ为ρ方向网格数，N_z为z方向网格数，标号1和2分别表示上层和下层；

步骤(1.2)、设置空间、时间步长：上层粗网格空间步长为δ1，其中，Δρ1＝Δz1＝δ1，下层细网格空间步长为δ2，其中，Δρ2＝Δz2＝δ2，其中Δρ为ρ方向网格大小，Δz为z方向网格大小，标号1和2分别表示上层和下层，时间步长为Δt，粗网格M(2，4)FDTD的时间步长设置为与细网格FDTD同样的时间步长；

步骤(1.3)、设置迭代次数为NT_Calc；

步骤(1.4)、设置传播路径电参数：ρ方向起始网格位置为ρ_Start、ρ方 向结束网格位置为ρ_Een、z方向起始网格位置为z_Start、z方向结束网格位置为z_Start、大地电导率为σ、相对介电常数为ε_r；

步骤(1.5)、设置吸收边界：CFS-PML层数为N_PML，相关参数为κ_ηmax，α_ηmax，σ_ηmax，其中η＝ρ,z；

步骤(1.6)、设置场源：源的个数N_Source、位置[ρ_Start，ρ_Een]和[z_Start，z_Een]，源的种类：有两种激励源供选择：单频正弦源、Loran_C源，激励场形式：有三种激励形式：E_z、E_ρ、源的类型：软源或硬源，幅度A，频率f，单频正弦源中的常数t₀，高斯脉宽T_p，时延/包周差τ；

步骤(1.7)、设置观测点：观测点个数N_VPoint，位置[ρ_Start，ρ_Een]和[z_Start，z_Een]，输出场量类型(E_z、E_ρ或H_φ)。

步骤2具体为：

步骤(2.1)、将整个区域的电磁场分量(E_z、E_ρ和H_φ)和电磁场分量系数(CA、CB)、中间变量中间变量系数(f_0z、f_1z、f_2z，f_0ρ、f_1ρ、f_2ρ，f_0φγ、f_1φγ、f_2φγ，)、辅助变量(ψ_φz、ψ_φρ、ψ_φγ、ψ_ρz、ψ_zρ)、观测量(传播时延t_w、采样点场强峰值场强)都初始化为零；

步骤(2.2)、初始化所有网格中的模型参数：将相对介电常数ε_r初始化为1，大地电导率σ初始化为0；

步骤(2.3)、根据步骤1所设置的模型文件中的路径信息，给相应网格的ε_r和σ赋值；

步骤(2.4)设置CFS-PML吸收边界参数k_η、σ_η、α_η，其中，k_η、σ_η、α_η具体按下式计算：

$\left(\begin{array}{c}{k}_{\eta }=1+(k_{\eta max}-1)\frac{|\eta -{\eta }_0{|}^m}{d^m}\\ {\sigma }_{\eta }=\frac{{\sigma }_{\eta max}|\eta -{\eta }_0{|}^m}{d^m}\\ {\alpha }_{\eta }={\alpha }_{\eta max}\frac{|\eta -{\eta }_0|}{d}\end{array}\right)$

式中，η＝ρ,z，η₀为PML层与非PML截面位置，d是PML吸收边界的厚度，κ_ηmax取整数，κ_ηmax取值范围为[1，60]，α_ηmax取值范围为[0，1)，σ_ηmax根据σ_opt设置，σ_ηmax/σ_opt取值范围为(0，12]，σ_opt＝(m+1)/150πΔη，m取值范围为[1，20]，其中m取值为4时边界的吸收效果最好，Δη取值范围λ为源的波长。

步骤3具体为：

添加的场源有两种类型：

(1)单频正弦源：

当辐射源采用单频正弦信号时，其电流源激励i_s(t)可以表示为：

$i_s\left(t\right)=\left(\begin{array}{cc}0& t<0\\ 0.5[1-cos\left(\frac{\pi t}{t_0}\right)]sin\left(2\pi ft\right)& 0\le t\le t_0\\ sin\left(2\pi ft\right)& t_0\le t\end{array}\right),$

其中t₀＝5×10^-6s。

(2)Loran-C源：

当辐射源采用正相位编码的Loran_C信号时，其电流波形前沿i_s(t)为：

$i_s\left(t\right)=\{\begin{array}{cc}0& t<\tau \\ A{(t-\tau )}^2{e}^{-\frac{2(t-\tau )}{65\times {10}^{-6}}}\sin \left(2\pi ft\right)& \tau \le t\le 65\times {10}^{-6}+\tau \end{array},$

其中，τ为包周差，单位为s。

步骤4具体为：

步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量E_ρ1进行更新，其中，上层区域不包含PML层，E_ρ1中的下标1表示上层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{\rho 1}{|}_{i1+1/2,k1}^{n+1}=CA\left(m\right)·E_{\rho 1}{|}_{i1+1/2,k1}^n\\ -CB\left(m\right)·\left(\begin{array}{c}{k}_1\left(\frac{H_{\phi 1}{|}_{i+1/2,k1+3/2}^{n+1/2}-H_{\phi 1}{|}_{i1+1/2,k1-3/2}^{n+1/2}}{3\delta 1}\right)\\ -(1-k_1)\left(\frac{H_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n+1/2}-H_{\phi 1}{|}_{i1+1/2,k1-1/2}^{n+1/2}}{\delta 1}\right)\end{array}\right)\end{array}\right)$

$\left(\begin{array}{c}CA\left(m\right)=\frac{\frac{{ϵ}_r\left(m\right)}{\Delta t}-\frac{\sigma \left(m\right)}{2{ϵ}_0}}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\\ CB\left(m\right)=\frac{c}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\end{array}\right)$

式中，n表示时间步，标号i1和k1分别表示上层区域中ρ向和z向的空间位置，ε₀为真空中介电常数，k₁为环路系数；

步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区域内的电磁场分量E_ρ2进行更新，其中，下层区域不包含PML层，E_ρ2中的下标2表示下层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{\rho 2}{|}_{i+1/2,k2}^{n+1}=CA\left(m\right)·E_{\rho 2}{|}_{i2+1/2,k2}^n\\ -CB\left(m\right)·\left(\frac{H_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n+1/2}-H_{\phi 2}{|}_{i2+1/2,k2-1/2}^{n+1/2}}{\delta 2}\right)\end{array}\right)$

$\left(\begin{array}{c}CA\left(m\right)=\frac{\frac{{ϵ}_r\left(m\right)}{\Delta t}-\frac{\sigma \left(m\right)}{2{ϵ}_0}}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\\ CB\left(m\right)=\frac{c}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\end{array}\right)$

式中，n表示时间步，标号i2和k2分别表示下层区域中ρ向和z向的空间位置；

步骤(4.3)、对边界区域的场量传递进行更新：

边界区域场量的传递分为如下几种情况：

a、当上层粗网格中E_ρ1的计算需用到下层细网格中的H_φ2时，由于上下粗细网格比为奇数，场量重叠，故直接取相应细网格的场量即可；

b、当下层细网格中边界网格处E_ρ2的计算需用到边界线上的H_φ时，直接采用边界线上细网格的H_φ2即可；

步骤(4.4)、对CFS-PML中的电场分量进行更新：

上方吸收边界中的网格空间大小与粗网格相同，右侧吸收边界中网格的空间大小与相邻左侧计算区域的网格空间大小相同，也分为上下两层，上层与粗网格一致，下层与细网格一致，

吸收边界中的场分量的差分公式为：

$\left(\begin{array}{c}{E}_{\rho }^{*}{|}_{i+1/2,k}^{n+1}={\mathrm{CA}}_{\rho }\left(m\right)E_{\rho }^{*}{|}_{i+1/2,k}^n\\ +{\mathrm{CB}}_{\rho }\left(m\right)[{\psi }_{\phi z}{|}_{i+1/2,k}^{n-1/2}+f_{1z}{|}_k\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\phi }{|}_{i+1/2,k-1/2}^{n+1/2}}{\Delta z}\right)]\end{array}\right)$

其中，

$\left(\begin{array}{c}{\psi }_{\phi z}{|}_{i+1/2,k}^{n-1/2}=f_{0z}{|}_k{\tilde{H}}_{\phi z}{|}_{i+1/2,k}^{n-1/2}+f_{2z}{|}_k\frac{\partial }{\partial z}{H}_{\phi }{|}^{n-1/2}\\ =f_{0z}{|}_k{\tilde{H}}_{\phi z}{|}_{i+1/2,k}^{n-1/2}+f_{2z}{|}_k\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n-1/2}-H_{\phi }{|}_{i+1/2,k-1/2}^{n-1/2}}{\Delta z}\right)\end{array}\right)$

$\left(\begin{array}{c}{f}_{0z}{|}_k=\frac{2{ϵ}_0{k}_z{|}_k-\Delta t(k_z{|}_k·{\alpha }_z{|}_k+{\sigma }_z{|}_k)}{2{ϵ}_0{k}_z{|}_k+\Delta t(k_z{|}_k·{\alpha }_z+{\sigma }_z{|}_k)}\\ f_{1z}{|}_k=\frac{{\alpha }_z{|}_k\Delta t+2{ϵ}_0}{2{ϵ}_0{k}_z{|}_k+\Delta t(k_z{|}_k·\alpha {|}_k+{\sigma }_z{|}_k)}\\ f_{2z}{|}_k=\frac{{\alpha }_z{|}_k\Delta t-2{ϵ}_0}{2{ϵ}_0{k}_z{|}_k+\Delta t(k_z{|}_k·{\alpha }_z{|}_k+{\sigma }_z{|}_k)}\end{array}\right),$

其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，k_z、σ_z和α_z为吸收边界参数。

步骤5具体为：

步骤(5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量E_z1进行更新，其中，上层区域不包含PML层，E_z1中的下标1表示上层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{z1}{|}_{i1,k1+1/2}^{n+1}=CA\left(m\right)·E_{z1}{|}_{i1,k1+1/2}^n\\ +CB\left(m\right)·\left(\begin{array}{c}{k}_1\frac{(i1+\frac{1}{2})H_{\phi 1}{|}_{i+13/2,k1+1/2}^{n+1/2}-(i1-\frac{1}{2})H_{\phi 1}{|}_{i1-3/2,k1+1/2}^{n+1/2}}{3i1\delta 1}\\ +(1+k_1)\frac{(i1+\frac{1}{2})H_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n+1/2}-(i1-\frac{1}{2})H_{\phi 1}{|}_{i1-1/2,k1+1/2}^{n+1/2}}{i1\delta 1}\end{array}\right)\end{array}\right)$

式中，n表示时间步，标号i1和k1分别表示上层区域中ρ向和z向的空间位置，ε₀为真空中介电常数，k₁为环路系数；

步骤(5.2)、根据步骤1中模型文件中定义的下层区域采用传统FDTD 方法，对该区域内的电磁场分量E_z2进行更新，其中，下层区域不包含PML层，E_z2中下标2表示下层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{z2}{|}_{i2,k2+1/2}^{n+1}=CA\left(m\right)·E_{z2}{|}_{i2,k2+1/2}^n\\ +CB\left(m\right)·\left[\frac{(i2+\frac{1}{2})H_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n+1/2}-(i2-\frac{1}{2})H_{\phi 2}{|}_{i2-1/2,k2+1/2}^{n+1/2}}{i2\delta 2}\right]\end{array}\right)$

式中，n表示时间步，标号i2和k2分别表示下层区域中ρ向和z向的空间位置；

步骤(5.3)、交界线上E_z均为细网格的E_z2，上层粗网格中E_z1的更新不包括交界线上的点；

步骤(5.4)、对CFS-PML中的电场分量进行更新：吸收边界中的场分量的差分公式为：

$\left(\begin{array}{c}{E}_z^{*}{|}_{i,k+1/2}^{n+1}=CA\left(m\right)E_z^{*}{|}_{i,k+1/2}^n\\ +CB\left(m\right)\left(\begin{array}{c}[{\psi }_{\phi \rho }{|}_{i,k+1/2}^{n-1/2}+f_{1\rho }{|}_i\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\phi }{|}_{i-1/2,k+1/2}^{n+1/2}}{\delta }\right)]\\ +[{\psi }_{\phi \gamma }{|}_{i,k+1/2}^{n-1/2}+f_{1\phi \gamma }{|}_i\frac{H_{\phi }}{i\delta }{|}_{i,k+1/2}^{n+1/2}]\end{array}\right)\end{array}\right)$

式中，标号m＝(i,k+1/2)，且

$\left(\begin{array}{c}{\psi }_{\phi \rho }{|}_{i,k+1/2}^{n-1/2}=f_{0\rho }{|}_i{\tilde{H}}_{\phi \rho }{|}_{i,k+1/2}^{n-1/2}+f_{2\rho }{|}_i\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n-1/2}-H_{\phi }{|}_{i-1/2,k+1/2}^{n-1/2}}{\Delta \rho }\right)\\ {\psi }_{\phi \gamma }{|}_{i,k+1/2}^{n-1/2}=f_{0\phi \gamma }{|}_i{\tilde{H}}_{\phi \gamma }{|}_{i,k+1/2}^{n-1/2}+f_{2\phi \gamma }{|}_i\frac{H_{\phi }}{i\delta }{|}_{i,k+1/2}^{n-1/2}\end{array}\right),$

$\left(\begin{array}{c}{f}_{0\rho }{|}_i=\frac{2{ϵ}_0{k}_{\rho }{|}_i-\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}{2{ϵ}_0{k}_{\rho }{|}_i+\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}\\ f_{1\rho }{|}_i=\frac{{\alpha }_{\rho }{|}_i\Delta t+2{ϵ}_0}{2{ϵ}_0{k}_{\rho }{|}_i+\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}\\ f_{2\rho }{|}_i=\frac{{\alpha }_{\rho }{|}_i\Delta t-2{ϵ}_0}{2{ϵ}_0{k}_{\rho }{|}_i+\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}\end{array}\right)$

$\left(\begin{array}{c}{f}_{0\phi \gamma }{|}_i=-\frac{\Delta t(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)-2{ϵ}_0{A}_{\rho }{|}_i}{\Delta t(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)+2{ϵ}_0{A}_{\rho }{|}_i}\\ f_{1\phi \gamma }{|}_i=\frac{({\mathrm{\Delta t\alpha }}_{\rho }{|}_i+2{ϵ}_0)}{\left(\Delta t\right(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)+2{ϵ}_0{A}_{\rho }{|}_i)}\\ f_{2\phi \gamma }{|}_i=\frac{({\mathrm{\Delta t\alpha }}_{\rho }{|}_i-2{ϵ}_0)}{\left(\Delta t\right(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)+2{ϵ}_0{A}_{\rho }{|}_i)}\end{array}\right)$

$\left(\begin{array}{c}{A}_{\rho }=1+(k_{\rho \max }-1)\frac{{(\rho -{\rho }_0)}^{m+1}}{(m+1)d^m}\frac{1}{\rho }\\ B_{\rho }=\frac{1}{{\alpha }_{\rho }+{\mathrm{j\omega ϵ}}_0}\frac{{\sigma }_{\rho \max }{(\rho -{\rho }_0)}^{m+1}}{(m+1)d^m}\frac{1}{\rho }\end{array}\right)$

其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，α_ρ、k_ρ、σ_ρ、k_ρmax和σ_ρmax为吸收边界参数，ρ＝i×δ表示ρ向长度，ρ₀是PML层与非PML层的分界位置。

步骤6具体为：

步骤(6.1)、首先根据步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量进行更新，其中，上层区域不包含PML层，具体更新公式为：

$\begin{array}{c}{H}_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n+1/2}=H_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n-1/2}\\ +\frac{c\Delta t}{{\mu }_r{|}_{i1+1/2,k1+1/2}^{n+1/2}}·\left(\begin{array}{c}{k}_1\frac{E_{z1}{|}_{i1+2,k1+1/2}^n-E_{z1}{|}_{i1-1,k1+1/2}^n}{3\delta 1}\\ +k_2\frac{E_{z1}{|}_{i1+2,k1+3/2}^n+E_{z1}{|}_{i1+2,k1-1/2}^n-E_{z1}{|}_{i1-1,k1+3/2}^n-E_{z1}{|}_{i1-1,k1-1/2}^n}{6\delta 1}\\ +(1-k_1-k_2)\frac{E_{z1}{|}_{i1+1,k1+1/2}^n-E_{z1}{|}_{i1,k1+1/2}^n}{\delta 1}\end{array}\right)\\ -\frac{c\Delta t}{{\mu }_r{|}_{i1+1/2,k1+1/2}^{n+1/2}}·\left(\begin{array}{c}{k}_1\frac{E_{\rho 1}{|}_{i1+1/2,k1+2}^n-E_{\rho 1}{|}_{i1+1/2k1-1}^n}{3\delta 1}\\ +k_2\frac{E_{\rho 1}{|}_{i1+3/2,k1+2}^n+E_{\rho 1}{|}_{i1-1/2,k1+2}^n-E_{\rho 1}{|}_{i1+3/2,k1-1}^n-E_{\rho 1}{|}_{i1-1/2,k1-1}^n}{6\delta 1}\\ +(1+k_1-k_2)\frac{E_{\rho 1}{|}_{i1+1/2,k1+1}^n-E_{\rho 1}{|}_{i1+1/2,k1}^n}{\delta 1}\end{array}\right)\end{array},$

其中，n表示时间步，i1和k1分别表示上层区域中ρ向和z向的空间位置，k₁和k₂为环路系数，c为光速，μ_r为磁导率；

步骤(6.2)、根据步骤1中模型文件中定义的下层区域(不包括PML层)采用传统FDTD方法，对该区域内的φ向磁场分量H_φ2进行更新，具体更新公式为：

$\left(\begin{array}{c}{H}_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n+1/2}=H_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n-1/2}\\ +\frac{c\Delta t}{{\mu }_r{|}_{i2+1/2,k2+1/2}^{n+1/2}}·\left(\begin{array}{c}\frac{E_{z2}{|}_{i2+1,k2+1/2}^n-E_{z2}{|}_{i2,k2+1/2}^n}{\delta 2}\\ -\frac{E_{\rho 2}{|}_{i2+1/2,k2+1}^n-E_{\rho 2}{|}_{i2+1/2,k2}^n}{\delta 2}\end{array}\right),\end{array}\right)$

其中，n表示时间步，i2和k2分别表示下层区域中ρ向和z向的空间位置，c为光速，μ_r为磁导率；

步骤(6.3)、边界区域的场量传递：

边界区域场量的传递分为如下2种情况：

a、当上层粗网格中H_φ1的计算需用到下层细网格中的E_ρ²时，由于上下粗细网格比为奇数，场量重叠，故直接取相应细网格的场量即可，同时粗网格中H_φ1的更新不包括交界线上的点；

b、当交界线上的H_φ均取细网格的H_φ2，且交界线上H_φ2的计算需用到上层细网格中的E_ρ1时，采用线性插值的方法获取，具体过程如下：

的差分公式为：

$\begin{array}{c}{H}_{\phi 2}{|}_{i2+1,k2+3}^{n+1/2}=H_{\phi 3}{|}_{i2+1,k2+3}^{n-1/2}\\ +\frac{c\Delta t}{{\mu }_r{|}_{i2+1,k2+3}}·\left(\begin{array}{c}\frac{E_{z2}{|}_{i2+3/2,k2+3}^n-E_{z2}{|}_{i2+1/2,k2+3}^n}{\delta 2}\\ -\frac{E_{\rho 2}{|}_{i2+1,k2+7/2}^n-E_{\rho 2}{|}_{i2+1,k2+5/2}^n}{\delta 2}\end{array}\right)\end{array},$

式中，根据粗网格中的电场分量和采用线性插值的方法得到即

$E_{\rho 2}{|}_{i2+1,k2+7/2}^n=\frac{1}{3}{E}_{\rho 1}{|}_{i1+1,k1+1/2}^n+\frac{2}{3}{E}_{\rho 1}{|}_{i1,k1+1/2}^n;$

步骤(6.4)、对CFS-PML中的电场分量H_φ进行更新：吸收边界中的场分量H_φ的差分公式为：

$\begin{array}{c}{H}_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}=H_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}\\ -\frac{c\Delta t}{{\mu }_r{|}_{i+1/2,k+1/2}^{n+1/2}}[{\psi }_{\rho z}{|}_{i+1/2,k}^{n-1/2}+f_{1z}\left(\frac{E_{\rho }^{*}{|}_{i+1/2,k+1/2}^{n+1/2}-E_{\rho }^{*}{|}_{i+1/2,k-1/2}^{n+1/2}}{\delta }\right)\\ -{\psi }_{z\rho }{|}_{i,k+1/2}^{n-1/2}+f_{1\rho }\left(\frac{E_z^{*}{|}_{i+1/2,k+1/2}^{n+1/2}-E_z^{*}{|}_{i-1/2,k+1/2}^{n+1/2}}{\delta }\right)]\end{array},$

式中，

$\left(\begin{array}{c}{\psi }_{\rho z}{|}_{i+1/2,k}^{n-1/2}=f_{0z}{\tilde{E}}_{\rho z}{|}_{i+1/2,k}^{n-1/2}+f_{2z}\left(\frac{E_{\rho }^{*}{|}_{i+1/2,k+1/2}^{n-1/2}-E_{\rho }^{*}{|}_{i+1/2,k-1/2}^{n-1/2}}{\delta }\right)\\ {\psi }_{z\rho }{|}_{i,k+1/2}^{n-1/2}=f_{0\rho }{\tilde{E}}_{z\rho }{|}_{i,k+1/2}^{n-1/2}+f_{2\rho }\left(\frac{E_z^{*}{|}_{i+1/2,k+1/2}^{n-1/2}-E_z^{*}{|}_{i-1/2,k+1/2}^{n-1/2}}{\delta }\right)\end{array}\right)$

$\left(\begin{array}{cc}\begin{array}{c}{f}_{0\eta }=\frac{2{ϵ}_0{k}_{\eta }-\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}{2{ϵ}_0{k}_z+\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}\\ f_{1\eta }=\frac{{\alpha }_{\eta }\Delta t+2{ϵ}_0}{2{ϵ}_0{k}_{\eta }+\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}\\ f_{2\eta }=\frac{{\alpha }_{\eta }\Delta t-2{ϵ}_0}{2{ϵ}_0{k}_{\eta }+\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}\end{array}& (\eta =\rho ,z)\end{array}\right),$

其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，α_η、k_η和σ_η为吸收边界参数。

步骤7具体为：对激励源的更新过程需根据步骤1中模型文件中所设置的源的种类、激励场形式、类型，对所有的源进行更新。

步骤8具体为：

判断是否达到总的迭代次数，若达到则进行步骤9；否则时间步加1，并返回步骤4。

步骤9具体为：

a、当激励源为单频正弦信号时：

根据峰值监测法提取观测点信号垂直方向电场幅值|E_z0|，则当天线辐射功率为P_t时，场强|E_z|为

$\left|E_z\right|=\left|E_{z0}\right|\sqrt{\frac{P_t}{160{\pi }^2}{\left(\frac{\lambda }{dl}\right)}^2},$

提取观测点处某一时域稳态周期内垂直电场波形过零点时刻T_s，则观测点的传播时延t_w为

$t_w=T_s-\frac{d}{c}-t_0+2.5\times {10}^{-6},$

式中，d为传播距离，c为光速，t₀为天线所加电流源中与所提取周期时刻相对应的参考时刻；

b、当激励源为Loran-C信号时：

提取观测点垂直方向电场波形中的第三载波周期负峰值场强，即采样点场强和峰值场强则当天线辐射功率为P_t时，采样点场强和峰值场强分别为：

$\left|E_z^c\right|=\frac{\left|E_{z0}^c\right|}{B}\sqrt{\frac{P_t}{160{\pi }^2}{\left(\frac{\lambda }{dl}\right)}^2},$

$\left|E_z^f\right|=\frac{\left|E_{z0}^f\right|}{B}\sqrt{\frac{P_t}{160{\pi }^2}{\left(\frac{\lambda }{dl}\right)}^2},$

式中，B＝A×(65×10^-6)²×e^-2≈5.717976×10^-10A，

提取第三载波周期正向过零点时刻T_s，观测点的传播时延t_w为

$t_W=T_s-\frac{d}{c}-3T+2.4\times {10}^{-6},$

式中，T为周期。

本发明的有益效果是，一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法，首先，“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法综合考虑了低频地波传播时延预测中整个电波传播区域中传播路径的特征：①包含空气层和地表层；②空气层对低频地波传播的影响主要为大气折射率的空时分布所产生的影响，而大气折射率的空时变化较缓慢，所产生的影响也较小；③地表层对低频地波传播的影响主要为地形起伏变化、地面电特性(大地电导率σ和相对介电常数ε_r)所产生的影响，而地形和地物的变化较复杂，所产生的影响也较大。因此，对整个电波传播区域进行分层处理，即分为上(空气层)下(地表层)两层；其次，“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法综合利用了M(2，4)FDTD方法和传统FDTD方法的特点：①传统FDTD方法可用于处理任何复杂结构/媒质的电波传播问题，但是网格剖分密度越小，数值色散误差越大；而网格剖分密度过大，则计算机消耗(内存占用、计算时长)越大；②M(2，4)FDTD通过修正系数的引进，可大大降低传统FDTD方法的数值色散误差，但不适用于复杂路径的电波传播问题。因此，将两种方法结合，发挥其各自优势，对上层空气层采用粗网格的M(2，4)FDTD进行传播时延预测，在保证预测精度的同时减少计算机内存占用、提高计算速度；对下层地表层采用细网格的传统FDTD方法进行传播时延预测，保证预测结果的可靠性和准确性；最后，巧妙利用亚网格技术，实现上下区域间场量的传递；利用CFS-PML技术，实现电磁波的有效吸收，保证预测结果的精度。

附图说明

图1是“本发明一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法原理示意图；

图2是本发明一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法中M(2，4)FDTD方法的网格示意图；

图3是本发明一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法中计算区域上下层网格排布及场量分布示意图；

图4是本发明一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法的流程图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法，理论基础及原理是：

针对低频地波传播时延预测问题，一方面希望达到较高的时延预测精度；另一方面要在有限的资源和一定的时间范围内实现长距离、复杂传播路径问题的求解。为了兼顾这两方面需求，解决网格剖分密度与计算机资源占用及耗时之间的矛盾，如图1所示，将M(2，4)FDTD方法与传统FDTD方法进行结合，以实现低频地波传播时延的高精度预测，同时，提高计算速度、减少计算机内存占用，其中，M(2，4)FDTD方法网格示意图如图2所示，包含C1、C2和C3三个环路，并引入环路系数k₁和k₂；计算区域上下层网格排布及场量分布示意图如图3所示，设定上下粗细网格比为奇数，便于场量的传递；流程图如图4所示，具体步骤如下：

步骤1、设置并生成模型文件、输入模型文件；

步骤2、参数设置及初始化；

步骤3、添加场源；

步骤4、更新整个计算区域的ρ向电场分量E_ρ；

步骤5、更新整个计算区域的z向电场分量E_z；

步骤6、更新整个计算区域的φ向磁场分量H_φ；

步骤7、更新激励源；

步骤8、判断结束条件，循环；

步骤9、观测点观测量计算并输出。

其中，每一步骤具体做法如下：

步骤1具体为：

步骤(1.1)、设模型文件的上层区域大小N_ρ1×N_z1，下层区域大小N_ρ2×N_z2，CFS-PML层数为N_PML，其中N_ρ为ρ方向网格数，N_z为z方向网格数，标号1和2分别表示上层和下层；

步骤(1.2)、设置空间、时间步长：上层粗网格空间步长为δ1，其中，Δρ1＝Δz1＝δ1，下层细网格空间步长为δ2，其中，Δρ2＝Δz2＝δ2，其中Δρ为ρ方向网格大小，Δz为z方向网格大小，标号1和2分别表示上层和下层，时间步长为Δt，粗网格M(2，4)FDTD的时间步长设置为与细网格FDTD同样的时间步长；

步骤(1.3)、设置迭代次数为NT_Calc；

步骤(1.4)、设置传播路径电参数：ρ方向起始网格位置为ρ_Start、ρ方向结束网格位置为ρ_Een、z方向起始网格位置为z_Start、z方向结束网格位置为z_Start、大地电导率为σ、相对介电常数为ε_r；

步骤(1.5)、设置吸收边界：CFS-PML层数为N_PML，相关参数为κ_ηmax，α_ηmax，σ_ηmax，其中η＝ρ,z；

步骤(1.6)、设置场源：源的个数N_Source、位置[ρ_Start，ρ_Een]和[z_Start，z_Een]，源的种类：有两种激励源供选择：单频正弦源、Loran_C源，激励场形式： 有三种激励形式：E_z、E_ρ、源的类型：软源或硬源，幅度A，频率f，单频正弦源中的常数t₀，高斯脉宽T_p，时延/包周差τ；

步骤(1.7)、设置观测点：观测点个数N_VPoint，位置[ρ_Start，ρ_Een]和[z_Start，z_Een]，输出场量类型(E_z、E_ρ或H_φ)。

步骤2具体为：

步骤(2.1)、将整个区域的电磁场分量(E_z、E_ρ和H_φ)和电磁场分量系数(CA、CB)、中间变量中间变量系数(f_0z、f_1z、f_2z，f_0ρ、f_1ρ、f_2ρ，f_0φγ、f_1φγ、f_2φγ，)、辅助变量(ψ_φz、ψ_φρ、ψ_φγ、ψ_ρz、ψ_zρ)、观测量(传播时延t_w、采样点场强峰值场强)都初始化为零；

步骤(2.2)、初始化所有网格中的模型参数：将相对介电常数ε_r初始化为1，大地电导率σ初始化为0；

步骤(2.3)、根据步骤1所设置的模型文件中的路径信息，给相应网格的ε_r和σ赋值；

步骤(2.4)设置CFS-PML吸收边界参数k_η、σ_η、α_η，其中，k_η、σ_η、α_η具体按下式计算：

$\left(\begin{array}{c}{k}_{\eta }=1+(k_{\eta max}-1)\frac{|\eta -{\eta }_0{|}^m}{d^m}\\ {\sigma }_{\eta }=\frac{{\sigma }_{\eta max}|\eta -{\eta }_0{|}^m}{d^m}\\ {\alpha }_{\eta }={\alpha }_{\eta max}\frac{|\eta -{\eta }_0|}{d}\end{array}\right)$

式中，η＝ρ,z，η₀为PML层与非PML截面位置，d是PML吸收边界的厚度，κ_ηmax取整数，κ_ηmax取值范围为[1，60]，α_ηmax取值范围为[0，1)， σ_ηmax根据σ_opt设置，σ_ηmax/σ_opt取值范围为(0，12]，σ_opt＝(m+1)/150πΔη，m取值范围为[1，20]，其中m取值为4时边界的吸收效果最好，Δη取值范围λ为源的波长。

步骤3具体为：

添加的场源有两种类型：

(1)单频正弦源：

当辐射源采用单频正弦信号时，其电流源激励i_s(t)可以表示为：

$i_s\left(t\right)=\left(\begin{array}{cc}0& t<0\\ 0.5[1-\cos \left(\frac{\pi t}{t_0}\right)]\sin \left(2\pi ft\right)& 0\le t\le t_0\\ \sin \left(2\pi ft\right)& t_0\le t\end{array}\right),$

其中t₀＝5×10^-6s。

(2)Loran-C源：

当辐射源采用正相位编码的Loran_C信号时，其电流波形前沿i_s(t)为：

$i_s\left(t\right)=\left(\begin{array}{cc}0& t<\tau \\ A{(t-\tau )}^2{e}^{-\frac{2(t-\tau )}{65\times {10}^{-6}}}sin\left(2\pi ft\right)& \tau \le t\le 65\times {10}^{-6}+\tau \end{array}\right),$

其中，τ为包周差，单位为s。

步骤4具体为：

步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量E_ρ1进行更新，其中，上层区域不包含PML层，E_ρ1中的下标1表示上层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{\rho 1}{|}_{i1+1/2,k1}^{n+1}=CA\left(m\right)·E_{\rho 1}{|}_{i1+1/2,k1}^n\\ -CB\left(m\right)·\left(\begin{array}{c}{k}_1\left(\frac{H_{\phi 1}{|}_{i+1/2,k1+3/2}^{n+1/2}-H_{\phi 1}{|}_{i1+1/2,k1-3/2}^{n+1/2}}{3\delta 1}\right)\\ -(1-k_1)\left(\frac{H_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n+1/2}-H_{\phi 1}{|}_{i1+1/2,k1-1/2}^{n+1/2}}{\delta 1}\right)\end{array}\right)\end{array}\right)$

$\left(\begin{array}{c}CA\left(m\right)=\frac{\frac{{ϵ}_r\left(m\right)}{\Delta t}-\frac{\sigma \left(m\right)}{2{ϵ}_0}}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\\ CB\left(m\right)=\frac{c}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\end{array}\right)$

式中，n表示时间步，标号i1和k1分别表示上层区域中ρ向和z向的空间位置，ε₀为真空中介电常数，k₁为环路系数；

步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区域内的电磁场分量E_ρ2进行更新，其中，下层区域不包含PML层，E_ρ2中的下标2表示下层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{\rho 2}{|}_{i+1/2,k2}^{n+1}=CA\left(m\right)·E_{\rho 2}{|}_{i2+1/2,k2}^n\\ -CB\left(m\right)·\left(\frac{H_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n+1/2}-H_{\phi 2}{|}_{i2+1/2,k2-1/2}^{n+1/2}}{\delta 2}\right)\end{array}\right)$

$\left(\begin{array}{c}CA\left(m\right)=\frac{\frac{{ϵ}_r\left(m\right)}{\Delta t}-\frac{\sigma \left(m\right)}{2{ϵ}_0}}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\\ CB\left(m\right)=\frac{c}{\frac{{ϵ}_r\left(m\right)}{\Delta t}+\frac{\sigma \left(m\right)}{2{ϵ}_0}}\end{array}\right)$

式中，n表示时间步，标号i2和k2分别表示下层区域 中ρ向和z向的空间位置；

步骤(4.3)、对边界区域的场量传递进行更新：

边界区域场量的传递分为如下几种情况，如图3所示：

a、当上层粗网格中E_ρ1的计算需用到下层细网格中的H_φ2时，由于上下粗细网格比为奇数，场量重叠，故直接取相应细网格的场量即可；

b、当下层细网格中边界网格处E_ρ2的计算需用到边界线上的H_φ时，直接采用边界线上细网格的H_φ2即可；

步骤(4.4)、对CFS-PML中的电场分量进行更新：

上方吸收边界中的网格空间大小与粗网格相同，右侧吸收边界中网格的空间大小与相邻左侧计算区域的网格空间大小相同，也分为上下两层，上层与粗网格一致，下层与细网格一致，

吸收边界中的场分量的差分公式为：

$\left(\begin{array}{c}{E}_{\rho }^{*}{|}_{i+1/2,k}^{n+1}={\mathrm{CA}}_{\rho }\left(m\right)E_{\rho }^{*}{|}_{i+1/2,k}^n\\ +{\mathrm{CB}}_{\rho }\left(m\right)[{\psi }_{\phi z}{|}_{i+1/2,k}^{n-1/2}+f_{1z}{|}_k\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\phi }{|}_{i+1/2,k-1/2}^{n+1/2}}{\Delta z}\right)]\end{array}\right)$

其中，

$\left(\begin{array}{c}{\psi }_{\phi z}{|}_{i+1/2,k}^{n-1/2}=f_{0z}{|}_k{\tilde{H}}_{\phi z}{|}_{i+1/2,k}^{n-1/2}+f_{2z}{|}_k\frac{\partial }{\partial z}{H}_{\phi }{|}^{n-1/2}\\ =f_{0z}{|}_k{\tilde{H}}_{\phi z}{|}_{i+1/2,k}^{n-1/2}+f_{2z}{|}_k\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n-1/2}-H_{\phi }{|}_{i+1/2,k-1/2}^{n-1/2}}{\Delta z}\right)\end{array}\right)$

$\left(\begin{array}{c}{f}_{0z}{|}_k=\frac{2{ϵ}_0{k}_z{|}_k-\Delta t(k_z{|}_k·{\alpha }_z{|}_k+{\sigma }_z{|}_k)}{2{ϵ}_0{k}_z{|}_k+\Delta t(k_z{|}_k·{\alpha }_z+{\sigma }_z{|}_k)}\\ f_{1z}{|}_k=\frac{{\alpha }_z{|}_k\Delta t+2{ϵ}_0}{2{ϵ}_0{k}_z{|}_k+\Delta t(k_z{|}_k·\alpha {|}_k+{\sigma }_z{|}_k)}\\ f_{2z}{|}_k=\frac{{\alpha }_z{|}_k\Delta t-2{ϵ}_0}{2{ϵ}_0{k}_z{|}_k+\Delta t(k_z{|}_k·{\alpha }_z{|}_k+{\sigma }_z{|}_k)}\end{array}\right),$

其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，k_z、σ_z和α_z为吸收边界参数。

步骤5具体为：

步骤(5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量E_z1进行更新，其中，上层区域不包含PML层，E_z1中的下标1表示上层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{z1}{|}_{i1,k1+1/2}^{n+1}=CA\left(m\right)·E_{z1}{|}_{i1,k1+1/2}^n\\ +CB\left(m\right)·\left(\begin{array}{c}{k}_1\frac{(i1+\frac{1}{2})H_{\phi 1}{|}_{i+13/2,k1+1/2}^{n+1/2}-(i1-\frac{1}{2})H_{\phi 1}{|}_{i1-3/2,k1+1/2}^{n+1/2}}{3i1\delta 1}\\ +(1+k_1)\frac{(i1+\frac{1}{2})H_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n+1/2}-(i1-\frac{1}{2})H_{\phi 1}{|}_{i1-1/2,k1+1/2}^{n+1/2}}{i1\delta 1}\end{array}\right)\end{array}\right)$

式中，n表示时间步，标号i1和k1分别表示上层区域中ρ向和z向的空间位置，ε₀为真空中介电常数，k₁为环路系数；

步骤(5.2)、根据步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区域内的电磁场分量E_z2进行更新，其中，下层区域不包含PML层，E_z2中下标2表示下层区域，具体更新公式为：

$\left(\begin{array}{c}{E}_{z2}{|}_{i2,k2+1/2}^{n+1}=CA\left(m\right)·E_{z2}{|}_{i2,k2+1/2}^n\\ +CB\left(m\right)·\left[\frac{(i2+\frac{1}{2})H_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n+1/2}-(i2-\frac{1}{2})H_{\phi 2}{|}_{i2-1/2,k2+1/2}^{n+1/2}}{i2\delta 2}\right]\end{array}\right)$

式中，n表示时间步，标号i2和k2分别表示下层区域中ρ向和z向的空间位置；

步骤(5.3)、如图3所示，交界线上E_z均为细网格的E_z2，上层粗网格中E_z1的更新不包括交界线上的点；

步骤(5.4)、对CFS-PML中的电场分量进行更新：吸收边界中的场分量的差分公式为：

$\left(\begin{array}{c}{E}_z^{*}{|}_{i,k+1/2}^{n+1}=CA\left(m\right)E_z^{*}{|}_{i,k+1/2}^n\\ +CB\left(m\right)\left(\begin{array}{c}[{\psi }_{\phi \rho }{|}_{i,k+1/2}^{n-1/2}+f_{1\rho }{|}_i\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}-H_{\phi }{|}_{i-1/2,k+1/2}^{n+1/2}}{\delta }\right)]\\ +[{\psi }_{\phi \gamma }{|}_{i,k+1/2}^{n-1/2}+f_{1\phi \gamma }{|}_i\frac{H_{\phi }}{i\delta }{|}_{i,k+1/2}^{n+1/2}]\end{array}\right)\end{array}\right)$

式中，标号m＝(i,k+1/2)，且

$\left(\begin{array}{c}{\psi }_{\phi \rho }{|}_{i,k+1/2}^{n-1/2}=f_{0\rho }{|}_i{\tilde{H}}_{\phi \rho }{|}_{i,k+1/2}^{n-1/2}+f_{2\rho }{|}_i\left(\frac{H_{\phi }{|}_{i+1/2,k+1/2}^{n-1/2}-H_{\phi }{|}_{i-1/2,k+1/2}^{n-1/2}}{\Delta \rho }\right)\\ {\psi }_{\phi \gamma }{|}_{i,k+1/2}^{n-1/2}=f_{0\phi \gamma }{|}_i{\tilde{H}}_{\phi \gamma }{|}_{i,k+1/2}^{n-1/2}+f_{2\phi \gamma }{|}_i\frac{H_{\phi }}{i\delta }{|}_{i,k+1/2}^{n-1/2}\end{array}\right),$

$\left(\begin{array}{c}{f}_{0\rho }{|}_i=\frac{2{ϵ}_0{k}_{\rho }{|}_i-\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}{2{ϵ}_0{k}_{\rho }{|}_i+\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}\\ f_{1\rho }{|}_i=\frac{{\alpha }_{\rho }{|}_i\Delta t+2{ϵ}_0}{2{ϵ}_0{k}_{\rho }{|}_i+\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}\\ f_{2\rho }{|}_i=\frac{{\alpha }_{\rho }{|}_i\Delta t-2{ϵ}_0}{2{ϵ}_0{k}_{\rho }{|}_i+\Delta t(k_{\rho }{|}_i{\alpha }_{\rho }{|}_i+{\sigma }_{\rho }{|}_i)}\end{array}\right)$

$\left(\begin{array}{c}{f}_{0\phi \gamma }{|}_i=-\frac{\Delta t(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)-2{ϵ}_0{A}_{\rho }{|}_i}{\Delta t(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)+2{ϵ}_0{A}_{\rho }{|}_i}\\ f_{1\phi \gamma }{|}_i=\frac{({\mathrm{\Delta t\alpha }}_{\rho }{|}_i+2{ϵ}_0)}{\left(\Delta t\right(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)+2{ϵ}_0{A}_{\rho }{|}_i)}\\ f_{2\phi \gamma }{|}_i=\frac{({\mathrm{\Delta t\alpha }}_{\rho }{|}_i-2{ϵ}_0)}{\left(\Delta t\right(A_{\rho }{|}_i·{\alpha }_{\rho }{|}_i+B_{\rho }{|}_i)+2{ϵ}_0{A}_{\rho }{|}_i)}\end{array}\right)$

$\left(\begin{array}{c}{A}_{\rho }=1+(k_{\rho \max }-1)\frac{{(\rho -{\rho }_0)}^{m+1}}{(m+1)d^m}\frac{1}{\rho }\\ B_{\rho }=\frac{1}{{\alpha }_{\rho }+{\mathrm{j\omega ϵ}}_0}\frac{{\sigma }_{\rho \max }{(\rho -{\rho }_0)}^{m+1}}{(m+1)d^m}\frac{1}{\rho }\end{array}\right)$

其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，α_ρ、k_ρ、σ_ρ、k_ρmax和σ_ρmax为吸收边界参数，ρ＝i×δ表示ρ向长度，ρ₀是PML层与非PML层的分界位置。

步骤6具体为：

步骤(6.1)、首先根据步骤1中模型文件中定义的上层区域采用粗网格M(2，4)FDTD方法，对该区域内的电磁场分量进行更新，其中，上层区域不包含PML层，具体更新公式为：

$\begin{array}{c}{H}_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n+1/2}=H_{\phi 1}{|}_{i1+1/2,k1+1/2}^{n-1/2}\\ +\frac{c\Delta t}{{\mu }_r{|}_{i1+1/2,k1+1/2}^{n+1/2}}·\left(\begin{array}{c}{k}_1\frac{E_{z1}{|}_{i1+2,k1+1/2}^n-E_{z1}{|}_{i1-1,k1+1/2}^n}{3\delta 1}\\ +k_2\frac{E_{z1}{|}_{i1+2,k1+3/2}^n+E_{z1}{|}_{i1+2,k1-1/2}^n-E_{z1}{|}_{i1-1,k1+3/2}^n-E_{z1}{|}_{i1-1,k1-1/2}^n}{6\delta 1}\\ +(1-k_1-k_2)\frac{E_{z1}{|}_{i1+1,k1+1/2}^n-E_{z1}{|}_{i1,k1+1/2}^n}{\delta 1}\end{array}\right)\\ -\frac{c\Delta t}{{\mu }_r{|}_{i1+1/2,k1+1/2}^{n+1/2}}·\left(\begin{array}{c}{k}_1\frac{E_{\rho 1}{|}_{i1+1/2,k1+2}^n-E_{\rho 1}{|}_{i1+1/2k1-1}^n}{3\delta 1}\\ +k_2\frac{E_{\rho 1}{|}_{i1+3/2,k1+2}^n+E_{\rho 1}{|}_{i1-1/2,k1+2}^n-E_{\rho 1}{|}_{i1+3/2,k1-1}^n-E_{\rho 1}{|}_{i1-1/2,k1-1}^n}{6\delta 1}\\ +(1+k_1-k_2)\frac{E_{\rho 1}{|}_{i1+1/2,k1+1}^n-E_{\rho 1}{|}_{i1+1/2,k1}^n}{\delta 1}\end{array}\right)\end{array},$

其中，n表示时间步，i1和k1分别表示上层区域中ρ向和z向的空间位置，k₁和k₂为环路系数，c为光速，μ_r为磁导率；

步骤(6.2)、根据步骤1中模型文件中定义的下层区域(不包括PML层)采用传统FDTD方法，对该区域内的φ向磁场分量H_φ2进行更新，具体更新公式为：

$\left(\begin{array}{c}{H}_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n+1/2}=H_{\phi 2}{|}_{i2+1/2,k2+1/2}^{n-1/2}\\ +\frac{c\Delta t}{{\mu }_r{|}_{i2+1/2,k2+1/2}^{n+1/2}}·\left(\begin{array}{c}\frac{E_{z2}{|}_{i2+1,k2+1/2}^n-E_{z2}{|}_{i2,k2+1/2}^n}{\delta 2}\\ -\frac{E_{\rho 2}{|}_{i2+1/2,k2+1}^n-E_{\rho 2}{|}_{i2+1/2,k2}^n}{\delta 2}\end{array}\right),\end{array}\right)$

其中，n表示时间步，i2和k2分别表示下层区域中ρ向和z向的空间位置，c为光速，μ_r为磁导率；

步骤(6.3)、边界区域的场量传递：

如图3所示，边界区域场量的传递分为如下2种情况：

a、当上层粗网格中H_φ1的计算需用到下层细网格中的E_ρ2时，由于上下粗细网格比为奇数，场量重叠，故直接取相应细网格的场量即可，同时粗网格中H_φ1的更新不包括交界线上的点；

b、当交界线上的H_φ均取细网格的H_φ2，且交界线上H_φ2的计算需用到上层细网格中的E_ρ1时，采用线性插值方法获取，具体过程如下：

的差分公式为：

$\begin{array}{c}{H}_{\phi 2}{|}_{i2+1,k2+3}^{n+1/2}=H_{\phi 3}{|}_{i2+1,k2+3}^{n-1/2}\\ +\frac{c\Delta t}{{\mu }_r{|}_{i2+1,k2+3}}·\left(\begin{array}{c}\frac{E_{z2}{|}_{i2+3/2,k2+3}^n-E_{z2}{|}_{i2+1/2,k2+3}^n}{\delta 2}\\ -\frac{E_{\rho 2}{|}_{i2+1,k2+7/2}^n-E_{\rho 2}{|}_{i2+1,k2+5/2}^n}{\delta 2}\end{array}\right)\end{array},$

式中，根据粗网格中的电场分量和采用线性插值的方法得到即

$E_{\rho 2}{|}_{i2+1,k2+7/2}^n=\frac{1}{3}{E}_{\rho 1}{|}_{i1+1,k1+1/2}^n+\frac{2}{3}{E}_{\rho 1}{|}_{i1,k1+1/2}^n;$

步骤(6.4)、对CFS-PML中的电场分量H_φ进行更新：吸收边界中的场分量H_φ的差分公式为：

$\begin{array}{c}{H}_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}=H_{\phi }{|}_{i+1/2,k+1/2}^{n+1/2}\\ -\frac{c\Delta t}{{\mu }_r{|}_{i+1/2,k+1/2}^{n+1/2}}[{\psi }_{\rho z}{|}_{i+1/2,k}^{n-1/2}+f_{1z}\left(\frac{E_{\rho }^{*}{|}_{i+1/2,k+1/2}^{n+1/2}-E_{\rho }^{*}{|}_{i+1/2,k-1/2}^{n+1/2}}{\delta }\right)\\ -{\psi }_{z\rho }{|}_{i,k+1/2}^{n-1/2}+f_{1\rho }\left(\frac{E_z^{*}{|}_{i+1/2,k+1/2}^{n+1/2}-E_z^{*}{|}_{i-1/2,k+1/2}^{n+1/2}}{\delta }\right)\end{array},$

式中，

$\left(\begin{array}{c}{\psi }_{\rho z}{|}_{i+1/2,k}^{n-1/2}=f_{0z}{\tilde{E}}_{\rho z}{|}_{i+1/2,k}^{n-1/2}+f_{2z}\left(\frac{E_{\rho }^{*}{|}_{i+1/2,k+1/2}^{n-1/2}-E_{\rho }^{*}{|}_{i+1/2,k-1/2}^{n-1/2}}{\delta }\right)\\ {\psi }_{z\rho }{|}_{i,k+1/2}^{n-1/2}=f_{0\rho }{\tilde{E}}_{z\rho }{|}_{i,k+1/2}^{n-1/2}+f_{2\rho }\left(\frac{E_z^{*}{|}_{i+1/2,k+1/2}^{n-1/2}-E_z^{*}{|}_{i-1/2,k+1/2}^{n-1/2}}{\delta }\right)\end{array}\right)$

$\left(\begin{array}{cc}\begin{array}{c}{f}_{0\eta }=\frac{2{ϵ}_0{k}_{\eta }-\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}{2{ϵ}_0{k}_z+\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}\\ f_{1\eta }=\frac{{\alpha }_{\eta }\Delta t+2{ϵ}_0}{2{ϵ}_0{k}_{\eta }+\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}\\ f_{2\eta }=\frac{{\alpha }_{\eta }\Delta t-2{ϵ}_0}{2{ϵ}_0{k}_{\eta }+\Delta t(k_{\eta }{\alpha }_{\eta }+{\sigma }_{\eta })}\end{array}& (\eta =\rho ,z)\end{array}\right),$

其中，n表示时间步，i和k分别表示计算区域中ρ向和z向的空间位置，α_η、k_η和σ_η为吸收边界参数。

步骤7具体为：对激励源的更新过程需根据步骤1中模型文件中所设置的源的种类、激励场形式、类型，对所有的源进行更新。

步骤8具体为：

判断是否达到总的迭代次数，若达到则进行步骤9；否则时间步加1，并返回步骤4。

步骤9具体为：

a、当激励源为单频正弦信号时：

根据峰值监测法提取观测点信号垂直方向电场幅值|E_z0|，则当天线辐射功率为P_t时，场强|E_z|为

$\left|E_z\right|=\left|E_{z0}\right|\sqrt{\frac{P_t}{160{\pi }^2}{\left(\frac{\lambda }{dl}\right)}^2},$

提取观测点处某一时域稳态周期内垂直电场波形过零点时刻T_s，则观测点的传播时延t_w为

$t_w=T_s-\frac{d}{c}-t_0+2.5\times {10}^{-6},$

式中，d为传播距离，c为光速，t₀为天线所加电流源中与所提取周期时刻相对应的参考时刻；

b、当激励源为Loran-C信号时：

提取观测点垂直方向电场波形中的第三载波周期负峰值场强，即采样点场强和峰值场强则当天线辐射功率为P_t时，采样点场强和峰值场强分别为：

$\left|E_z^c\right|=\frac{\left|E_{z0}^c\right|}{B}\sqrt{\frac{P_t}{160{\pi }^2}{\left(\frac{\lambda }{dl}\right)}^2},$

$\left|E_z^f\right|=\frac{\left|E_{z0}^f\right|}{B}\sqrt{\frac{P_t}{160{\pi }^2}{\left(\frac{\lambda }{dl}\right)}^2},$

式中，B＝A×(65×10^-6)²×e^-2≈5.717976×10^-10A，

提取第三载波周期正向过零点时刻T_s，观测点的传播时延t_w为

$t_W=T_s-\frac{d}{c}-3T+2.4\times {10}^{-6},$

式中，T为周期。