CN101937026B - 高精度预测地波传播衰减因子的方法 - Google Patents

Abstract

本发明公开了一种高精度预测地波传播衰减因子的积分方法，该方法按照下列步骤实施：利用积分方程方法计算电波传播路径的衰减因子；在地形变化剧烈及地物突变的局部区域选取数值计算区域，在该区域的电波传播入射面上取积分方法计算的电场分布作为激励源，用电磁场数值计算方法重新计算该区域的传播衰减因子；用数值计算的衰减因子代替采用积分方法在该区域的计算结果，得到整个传播区域的衰减因子。本发明的方法能够在计算资源少量增加的条件下，有效减小积分方程方法在地形和地物突变区域衰减因子预测的误差。

Description

高精度预测地波传播衰减因子的方法

技术领域

本发明属于电波传播理论计算领域，具体涉及一种高精度预测地波传播衰减因子的积分方法。

背景技术

现有地波传播衰减计算采用的主要计算模型和方法有：均匀光滑球面地模型（如：Fock的地波绕射计算方法）、分段均匀光滑球面地模型（如：Millington的经验公式方法、Wait积分方法、波模转换法）、不均匀不光滑球面地模型（如：积分方程方法、抛物线方法）。

早期研究主要采用均匀和分段均匀光滑球面地模型，但预测误差较大。近些年，将积分方程方法用于解决复杂路径地波传播问题具有更好的预测精度，得到了广泛重视和应用。但是用于复杂路径地波传播衰减预测的积分方程方法本身存在如下问题：①在积分方程推导过程中假定地面是平滑的，赫兹矢量只有单一的垂直于地表的分量，当沿传播方向（纵向）地形起伏坡度较大时，预测结果误差较大；②实用的积分方程方法是二维积分方程的一维近似，认为在第一夫累涅尔带内地形横向均匀，从而导致当地形沿横向区域在第一夫累涅尔带内起伏变化较大时，计算结果误差会增大；③积分方程方法是沿传播路径积分到观测点为止，没有考虑积分路径后向的反射场影响，在传播路径后向地形起伏剧烈（陡峭的山脉）或地物突变（如：陆海交界的电导率突变等）时，会引起较大误差。

因此采用积分方程方法在预测地波传播的衰减和延迟时，在地形和地物突变区域的预测误差较大。

电磁场的数值计算可用于对复杂区域电磁场分布的较精确计算，但其占用计算资源较大，在一般计算设备上无法完成地波长距离传播时衰减因子的预测。

发明内容

本发明的目的是提供一种高精度预测地波传播衰减因子的积分方法，解决了现有技术存在的利用积分方程方法在地形和地物突变区域预测衰减因子的误差较大，而利用电磁场的数值计算所占用计算资源过大的问题。

本发明所采用的技术方案是，高精度预测地波传播衰减因子的积分方法，包括以下操作步骤：

步骤1，利用积分方程方法计算入射馈源面上E_z和整个地波传播路径衰减因子的幅度和相位；

步骤2，在地形变化剧烈及地物突变的局部区域选取电磁场数值计算区域，在该区域的电波传播入射馈源面上取积分方程方法计算出的电场E_z分布作为激励源，用电磁场数值计算方法重新计算该区域地波传播衰减因子的幅度和相位；

步骤3，用步骤2得到在地形变化剧烈及地物突变的局部区域的衰减因子幅度和相位代替步骤1得到的积分结果在该区域的衰减因子幅度和相位，得到整个传播路径的地波衰减因子幅度和相位。

其中，在步骤1中利用积分方程方法计算入射馈源面上E_z和整个地波传播路径的衰减因子W_g的公式如下：

$W_g=A\{1-e^{-i\frac{\mathrm{\π}}{4}}\frac{r_0}{x}\sqrt{\frac{k_0}{2\mathrm{\π}}}{\∫}_0^x{W}_g\left(Q\right)e^{i{k}_0(r_1+r_2-r_0)}[{\mathrm{\Δ}}_g(1+\frac{i}{k_0{r}_2})\frac{{\∂r}_2}{\∂n}\left]\sqrt{\frac{x}{L(x-L)}}\mathrm{dL}\right\}---\left(1\right)$

