一种高精度预测ASF的窄带离散分布抛物方程方法与流程

本发明属于电波传播技术领域，具体涉及一种高精度预测ASF的窄带离散分布抛物方程方法。

背景技术：

目前获取罗兰-C附加二次时延(Additional Secondary Factor，ASF)的方法可以分为两类：时域方法和频域方法。其中时域方法都是基于FDTD方法展开的，这类方法预测精度高，但是内存消耗大、计算时间过长，因此其不适用于长距离电波传播问题。频域方法主要包括：平地面公式、Fock绕射方法、Wait积分、Millington公式、积分方程(Integral Equation，IE)方法、抛物方程(Parabolic Equation，PE)方法等。频域方法仅适用于单频信号电波传播预测，其对ASF的求解均是依据对100kHz单频信号展开的，而实际的罗兰-C信号为100kHz载波调制的高斯脉冲，周围环境对罗兰-C信号和单频信号的影响不同。因此，有必要对实际长距离罗兰源的ASF分布进行仿真预测。

技术实现要素：

本发明的目的是提供一种高精度预测ASF的窄带离散分布抛物方程方法，解决现有时域方法不适用于长距离电波传播，而现有频域方法预测精度低的问题。

本发明所采用的技术方案是，一种高精度预测ASF的窄带离散分布抛物方程方法，具体按照以下步骤实施：

步骤1：对实测罗兰-C时域电流时域信号i(t)进行采样，采样后得到信号i(n)，其中，n＝0,1,2,…,N_t，将采样后的信号作离散傅里叶变换，分解为多个频率电流分量I(m)，其中，m＝0,1,2,…,N_f；

步骤2：当m＝0时，利用平地面公式计算频率f(m)对应的电流分量I(m)在近区即区域I(ρ≤ρ₀)产生的磁场分布

步骤3：以步骤2计算的近区最大处即近远区边界ρ＝ρ₀处的磁场作为远区即区域II(ρ≥ρ₀)的激励源，用非均匀网格剖分的离散分布抛物方程方法用于求解远区的磁场分布

步骤4：通过步骤2和步骤3可得到频率为f(m)的电流分量I(m)在地面接收点辐射的磁场强度；更新步骤2和步骤3中的m，对其赋值为m＝m+1，当f(m)≤f_max时，重复步骤2和步骤3即可求得罗兰-C带宽中每个频率处电流分量I(m)在地面接收点辐射的磁场强度

步骤5：基于滑动窗的思想，对步骤4中计算的磁场分布，采用离散傅里叶逆变换求解地表任意距离处接收到的时域罗兰-C信号和滑动时间t_window(ρ)，并提取ASF(ρ)分布。

本发明的特点还在于：

步骤1中得到电流分量I(m)具体为：

对采样后的信号i(n)作离散傅里叶变换，公式为：

其中，N_t和Δt分别为离散时域电流信号的样本长度和采样步长，N_f和Δf分别为离散频域电流信号的样本长度和采样步长，记I_max为I(m)的最大值，f_max和f_min为频率上下边界且均对应于0.1I_max处的频率，为了简便，记f(m)＝f_min+mΔf，且有f(m)≤f_max。

步骤2具体为：

设时谐因子为e^jωt，采用二维柱坐标系(ρ,z)，其中ρ和z分别表示为距离和高度坐标，根据实际发射天线尺寸，通过测量得到垂直电偶极子的电荷间距dl，放置在距离地面高度为d的位置，利用平地面公式计算近区即ρ≤ρ₀频率为f(m)的电流分量所产生的辐射地波磁场采用如下公式计算：

其中，ρ₀为近区与远区的分界处，k₀和k₁分别对应频率为f(m)真空和地面的波数，r₁表示从源点到观测点的直线距离，r₀表示从源的镜像点到观测点的直线距离，P为中间参量为：

F(z)是Fresnel积分，其定义为：

步骤3中用非均匀网格剖分的离散分布抛物方程方法用于求解远区的磁场分布具体为：

首先将近远场边界处的磁场作为远区的激励源，即

定义沿ρ轴正方向传播的波函数为：

离散分布抛物方程方法的求解过程为步进迭代算法：

其中，和为离散混合傅里叶变换对，为大气折射率，Δρ为ρ方向的网格大小，对于实际的复杂地形，Δρ为非均匀网格，在均匀平地面Δρ为大网格，在起伏地形处，Δρ取小网格；

由步骤2计算可得到频率为f(m)的电流分量在近远区边界产生的磁场将其代入式(4)和(5)中求得u(m,ρ₀,z)，通过式(6)求出每一步进处的u(m,ρ,z)，从而进一步确定该频率成分电流分量在远区辐射的磁场

步骤4具体为：

步骤2和步骤3可得到频率为f(m)的电流分量I(m)在整个计算区域产生的磁场强度包含区域I和区域II产生的磁场和为了求解地面ASF分布，提取地面接收点z_x＝T(x)辐射的磁场强度其中T(x)是地形高度；

更新步骤2和步骤3中的m，对其赋值为m＝m+1，当f(m)≤f_max时，重复步骤2和步骤3求得罗兰-C带宽中每个频率处电流分量在地表面的磁场强度

步骤5具体为：

定义滑动窗IDFT：

其中，N_w是滑动窗的网格数且N_w＝4T/dt，T为载波周期，4T既表示罗兰-C信号的前四载波周期，也表示滑动窗的持续时间，t_window(ρ)表示每一距离窗口的滑动时间；

下面给出ASF(ρ)提取公式：

ASF(ρ)＝T_AT-T_AT0 (8)

其中

T_AT(ρ)＝t_window(ρ)+3T (9)

T_AT0(ρ)＝t_window0(ρ)+3T (10)

t_window(ρ)和t_window0(ρ)分别表示实际路径和良导体路径的滑动时间，T_AT(ρ)和T_AT0(ρ)分别表示实际路径和良导体路径的信号跟踪点的到达时刻；

结合式(8)、(9)和(10)求出地表的ASF(ρ)分布。

t_window(ρ)具体为：

地波的传播速度小于光速C，所以在应用滑动窗的时候，要随时监测相位跟踪点在窗中的位置N_s有没有发生变化，已知初始时刻相位跟踪点的位置为3N_w/4，开始窗口以光速滑动，每当窗口滑动一次，就需要对此刻相位跟踪点的位置N_s进行校对，当N_s等于3N_w/4，窗口不动且将ρ赋值为ρ+Δρ，否则，向后滑动N_s-3N_w/4个时间网格并且将ρ保持不变，校正后的t_window(ρ)可用公式表示为：

通过步骤4已求出各个频率成分电流分量在地面产生的磁场的将代入式(7)中的H(m)中；结合式(11)，即可求出地表任意距离处接收到的t_window(ρ)。

本发明的有益效果是：

①本发明克服了传统方法难于预测实际罗兰-C信号源在长距离复杂环境下的ASF分布，提出一种高精度预测ASF的窄带离散分布抛物方程方法。与传统SSPE方法相比，该方法直接从罗兰-C信号源入手，对其进行分析，而非100kHz的单频信号，这样避免了对源的近似带来的误差；并且该方法的到达时刻是从时域恢复的罗兰-C信号中提取第三时刻载波周期的过零点，该操作与实际相符，而非100kHz单频信号是从频域对衰减因子相位进行预测。从而，该方法的ASF预测结果比传统SSPE方法与实际吻合的更好。

②本发明一种高精度预测ASF的窄带离散分布抛物方程方法基于传统SSPE方法对远区场分布进行求解，为了提高分布离散傅里叶变换(Split-Step Fourier Transform，SSFT)算法的总体效率，采用非均匀剖分网格思想，在均匀平缓路径，使用大网格，在复杂地形路径，使用小网格。

③本发明一种高精度预测ASF的窄带离散分布抛物方程方法将滑动窗思想和IDFT技术结合，提出滑动窗IDFT，从而能够高效的从频域信号中恢复出有用的截断时域信号，而不产生任何冗余，大大提高了计算的效率。

附图说明

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

图2是本发明方法的结构示意图；

图3是本发明实施例1，在两条200km长平地面路径中，本发明方法与FDTD方法和SSPE方法得到的ASF分布比较；

