高精度预测低频电波传播特性的双向抛物方程方法与流程

本发明属于电波传播技术领域，具体涉及一种高精度预测低频电波传播特性的双向抛物方程方法。

背景技术：

低频电波以其波长较长，信号的传播损耗小，信号的幅度与相位稳定，被广泛应用于授时、导航、通信等领域。为提高低频无线电工程的使用性能和精度，需要深入研究低频电波在各种地域、时域、频域范围内的变化规律及预测技术。目前，预测复杂环境下低频电波传播已有的理论方法有：积分方程方法、抛物方程方法和时域有限差分(finite-differencetime-domain，fdtd)方法。积分方程方法和抛物方程方法都可以考虑地形变换，对复杂路径具有较高精度，但是，这两种方法都忽略了沿路径后向传播的反射场影响，在地形起伏剧烈的传播路径上会引起较大误差。fdtd方法可方便的考虑任何地质和地形路径的电波传播问题，并且精度较高，但其存在过长的计算时间、较多的内存资源占用等问题，不适合大区域的工程化推广。

技术实现要素：

本发明的目的是提供一种高精度预测低频电波传播特性的双向抛物方程方法，既能解决fdtd方法计算量大、耗时长的缺点，又能解决积分方程方法和抛物方程方法在地形起伏剧烈路径中由于忽略后向电波传播影响引起的误差大的问题。

本发明所采用的技术方案是：高精度预测低频电波传播特性的双向抛物方程方法，具体按照以下步骤实施：

步骤1：输入模型文件；

步骤2：利用平地面公式计算初始位置ρ＝ρ0处纵向横切面上的磁场hf(ρ0,z)，并通过其求解双向抛物方程的前向初始场uf(ρ0,z)；

步骤3：结合步骤2计算出来的前向初始场uf(ρ0,z)，基于坐标变换模型，采用分布离散傅里叶变换算法，求解计算区域任意位置电波传播的前向场uf(ρ,z)；

步骤4：结合步骤2计算出来的前向初始场uf(ρ0,z)，基于阶梯近似模型，采用ssft算法，递归的求解由于地形影响对电波传播产生的后向场和由多次反射产生的前向场τ＝1,2,3,…；

步骤5：结合步骤3和步骤4中的传播场结果uf(ρ,z)、τ＝1,2,3,…和求出电波传播的总磁场

本发明的特点还在于：

步骤1具体为：

计算区域大小nρ×nz，其中nρ为ρ方向的网格数，nz为z方向的网格数；空间网格步长分别为δρ和δz，ρ和z分别表示横向和纵向坐标；ρ0为初始场位置；源的参数；地面相对介电常数εr和电导率σ；地形模型参数。

步骤2具体为：

选取二维柱坐标系(ρ,z)，其中ρ和z分别表示为横向和纵向坐标，根据实际发射天线尺寸，通过测量得到垂直电偶极子的电荷间距dl，放置在距离地面高度为d的位置，假设时谐因子为e^-iωt，利用平地面公式计算初始距离ρ＝ρ0处纵向横切面的磁场hf(ρ0,z)：

其中，i为电流大小，k0和kg分别为真空和地面的波数，r0表示从源点到观测点的直线距离，r1表示从源的镜像到观测点的直线距离，p2为中间参量：

f(z)是fresnel积分并定义为

定义沿ρ轴正方向传播的波函数，即前向场uf(ρ0,z)为：

步骤3具体为：

在原坐标系(ρ,z)中，地形函数为t(ρ)，以地形表面上任意一点为新的坐标原点o′建立新的坐标系则

在新坐标系中，波函数已经修正为：

其中，

采用起伏地形ssft算法求解下一步进处的前向场：

其中，n为大气折射率，和为离散混合傅里叶变换对，p＝k0sinθ，θ为电波传播的仰角；

步骤2中得到初始距离ρ＝ρ0处电波传播的前向初始场为uf(ρ0,z)，结合式(4)和(5)，即可得到新坐标系下初始场

将代入式(6)中步进计算，得到新坐标系下任意位置处的前向场

再考虑原坐标系和新坐标系前向场之间的关系，即式(4)，求出原坐标系中前向场uf(ρ,z)。

步骤4具体为：

沿ρ轴正方向传播的各个前向场为：

其中，上标τ表示前向场的序号，τ＝0表示电波传播的主传播场，而τ＝1,2,3…则表示由后向场多次反射所产生的前向场，为前向磁场分量；

当电波沿ρ轴正方向传播时，采用ssft算法即可求得下一位置处的前向场：