其中，A为常量（观测点在空中时取值为0.5，观测点在地面时取值为1）；i为虚数单位；π为圆周率；k₀为自由空间波数；r₀表示从源点到观测点的直线距离；r₁表示从源点到地面上积分动点Q之间的直线距离；r₂表示从积分动点Q到观测点的距离；W_g(Q)为积分动点Q的地波衰减因子；

表示r₂在表面上的法向导数，n为地面的外法线方向；x表示观测点至源点大圆距离的垂直投影；L是积分动点Q与源点之间的大圆距离的垂直投影；e代表自然对数的底；Δ_g为归一化地表面阻抗，由下式表示：

${\mathrm{\Δ}}_g=i{\left(\frac{k_{0}a}{2}\right)}^{\frac{1}{3}}\frac{\sqrt{{\mathrm{\ϵ}}_r-1+i60\mathrm{\λ\σ}}}{{\mathrm{\ϵ}}_r+i60\mathrm{\λ\σ}}---\left(2\right)$

其中ε_r和σ分别为地面相对介电常数和地面电导率，a为地球半径，λ为波长，k₀为自由空间波数，i为虚数单位；

任意观测点的垂直电场E_z与地波衰减因子W_g满足下列关系式：

E_z＝E₀W_g（3）

其中E₀代表理想导电平面上相同接收距离的电场，

I为电流大小，l电偶极子的电荷间的距离，θ为观测点至场源的大圆角距离，λ为波长，a为地球半径，k₀为自由空间波数，e代表自然对数的底，i为虚数单位，π为圆周率。

其中，在步骤2中的电磁场数值方法采用2维柱坐标系(ρ，φ，z)下时域有限差分方法，该方法中包含ρ方向电场E_ρ、z方向电场E_z以及φ方向磁场H_φ三个分量；麦克斯韦方程组在柱坐标系下FDTD差分表达式如下：

$E_{\mathrm{\ρ}}^{n+1}(i+\frac{1}{2},k)=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{E}_{\mathrm{\ρ}}^n(i+\frac{1}{2},k)-\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}$

$\·\frac{H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k-\frac{1}{2})}{\mathrm{\Δz}}---\left(4\right)$

$E_z^{n+1}(i,k+\frac{1}{2})=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{E}_z^n(i,k+\frac{1}{2})+\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}$

$\·\frac{(i+\frac{1}{2})H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-(i-\frac{1}{2})H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i-\frac{1}{2},k+\frac{1}{2})}{\mathrm{i\Δ\ρ}}---\left(5\right)$

$H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})=H_{\mathrm{\φ}}^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-\frac{\mathrm{\Δt}}{\mathrm{\μ}}$

$\·[\frac{E_{\mathrm{\ρ}}^n(i+\frac{1}{2},k+1)-E_{\mathrm{\ρ}}^n(i+\frac{1}{2},k)}{\mathrm{\Δz}}-\frac{E_z^n(i+1,k+\frac{1}{2})-E_z^n(i,k+\frac{1}{2})}{\mathrm{\Δ\ρ}}]---\left(6\right)$

其中n表示迭代时刻；i和k分别表示ρ方向和z方向的空间网格节点；ε为介电常数；σ为电导率；μ为磁导率；Δt为迭代时间步长；Δρ和Δz分别为ρ方向和z方向的网格剖分大小；

在ρ=0的轴对称边界上，E_z需要通过以下差分形式迭代：

$E_z^{n+1}(0,k+\frac{1}{2})=E_z^n(0,k+\frac{1}{2})+\frac{4\mathrm{\Δt}}{\mathrm{\ϵ\Δ\ρ}}{H}_{\mathrm{\φ}}^{n+\frac{1}{2}}(\frac{1}{2},k+\frac{1}{2})+\frac{4\mathrm{I\Δt}}{{\mathrm{\ϵ\π\Δ\ρ}}^2}---\left(7\right)$

其中，I为电流大小，π为圆周率；根据公式（4）-（7）可得出任意观测点每个迭代时刻的垂直电场值E_Z，再根据垂直电场E_z与地波衰减因子W_g满足的关系式：

