一种“m(2，4)fdtd+fdtd”的低色散低频地波传播时延预测方法
【专利摘要】本发明公开了一种“M(2，4)FDTD+FDTD”的低色散低频地波传播时延预测方法，该方法通过对整个计算区域进行分层处理，将粗网格M(2，4)FDTD算法(用于上层区域)与细网格传统FDTD算法(用于下层区域)相结合，利用修正因子以及亚网格技术，进行低频地波传播时延的高精度快速预测，本发明解决了传统全区域FDTD算法在减小数值色散误差与降低计算机消耗之间的矛盾，提高了低频地波传播时延预测精度的同时，减少了计算机内存资源占用，提高了计算速度。
【专利说明】
一种"Μ (2，4) FDTD+FDTD"的低色散低频地波传播时延预测 方法
技术领域
[0001] 本发明属于电波传播理论和电磁场数值计算技术领域，具体涉及一种"M(2,4) FDTD+FDTD"的低色散低频地波传播时延预测方法。
【背景技术】
[0002] 低频地波传播时延预测精度是提高陆基无线电定位/导航/授时系统精度的关键。 现有的经典理论预测方法主要有:基于均匀光滑地面模型的平地面公式和Fock地波绕射公 式等；基于分段不均匀光滑地面模型的Millington经验公式、Wait积分公式、波模转换法 等;基于不均匀不光滑地面模型的积分方程方法、抛物方程方法等。然而，一方面，经典的理 论预测方法均是在一定的模型近似及理论假设条件下得到的，使用时无法考虑实际传播中 的各种因素的复杂变化(如复杂路径的后向反射波的影响、地形横向剧烈变化的影响、大气 折射率空时变化的影响等）；另一方面，经典理论预测方法多是针对单一频率的频域预测结 果，而对经过调制的脉冲信号(如罗兰C信号），计算距离越长、结果误差越大。
[0003] 近些年，随着计算机技术的不断提高以及数值计算理论和方法的进一步发展，时 域有限差分方法(Finite-Difference Time-Domain Method，FDTD)被用于该领域来进行低 频地波传播时延的高精度预测。它从Maxwell旋度方程出发，利用二阶精度的中心差分近 似，将旋度方程中的时间及空间微分算符直接转换为差分形式，从而在有限的空间和时间 上对连续的电磁场数据进行抽样。与经典的地波传播时延理论预测方法相比较，FDTD方法 具有较高的预测精度，但是，随着传播距离的增长以及路径复杂度的增大，FDTD方法的预测 精度受到显著影响，究其原因主要是H)TD方法存在数值色散误差。理论上，通过进一步增加 网格剖分密度或采用更高阶的H)TD算法等措施可以减小其数值色散误差，但这就意味着需 要占用更大的计算机内存资源、消耗更长的计算时间，难以用于长距离、大区域甚至三维 FDTD在低频地波传播时延预测中的工程推广。
[0004] 为了解决传统roTD方法数值色散问题，国内外学者分别从时间离散方式和空间差 分格式出发对传统roTD方法进行了一系列的改进，其中M(2,4)roTD是Hadi和Piket-May从 Maxwe 11方程的积分形式出发，结合Fang的H)TD (2，4)算法提出的一种高阶H)TD方法，通过 引入两个修正系数k#Pk2，可降低特定频率所有方向的数值色散。但是，该方法无法直接用 于非均匀媒质。
[0005] 目前，尚未见M(2,4)H)TD与FDTD相结合方法在低频地波传播时延预测方面的研究 报道和专利。
【发明内容】
[0006] 本发明的目的是提供一种"M(2,4)FDTD+FDTD"的低色散低频地波传播时延预测方 法，解决了现有技术中存在的roTD数值方法在进行低频地波传播时延预测时，随传播距离 增长、路径复杂度提高而存在的"数值计算误差增大、计算速度缓慢、内存资源占用较大"的 问题。
[0007] 本发明所采用的技术方案是，一种"M(2,4)FDTD+roTD"的低色散低频地波传播时 延预测方法，通过对低色散低频地波计算区域进行分层处理，将粗网格M(2,4)FDTD与细网 格传统H)TD方法相结合，进行低频地波传播时延的预测，从而达到在保证预测精度的同时， 提高计算速度、减少计算机内存占用，具体步骤如下：
[0008] 步骤1、设置并生成模型文件、输入模型文件；
[0009] 步骤2、参数设置及初始化；
[0010] 步骤3、添加场源；
[0011] 步骤4、更新整个计算区域的P向电场分量Ep;
[0012] 步骤5、更新整个计算区域的z向电场分量Ez;
[0013] 步骤6、更新整个计算区域的Φ向磁场分量ΗΦ;
[0014] 步骤7、更新激励源；
[0015]步骤8、判断结束条件，循环；
[0016] 步骤9、观测点观测量计算并输出。
[0017]本发明的特点还在于，
[0018] 步骤1具体为：
[0019] 步骤（1.1)、设模型文件的上层区域大小NPl X Νζ1，下层区域大小Νρ2 X Nz2，CFS-PML 层数为Npml，其中ΝΡ为Ρ方向网格数具为ζ方向网格数，标号1和2分别表示上层和下层；
[0020] 步骤(1.2)、设置空间、时间步长:上层粗网格空间步长为δ 1，其中，Δ ρ 1 = Δ ζ 1 = δ 1，下层细网格空间步长为处，其中，Αρ2= Δζ2 = δ2,其中Δρ为ρ方向网格大小，ΔΖ为ζ方 向网格大小，标号1和2分别表示上层和下层，时间步长为Δ t，粗网格M( 2，4)H)TD的时间步 长设置为与细网格FDTD同样的时间步长；
[0021]步骤(1.3)、设置迭代次数为NTCalc;
[0022] 步骤（1.4)、设置传播路径电参数:P方向起始网格位置为pstart、P方向结束网格位 置为Pe?、ζ方向起始网格位置为zst art、ζ方向结束网格位置为zstart、大地电导率为σ、相对介 电常数为
[0023] 步骤(1.5)、设置吸收边界:CFS-PML层数为Npml，相关参数为icnmax，anmax，o nmax，其中η = P,z；
[0024] 步骤（1.6 )、设置场源:源的个数Nsouixe、位置[Pstart，PEen ]和[ZStart，ZEen]，源的种 类:有两种激励源供选择：单频正弦源、1^^811_(：源，激励场形式：有三种激励形式:E z、Ep、 源的类型:软源或硬源，幅度A，频率f，单频正弦源中的常数t〇,高斯脉宽!^，时延/包周 差丁；
[0025] 步骤（1 · 7)、设置观测点：观测点个数NVPciint，位置[pstart，pEen]和[z Start，zEen]，输出 场量类型(EZ、EP或Η Φ)。
[0026] 步骤2具体为：
[0027]步骤(2.1)、将整个区域的电磁场分量(Εζ、ΕΡ和ΗΦ)和电磁场分量系数(CA、CB)、中
、中间变量系数(fQz、flz、f 2z，fQp、flp、f2p，?Οφγ、 fl4)Y、f24)Y，）、辅助变直(Φφζ、Φφρ、Φφγ、Φρζ、Φζρ)、观测直(传播时延tw、米样点场强五；；。峰值 s 场强I)都初始化为零；
[0028] 步骤(2.2 )、初始化所有网格中的模型参数:将相对介电常数ε4刀始化为1，大地电 导率〇初始化为〇;
[0029] 步骤（2.3)、根据步骤1所设置的模型文件中的路径信息，给相应网格的^和〇赋 值；
[0030] 步骤(2.4)设置CFS-PML吸收边界参数kn、〇n、an，其中，knWi^ai^ic体按下式计算：
[0032]式中，n = p，Z，n〇为PML层与非PML截面位置，d是PML吸收边界的厚度，1^_取整数， Knmax取值范围为[1，60]，anmax取值范围为[0，1)，onmax根据σ。^设置，στ^χ/σ。^取值范围为(0， 12]，〇。的=( 111+1)/15031八11，111取值范围为[1，20]，其中111取值为4时边界的吸收效果最好，八11
[0033] 步骤3具体为：
[0034] 添加的场源有两种类型：
[0035] (1)单频正弦源：
[0036]当辐射源采用单频正弦信号时，其电流源激励。(〇可以表示为：
[0038] 其中 to = 5X10-6s。
[0039] (2)Loran_C 源：
[0040] 当辐射源采用正相位编码的Loran_C信号时，其电流波形前沿is(t)为：
[0042] 其中，τ为包周差，单位为s。
[0043] 步骤4具体为：
[0044]步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4) H)TD方法，对该区域内的电磁场分量EP^行更新，其中，上层区域不包含PML层，EP冲的下 标1表示上层区域，具体更新公式为：
[0047] 式中，η表示时间步
.，.il和kl分别表示上层区域中Ρ向和ζ向 的空间位置，为真空中介电常数，lu为环路系数；
[0048]步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法，对 该区域内的电磁场分量Ep2进行更新，其中，下层区域不包含PML层，EP冲的下标2表示下层 区域，具体更新公式为：
[0051] 式中，η表示时间步，
.i2和k2分别表示下层区域中Ρ向和ζ 向的空间位置；
[0052] 步骤(4.3)、对边界区域的场量传递进行更新：
[0053]边界区域场量的传递分为如下几种情况：
[0054] a、当上层粗网格中Epi的计算需用到下层细网格中的Ηφ2时，由于上下粗细网格比 为奇数，场量重叠，故直接取相应细网格的场量即可；
[0055] b、当下层细网格中边界网格处Ερ2的计算需用到边界线上的ΗΦ时，直接采用边界线 上细网格的Η Φ2即可；
[0056] 步骤(4.4)、对CFS-PML中的电场分量进行更新：
[0057] 上方吸收边界中的网格空间大小与粗网格相同，右侧吸收边界中网格的空间大小 与相邻左侧计算区域的网格空间大小相同，也分为上下两层，上层与粗网格一致，下层与细 网格一致，
[0058] 吸收边界中的场分量的差分公式为：
[0063] 其中，η表示时间步，i和k分别表示计算区域中Ρ向和ζ向的空间位置，kz、σ ζ和<^为 吸收边界参数。
[0064] 步骤5具体为：
[0065]步骤（5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4) H)TD方法，对该区域内的电磁场分量Ezl进行更新，其中，上层区域不包含PML层，Ez冲的下 标1表示上层区域，具体更新公式为：
[0067] 式中，η表示时间步，
，，il和kl分别表示上层区域中Ρ向和ζ向 的空间位置，为真空中介电常数，lu为环路系数；
[0068]步骤(5.2)、根据步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区 域内的电磁场分量Ez2进行更新，其中，下层区域不包含PML层，Ez2中下标2表示下层区域，具 体更新公式为：
[0070] 式中，η表示时间步 i2和k2分别表示下层区域中Ρ向和ζ向 ，. 的空间位置；
[0071] 步骤(5.3)、交界线上Ez均为细网格的Ez2,上层粗网格中E zl的更新不包括交界线上 的点；
[0072] 步骤(5.4)、对CFS-PML中的电场分量<进行更新:吸收边界中的场分量五]的差分 公式为：
[0079] 其中，η表示时间步，i和k分别表示计算区域中P向和z向的空间位置，aP、kP、〇P、 kp^x和〇Pmax为吸收边界参数，p = i X δ表示p向长度，pQ是PML层与非PML层的分界位置。
[0080] 步骤6具体为：
[0081]步骤(6.1)、首先根据步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD 方法，对该区域内的电磁场分量好^进行更新，其中，上层区域不包含PML层，具体更新公式 为：
[0083] 其中，η表示时间步，il和kl分别表示上层区域中P向和z向的空间位置，ki和k2为环 路系数，c为光速，为磁导率；
[0084]步骤（6.2)、根据步骤1中模型文件中定义的下层区域（不包括PML层）采用传统 FDTD方法，对该区域内的Φ向磁场分量ΗΦ2进行更新，具体更新公式为：
[0086] 其中，η表示时间步，i2和k2分别表示下层区域中P向和z向的空间位置，c为光速， 为磁导率；
[0087] 步骤(6.3)、边界区域的场量传递：
[0088]边界区域场量的传递分为如下2种情况：
[0089] a、当上层粗网格中山^的计算需用到下层细网格中的Ερ2时，由于上下粗细网格比 为奇数，场量重叠，故直接取相应细网格的场量即可，同时粗网格中ΗΦ1的更新不包括交界 线上的点；
[0090] b、当交界线上的Ηφ均取细网格的ΗΦ2,且交界线上ΗΦ2的计算需用到上层细网格中 的，采用线性插值的方法获取，具体过程如下：
_3]式中，~1^_2根据粗嘯中的电场分量&11_和~1^_采用线性插 值的方法得到^_4U2+7&即
[0095] 步骤(6.4)、对CFS-PML中的电场分量ΗΦ进行更新：吸收边界中的场分量ΗΦ的差分 公式为：
[0100] 其中，η表示时间步，i和k分别表示计算区域中P向和z向的空间位置，αη、kn和ση为 吸收边界参数。
[0101] 步骤7具体为：对激励源的更新过程需根据步骤1中模型文件中所设置的源的种 类、激励场形式、类型，对所有的源进行更新。
[0102] 步骤8具体为：
[0103] 判断是否达到总的迭代次数，若达到则进行步骤9;否则时间步加1，并返回步骤4。
[0104] 步骤9具体为：
[0105] a、当激励源为单频正弦信号时：
[0106]根据峰值监测法提取观测点信号垂直方向电场幅值|EzQ|，则当天线辐射功率为Pt 时，场强|EZ|为
[0108]提取观测点处某一时域稳态周期内垂直电场波形过零点时刻Ts，则观测点的传播 时延^为
[0110]式中，d为传播距离，c为光速，to为天线所加电流源中与所提取周期时刻相对应的 参考时刻；
[0111] b、当激励源为Loran-C信号时：
[0112]提取观测点垂直方向电场波形中的第三载波周期负峰值场强，即采样点场强 和峰值场强pil，则当天线辐射功率为pt时，采样点场强|五_〔 |和峰值场强|分别为：
[0115] 式中，B=AX(65X10-6)2Xe-2 ? 5.717976X 10-10Α，
[0116] 提取第三载波周期正向过零点时刻Ts，观测点的传播时延tw为
[0118] 式中，T为周期。
[0119] 本发明的有益效果是，一种"M(2，4)FDTD+FDTD"的低色散低频地波传播时延预测 方法，首先，"M(2，4)FDTD+FDTD"的低色散低频地波传播时延预测方法综合考虑了低频地波 传播时延预测中整个电波传播区域中传播路径的特征:①包含空气层和地表层;②空气层 对低频地波传播的影响主要为大气折射率的空时分布所产生的影响，而大气折射率的空时 变化较缓慢，所产生的影响也较小;③地表层对低频地波传播的影响主要为地形起伏变化、 地面电特性(大地电导率σ和相对介电常数所产生的影响，而地形和地物的变化较复杂， 所产生的影响也较大。因此，对整个电波传播区域进行分层处理，即分为上（空气层）下（地 表层）两层；其次，"M(2,4)FDTD+FDTD"的低色散低频地波传播时延预测方法综合利用了Μ (2，4 )FDTD方法和传统H)TD方法的特点：①传统H)TD方法可用于处理任何复杂结构/媒质的 电波传播问题，但是网格剖分密度越小，数值色散误差越大;而网格剖分密度过大，则计算 机消耗（内存占用、计算时长)越大；②M(2,4)FDTD通过修正系数的引进，可大大降低传统 H)TD方法的数值色散误差，但不适用于复杂路径的电波传播问题。因此，将两种方法结合， 发挥其各自优势，对上层空气层采用粗网格的M( 2，4 )H)TD进行传播时延预测，在保证预测 精度的同时减少计算机内存占用、提高计算速度;对下层地表层采用细网格的传统H)TD方 法进行传播时延预测，保证预测结果的可靠性和准确性;最后，巧妙利用亚网格技术，实现 上下区域间场量的传递;利用CFS-PML技术，实现电磁波的有效吸收，保证预测结果的精度。
【附图说明】
[0120]图1是"本发明一种"M(2,4)FDTD+FDTD"的低色散低频地波传播时延预测方法原理 示意图；
[0121] 图2是本发明一种"M(2，4)FDTD+FDTD"的低色散低频地波传播时延预测方法中Μ (2，4) FDTD方法的网格示意图；
[0122] 图3是本发明一种"M(2，4)FDTD+H)TD"的低色散低频地波传播时延预测方法中计 算区域上下层网格排布及场量分布示意图；
[0123] 图4是本发明一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方法的流 程图。
【具体实施方式】
[0124] 下面结合附图和【具体实施方式】对本发明进行详细说明。
[0125] 本发明一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方法，理论基础 及原理是：
[0126] 针对低频地波传播时延预测问题，一方面希望达到较高的时延预测精度；另一方 面要在有限的资源和一定的时间范围内实现长距离、复杂传播路径问题的求解。为了兼顾 这两方面需求，解决网格剖分密度与计算机资源占用及耗时之间的矛盾，如图1所示，将Μ (2,4)H)TD方法与传统H)TD方法进行结合，以实现低频地波传播时延的高精度预测，同时， 提高计算速度、减少计算机内存占用，其中，M(2,4)FDTD方法网格示意图如图2所示，包含 C1、C2和C3三个环路，并引入环路系数k#Pk2;计算区域上下层网格排布及场量分布示意图 如图3所示，设定上下粗细网格比为奇数，便于场量的传递;流程图如图4所示，具体步骤如 下：
[0127] 步骤1、设置并生成模型文件、输入模型文件；
[0128] 步骤2、参数设置及初始化；
[0129] 步骤3、添加场源；
[0130] 步骤4、更新整个计算区域的P向电场分量EP;
[0131 ]步骤5、更新整个计算区域的z向电场分量Ez;
[0132] 步骤6、更新整个计算区域的Φ向磁场分量ΗΦ;
[0133] 步骤7、更新激励源；
[0134]步骤8、判断结束条件，循环；
[0135] 步骤9、观测点观测量计算并输出。
[0136] 其中，每一步骤具体做法如下：
[0137] 步骤1具体为：
[0138] 步骤（1.1)、设模型文件的上层区域大小NPl ΧΝζ1，下层区域大小Νρ2 X Nz2，CFS_PML 层数为NPML，其中NP为P方向网格数具为z方向网格数，标号1和2分别表示上层和下层；
[0139] 步骤（1.2)、设置空间、时间步长:上层粗网格空间步长为δ?，其中，Δρ1=Δζ1 = δ 1，下层细网格空间步长为处，其中，Αρ2= Δζ2 = δ2,其中Δρ为ρ方向网格大小，ΔΖ为ζ方 向网格大小，标号1和2分别表示上层和下层，时间步长为Δ t，粗网格M( 2，4)H)TD的时间步 长设置为与细网格FDTD同样的时间步长；
[0140]步骤(1.3)、设置迭代次数为NTCalc;
[0141]步骤（1.4)、设置传播路径电参数:P方向起始网格位置为pstart、P方向结束网格位 置为Pe?、z方向起始网格位置为zstart、z方向结束网格位置为zstart、大地电导率为σ、相对介 电常数为
[0142] 步骤(1.5)、设置吸收边界:CFS-PML层数为Npml，相关参数为Knmax，a nmax，onmax，其中η = Ρ,ζ；
[0143] 步骤（1.6 )、设置场源:源的个数Nsouixe、位置[Pstart，PEen ]和[ZStart，ZEen]，源的种 类:有两种激励源供选择：单频正弦源、1^^811_(：源，激励场形式：有三种激励形式:E z、Ep、 源的类型：软源或硬源，幅度A，频率f，单频正弦源中的常数t〇,高斯脉宽!^，时延/包周 差丁；
[0144] 步骤（1.7)、设置观测点：观测点个数NvPoint，位置[PStart，PEen]和[ZStart，ZEen]，输出 场量类型(E Z、EP或Ηφ)。
[0145] 步骤2具体为：
[0146]步骤(2.1)、将整个区域的电磁场分量(Εζ、ΕΡ和ΗΦ)和电磁场分量系数(CA、CB)、中 间变好#、/./&,、￡W、馬Ρ)、中间变Μ系数（f〇z、flz、f2z，f〇P、flP、f2p，f〇4>Y、 ?^φγ、?2Φγ，）、辅助变莖(Φφζ、Φφρ、Φφγ、Φρζ、Φζρ)、观测莖(传播时延tw、米样点场强￡:〇峰值 场强|$^〇|)都初始化为零；
[0147]步骤(2.2 )、初始化所有网格中的模型参数:将相对介电常数ε4刀始化为1，大地电 导率〇初始化为〇;
[0148] 步骤（2.3)、根据步骤1所设置的模型文件中的路径信息，给相应网格的^和〇赋 值；
[0149] 步骤(2.4)设置CFS-PML吸收边界参数kn、〇n、an，其中，体按下式计算：
[0151] 式中，n = P，Z，no为PML层与非PML截面位置，d是PML吸收边界的厚度，1^_取整数， Knmax取值范围为[1，60]，anmax取值范围为[0，1)，onmax根据σ。^设置，στ^χ/σ。^取值范围为(0， 12]，〇。的=( 111+1)/15031八11，111取值范围为[1，20]，其中111取值为4时边界的吸收效果最好，八11
[0152] 步骤3具体为：
[0153] 添加的场源有两种类型：
[0154] (1)单频正弦源：
[0155] 当辐射源采用单频正弦信号时，其电流源激励。(〇可以表示为：
[0157] 其中 to = 5X10-6s。
[0158] (2)Loran_C 源：
[0159] 当辐射源采用正相位编码的Loran_C信号时，其电流波形前沿is(t)为：
[0161] 其中，τ为包周差，单位为s。
[0162] 步骤4具体为：
[0163] 步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格Μ(2,4) H)TD方法，对该区域内的电磁场分量ΕΡ^行更新，其中，上层区域不包含PML层，ΕΡ冲的下 标1表示上层区域，具体更新公式为：
[0166] 式中，η表示时间步
：，il和kl分别表示上层区域中Ρ向和ζ向 的空间位置，为真空中介电常数，lu为环路系数；
[0167] 步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法，对 该区域内的电磁场分量Ep2进行更新，其中，下层区域不包含PML层，E P冲的下标2表示下层 区域，具体更新公式为：
[0170] 式中，η表示时间步，
、i2和k2分别表示下层区域中Ρ向和ζ 向的空间位置；
[0171] 步骤(4.3)、对边界区域的场量传递进行更新：
[0172] 边界区域场量的传递分为如下几种情况，如图3所示：
[0173] a、当上层粗网格中Epi的计算需用到下层细网格中的Ηφ2时，由于上下粗细网格比 为奇数，场量重叠，故直接取相应细网格的场量即可；
[0174] b、当下层细网格中边界网格处Ερ2的计算需用到边界线上的ΗΦ时，直接采用边界线 上细网格的ΗΦ2即可；
[0175] 步骤(4.4)、对CFS-PML中的电场分量进行更新：
[0176] 上方吸收边界中的网格空间大小与粗网格相同，右侧吸收边界中网格的空间大小 与相邻左侧计算区域的网格空间大小相同，也分为上下两层，上层与粗网格一致，下层与细 网格一致，
[0177] 吸收边界中的场分量g的差分公式为：
[0182] 其中，η表示时间步，i和k分别表示计算区域中P向和z向的空间位置，kz、σ ζ和<^为 吸收边界参数。
[0183] 步骤5具体为：
[0184] 步骤（5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格Μ(2,4) H)TD方法，对该区域内的电磁场分量Εζ1进行更新，其中，上层区域不包含PML层，Εζ冲的下 标1表示上层区域，具体更新公式为：
[0186] 式中，η表示时间步
：，:il和kl分别表示上层区域中Ρ向和ζ向 的空间位置，为真空中介电常数，lu为环路系数；
[0187] 步骤(5.2)、根据步骤1中模型文件中定义的下层区域采用传统FDTD方法，对该区 域内的电磁场分量Ez2进行更新，其中，下层区域不包含PML层，E z2中下标2表示下层区域，具 体更新公式为：
[0189] 式中，η表示时间步，
i2和k2分别表示下层区域中Ρ向和ζ向 的空间位置；
[0190] 步骤(5.3)、如图3所示，交界线上Ez均为细网格的Ez2,上层粗网格中E zl的更新不包 括交界线上的点；
[0191] 步骤(5.4)、对CFS-PML中的电场分量 < 进行更新:吸收边界中的场分量 <的差分 公式为：
[0198] 其中，η表示时间步，i和k分别表示计算区域中P向和z向的空间位置，aP、k P、〇P、 kp^x和〇Pmax为吸收边界参数，p = i X δ表示p向长度，pQ是PML层与非PML层的分界位置。
[0199] 步骤6具体为：
[0200]步骤(6.1)、首先根据步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD 方法，对该区域内的电磁场分量/_^进行更新，其中，上层区域不包含PML层，具体更新公式 为：
[0201]
[0202] 其中，η表示时间步，il和kl分别表示上层区域中P向和z向的空间位置，ki和k2为环 路系数，c为光速，为磁导率；
[0203] 步骤（6.2)、根据步骤1中模型文件中定义的下层区域（不包括PML层）采用传统 FDTD方法，对该区域内的Φ向磁场分量ΗΦ2进行更新，具体更新公式为：
[0205] 其中，η表示时间步，i2和k2分别表示下层区域中Ρ向和ζ向的空间位置，c为光速， 为磁导率；
[0206] 步骤(6.3)、边界区域的场量传递：
[0207] 如图3所示，边界区域场量的传递分为如下2种情况：
[0208] a、当上层粗网格中Ηφι的计算需用到下层细网格中的Ep2时，由于上下粗细网格比 为奇数，场量重叠，故直接取相应细网格的场量即可，同时粗网格中Η Φ1的更新不包括交界 线上的点；
[0209] b、当交界线上的ΗΦ均取细网格的ΗΦ2,且交界线上ΗΦ2的计算需用到上层细网格中 的，采用线性插值方法获取，具体过程如下：
[0212]式中，￡；2[^2+7/2根据粗网格中的电场分量￡? Η+ι/，￡ρ1|1ι+ι/2采用线性插 值的方法得到^：4砂'即
[0214] 步骤(6.4)、对CFS-PML中的电场分量ΗΦ进行更新：吸收边界中的场分量Η Φ的差分 公式为：
7
[0219] 其中，n表示时间步，i和k分别表示计算区域中p向和Z向的空间位置，α η、kn和ση为 吸收边界参数。
[0220] 步骤7具体为：对激励源的更新过程需根据步骤1中模型文件中所设置的源的种 类、激励场形式、类型，对所有的源进行更新。
[0221] 步骤8具体为：
[0222] 判断是否达到总的迭代次数，若达到则进行步骤9;否则时间步加1，并返回步骤4。
[0223] 步骤9具体为：
[0224] a、当激励源为单频正弦信号时：
[0225]根据峰值监测法提取观测点信号垂直方向电场幅值|EzQ|，则当天线辐射功率为Pt 时，场强|EZ|为
[0227]提取观测点处某一时域稳态周期内垂直电场波形过零点时刻Ts，则观测点的传播 时延^为
[0229]式中，d为传播距离，c为光速，to为天线所加电流源中与所提取周期时刻相对应的 参考时刻；
[0230] b、当激励源为Loran-c信号时：
[0231]提取观测点垂直方向电场波形中的第三载波周期负峰值场强，即采样点场强|五^ 和峰值场强|￡^|，则当天线福射功率为Pt时，采样点场强|^:|和峰值场强|纪|分别为：
[0234] 式中，B=AX (65X10-6)2Xe-2 ? 5.717976X 10-10A，
[0235] 提取第三载波周期正向过零点时刻Ts，观测点的传播时延tw为
【主权项】
1. 一种"M(2，4)FDTD+roTD"的低色散低频地波传播时延预测方法，其特征在于，通过对 低色散低频地波计算区域进行分层处理，将粗网格M(2,4)H)TD与细网格传统FDTD方法相结 合，进行低频地波传播时延的预测，从而达到在保证预测精度的同时，提高计算速度、减少 计算机内存占用，具体步骤如下： 步骤1、设置并生成模型文件、输入模型文件； 步骤2、参数设置及初始化； 步骤3、添加场源； 步骤4、更新整个计算区域的P向电场分量EP; 步骤5、更新整个计算区域的z向电场分量Ez; 步骤6、更新整个计算区域的Φ向磁场分量Ηφ; 步骤7、更新激励源； 步骤8、判断结束条件，循环； 步骤9、观测点观测量计算并输出。2. 根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤1具体为： 步骤（1.1)、设模型文件的上层区域大小NPl XNzl，下层区域大小Np2 XNz2，CFS-PML层数 为Npml，其中Np为P方向网格数，Nz为z方向网格数，标号1和2分别表不上层和下层； 步骤（1.2)、设置空间、时间步长:上层粗网格空间步长为δ?，其中，Δρ1=Δζ1 = δ1，Τ 层细网格空间步长为δ2,其中，Δρ2= Δζ2 = δ2,其中Δρ为ρ方向网格大小，ΔΖ为ζ方向网 格大小，标号1和2分别表示上层和下层，时间步长为△ t，粗网格M(2，4)H)TD的时间步长设 置为与细网格FDTD同样的时间步长； 步骤（1.3)、设置迭代次数为NTCalc; 步骤（1.4)、设置传播路径电参数:p方向起始网格位置为pstart、p方向结束网格位置为 Pe?、ζ方向起始网格位置为zstart、ζ方向结束网格位置为zstart、大地电导率为σ、相对介电常 数为 步骤(1.5)、设置吸收边界:CFS-PML层数为Npml，相关参数为Knmax，anmax，o nmax，其中η = ρ， ζ； 步骤（1.6 )、设置场源:源的个数Nsource、位置[PStart，PEen]和[ZStart，ZEen]，源的种类:有两 种激励源供选择:单频正弦源、1^〇瓜11_(：源，激励场形式:有三种激励形式:Ez、E P、％，源的类 型:软源或硬源，幅度A，频率f，单频正弦源中的常数to，高斯脉宽!^，时延/包周差τ; 步骤（1.7)、设置观测点:观测点个数NvPoint，位置[Pstart，PEen]和[ZStart，ZEen]，输出场量 类型(EZ、EP或Ηφ)。3. 根据权利要求1所述的一种"M(2,4)FDTD+roTD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤2具体为： 步骤(2.1)、将整个区域的电磁场分量(Ez、EP和ΗΦ)和电磁场分量系数(CA、CB)、中间变 量）、中间变量系数(1^、；1^、；1^，；1^、；1^、；1^，；^()4^、；^14^、 5<^，）、辅助变量(^、1^、1^、1^、'）、观测量(传播时延^、采样点场强€()峰值场强 都初始化为零； 步骤(2.2)、初始化所有网格中的模型参数:将相对介电常数^初始化为1，大地电导率σ 初始化为〇; 步骤(2.3)、根据步骤1所设置的模型文件中的路径信息，给相应网格的^和〇赋值； 步骤(2.4)设置CFS-PML吸收边界参数kn、ση、αη，其中，k n、ση、αη具体按下式计算：式中，η = Ρ，ζ，％为PML层与非PML截面位置，d是PML吸收边界的厚度，Knmax取整数，Knmax 取值范围为[1，60 ]，anmax取值范围为[〇，1)，onmax根据σ。^设置，στ^χ/σ。^取值范围为(〇，12 ]， 〇opt= (m+l)/15〇Ji Δ ri，m取值范围为[1，20]，其中m取值为4时边界的吸收效果最好，Δ τι取值 范围人为源的波长。 30 10 ，4. 根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤3具体为： 添加的场源有两种类型： (1) 单频正弦源： 当辐射源采用单频正弦信号时，其电流源激励is(t)可以表示为： 、 ' -， 其中 t〇 = 5X10-6s， (2) Loran_C 源： 当辐射源采用正相位编码的Loran_C信号时，其电流波形前沿is(t)为：7 其中，τ为包周差，单位为s。5. 根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤4具体为： 步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)H)TD 方法，对该区域内的电磁场分量EPl进行更新，其中，上层区域不包含PML层，EPl*的下标1表 示上层区域，具体更新公式为： L. -一 〇 式中，η表示时间步，标,il和kl分别表示上层区域中Ρ向和ζ向的空 间位置，ε〇为真空中介电常数，ki为坏路糸数； 步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统H)TD方法，对该区 域内的电磁场分量Ep2进行更新，其中，下层区域不包含PML层，Ep2中的下标2表示下层区域， 具体更新公式为：\式中，η表示时间步，标号i _、i 2和k2分别表示下层区域中P向和z向的空 V " ) 间位置； 步骤(4.3)、对边界区域的场量传递进行更新： 边界区域场量的传递分为如下几种情况： a、 当上层粗网格中计算需用到下层细网格中的ΗΦ2时，由于上下粗细网格比为奇 数，场量重叠，故直接取相应细网格的场量即可； b、 当下层细网格中边界网格处Ερ2的计算需用到边界线上的ΗΦ时，直接采用边界线上细 网格的Η Φ2即可； 步骤(4.4)、对CFS-PML中的电场分量进行更新： 上方吸收边界中的网格空间大小与粗网格相同，右侧吸收边界中网格的空间大小与相 邻左侧计算区域的网格空间大小相同，也分为上下两层，上层与粗网格一致，下层与细网格 一致， 吸收边界中的场分量E；；的差分公式为：其中，η表示时间步，i和k分别表示计算区域中P向和z向的空间位置，kz、〇4Paz为吸收 边界参数。6.根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤5具体为： 步骤(5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)H)TD方 法，对该区域内的电磁场分量Ezl进行更新，其中，上层区域不包含PML层，Ezl中的下标1表示 上层区域，具体更新公式为：式中，η表示时间步，标号，_il和kl分别表示上层区域中Ρ向和ζ向的空 间位置，ε〇为真空中介电常数，ki为环路系数； 步骤(5.2)、根据所述步骤1中模型文件中定义的下层区域采用传统H)TD方法，对该区 域内的电磁场分量Ez2进行更新，其中，下层区域不包含PML层，Ez2中下标2表示下层区域，具 体更新公式为：式中，η表示时间步，标,i2和k2分别表示下层区域中Ρ向和ζ向的空 间位置； 步骤（5.3)、交界线上Ez均为细网格的Ez2,上层粗网格中Ezl的更新不包括交界线上的 占 . 步骤(5.4)、对CFS-PML中的电场分量 < 进行更新:吸收边界中的场分量 <的差分公式 为： V ·--1 J 式中，标号m= (i，k+l/2)，且其中，η表示时间步，i和k分别表示计算区域中P向和z向的空间位置，αρ、kp、〇p、kpmax和 〇Pmax为吸收边界参数，p = i X δ表示p向长度，pQ是PML层与非PML层的分界位置。7.根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤6具体为： 步骤(6.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)H)TD 方法，对该区域内的电磁场分量进行更新，其中，上层区域不包含PML层，具体更新公式 为：其中，η表示时间步，i 1和kl分别表示上层区域中P向和z向的空间位置，1^和1?为环路系 数，c为光速，yr为磁导率； 步骤（6.2)、根据所述步骤1中模型文件中定义的下层区域（不包括PML层）采用传统 FDTD方法，对该区域内的Φ向磁场分量ΗΦ2进行更新，具体更新公式为：其中，η表示时间步，i 2和k2分别表示下层区域中P向和z向的空间位置，c为光速，yr为磁 导率； 步骤(6.3)、边界区域的场量传递： 边界区域场量的传递分为如下2种情况： a、 当上层粗网格中山^的计算需用到下层细网格中的Ep2时，由于上下粗细网格比为奇 数，场量重叠，故直接取相应细网格的场量即可，同时粗网格中取^的更新不包括交界线上 的点； b、 当交界线上的ΗΦ均取细网格的ΗΦ2,且交界线上ΗΦ2的计算需用到上层细网格中的E Pl 时，采用线性插值的方法获取，具体过程如下：L U-Δ,· 」 式中，44-,根据粗网格中的电场分量 方法得即步骤(6.4)、对CFS-PML中的电场分量Ηφ进行更新：吸收边界中的场分量Ηφ的差分公式 为：其中，η表示时间步，i和k分别表示计算区域中P向和z向的空间位置，αη、kn和〇n为吸收 边界参数。8. 根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤7具体为:对激励源的更新过程需根据所述步骤1中模型文件中所 设置的源的种类、激励场形式、类型，对所有的源进行更新。9. 根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤8具体为： 判断是否达到总的迭代次数，若达到则进行步骤9;否则时间步加1，并返回步骤4。10. 根据权利要求1所述的一种"M(2,4)FDTD+H)TD"的低色散低频地波传播时延预测方 法，其特征在于，所述步骤9具体为： a、当激励源为单频正弦信号时： 根据峰值监测法提取观测点信号垂直方向电场幅值I EzQ |，则当天线辐射功率为Pt时， 场强|EZ|为提取观测点处某一时域稳态周期内垂直电场波形过零点时刻Ts，则观测点的传播时延 tw为式中，d为传播距离，C为光速，to为天线所加电流源中与所提取周期时刻相对应的参考 时刻； b、当激励源为Loran-C信号时： 提取观测点垂直方向电场波形中的第三载波周期负峰值场强，即采样点场强|五^|和峰 值场强，则当天线辐射功率为Pt时，采样点场强和峰值场强|￡「/|分别为：式中，B=AX(65X10-6)2Xe-2 ? 5.717976X10-10A， 提取第三载波周期正向过零点时刻Ts，观测点的传播时延U为 式中，T为周期。
【文档编号】G06F19/00GK105868571SQ201610251307
【公开日】2016年8月17日
【申请日】2016年4月21日
【发明人】蒲玉蓉, 周丽丽, 张金生, 席晓莉, 顾妍
【申请人】西安理工大学