将步骤2中求得的前向初始场uf(ρ0,z)赋值给若传播路径不含起伏地形，则计算区域仅存在τ＝0时对应的将代入式(8)依次迭代计算出各个位置的场若存在起伏地形，需要考虑后向场影响；

定义沿ρ轴负方向传播的各个后向波函数，即后向场：

其中，表示后向场的序号，每当前向场遇到一个相对于其传播方向的上升阶梯面，就会产生一个后向磁场分量

对应于式(9)的ssft算法为：

假设前向波沿ρ轴正方向传播到ρ＝ρe处存在相对于传播方向的上升阶梯面，由于该上升阶梯面对电波传播的反射作用，前向波分为两部分：继续向前传播的前向波和反射作用产生的后向波；

在上升阶梯面上，根据磁场所满足切向边界条件可以得到后向初始场为：

修正前向场：

其中，τ和都是具体的某一个标号，具体过程为标记为τ的前向波，每遇到一个相对于ρ正方向的上升沿就会产生一个标记为的后向波，并且更新为

假设电波沿ρ轴负方向传播到ρ＝ρt处存在相对于传播方向的上升阶梯面，后向波分为两部分：继续向后传播的后向波和由后向波反射产生的前向波；

在上升阶梯面上，根据磁场所满足切向边界条件可以得到后向场多次反射产生的前向场为：

修正后向场：

其中，τ和都是具体的某一个标号。具体过程为标记为的后向波，每遇到一个相对于ρ负方向的上升沿就会产生一个标记为τ的前向波，并且τ更新为τ＝τ+1；

转到式(7)和式(9)继续循环对所有前向场和后向场迭代计算下一步进处的场强，直至步进到计算区域边界；其中，前向场的边界为ρmax＝nρδρ，后向场的边界为ρ0。

步骤4中循环时设置门限终止计算，循环计算的判定依据为：

其中，u(ρ,z)表示τ＝1,2,3…或者判决系数ξ设置为0.01，||||为矢量求模。

步骤5具体为：

将步骤3和步骤4中的前向场uf(ρ,z)和多次反射产生的前向场τ＝1,2,3,…分别代入式(2)和式(7)，分别求出所有存在的前向场hf(ρ,z)和

将步骤4计算的代入式(9)中，求出所有存在的后向磁场分量

电波传播的总磁场为所有前向磁场分量和后向磁场分量的叠加：

本发明的有益效果是：

①本发明高精度预测低频电波传播特性的双向抛物方程方法，是基于柱坐标系实现的，因此，能够更为方便的处理圆柱对称的电波传播问题，对工程指导比较有意义；

②本发明高精度预测低频电波传播特性的双向抛物方程方法，能够准确的预测复杂环境中电波传播的前向和后向传播影响，克服了积分方程方法和传统的抛物方程方法忽略后向波影响引起的误差。此外，与fdtd方法相比，在相同的计算精度前提下，该方法基于高效的ssft算法步进求解场强，大大的缩短了计算时间，更适合工程推广。

附图说明

图1是本发明方法的流程图；

图2(a)是高度为0.5km时本发明方法与fdtd方法的地面场强比较；

图2(b)是高度为1.5km时本发明方法与fdtd方法的地面场强比较；

图3是在多个高斯山脉地形情况下，本发明方法与fdtd方法的地面场强比较图；

图4(a)是fdtd方法的空间场强伪彩色图比较；

图4(b)是本发明方法的空间场强伪彩色图比较。

具体实施方式

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

本发明提供了高精度预测低频电波传播特性的双向抛物方程方法，改进了传统抛物方程方法在复杂地形传播路径上由于忽略后向传播影响引起的预测误差，具体实现过程是：首先利用平地面公式计算出前向初始场；接着，基于坐标变化模型，采用抛物方程方法预测整个区域在忽略后向传播影响情况下的前向场的场值；然后，基于阶梯近似模型，通过双向抛物方程方法预测由于电波在复杂地形作用下引起的反射场(包括后向场和多次反射产生的前向场)影响。最后，通过抛物方程方法计算的前向场和双向抛物方程方法计算的后向场和多次反射产生的前向场叠加求解总磁场，具体实施步骤流程图见图1。

在采用双向抛物方程方法求圆柱对称模型问题中的前向场和后向场时，首先需要推导柱坐标系下的双向抛物方程及其对应的ssft算法：

假设电磁波的时谐因子是e^-iωt，取二维的柱坐标系(ρ,z)，其中ρ和z分别表示为横向和纵向坐标，并以标量φ表示任意标量场分量，假设φ与方位角无关，则对于水平极化波，而对于垂直极化波，φ满足二维波动方程

其中k0是真空中的波数，n为大气折射率。首先，作如下变量代换：