E_z＝E₀W_g（8）

即可求出地波衰减因子W_g；其中E₀代表理想导电平面上相同接收距离的电场，I为电流大小，l电偶极子的电荷间的距离，θ为观测点至场源的大圆角距离，λ为波长，a为地球半径，k₀为自由空间波数，e代表自然对数的底，i为虚数单位，π为圆周率。

其中，地波传播衰减因子频率范围为30Hz-300MHz。

本发明的有益效果是，采用电磁场数值计算方法与积分方程方法相结合，在积分方程方法预测误差较大的区域采用电磁场数值计算方法进行修正，这样可在计算资源少量增加的条件下，减小积分方程方法在地形和地物突变区域衰减因子预测的误差，提高积分方程方法的精度。

附图说明

图1是本发明实施例1短距离传播应用图；

图2是本发明实施例1中三种不同算法得到的衰减因子幅度比较图（a为传统FDTD方法测得的衰减因子幅度图，b为本发明方法测得的衰减因子幅度图，c为传统积分方程方法测得的衰减因子幅度图）；

图3是本发明实施例1中三种不同算法得到的衰减因子相位比较图（a为传统FDTD方法测得的衰减因子相位图，b为本发明方法测得的衰减因子相位图，c为传统积分方程方法测得的衰减因子相位图）；

图4是本发明实施例2长距离传播应用图；

图5是本发明实施例2方法得到的衰减因子幅度图；

图6是本发明实施例2方法得到的衰减因子相位图；

图7是本发明实施例3陆海分界区域应用示意图；

图8是本发明实施例3方法得到的衰减因子幅度图；

图9是本发明实施例3方法得到的衰减因子相位图。

具体实施方式

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

本发明提供了一种高精度预测地波传播衰减因子的方法，采用电磁场数值计算的方法改进积分方程方法在复杂传播路径上的地波传播衰减因子的预测精度，具体实现过程是：首先用积分方程方法预测整个传播路径的衰减因子，然后在地形变化剧烈和地物突变区域的局部区域采用电磁场的数值计算方法进行更高精度的预测（电磁场数值计算的源采用积分方程方法的计算结果给出），最后用电磁场数值计算结果替代相同位置的积分结果，得到整个传播路径的衰减因子。该方法在占用计算资源有限的条件下，能够完成了积分方程方法在局部区域误差的修正。

其中的地波传播衰减因子包括地波的强度衰减因子和相位延迟量，相位延迟量包括二次相位因子、二次相位、二次相位延迟、二次相位延时、二次相位时延、附加二次相位因子、附加二次因子、附加二次相位、附加二次相位延迟、附加二次相位时延、附加二次相位延时。

其中所指的地波传播衰减因子频率范围为30Hz-300MHz。

地形变化剧烈区域是指地形沿电波传播方向（纵向）起伏剧烈的区域及沿垂直于电波传播方向（横向）起伏剧烈的区域；地物突变区域是地物从平原、山脉、城市、湖泊、海水、沙漠等不同地物任意二者间的过渡区域。

电磁场数值计算区域选取原则为：包含整个突变区域，并且数值计算区域边界距突变区域距离不小于3个波长。

其中的电磁场数值计算方法是指时域有限差分方法、二维数值计算方法。

下面给出本发明预测地波传播衰减因子精度的方法所涉及到的相关公式：

传统积分方程方法求解地波衰减因子W_g的公式如下：

$W_g=A\{1-e^{-i\frac{\mathrm{\π}}{4}}\frac{r_0}{x}\sqrt{\frac{k_0}{2\mathrm{\π}}}{\∫}_0^x{W}_g\left(Q\right)e^{i{k}_0(r_1+r_2-r_0)}[{\mathrm{\Δ}}_g(1+\frac{i}{k_0{r}_2})\frac{{\∂r}_2}{\∂n}\left]\sqrt{\frac{x}{L(x-L)}}\mathrm{dL}\right\}---\left(1\right)$

其中，A为常量（观测点在空中时取值为0.5，观测点在地面时取值为1）；i为虚数单位；π为圆周率；k₀为自由空间波数；r₀表示从源点到观测点的直线距离；r₁表示从源点到地面上积分动点Q之间的直线距离；r₂表示从积分动点Q到观测点的距离；W_g(Q)为积分动点Q的地波衰减因子；

表示r₂在表面上的法向导数，n为地面的外法线方向；x表示观测点至源点大圆距离的垂直投影；L是积分动点Q与源点之间的大圆距离的垂直投影；e代表自然对数的底；Δ_g为归一化地表面阻抗，由下式表示：

${\mathrm{\Δ}}_g=i{\left(\frac{k_{0}a}{2}\right)}^{\frac{1}{3}}\frac{\sqrt{{\mathrm{\ϵ}}_r-1+i60\mathrm{\λ\σ}}}{{\mathrm{\ϵ}}_r+i60\mathrm{\λ\σ}}---\left(2\right)$

其中ε_r和σ分别为地面相对介电常数和地面电导率，a为地球半径，λ为波长，k₀为自由空间波数，i为虚数单位。

任意观测点的垂直电场E_z与地波衰减因子W_g满足下列关系式：

E_z＝E₀W_g（3）

其中，E₀代表理想导电平面上相同接收距离的电场，

I为电流大小，l电偶极子的电荷间的距离，θ为观测点至场源的大圆角距离，λ为波长，a为地球半径，k₀为自由空间波数，e代表自然对数的底，i为虚数单位，π为圆周率。

有关公式（1）的推导及数值实现参见文献：潘威炎.长波超长波极长波传播[M].成都,电子科技大学出版社.2004。

电磁场数值方法采用2维柱坐标系(ρ，φ，z)下时域有限差分方法（FDTD）。该方法中包含ρ方向电场E_ρ、z方向电场E_z以及φ方向磁场H_φ三个分量。麦克斯韦方程组在柱坐标系下FDTD差分表达式如下：

$E_{\mathrm{\ρ}}^{n+1}(i+\frac{1}{2},k)=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{E}_{\mathrm{\ρ}}^n(i+\frac{1}{2},k)-\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}$

$\·\frac{H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k-\frac{1}{2})}{\mathrm{\Δz}}---\left(4\right)$

$E_z^{n+1}(i,k+\frac{1}{2})=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{E}_z^n(i,k+\frac{1}{2})+\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}$

$\·\frac{(i+\frac{1}{2})H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-(i-\frac{1}{2})H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i-\frac{1}{2},k+\frac{1}{2})}{\mathrm{i\Δ\ρ}}---\left(5\right)$

$H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})=H_{\mathrm{\φ}}^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-\frac{\mathrm{\Δt}}{\mathrm{\μ}}$

$\·[\frac{E_{\mathrm{\ρ}}^n(i+\frac{1}{2},k+1)-E_{\mathrm{\ρ}}^n(i+\frac{1}{2},k)}{\mathrm{\Δz}}-\frac{E_z^n(i+1,k+\frac{1}{2})-E_z^n(i,k+\frac{1}{2})}{\mathrm{\Δ\ρ}}]---\left(6\right)$

其中，n表示迭代时刻；i和k分别表示ρ方向和z方向的空间网格节点；ε为介电常数；σ为电导率；μ为磁导率；Δt为迭代时间步长；Δρ和Δz分别为ρ方向和z方向的网格剖分大小。

在ρ=0的轴对称边界上，E_z需要通过以下差分形式迭代：

$E_z^{n+1}(0,k+\frac{1}{2})=E_z^n(0,k+\frac{1}{2})+\frac{4\mathrm{\Δt}}{\mathrm{\ϵ\Δ\ρ}}{H}_{\mathrm{\φ}}^{n+\frac{1}{2}}(\frac{1}{2},k+\frac{1}{2})+\frac{4\mathrm{I\Δt}}{{\mathrm{\ϵ\π\Δ\ρ}}^2}---\left(7\right)$

其中，I为电流大小，π为圆周率。根据公式（4）-（7）可得出任意观测点每个迭代时刻的垂直电场值E_Z，再根据垂直电场E_z与地波衰减因子W_g满足的关系式：