图4是本发明实施例2，在200km长路径中包含有高斯山的情况下，本发明方法在150km的观测点处得到的归一化罗兰-C信号；

图5是本发明实施例2，在200km长路径中包含有高斯山的情况下，本发明方法与FDTD方法和SSPE方法得到的ASF分布比较；

图6是本发明实施例3，在1000km长路径中包含有高斯山的情况下，本发明方法与FDTD方法和SSPE方法得到的ASF分布比较。

具体实施方式

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

本发明一种高精度预测ASF的窄带离散分布抛物方程方法，所依据的原理为：在传统SSPE的基础上，加入离散傅里叶变换技术对时域信号进行处理，使处理后的频域信号能采用改进的非均匀网格SSPE方法求解，最后通过基于滑动窗的离散傅里叶逆变换还原时域信号并在时域提取ASF，具体实施步骤流程图见图1。

具体按照以下步骤实施：

步骤1：如图2，在区域I中，天线激励为实测罗兰-C时域信号。设罗兰-C电流时域信号为i(t)，起始时刻为0，即t＝0，对实测罗兰-C时域电流时域信号i(t)进行采样，采样后得到信号i(n)，其中，n＝0,1,2,…,N_t，将采样后的信号作离散傅里叶变换(Discrete Fourier Transform，DFT)，分解为多个频率电流分量I(m)，其中，m＝0,1,2,…,N_f；

步骤1中为得到罗兰-C频域离散电流信号I(m)，m＝0,1,2,…,N_f，对采样后的罗兰-C时域电流信号i(n)，n＝0,1,2,…,N_t，作DFT，其公式如下：

其中，N_t和Δt分别为离散时域电流信号的样本长度和采样步长，N_f和Δf分别为离散频域电流信号的样本长度和采样步长，记I_max为I(m)的最大值，f_max和f_min为频率上下边界且均对应于0.1I_max处的频率，为了简便，记f(m)＝f_min+mΔf，且有f(m)≤f_max。

当时域采样步长Δt给定后，N_t根据i(t)的信号周期及持续时间选取。根据奈奎斯特采样定理，频域采样步长则可得到频域网格数为

本发明采用的正相位编码的罗兰-C信号电流波形为：

其中A是与天线峰值有关的常数，τ为包络差，单位s。这里设置A为1，τ为0。

步骤2：当m＝0时，利用平地面公式(Flat Earth Formula，FEF)计算频率f(m)对应的电流分量I(m)在近区即区域I(ρ≤ρ₀)产生的磁场分布具体为：

图2为NB-SSPE方法的传播模型。整个计算区域包含分为两部分，区域I为FEF计算域，在区域I中，采用FEF计算每一频率成分电流分量在近区场(即ρ≤ρ₀)产生的场值(注：区域I中地面必须为光滑平地面)。

设时谐因子为e^jωt，采用二维柱坐标系(ρ,z)，其中ρ和z分别表示为距离和高度坐标，根据实际发射天线尺寸，通过测量得到垂直电偶极子的电荷间距dl，放置在距离地面高度为d的位置，利用平地面公式计算近区即ρ≤ρ₀频率为f(m)的电流分量所产生的辐射地波磁场采用如下公式计算：

其中，ρ₀为近区与远区的分界处，k₀和k₁分别对应频率为f(m)真空和地面的波数，表示从源点到观测点的直线距离，表示从源的镜像点到观测点的直线距离，P为中间参量为：

F(z)是Fresnel积分，其定义为：

步骤3：如图2所示，在区域I和区域II计算的边界ρ＝ρ₀处，将区域I(ρ≤ρ₀)的计算出来的磁场作为区域II(ρ≥ρ₀)的初始边界条件(激励源)；在区域II中，通过处理PE的非均匀网格SSFT算法对场值进行步进计算，并且其对应的上边界采用Cosine窗函数，下边界采用阻抗边界条件。使用非均匀网格剖分的离散分布抛物方程(Spilt-Step PE，SSPE)方法计算频率为f(m)的电流分量在远区产生的磁场分布具体为：

首先将近远场边界处的磁场作为远区的激励源，即

定义沿ρ轴正方向传播的波函数为：

SSPE方法的求解过程为步进迭代算法：

其中，和为离散混合傅里叶变换(Discrete Mixed Fourier Transform，DMFT)对，为大气折射率，Δρ为ρ方向的网格大小，对于实际的复杂地形，Δρ为非均匀网格，在均匀平地面Δρ为大网格，取0.5倍的波长，在起伏地形处，Δρ取小网格，取0.1倍的波长，从而确保计算精度和速度。p＝k₀sinθ，θ为电波传播的仰角，根据Nyquist采样定理，p_max满足z_maxp_max＝πN_z，N_z和Δz分别为z方向的网格个数和大小，且Z_max＝N_zΔz，则Δp·Δz＝π/N_z。

由步骤2计算可得到频率为f(m)的电流分量在近远区边界产生的磁场将其代入式(4)和(5)中求得u(m,ρ₀,z)，通过式(6)求出每一步进处的u(m,ρ,z)，从而进一步确定该频率成分电流分量在远区辐射的磁场

步骤4：通过步骤2和步骤3可得到频率为f(m)的电流分量I(m)在地面接收点辐射的磁场强度包含区域I和区域II产生的磁场和为了求解地面ASF分布，提取地面接收点z_x＝T(x)辐射的磁场强度其中T(x)是地形高度；

更新步骤2和步骤3中的m，对其赋值为m＝m+1，当f(m)≤f_max时，重复步骤2和步骤3即可求得罗兰-C带宽中每个频率处电流分量I(m)在地面接收点辐射的磁场强度

步骤5：

整个罗兰-C信号的持续时间大概为500μs，若采用传统的IDFT方法将频域信号转换为时域罗兰-C信号，会消耗大量计算内存和时间。罗兰-C信号的相位跟踪点为载波周期第三正向过零点，与第三周期后的波形无关，因此为了提高效率，采用IDFT恢复时域信号时，仅需要恢复截断的罗兰-C时域信号，这里取罗兰-C信号的前四载波周期。为了确保罗兰-C信号前四载波周期刚好被恢复出来，引入滑动窗的思想，并将其与IDFT相结合。为了使信号时刻处于窗口的固定位置，在窗口滑动的过程中要保证窗口滑动的速度与信号同步。

基于滑动窗的思想，对步骤4中计算的磁场分布，采用离散傅里叶逆变换(Inverse DFT，IDFT)求解地表任意距离处接收到的时域罗兰-C信号和滑动时间t_window(ρ)，并提取ASF(ρ)分布，具体为：

定义滑动窗IDFT：

其中，N_w是滑动窗的网格数且N_w＝4T/dt，T＝10μs为载波周期，4T既表示罗兰-C信号的前四载波周期，也表示滑动窗的持续时间，t_window(ρ)表示每一距离窗口的滑动时间；

地波的传播速度小于光速C，所以在应用滑动窗的时候，要随时监测相位跟踪点在窗中的位置N_s有没有发生变化，已知初始时刻相位跟踪点的位置为3N_w/4，开始窗口以光速滑动，每当窗口滑动一次，就需要对此刻相位跟踪点的位置N_s进行校对，当N_s等于3N_w/4，窗口不动且将ρ赋值为ρ+Δρ，否则，向后滑动N_s-3N_w/4个时间网格并且将ρ保持不变，校正后的t_window(ρ)可用公式表示为：

通过步骤4已求出各个频率成分电流分量在地面产生的磁场的将代入式(7)中的H(m)中；结合式(11)，即可求出地表任意距离处接收到的t_window(ρ)。

下面给出ASF(ρ)提取公式：

ASF(ρ)＝T_AT-T_AT0 (8)

其中

T_AT(ρ)＝t_window(ρ)+3T (9)

T_AT0(ρ)＝t_window0(ρ)+3T (10)

t_window(ρ)和t_window0(ρ)分别表示实际路径和良导体路径的滑动时间，T_AT(ρ)和T_AT0(ρ)分别表示实际路径和良导体路径的信号跟踪点的到达时刻；

结合式(8)、(9)和(10)求出地表的ASF(ρ)分布。

实施例1

长距离平地面路径中的地表ASF预测