其中，ψ(ρ,z)仅考虑了电磁波的柱面传播影响，没有涉及相位修正。将式(18)代入(17)中，可以得到：

然后，对上式作远场近似，即假设ρ>>1且k0ρ>>1，式(19)变为：

当折射指数n随ρ的变化很小时，式(20)可分解为

由于式(21)中的两个因式互不相关，故可以得到关于ρ的两个独立的式子

ψf和ψb分别对应于前向波和后向波。接下来，考虑相位修正，引入合适的相位因子，对ψf和ψb分别作如下修正：

其中，φf和φb分别代表前向和后向场分量。将式(18)代入式(24)和(25)中，可以分别得到

其中，φf和φb分别为前向和后向磁场分量。

下面将式(24)和(25)分别代入(22)和(23)，可以得到双向抛物方程：

上面两式分别称为前向和后向抛物方程。采用ssft算法求解上面两式，可得：

式中和为离散混合傅里叶变换对，为大气折射率，δρ为ρ方向的网格大小。实际计算中，φf和φb可能分别为一系列的和其中τ＝0,1,2,3,…和

本发明高精度预测低频电波传播特性的双向抛物方程方法，具体按照以下步骤实施：

步骤1：输入模型文件，具体为：

计算区域大小nρ×nz，其中nρ为ρ方向的网格数，nz为z方向的网格数；空间网格步长分别为δρ和δz，ρ和z分别表示横向和纵向坐标；ρ0为初始场位置；源的参数；地面相对介电常数εr和电导率σ；地形模型参数。

步骤2：利用平地面公式计算初始位置ρ＝ρ0处纵向横切面上的磁场hf(ρ0,z)，并通过其求解双向抛物方程的前向初始场uf(ρ0,z)，具体为：

为了更好地解决圆柱对称模型问题，在此选取二维柱坐标系(ρ,z)，其中ρ和z分别表示为横向和纵向坐标，根据实际发射天线尺寸，通过测量得到垂直电偶极子的电荷间距dl，放置在距离地面高度为d的位置，假设时谐因子为e^-iωt，利用平地面公式计算初始距离ρ＝ρ0处纵向横切面的磁场hf(ρ0,z)：

其中，i为电流大小，k0和kg分别为真空和地面的波数，r0表示从源点到观测点的直线距离，r1表示从源的镜像到观测点的直线距离，p2为中间参量：

f(z)是fresnel积分并定义为

定义沿ρ轴正方向传播的波函数，即前向场uf(ρ0,z)为：

步骤3：结合步骤2计算出来的前向初始场uf(ρ0,z)，基于坐标变换模型，采用分布离散傅里叶变换(split-stepfouriertransform，ssft)算法，求解计算区域任意位置电波传播的前向场uf(ρ,z)，具体为：

在原坐标系(ρ,z)中，地形函数为t(ρ)，以地形表面上任意一点为新的坐标原点o′建立新的坐标系则

在新坐标系中，波函数已经修正为：

其中，

采用起伏地形ssft算法求解下一步进处的前向场：

其中，n为大气折射率，和为离散混合傅里叶变换对，p＝k0sinθ，θ为电波传播的仰角；

步骤2中得到初始距离ρ＝ρ0处电波传播的前向初始场为uf(ρ0,z)，结合式(4)和(5)，即可得到新坐标系下初始场

将代入式(6)中步进计算，得到新坐标系下任意位置处的前向场

再考虑原坐标系和新坐标系前向场之间的关系，即式(4)，求出原坐标系中前向场uf(ρ,z)。

步骤4：结合步骤2计算出来的前向场uf(ρ0,z)，基于阶梯近似模型，采用ssft算法，递归的求解由于地形影响对电波传播产生的后向场和由多次反射产生的前向场τ＝1,2,3,…，具体为：

沿ρ轴正方向传播的各个前向场为：

其中，上标τ表示前向场的序号，τ＝0表示电波传播的主传播场，而τ＝1,2,3…则表示由后向场多次反射所产生的前向场，为前向磁场分量；与式(2)相比，可以看出，前向场的磁场求解具有相同的表达形式。

当电波沿ρ轴正方向传播时，采用ssft算法即可求得下一位置处的前向场：

将步骤2中求得的前向初始场uf(ρ0,z)赋值给若传播路径不含起伏地形，则计算区域仅存在τ＝0时对应的将代入式(8)依次迭代计算出各个位置的场若存在起伏地形，就需要考虑后向场影响；

下面定义沿ρ轴负方向传播的各个后向波函数，即后向场：

其中，表示后向场的序号，每当前向场遇到一个相对于其传播方向的上升阶梯面，就会产生一个后向磁场分量

对应于式(9)的ssft算法为：