E_z＝E₀W_g（8）

即可求出地波衰减因子；其中E₀代表理想导电平面上相同接收距离的电场，

I为电流大小，l电偶极子的电荷间的距离，θ为观测点至场源的大圆角距离，λ为波长，a为地球半径，k₀为自由空间波数，e代表自然对数的底，i为虚数单位，π为圆周率。

有关公式（4）-（7）的推导参见文献：葛德彪，闫玉波.电磁波时域有限差分法[M].西安,西安电子科技大学出版社.2004。

实施例1

100kHz地波传播衰减因子短距离传播的预测

图1为100kHz信号从发射台A点到接收机B点传播路径的地形分布，传播距离总长为100km，其中C为电磁场数值计算区域（距离发射点40km，范围20km），包含的地形变化剧烈区域为一较陡的高斯型山脉。地面相对介电常数ε_r为13，电导率σ为0.003S/m，地形函数z(ρ)表达式如下：

$z\left(\mathrm{\ρ}\right)={1500e}^{-9{\left(\frac{\mathrm{\ρ}-50}{2}\right)}^2}---\left(9\right)$

其中ρ表示传播路径上点到发射台A点距离，单位为km；e代表自然对数的底。

对该传播路径分别采用传统的积分方程方法（公式（1））、传统的FDTD方法（公式（4）-（7））和本发明方法计算的地面接收的传播衰减因子的幅度和相位比较分别参见图2和图3。由图2和图3可见，传统的FDTD方法和本发明的改进积分方程方法结果在数值计算区域一致，验证了本发明方法的正确性，传统的积分方程方法与FDTD方法的差异，说明在地形变化剧烈区域传统的积分方程方法存在较大误差，无法在实践中使用。

采用本发明方法预测地波传播衰减因子的具体步骤是：

1）采用积分方程方法（公式（1）和（3））计算图1所示的入射馈源面上E_z和整个传播路径A点到B点的衰减因子的幅度和相位；

2）在山脉前后取电磁场数值计算区域，即为40-60km的区段，在该区域的电波传播入射馈源面上取积分方程方法计算出的电场E_z分布作为电磁场数值计算的激励源，采用公式（4）-（6）重新计算地形变化剧烈区域的电场分布，并采用公式（8）计算提取衰减因子幅度和相位，得到该地形变化剧烈区域地波传播衰减因子的幅度和相位如图2（b）和图3（b）所示；

3）用步骤2得到在地形变化剧烈区域（40-60km的区段）的衰减因子幅度和相位代替步骤1得到的积分结果在地形变化剧烈区域（40-60km的区段）的衰减因子幅度和相位，得到整个传播路径的衰减因子幅度和相位。

实施例2

100kHz地波传播衰减因子长距离传播的预测

图4为100kHz信号从发射台A点到接收机B点传播路径的地形分布，其中C为电磁场数值计算区域（距离发射点300km-320km范围），包含的地形变化剧烈区域同实施例1，地面相对介电常数ε_r为13，电导率σ为0.003S/m，地形函数z(ρ)表达式如下：

$z\left(\mathrm{\ρ}\right)={1500e}^{-9{\left(\frac{\mathrm{\ρ}-210}{2}\right)}^2}---\left(10\right)$

其中ρ表示传播路径上点到发射台A点距离，单位为km；e代表自然对数的底。

采用本发明方法进行传播衰减因子幅度和相位计算的具体步骤是：

1）采用积分方程方法（公式（1）和（3））计算图4所示的入射馈源面上E_z和整个传播路径A点到B点的衰减因子的幅度和相位；

2）在山脉前后取电磁场数值计算区域（距离发射点300km-320km范围），在该区域的电波传播入射馈源面上取积分方程方法计算出的电场E_z分布作为电磁场数值计算的激励源，采用公式（4）-（6）重新计算地形变化剧烈区域的电场分布，并采用公式（8）计算提取衰减因子的幅度和相位，得到该地形变化剧烈区域地波传播衰减因子的幅度和相位如图5和图6所示；