垂直电偶极子由罗兰-C形式电流信号i_s(t)激励，计算区域总大小为ρ_max：200km×z_max：102.4km，大地电参数取两组不同情况：ε_r＝8，σ＝6×10^-4S/m(森林)和ε_r＝¹³，σ＝3×10^-3S/m(陆地)。FEF计算区域为前20km，SSPE计算区域为后180km，网格剖分大小均分别为dρ为500m、dz为100m。时域网格数N_t＝23333，步长Δt＝1.668×10^-8s；频域积分范围75kHz～250kHz，即f_min和f_max分别为75kHz和125kHz，网格数N_f＝51，步长Δf＝1kHz。滑动窗口网格数N_w＝2400。图3为两条200km长平地面路径中，本发明的方法与FDTD方法和SSPE方法得到的ASF分布比较。由图3可见，本发明的方法与FDTD方法的预测结果具有很好的一致性，验证了本发明方法的正确性。此外，与SSPE方法相比，随着距离的增加，本发明的方法和FDTD方法的ASF结果更为接近。

实施例2

长距离复杂地形路径中的电波时域特性和地表ASF预测

垂直电偶极子由罗兰-C形式电流信号i_s(t)激励，计算区域总大小为ρ_max：200km×z_max：102.4km，大地电参数取：ε_r＝8，σ＝6×10^-4S/m(森林)。地形函数为：

其中ρ表示传播路径上的距离，单位km。FEF计算区域为前20km，SSPE计算区域为后180km，平地面的网格剖分大小均分别为dρ为500m，dz为100m，而起伏地形处dρ缩小为100m。其他参数设置不变。图4给出了距离为150km处，归一化的时域磁场波形，从图中可以看出本发明方法与FDTD方法得到的罗兰-C磁场信号曲线几乎完全重合。图5为200km长的路径中包含有高斯山的情况下，三种方法的ASF预测结果比较。通过比较发现，尤其是在远区，本发明的方法与FDTD方法的预测结果更为吻合。

实施例3

超长距离复杂地形路径中的地表ASF预测

将实施例2中的计算区域总大小扩大为ρ_max：1000km×z_max：409.6km，FEF计算区域为前20km，SSPE计算区域为后980km。其他参数设置不变。垂直电偶极子由罗兰-C形式电流信号i_s(t)激励，大地电参数取：ε_r＝8，σ＝6×10^-4S/m(森林)。地形函数为：

其中ρ表示传播路径上的距离，单位km。FEF计算区域为前20km，SSPE计算区域为180km，平地面的网格剖分大小均分别为dρ为500m，dz为100m，而起伏地形处dρ缩小为100m。时域和频域网格剖分规则不变。图6是在1000km长路径包含高斯山的情况下，本发明方法与SSPE方法得到的ASF分布比较。由于FDTD方法的自身缺点，不适用于长距离电波传播预测问题。因此采用本发明的方法作为参考结果，分析SSPE方法的误差是否随距离变化。由图6可见，随着距离的增加，SSPE方法的误差从开始变大慢慢转为稳定。该分析和结果对地波传播特性工程测量与分析具有一定的指导意义。

本发明克服了传统方法难于预测实际罗兰-C信号源在长距离复杂环境下的ASF分布，提出一种高精度预测ASF的窄带离散分布抛物方程方法。与传统SSPE方法相比，该方法直接从罗兰-C信号源入手，对其进行分析，而非100kHz的单频信号，这样避免了对源的近似带来的误差；并且该方法的到达时刻是从时域恢复的罗兰-C信号中提取第三时刻载波周期的过零点，该操作与实际相符，而非100kHz单频信号是从频域对衰减因子相位进行预测。从而，该方法的ASF预测结果比传统SSPE方法与实际吻合的更好。

本发明基于传统SSPE方法对远区场分布进行求解，为了提高SSFT算法的总体效率，采用非均匀剖分网格思想，在均匀平缓路径，使用大网格，在复杂地形路径，使用小网格。

本发明将滑动窗思想和IDFT技术结合，提出滑动窗IDFT，从而能够高效的从频域信号中恢复出有用的截断时域信号，而不产生任何冗余，大大提高了计算的效率。