比较式(8)和(10)可以得出，前向场和后向场的ssft算法迭代形式基本一致。

假设前向波沿ρ轴正方向传播到ρ＝ρe处存在相对于传播方向的上升阶梯面，由于该上升阶梯面对电波传播的反射作用，前向波分为两部分：继续向前传播的前向波和反射作用产生的后向波。在上升阶梯面上，根据磁场所满足切向边界条件可以得到后向初始场为：

修正前向场：

这里的τ和都是具体的某一个标号。具体过程为标记为τ的前向波，每遇到一个相对于ρ正方向的上升沿就会产生一个标记为的后向波，并且更新为

同理，假设电波沿ρ轴负方向传播到ρ＝ρt处存在相对于传播方向的上升阶梯面，后向波分为两部分：继续向后传播的后向波和由后向波反射产生的前向波。在上升阶梯面上，根据磁场所满足切向边界条件可以得到后向场多次反射产生的前向场为：

修正后向场：

这里的τ和都是具体的某一个标号。具体过程为标记为的后向波，每遇到一个相对于ρ负方向的上升沿就会产生一个标记为τ的前向波，并且τ更新为τ＝τ+1。

然后，采用式(7)和(9)继续对所有前向场和后向场迭代计算下一步进处的场强，接着判断是否存相对于传播方向上的上升阶梯面，具体执行过程见上述假设，循环此过程直至步进到计算区域边界(对于前向场，边界为ρmax＝nρδρ，对于后向场，边界为ρ0)。

由此可以看出对于孤立山峰，电波向前传播遇到山峰反射会产生许多后向波，而后向波不存在反射，即不会产生前向波；而对于多个山峰，当电波向前传播遇到第一个山峰时，一部分电波会由于反射作用而产生后向波，另一部分电波则会继续向前传播，当遇到第二个山峰时，又会出现后向波，并且，后向波会在两个山峰之间往复反射，以此类推，从而产生一系列的反射波(其中包括多次反射产生的前向波和后向波)。由于受空间散射及介质的吸收影响，循环计算的反射波场强逐渐减小，继而不会对全域场产生影响。因此，有必要设置门限来终止计算。为了方便计算，给出递归计算的判定依据为：

其中，u(ρ,z)表示τ＝1,2,3…或者判决系数ξ设置为0.01，||||为矢量求模。

步骤5：结合步骤3和步骤4中的传播场结果uf(ρ,z)、和τ＝1,2,3,…，求出电波传播的总磁场具体为：

将步骤3和步骤4中的前向场uf(ρ,z)和多次反射产生的前向场τ＝1,2,3,…分别代入式(2)和式(7)，分别求出所有存在的前向场hf(ρ,z)和注意：这里没有采用而是用步骤3求出来的hf(ρ,z)替代，这是因为采用阶梯近似模型求出来的前向场分量精度在山后不够。

将步骤4计算的代入式(9)中，求出所有存在的后向磁场分量

电波传播的总磁场为所有前向磁场分量和后向磁场分量的叠加：

实施例1

单个高斯山脉地形地面场强预测

垂直电偶极子天线的辐射功率均采用1kw，信号源频率为100khz。计算区域总大小为ρmax：100km×zmax：102.4km，网格剖分大小均分别为dρ＝200m和dz＝100m，初始距离ρ0＝10km；地面电参数为εr＝13，σ＝3×10^-3s/m(陆地)；在中心位置ρc＝50km处有一孤立高斯山峰，其高度函数为

l为山峰宽度取2km，h为山峰高度且分别取0.5km和1.5km。图2(a)和图2(b)分别为在相同宽度不同高度的单个高斯山脉地形情况下，本发明方法与fdtd方法的地面场强比较。由图2(a)和图2(b)可以看出，两种方法的计算结果吻合一致，并且都可以预测到山前电波反射影响的振荡存在。但是，相比于fdtd方法，双向抛物方程方法计算更快。

实施例2

多个高斯山脉地形地面场强预测

将实施例1中的山峰改为多个高斯山脉地形，第一个山峰的高度为1km，宽度为4km，中心位置为40km；第二个山峰的高度为1.5km，宽度为2km，中心位置为60km，其他参数不变。图3给出了在多个高斯山脉地形情况下，本发明方法与fdtd方法的地面场强比较。由图3可以看出，对于多个高斯山脉地形，双向抛物方程方法也能精确的预测山峰之间波的往复反射影响。并且经统计，fdtd方法的计算时间是双向抛物方程方法的21倍。图4(a)和图4(b)分别为fdtd方法与本发明方法的空间场强伪彩色图比较，从图中可以直观的看到由于山峰的影响，电波在山前和山峰之间的反射影响。