3）用步骤2得到在地形变化剧烈区域（距离发射点300km-320km范围）的衰减因子幅度和相位代替步骤1得到的积分结果在地形变化剧烈区域（距离发射点300km-320km范围）的衰减因子幅度和相位，得到整个传播路径的衰减因子幅度和相位。

实施例3

100kHz地波传播衰减因子陆海分界区域传播的预测

图7为100kHz信号从发射台A点到接收机B点传播路径的地形、地物分布，其中C为电磁场数值计算区域（距离发射点300km-320km范围），包含的地物突变区域为不同地物陆海交界区域（分界点为310km处），陆地相对介电常数ε_r为13，电导率为0.003S/m，海水相对介电常数ε_r为80，电导率σ为5S/m。

采用本发明方法传播衰减因子幅度和相位计算的具体步骤是：

1）采用积分方程方法（公式（1）和（3））计算图7所示的入射馈源面上E_z和整个传播路径A点到B点的衰减因子的幅度和相位；

2）在陆海交界区域前后取电磁场数值计算区域（距离发射点300km-320km范围），在该区域的电波传播入射馈源面上取积分方程方法计算出的电场E_z分布作为电磁场数值计算的激励源，采用公式（4）-（6）重新计算地物突变区域的电场分布，并采用公式（8）计算提取衰减因子的幅度和相位，得到该突变区域地波传播衰减因子的幅度和相位如图8和图9所示；

3）用步骤2得到在突变区域（距离发射点300km-320km范围）的衰减因子幅度和相位代替步骤1得到的积分结果在突变区域（距离发射点300km-320km范围）的衰减因子幅度和相位，得到整个传播路径的衰减因子幅度和相位。

Claims (4)

1.高精度预测地波传播衰减因子的方法，其特征在于：包括以下操作步骤：

步骤1，利用积分方程方法计算入射馈源面上的垂直电场E_z和整个路径上地波传播衰减因子的幅度和相位；

步骤2，在地形变化剧烈及地物突变的局部区域选取电磁场数值计算区域，在电磁场数值计算区域的电波传播入射馈源面上取积分方程方法计算出的垂直电场E_z分布作为激励源，用电磁场数值计算方法重新计算地形变化剧烈及地物突变区域地波传播衰减因子的幅度和相位；

步骤3，用步骤2得到在地形变化剧烈及地物突变的局部区域的地波传播衰减因子幅度和相位代替步骤1得到的积分结果在地形变化剧烈及地物突变的局部区域的地波传播衰减因子幅度和相位，得到整个路径的地波传播衰减因子幅度和相位。

2.根据权利要求1所述高精度预测地波传播衰减因子的方法，其特征在于：在步骤1中利用积分方程方法计算入射馈源面上的垂直电场E_z和整个路径上的地波传播衰减因子W_g的公式如下：

$W_g=A\{1-e^{-i\frac{\mathrm{\π}}{4}}\frac{r_0}{x}\sqrt{\frac{k_0}{2\mathrm{\π}}}{\∫}_0^x{W}_g\left(Q\right)e^{i{k}_0(r_1+r_2-r_0)}[{\mathrm{\Δ}}_g(1+\frac{i}{k_0{r}_2})\frac{{\∂r}_2}{\∂n}\left]\sqrt{\frac{x}{L(x-L)}}\mathrm{dL}\right\}---\left(1\right)$

其中，A为常量，当观测点在空中时取值为0.5，当观测点在地面时取值为1；i为虚数单位；π为圆周率；k₀为自由空间波数；r₀表示从源点到观测点的直线距离；r₁表示从源点到地面上积分动点Q之间的直线距离；r₂表示从积分动点Q到观测点的距离；W_g(Q)为积分动点Q的地波传播衰减因子；

表示r₂在表面上的法向导数，n为地面的外法线方向；x表示观测点至源点大圆距离的垂直投影；L是积分动点Q与源点之间的大圆距离的垂直投影；e代表自然对数的底；Δ_g为归一化地表面阻抗，由下式表示：

${\mathrm{\Δ}}_g=i{\left(\frac{k_{0}a}{2}\right)}^{\frac{1}{3}}\frac{\sqrt{{\mathrm{\ϵ}}_r-1+i60\mathrm{\λ\σ}}}{{\mathrm{\ϵ}}_r+i60\mathrm{\λ\σ}}---\left(2\right)$

