CN105868571B - 一种“m(2,4)fdtd+fdtd”的低色散低频地波传播时延预测方法 - Google Patents
一种“m(2,4)fdtd+fdtd”的低色散低频地波传播时延预测方法 Download PDFInfo
- Publication number
- CN105868571B CN105868571B CN201610251307.8A CN201610251307A CN105868571B CN 105868571 B CN105868571 B CN 105868571B CN 201610251307 A CN201610251307 A CN 201610251307A CN 105868571 B CN105868571 B CN 105868571B
- Authority
- CN
- China
- Prior art keywords
- fdtd
- source
- formula
- pml
- grid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16Z—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
- G16Z99/00—Subject matter not provided for in other main groups of this subclass
Landscapes
- Geophysics And Detection Of Objects (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种“M(2,4)FDTD+FDTD”的低色散低频地波传播时延预测方法,该方法通过对整个计算区域进行分层处理,将粗网格M(2,4)FDTD算法(用于上层区域)与细网格传统FDTD算法(用于下层区域)相结合,利用修正因子以及亚网格技术,进行低频地波传播时延的高精度快速预测,本发明解决了传统全区域FDTD算法在减小数值色散误差与降低计算机消耗之间的矛盾,提高了低频地波传播时延预测精度的同时,减少了计算机内存资源占用,提高了计算速度。
Description
技术领域
本发明属于电波传播理论和电磁场数值计算技术领域,具体涉及一种“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方法,通过引入两个修正系数k1和k2,可降低特定频率所有方向的数值色散。但是,该方法无法直接用于非均匀媒质。
目前,尚未见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向电场分量Ez;
步骤6、更新整个计算区域的φ向磁场分量Hφ;
步骤7、更新激励源;
步骤8、判断结束条件,循环;
步骤9、观测点观测量计算并输出。
本发明的特点还在于,
步骤1具体为:
步骤(1.1)、设模型文件的上层区域大小Nρ1×Nz1,下层区域大小Nρ2×Nz2,CFS-PML层数为NPML,其中Nρ为ρ方向网格数,Nz为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)、设置迭代次数为NTCalc;
步骤(1.4)、设置传播路径电参数:ρ方向起始网格位置为ρStart、ρ方向结束网格位置为ρEen、z方向起始网格位置为zStart、z方向结束网格位置为zStart、大地电导率为σ、相对介电常数为εr;
步骤(1.5)、设置吸收边界:CFS-PML层数为NPML,相关参数为κηmax,αηmax,σηmax,其中η=ρ,z;
步骤(1.6)、设置场源:源的个数NSource、位置[ρStart,ρEen]和[zStart,zEen],源的种类:有两种激励源供选择:单频正弦源、Loran_C源,激励场形式:有三种激励形式:Ez、Eρ、源的类型:软源或硬源,幅度A,频率f,单频正弦源中的常数t0,高斯脉宽Tp,时延/包周差τ;
步骤(1.7)、设置观测点:观测点个数NVPoint,位置[ρStart,ρEen]和[zStart,zEen],输出场量类型(Ez、Eρ或Hφ)。
步骤2具体为:
步骤(2.1)、将整个区域的电磁场分量(Ez、Eρ和Hφ)和电磁场分量系数(CA、CB)、中间变量中间变量系数(f0z、f1z、f2z,f0ρ、f1ρ、f2ρ,f0φγ、f1φγ、f2φγ,)、辅助变量(ψφz、ψφρ、ψφγ、ψρz、ψzρ)、观测量(传播时延tw、采样点场强峰值场强)都初始化为零;
步骤(2.2)、初始化所有网格中的模型参数:将相对介电常数εr初始化为1,大地电导率σ初始化为0;
步骤(2.3)、根据步骤1所设置的模型文件中的路径信息,给相应网格的εr和σ赋值;
步骤(2.4)设置CFS-PML吸收边界参数kη、ση、αη,其中,kη、ση、αη具体按下式计算:
式中,η=ρ,z,η0为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)单频正弦源:
当辐射源采用单频正弦信号时,其电流源激励is(t)可以表示为:
其中t0=5×10-6s。
(2)Loran-C源:
当辐射源采用正相位编码的Loran_C信号时,其电流波形前沿is(t)为:
其中,τ为包周差,单位为s。
步骤4具体为:
步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量Eρ1进行更新,其中,上层区域不包含PML层,Eρ1中的下标1表示上层区域,具体更新公式为:
式中,n表示时间步,标号i1和k1分别表示上层区域中ρ向和z向的空间位置,ε0为真空中介电常数,k1为环路系数;
步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法,对该区域内的电磁场分量Eρ2进行更新,其中,下层区域不包含PML层,Eρ2中的下标2表示下层区域,具体更新公式为:
式中,n表示时间步,标号i2和k2分别表示下层区域中ρ向和z向的空间位置;
步骤(4.3)、对边界区域的场量传递进行更新:
边界区域场量的传递分为如下几种情况:
a、当上层粗网格中Eρ1的计算需用到下层细网格中的Hφ2时,由于上下粗细网格比为奇数,场量重叠,故直接取相应细网格的场量即可;
b、当下层细网格中边界网格处Eρ2的计算需用到边界线上的Hφ时,直接采用边界线上细网格的Hφ2即可;
步骤(4.4)、对CFS-PML中的电场分量进行更新:
上方吸收边界中的网格空间大小与粗网格相同,右侧吸收边界中网格的空间大小与相邻左侧计算区域的网格空间大小相同,也分为上下两层,上层与粗网格一致,下层与细网格一致,
吸收边界中的场分量的差分公式为:
其中,
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,kz、σz和αz为吸收边界参数。
步骤5具体为:
步骤(5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量Ez1进行更新,其中,上层区域不包含PML层,Ez1中的下标1表示上层区域,具体更新公式为:
式中,n表示时间步,标号i1和k1分别表示上层区域中ρ向和z向的空间位置,ε0为真空中介电常数,k1为环路系数;
步骤(5.2)、根据步骤1中模型文件中定义的下层区域采用传统FDTD方法,对该区域内的电磁场分量Ez2进行更新,其中,下层区域不包含PML层,Ez2中下标2表示下层区域,具体更新公式为:
式中,n表示时间步,标号i2和k2分别表示下层区域中ρ向和z向的空间位置;
步骤(5.3)、交界线上Ez均为细网格的Ez2,上层粗网格中Ez1的更新不包括交界线上的点;
步骤(5.4)、对CFS-PML中的电场分量进行更新:吸收边界中的场分量的差分公式为:
式中,标号m=(i,k+1/2),且
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,αρ、kρ、σρ、kρmax和σρmax为吸收边界参数,ρ=i×δ表示ρ向长度,ρ0是PML层与非PML层的分界位置。
步骤6具体为:
步骤(6.1)、首先根据步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量进行更新,其中,上层区域不包含PML层,具体更新公式为:
其中,n表示时间步,i1和k1分别表示上层区域中ρ向和z向的空间位置,k1和k2为环路系数,c为光速,μr为磁导率;
步骤(6.2)、根据步骤1中模型文件中定义的下层区域(不包括PML层)采用传统FDTD方法,对该区域内的φ向磁场分量Hφ2进行更新,具体更新公式为:
其中,n表示时间步,i2和k2分别表示下层区域中ρ向和z向的空间位置,c为光速,μr为磁导率;
步骤(6.3)、边界区域的场量传递:
边界区域场量的传递分为如下2种情况:
a、当上层粗网格中Hφ1的计算需用到下层细网格中的Eρ 2时,由于上下粗细网格比为奇数,场量重叠,故直接取相应细网格的场量即可,同时粗网格中Hφ1的更新不包括交界线上的点;
b、当交界线上的Hφ均取细网格的Hφ2,且交界线上Hφ2的计算需用到上层细网格中的Eρ1时,采用线性插值的方法获取,具体过程如下:
的差分公式为:
式中,根据粗网格中的电场分量和采用线性插值的方法得到即
步骤(6.4)、对CFS-PML中的电场分量Hφ进行更新:吸收边界中的场分量Hφ的差分公式为:
式中,
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,αη、kη和ση为吸收边界参数。
步骤7具体为:对激励源的更新过程需根据步骤1中模型文件中所设置的源的种类、激励场形式、类型,对所有的源进行更新。
步骤8具体为:
判断是否达到总的迭代次数,若达到则进行步骤9;否则时间步加1,并返回步骤4。
步骤9具体为:
a、当激励源为单频正弦信号时:
根据峰值监测法提取观测点信号垂直方向电场幅值|Ez0|,则当天线辐射功率为Pt时,场强|Ez|为
提取观测点处某一时域稳态周期内垂直电场波形过零点时刻Ts,则观测点的传播时延tw为
式中,d为传播距离,c为光速,t0为天线所加电流源中与所提取周期时刻相对应的参考时刻;
b、当激励源为Loran-C信号时:
提取观测点垂直方向电场波形中的第三载波周期负峰值场强,即采样点场强和峰值场强则当天线辐射功率为Pt时,采样点场强和峰值场强分别为:
式中,B=A×(65×10-6)2×e-2≈5.717976×10-10A,
提取第三载波周期正向过零点时刻Ts,观测点的传播时延tw为
式中,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三个环路,并引入环路系数k1和k2;计算区域上下层网格排布及场量分布示意图如图3所示,设定上下粗细网格比为奇数,便于场量的传递;流程图如图4所示,具体步骤如下:
步骤1、设置并生成模型文件、输入模型文件;
步骤2、参数设置及初始化;
步骤3、添加场源;
步骤4、更新整个计算区域的ρ向电场分量Eρ;
步骤5、更新整个计算区域的z向电场分量Ez;
步骤6、更新整个计算区域的φ向磁场分量Hφ;
步骤7、更新激励源;
步骤8、判断结束条件,循环;
步骤9、观测点观测量计算并输出。
其中,每一步骤具体做法如下:
步骤1具体为:
步骤(1.1)、设模型文件的上层区域大小Nρ1×Nz1,下层区域大小Nρ2×Nz2,CFS-PML层数为NPML,其中Nρ为ρ方向网格数,Nz为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)、设置迭代次数为NTCalc;
步骤(1.4)、设置传播路径电参数:ρ方向起始网格位置为ρStart、ρ方向结束网格位置为ρEen、z方向起始网格位置为zStart、z方向结束网格位置为zStart、大地电导率为σ、相对介电常数为εr;
步骤(1.5)、设置吸收边界:CFS-PML层数为NPML,相关参数为κηmax,αηmax,σηmax,其中η=ρ,z;
步骤(1.6)、设置场源:源的个数NSource、位置[ρStart,ρEen]和[zStart,zEen],源的种类:有两种激励源供选择:单频正弦源、Loran_C源,激励场形式:有三种激励形式:Ez、Eρ、源的类型:软源或硬源,幅度A,频率f,单频正弦源中的常数t0,高斯脉宽Tp,时延/包周差τ;
步骤(1.7)、设置观测点:观测点个数NVPoint,位置[ρStart,ρEen]和[zStart,zEen],输出场量类型(Ez、Eρ或Hφ)。
步骤2具体为:
步骤(2.1)、将整个区域的电磁场分量(Ez、Eρ和Hφ)和电磁场分量系数(CA、CB)、中间变量中间变量系数(f0z、f1z、f2z,f0ρ、f1ρ、f2ρ,f0φγ、f1φγ、f2φγ,)、辅助变量(ψφz、ψφρ、ψφγ、ψρz、ψzρ)、观测量(传播时延tw、采样点场强峰值场强)都初始化为零;
步骤(2.2)、初始化所有网格中的模型参数:将相对介电常数εr初始化为1,大地电导率σ初始化为0;
步骤(2.3)、根据步骤1所设置的模型文件中的路径信息,给相应网格的εr和σ赋值;
步骤(2.4)设置CFS-PML吸收边界参数kη、ση、αη,其中,kη、ση、αη具体按下式计算:
式中,η=ρ,z,η0为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)单频正弦源:
当辐射源采用单频正弦信号时,其电流源激励is(t)可以表示为:
其中t0=5×10-6s。
(2)Loran-C源:
当辐射源采用正相位编码的Loran_C信号时,其电流波形前沿is(t)为:
其中,τ为包周差,单位为s。
步骤4具体为:
步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量Eρ1进行更新,其中,上层区域不包含PML层,Eρ1中的下标1表示上层区域,具体更新公式为:
式中,n表示时间步,标号i1和k1分别表示上层区域中ρ向和z向的空间位置,ε0为真空中介电常数,k1为环路系数;
步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法,对该区域内的电磁场分量Eρ2进行更新,其中,下层区域不包含PML层,Eρ2中的下标2表示下层区域,具体更新公式为:
式中,n表示时间步,标号i2和k2分别表示下层区域中ρ向和z向的空间位置;
步骤(4.3)、对边界区域的场量传递进行更新:
边界区域场量的传递分为如下几种情况,如图3所示:
a、当上层粗网格中Eρ1的计算需用到下层细网格中的Hφ2时,由于上下粗细网格比为奇数,场量重叠,故直接取相应细网格的场量即可;
b、当下层细网格中边界网格处Eρ2的计算需用到边界线上的Hφ时,直接采用边界线上细网格的Hφ2即可;
步骤(4.4)、对CFS-PML中的电场分量进行更新:
上方吸收边界中的网格空间大小与粗网格相同,右侧吸收边界中网格的空间大小与相邻左侧计算区域的网格空间大小相同,也分为上下两层,上层与粗网格一致,下层与细网格一致,
吸收边界中的场分量的差分公式为:
其中,
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,kz、σz和αz为吸收边界参数。
步骤5具体为:
步骤(5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量Ez1进行更新,其中,上层区域不包含PML层,Ez1中的下标1表示上层区域,具体更新公式为:
式中,n表示时间步,标号i1和k1分别表示上层区域中ρ向和z向的空间位置,ε0为真空中介电常数,k1为环路系数;
步骤(5.2)、根据步骤1中模型文件中定义的下层区域采用传统FDTD方法,对该区域内的电磁场分量Ez2进行更新,其中,下层区域不包含PML层,Ez2中下标2表示下层区域,具体更新公式为:
式中,n表示时间步,标号i2和k2分别表示下层区域中ρ向和z向的空间位置;
步骤(5.3)、如图3所示,交界线上Ez均为细网格的Ez2,上层粗网格中Ez1的更新不包括交界线上的点;
步骤(5.4)、对CFS-PML中的电场分量进行更新:吸收边界中的场分量的差分公式为:
式中,标号m=(i,k+1/2),且
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,αρ、kρ、σρ、kρmax和σρmax为吸收边界参数,ρ=i×δ表示ρ向长度,ρ0是PML层与非PML层的分界位置。
步骤6具体为:
步骤(6.1)、首先根据步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量进行更新,其中,上层区域不包含PML层,具体更新公式为:
其中,n表示时间步,i1和k1分别表示上层区域中ρ向和z向的空间位置,k1和k2为环路系数,c为光速,μr为磁导率;
步骤(6.2)、根据步骤1中模型文件中定义的下层区域(不包括PML层)采用传统FDTD方法,对该区域内的φ向磁场分量Hφ2进行更新,具体更新公式为:
其中,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时,采用线性插值方法获取,具体过程如下:
的差分公式为:
式中,根据粗网格中的电场分量和采用线性插值的方法得到即
步骤(6.4)、对CFS-PML中的电场分量Hφ进行更新:吸收边界中的场分量Hφ的差分公式为:
式中,
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,αη、kη和ση为吸收边界参数。
步骤7具体为:对激励源的更新过程需根据步骤1中模型文件中所设置的源的种类、激励场形式、类型,对所有的源进行更新。
步骤8具体为:
判断是否达到总的迭代次数,若达到则进行步骤9;否则时间步加1,并返回步骤4。
步骤9具体为:
a、当激励源为单频正弦信号时:
根据峰值监测法提取观测点信号垂直方向电场幅值|Ez0|,则当天线辐射功率为Pt时,场强|Ez|为
提取观测点处某一时域稳态周期内垂直电场波形过零点时刻Ts,则观测点的传播时延tw为
式中,d为传播距离,c为光速,t0为天线所加电流源中与所提取周期时刻相对应的参考时刻;
b、当激励源为Loran-C信号时:
提取观测点垂直方向电场波形中的第三载波周期负峰值场强,即采样点场强和峰值场强则当天线辐射功率为Pt时,采样点场强和峰值场强分别为:
式中,B=A×(65×10-6)2×e-2≈5.717976×10-10A,
提取第三载波周期正向过零点时刻Ts,观测点的传播时延tw为
式中,T为周期。
Claims (6)
1.一种M(2,4)FDTD+FDTD的低色散低频地波传播时延预测方法,其特征在于,通过对低色散低频地波计算区域进行分层处理,将粗网格M(2,4)FDTD与细网格传统FDTD方法相结合,进行低频地波传播时延的预测,从而达到在保证预测精度的同时,提高计算速度、减少计算机内存占用,具体步骤如下:
步骤1、设置并生成模型文件、输入模型文件;
步骤2、参数设置及初始化;
步骤3、添加场源;
步骤4、更新整个计算区域的ρ向电场分量Eρ;
步骤5、更新整个计算区域的z向电场分量Ez;
步骤6、更新整个计算区域的φ向磁场分量Hφ;
步骤7、更新激励源;
步骤8、判断结束条件,循环;
步骤9、观测点观测量计算并输出;
步骤1具体为:
步骤(1.1)、设模型文件的上层区域大小Nρ1×Nz1,下层区域大小Nρ2×Nz2,CFS-PML层数为NPML,其中Nρ为ρ方向网格数,Nz为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)、设置迭代次数为NTCalc;
步骤(1.4)、设置传播路径电参数:ρ方向起始网格位置为ρStart、ρ方向结束网格位置为ρEen、z方向起始网格位置为zStart、z方向结束网格位置为zEen、大地电导率为σ、相对介电常数为εr;
步骤(1.5)、设置吸收边界:CFS-PML层数为NPML,相关参数为κηmax,αηmax,σηmax,其中η=ρ,z;
步骤(1.6)、设置场源:源的个数NSource、位置[ρStart,ρEen]和[zStart,zEen],源的种类:有两种激励源供选择:单频正弦源、Loran_C源,激励场形式:有三种激励形式:Ez、Eρ、Hφ,源的类型:软源或硬源,幅度A,频率f,单频正弦源中的常数t0,高斯脉宽Tp,时延/包周差τ;
步骤(1.7)、设置观测点:观测点个数NVPoint,位置[ρStart,ρEen]和[zStart,zEen],输出场量类型Ez、Eρ或Hφ;
步骤4具体为:
步骤(4.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量Eρ1进行更新,其中,上层区域不包含PML层,Eρ1中的下标1表示上层区域,具体更新公式为:
式中,n表示时间步,标号i1和k1分别表示上层区域中ρ向和z向的空间位置,ε0为真空中介电常数,k1为环路系数;
步骤(4.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法,对该区域内的电磁场分量Eρ2进行更新,其中,下层区域不包含PML层,Eρ2中的下标2表示下层区域,具体更新公式为:
式中,n表示时间步,标号i2和k2分别表示下层区域中ρ向和z向的空间位置;
步骤(4.3)、对边界区域的场量传递进行更新:
边界区域场量的传递分为如下几种情况:
a、当上层粗网格中Eρ1的计算需用到下层细网格中的Hφ2时,由于上下粗细网格比为奇数,场量重叠,故直接取相应细网格的场量即可;
b、当下层细网格中边界网格处Eρ2的计算需用到边界线上的Hφ时,直接采用边界线上细网格的Hφ2即可;
步骤(4.4)、对CFS-PML中的电场分量进行更新:
上方吸收边界中的网格空间大小与粗网格相同,右侧吸收边界中网格的空间大小与相邻左侧计算区域的网格空间大小相同,也分为上下两层,上层与粗网格一致,下层与细网格一致,
吸收边界中的场分量的差分公式为:
其中,
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,kz、σz和αz为吸收边界参数;
所述步骤5具体为:
步骤(5.1)、先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量Ez1进行更新,其中,上层区域不包含PML层,Ez1中的下标1表示上层区域,具体更新公式为:
式中,n表示时间步,标号i1和k1分别表示上层区域中ρ向和z向的空间位置,ε0为真空中介电常数,k1为环路系数;
步骤(5.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法,对该区域内的电磁场分量Ez2进行更新,其中,下层区域不包含PML层,Ez2中下标2表示下层区域,具体更新公式为:
式中,n表示时间步,标号i2和k2分别表示下层区域中ρ向和z向的空间位置;
步骤(5.3)、交界线上Ez均为细网格的Ez2,上层粗网格中Ez1的更新不包括交界线上的点;
步骤(5.4)、对CFS-PML中的电场分量进行更新:吸收边界中的场分量的差分公式为:
式中,标号m=(i,k+1/2),且
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,αρ、kρ、σρ、kρmax和σρmax为吸收边界参数,ρ=i×δ表示ρ向长度,ρ0是PML层与非PML层的分界位置;
所述步骤6具体为:
步骤(6.1)、首先根据所述步骤1中模型文件中定义的上层区域采用粗网格M(2,4)FDTD方法,对该区域内的电磁场分量Hφ1进行更新,其中,上层区域不包含PML层,具体更新公式为:
其中,n表示时间步,i1和k1分别表示上层区域中ρ向和z向的空间位置,k1和k2为环路系数,c为光速,μr为磁导率;
步骤(6.2)、根据所述步骤1中模型文件中定义的下层区域采用传统FDTD方法,对该区域内的φ向磁场分量Hφ2进行更新,不包括PML层,具体更新公式为:
其中,n表示时间步,i2和k2分别表示下层区域中ρ向和z向的空间位置,c为光速,μr为磁导率;
步骤(6.3)、边界区域的场量传递:
边界区域场量的传递分为如下2种情况:
a、当上层粗网格中Hφ1的计算需用到下层细网格中的Eρ2时,由于上下粗细网格比为奇数,场量重叠,故直接取相应细网格的场量即可,同时粗网格中Hφ1的更新不包括交界线上的点;
b、当交界线上的Hφ均取细网格的Hφ2,且交界线上Hφ2的计算需用到上层细网格中的Eρ1时,采用线性插值的方法获取,具体过程如下:
的差分公式为:
式中,根据粗网格中的电场分量和采用线性插值的方法得到即
步骤(6.4)、对CFS-PML中的电场分量Hφ进行更新:吸收边界中的场分量Hφ的差分公式为:
式中,
其中,n表示时间步,i和k分别表示计算区域中ρ向和z向的空间位置,αη、kη和ση为吸收边界参数。
2.根据权利要求1所述的一种M(2,4)FDTD+FDTD的低色散低频地波传播时延预测方法,其特征在于,所述步骤2具体为:
步骤(2.1)、将整个区域的电磁场分量Ez、Eρ和Hφ和电磁场分量系数CA、CB、中间变量中间变量系数f0z、f1z、f2z,f0ρ、f1ρ、f2ρ,f0φγ、f1φγ、f2φγ、辅助变量ψφz、ψφρ、ψφγ、ψρz、ψzρ、观测量:传播时延tw、采样点场强峰值场强都初始化为零;
步骤(2.2)、初始化所有网格中的模型参数:将相对介电常数εr初始化为1,大地电导率σ初始化为0;
步骤(2.3)、根据步骤1所设置的模型文件中的路径信息,给相应网格的εr和σ赋值;
步骤(2.4)设置CFS-PML吸收边界参数kη、ση、αη,其中,kη、ση、αη具体按下式计算:
式中,η=ρ,z,η0为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所述的一种M(2,4)FDTD+FDTD的低色散低频地波传播时延预测方法,其特征在于,所述步骤3具体为:
添加的场源有两种类型:
(1)单频正弦源:
当辐射源采用单频正弦信号时,其电流源激励is(t)可以表示为:
其中t0=5×10-6s,
(2)Loran-C源:
当辐射源采用正相位编码的Loran_C信号时,其电流波形前沿is(t)为:
其中,τ为包周差,单位为s。
4.根据权利要求1所述的一种M(2,4)FDTD+FDTD的低色散低频地波传播时延预测方法,其特征在于,所述步骤7具体为:对激励源的更新过程需根据所述步骤1中模型文件中所设置的源的种类、激励场形式、类型,对所有的源进行更新。
5.根据权利要求1所述的一种M(2,4)FDTD+FDTD的低色散低频地波传播时延预测方法,其特征在于,所述步骤8具体为:
判断是否达到总的迭代次数,若达到则进行步骤9;否则时间步加1,并返回步骤4。
6.根据权利要求1所述的一种M(2,4)FDTD+FDTD的低色散低频地波传播时延预测方法,其特征在于,所述步骤9具体为:
a、当激励源为单频正弦信号时:
根据峰值监测法提取观测点信号垂直方向电场幅值|Ez0|,则当天线辐射功率为Pt时,场强|Ez|为
提取观测点处某一时域稳态周期内垂直电场波形过零点时刻Ts,则观测点的传播时延tw为
式中,d为传播距离,c为光速,t0为天线所加电流源中与所提取周期时刻相对应的参考时刻;
b、当激励源为Loran-C信号时:
提取观测点垂直方向电场波形中的第三载波周期负峰值场强,即采样点场强和峰值场强则当天线辐射功率为Pt时,采样点场强和峰值场强分别为:
式中,B=A×(65×10-6)2×e-2≈5.717976×10-10A,
提取第三载波周期正向过零点时刻Ts,观测点的传播时延tw为
式中,T为周期。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610251307.8A CN105868571B (zh) | 2016-04-21 | 2016-04-21 | 一种“m(2,4)fdtd+fdtd”的低色散低频地波传播时延预测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610251307.8A CN105868571B (zh) | 2016-04-21 | 2016-04-21 | 一种“m(2,4)fdtd+fdtd”的低色散低频地波传播时延预测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105868571A CN105868571A (zh) | 2016-08-17 |
CN105868571B true CN105868571B (zh) | 2018-09-14 |
Family
ID=56633696
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610251307.8A Active CN105868571B (zh) | 2016-04-21 | 2016-04-21 | 一种“m(2,4)fdtd+fdtd”的低色散低频地波传播时延预测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105868571B (zh) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106874549B (zh) * | 2017-01-10 | 2020-02-18 | 西安理工大学 | 一种高精度预测asf的窄带离散分布抛物方程方法 |
US10554299B2 (en) * | 2017-05-09 | 2020-02-04 | Huawei Technologies Co., Ltd. | Method and apparatus for characterizing a dispersion of an optical medium |
CN110837688B (zh) * | 2019-09-30 | 2021-11-16 | 西安电子科技大学 | 等离子体鞘套3d-fdtd建模中总场/散射场平面波源产生方法 |
CN111783339B (zh) * | 2020-06-30 | 2024-04-16 | 西安理工大学 | 电磁波在随机色散介质中传播的pce-fdtd方法 |
CN113627074B (zh) * | 2021-07-13 | 2024-04-19 | 西安理工大学 | 一种基于迁移学习的地波传播时延预测方法 |
CN113468759A (zh) * | 2021-07-21 | 2021-10-01 | 安徽大学 | 一种电磁问题确定方法、系统及存储介质 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104794289A (zh) * | 2015-04-23 | 2015-07-22 | 西安理工大学 | 一种扩展直角坐标系下完全匹配吸收边界的实现方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7194497B2 (en) * | 2001-11-19 | 2007-03-20 | Em Photonics, Inc. | Hardware-based accelerator for time-domain scientific computing |
-
2016
- 2016-04-21 CN CN201610251307.8A patent/CN105868571B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104794289A (zh) * | 2015-04-23 | 2015-07-22 | 西安理工大学 | 一种扩展直角坐标系下完全匹配吸收边界的实现方法 |
Non-Patent Citations (1)
Title |
---|
CFS-PML边界条件在PSTD算法中的实现与性能分析;姜永金等;《微波学报》;20041231;第20卷(第4期);第36-39页 * |
Also Published As
Publication number | Publication date |
---|---|
CN105868571A (zh) | 2016-08-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105868571B (zh) | 一种“m(2,4)fdtd+fdtd”的低色散低频地波传播时延预测方法 | |
CN103258246B (zh) | 一种基于元胞历史数据的负荷密度指标的求取方法 | |
CN102770009B (zh) | 一种吸波超材料 | |
CN101621825B (zh) | 一种充分利用实测数据的无线信号场强预测及网络覆盖地图绘制的方法 | |
CN104809343B (zh) | 一种等离子体中使用电流密度卷积完全匹配层的实现方法 | |
CN108494112B (zh) | 一种用于无线电能传输系统的超材料等效电路的分析方法 | |
CN101673869A (zh) | 陶瓷小球和铜线构成的左手材料的制备方法 | |
CN107341284B (zh) | 高精度预测低频电波传播特性的双向抛物方程方法 | |
CN105354635A (zh) | 一种配电网不平衡谐波负荷指数预测方法 | |
CN104573240B (zh) | 周期性非均匀介质波导特征模分析的七点频域有限差分方法 | |
CN106874549A (zh) | 一种高精度预测asf的窄带离散分布抛物方程方法 | |
CN104809286B (zh) | 一种等离子体中扩展坐标的完全匹配吸收边界的实现方法 | |
CN105356500A (zh) | 一种多光伏区域电网电压稳定性裕度预测方法 | |
CN104731996B (zh) | 一种快速提取电大尺寸金属腔体目标瞬态散射信号的仿真方法 | |
CN102572859B (zh) | 基于自适应邻域csa算法的无线传感器网络节点布局方法 | |
CN106777472A (zh) | 基于拉盖尔多项式的减少分裂误差的完全匹配层实现方法 | |
Wang et al. | An efficient integral equation/modified surface integration method for analysis of antenna-radome structures in receiving mode | |
Chai et al. | Indoor radio propagation models and wireless network planning | |
CN104778286B (zh) | 掠海飞行器电磁散射特性快速仿真方法 | |
CN104852154B (zh) | 一种确定网状反射面天线无源互调量的方法 | |
CN102332055B (zh) | 一种极低频电磁波的仿真计算方法 | |
CN107016184A (zh) | 一种二维高精度迭代的非磁化等离子体中的实现方法 | |
Li et al. | Research on TD-LTE wireless communication network propagation model optimization based on visual simulation and GIS | |
Li et al. | An efficient algorithm based on the equivalence principle and FDTD for wave propagation prediction in tunnels | |
CN101937026B (zh) | 高精度预测地波传播衰减因子的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |