1.一种分区混合路径低频地波传播特性预测方法,其特征在于,包括如下步骤:步骤1.输入模型文件,输入的模型文件的内容包括计算区域的网格参数、激励源参数、电波传播路径的电参数;
步骤2.利用平地面公式计算由激励源在电波传播路径x=x0处的磁场Hy(x0,y,z),通过磁场Hy(x0,y,z)求解三维抛物方程的初始场分布u(x0,y,z);
步骤3.采用离散混合傅里叶变换和快速傅里叶变换相结合的分步傅里叶变换算法步进求解三维抛物方程,得到整个计算区域内的场分布结果;
步骤4.利用整个计算区域内的场分布结果,求出低频地波在分区混合路径上传播的磁场强度Hy以及产生的二次时延tw;
所述步骤3具体为:
步骤3.1.将初始场分布u(x0,y,z)写成离散化的表达形式u(x0,lΔy,mΔz);
其中,l=1,2,...,Ny,m=0,1,...,Nz,lΔy、mΔz分别表示沿y、z方向的每一个网格点;
在三维直角坐标系(x,y,z)中,Nx为横向坐标x方向的网格数、Ny为纵向坐标y方向的网格数,Nz为高度坐标z方向的网格数;Δx、Δy和Δz为空间网格步长;
步骤3.2.沿y方向求初始场分布u(x0,lΔy,mΔz)的快速傅里叶变换FFT,计算得到场分布其中,Δky表示在y方向的离散波数,i表示虚数单位;
步骤3.3.沿z方向求场分布 的离散混合傅里叶变换DMFT,计算区域在高度z方向被划分为Nz个网格,结合Nyquist准则,得到场分布 的表达式为:其中,∑'表示在求和计算时给m=0和m=N两项均乘以1/2,j=0,1,…,Nz,Δkz表示在z方向的离散波数;
不考虑折射指数项,对 乘以z方向的绕射指数项 得到:其中,kx表示x方向的波数;
步骤3.4.对场分布U(x0+Δx,lΔky,jΔkz)沿z方向求逆离散傅里叶变换IDMFT,得到场分布步骤3.5.对求得的 沿y方向进行快速傅里叶逆变换IFFT,得到场分布u(x0+Δx,lΔy,mΔz);
步骤3.6.考虑折射指数项,用 乘以u(x0+Δx,lΔy,mΔz),得到下一步进面上的场分布;
步骤3.7.重复步骤3.1至步骤3.6,直至求得整个计算区域内的场分布结果。
2.根据权利要求1所述的分区混合路径低频地波传播特性预测方法,其特征在于,所述步骤1具体为:在三维直角坐标系(x,y,z)中,定义计算区域大小为Nx×Ny×Nz;
引入替代变量ξ代表每个方向的最大距离,得到:
ξ=NξΔξ (1)
其中,ξ=x,y,z;
激励源为垂直电偶极子,预设垂直磁偶极子的电流I和电荷间距dl;
电波传播路径即平地面、海陆混合路径的电参数包括介电常数εr和电导率σ。
3.根据权利要求2所述的分区混合路径低频地波传播特性预测方法,其特征在于,所述步骤2具体为:选择解析解法中的平地面公式作为三维抛物方程3DPE的初始场,预设的初始距离为x=x0,x=x0处的地波磁场分量Hy(x0,y,z)表示为:其中,垂直电偶极子位于z轴上离地h的高度处,k0为真空中的传播常数,kg为地面波数,r1表示从计算区域中放置垂直电偶极子的源点到观测初始场位置x0处的观测点的直线距离,r2表示从源点的镜像点到观测点的直线距离,P2为中间参量;有:对于初始距离x=x0,根据公式(2)得到低频地波的初始场分布u(x0,y,z)为:
4.根据权利要求3所述的分区混合路径低频地波传播特性预测方法,其特征在于,所述步骤3中,三维抛物方程3DPE的计算区域由上、下边界以及初始边界组成;
上边界为吸收边界,吸收边界设置为加在y方向和z方向的窗函数W(y)和W(z):其中,ymax为计算区域沿y方向的最大距离,zmax为计算区域沿z方向的最大高度;
从初始距离开始,在每一步进上抛物方程PE计算出来的场分布都要乘以W(y)和W(z);
下边界条件为地表边界,采用阻抗边界条件,即:
其中,u(x,y,z)表示(x,y,z)处的场分布,α为边界面上的阻抗特性参数;
初始边界为初始场分布,即步骤2中通过公式(2)和公式(6)得到的3DPE的初始场分布u(x0,y,z)。
5.根据权利要求4所述的分区混合路径低频地波传播特性预测方法,其特征在于,所述步骤3.4具体为:将公式(9)的阻抗边界条件转为差分方程:
其中,u(x,y,mΔz)表示(x,y,mΔz)处的场值分布,u(x,y,(m+1)Δz)表示u(x,y,mΔz)下一个网格点的场值,u(x,y,(m‑1)Δz)表示u(x,y,mΔz)上一个网格点的场值;
公式(13)的特征方程为:
2
λ+2λαΔz‑1=0 (14)
其中,λ表示特征方程的特征值,解得:
为确保最终求出的解的稳定性,在λ1和λ2中选取模小于1的一项作为λ;设:其中,C1和C2为中间变量计算系数,
则场分布U(x0+Δx,lΔky,jΔkz)沿z方向的逆离散傅里叶变换定义为:其中, 表示u(x,y,mΔz)的离散傅里叶变换形式;
将阻抗边界条件的差分形式作为辅助函数g(x,y,mΔz):定义场u(x,y,z)沿z方向的离散正弦变换DST及其逆变换IDST分别为:其中,DST[u]表示u(x,y,z)沿z方向的离散正弦变换得到的场分布;
对公式(17)进行DST计算,得到 即:
其中, 表示离散混合傅里叶变换,g(x,y,mΔz)由 的IDST解出;
根据辅助函数g(x,y,mΔz),由公式(17)求出场分布u(x,y,mΔz):其中,up(x,y,mΔz)为公式(20)的特解,B1、B2为系数;
up(x,y,mΔz)通过后向递归算法求解,由于up(x,y,mΔz)满足公式(16),则关于up(x,y,mΔz)的差分方程分解为:公式(21)中有λ‑1/λ=‑2αΔz;若令:其中,up(x,y,(m+1)Δz)表示在z方向m+1网格点处的特解值;则公式(21)变为:s(x,y,mΔz)‑λ·s(x,y,(m‑1)Δz)=2Δz·g(x,y,mΔz) (23)根据公式(22)和公式(23),设s(0)=0、up(Nz)=0,s(0)为s(x,y,0Δz),s(x,y,0Δz)表示在z方向初始网格点处的特解值,up(Nz)为up(x,y,NzΔz),up(x,y,NzΔz)表示在z方向Nz网格点处的特解值,通过后向的递归算法求出up(x,y,mΔz);
B1和B2的计算步骤为:
对公式(20)两端分别以 和 形式求和计算,得到:其中,C1和C2的初始值C1(x)和C2(x)由公式(15)计算得到,步进值C1(x+Δx)和C2(x+Δx)为:利用公式(20)至公式(28)对场分布U(x0+Δx,lΔky,jΔkz)沿z方向求逆离散傅里叶变换,得到场分布
6.根据权利要求5所述的分区混合路径低频地波传播特性预测方法,其特征在于,所述步骤4具体为:将步骤3中计算得到的每一步进处的场分布u(x,y,z)代入公式(29)中求出计算区域内的所有磁场Hy(x,y,z):对计算区域内每个位置的磁场Hy(x,y,z)提取实部和虚部,并计算相位其中,Re(·)表示提取实部操作,Im(·)表示提取虚部操作; 表示电磁波在传播路径中的累积相位延迟;
当电波传播路径设置为分区混合路径时,通过公式(30)得到低频地波在分区混合路径上传播的相位 当电波传播路径设置为良导体路径时,通过公式(30)得到低频地波在良导体路径上传播的相位计算低频地波在分区混合路径上传播产生的二次时延tw的公式为:
7.一种用于实现如权利要求1所述的分区混合路径低频地波传播特性预测方法的分区混合路径低频地波传播特性预测系统,其特征在于,所述分区混合路径低频地波传播特性预测系统包括:模型参数初始化模块,用于输入模型文件,输入的模型文件的内容包括计算区域的网格参数、激励源参数、电波传播路径的电参数;
初始场分布求解模块,用于利用平地面公式计算由激励源在电波传播路径x=x0处的磁场Hy(x0,y,z),通过磁场Hy(x0,y,z)求解三维抛物方程的初始场分布u(x0,y,z);
抛物方程求解模块,用于采用离散混合傅里叶变换和快速傅里叶变换相结合的分步傅里叶变换算法步进求解三维抛物方程,得到整个计算区域内的场分布结果;
以及传播特性预测模块,用于利用整个计算区域内的场分布结果,求出低频地波在分区混合路径上传播的磁场强度Hy以及产生的二次时延tw。
8.一种计算机设备,包括存储器和一个或多个处理器,所述存储器中存储有可执行代码,其特征在于,所述处理器执行所述可执行代码时,实现如权利要求1至6任一项所述的分区混合路径低频地波传播特性预测方法的步骤。
9.一种计算机可读存储介质,其上存储有程序,其特征在于,程序被处理器执行时,实现如权利要求1至6任一项所述的分区混合路径低频地波传播特性预测方法的步骤。