其中ε_r和σ分别为地面相对介电常数和地面电导率，a为地球半径，λ为波长，k₀为自由空间波数，i为虚数单位；

任意观测点的垂直电场E_z与地波传播衰减因子W_g满足下列关系式：

E_z＝E₀W_g（3）

其中E₀代表理想导电平面上相同接收距离的电场，

I为电流大小，l为电偶极子的电荷间的距离，θ为观测点至场源的大圆角距离，λ为波长，a为地球半径，k₀为自由空间波数，e代表自然对数的底，i为虚数单位，π为圆周率。

3.根据权利要求1所述高精度预测地波传播衰减因子的方法，其特征在于：步骤2中的电磁场数值方法采用2维柱坐标系(ρ，φ，z)下时域有限差分方法，该方法中包含ρ方向电场E_ρ、z方向电场E_z以及φ方向磁场H_φ三个分量；麦克斯韦方程组在柱坐标系下FDTD差分表达式如下：

$E_{\mathrm{\ρ}}^{n+1}(i+\frac{1}{2},k)=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{E}_{\mathrm{\ρ}}^n(i+\frac{1}{2},k)-\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}$

$\·\frac{H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k-\frac{1}{2})}{\mathrm{\Δz}}---\left(4\right)$

$E_z^{n+1}(i,k+\frac{1}{2})=\frac{2\mathrm{\ϵ}-\mathrm{\σ\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}{E}_z^n(i,k+\frac{1}{2})+\frac{2\mathrm{\Δt}}{2\mathrm{\ϵ}+\mathrm{\σ\Δt}}$

$\·\frac{(i+\frac{1}{2})H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-(i-\frac{1}{2})H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i-\frac{1}{2},k+\frac{1}{2})}{\mathrm{i\Δ\ρ}}---\left(5\right)$

$H_{\mathrm{\φ}}^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})=H_{\mathrm{\φ}}^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2})-\frac{\mathrm{\Δt}}{\mathrm{\μ}}$

$\·[\frac{E_{\mathrm{\ρ}}^n(i+\frac{1}{2},k+1)-E_{\mathrm{\ρ}}^n(i+\frac{1}{2},k)}{\mathrm{\Δz}}-\frac{E_z^n(i+1,k+\frac{1}{2})-E_z^n(i,k+\frac{1}{2})}{\mathrm{\Δ\ρ}}]---\left(6\right)$

其中n表示迭代时刻；i和k分别表示ρ方向和z方向的空间网格节点；ε为介电常数；σ为电导率；μ为磁导率；Δt为迭代时间步长；Δρ和Δz分别为ρ方向和z方向的网格剖分大小；

在ρ=0的轴对称边界上，E_z需要通过以下差分形式迭代：

$E_z^{n+1}(0,k+\frac{1}{2})=E_z^n(0,k+\frac{1}{2})+\frac{4\mathrm{\Δt}}{\mathrm{\ϵ\Δ\ρ}}{H}_{\mathrm{\φ}}^{n+\frac{1}{2}}(\frac{1}{2},k+\frac{1}{2})+\frac{4\mathrm{I\Δt}}{{\mathrm{\ϵ\π\Δ\ρ}}^2}---\left(7\right)$

其中，I为电流大小，π为圆周率；根据公式（4）-（7）可得出任意观测点每个迭代时刻的垂直电场E_Z，再根据垂直电场E_z与地波传播衰减因子W_g满足的关系式：

E_z＝E₀W_g（8）

即可求出地波传播衰减因子W_g；其中E₀代表理想导电平面上相同接收距离的电场，

I为电流大小，l为电偶极子的电荷间的距离，θ为观测点至场源的大圆角距离，λ为波长，a为地球半径，k₀为自由空间波数，e代表自然对数的底，i为虚数单位，π为圆周率。

4.根据权利要求1所述高精度预测地波传播衰减因子的方法，其特征在于：所述的地波传播衰减因子频率范围为30Hz-300MHz。

Images