Terrain optimization method for low-frequency ground-wave propagation delay calculation
-
摘要: 附加二次相位因子(ASPF)是影响长波信号传播时间延迟的主要因素. 积分方程方法计算ASPF适用于地形复杂的路径，是目前广泛使用的计算方法. 它在推导过程中使用了稳定相位积分的近似方法，在地形起伏剧烈地区不再适用. 实际传播路径地形起伏变化巨大，通过数学形态法对实际传播路径修正，可以使路径保留基本几何轮廓并且光滑. 积分方程方法计算结果表明: 利用数学形态法可以有效优化处理实际传播路径.
-
关键词:
- 积分方程方法 /
- 附加二次相位因子(ASPF) /
- 数学形态学 /
- 长波时延预测 /
- 地形优化处理
Abstract: The additional secondary phase factor (ASPF) plays an important role in propagation delay of the low-frequency ground-wave. A method of calculating attenuation for ground-wave propagating over irregular terrain called integral function has been developed recently. It is derived by means of a stationary-phase integration that reduces the dimensionality of the general version, but such an approximation is not valid for all terrain types. The terrain of the actual propagation path changes greatly, and the actual propagation path is corrected by the mathematical morphology method to keep the basic geometric contour and smoothness of the path. The calculation results of the integral equation method show that the actual propagation path can be effectively optimized through the mathematical morphology method. -
0. 引 言
近年来国际上越来越重视罗兰C (Loran-C)陆基远程长波定位导航和授时系统(PNT)，它与卫星导航系统独立，同时也是卫星导航系统的备份. Loran-C系统发射峰值功率达到2 MW，信号传播距离达到几千公里，具有很强的抗干扰能力和很广泛的覆盖范围. Loran-C信号是波长为3 km的低频长波电磁波，因其沿地球表面传播又称为低频地波.
附加二次相位因子(ASPF)定义为低频地波传播路径，是陆地时对传播时延产生的影响，ASPF已经成为影响Loran-C系统PNT精度的主要因素. 计算ASPF要通过求解麦克斯韦方程组在满足边界条件的解，这是非常复杂的计算过程. 对低频地波的传播研究方兴未艾，最早由熊皓[1]研究得出均匀光滑平地面模型下地波场的积分表达式；Fock[2]得出均匀光滑球面模型下地波场的级数表达式；Wait[3]、Millington[4]等得出分段均匀光滑地面模型下的计算公式，Millington公式在工程中被广泛应用. 《中华人民共和国电子行业军用标准》中长波地波传输信道计算方法采用的就是以上三种模型与计算方法.
低频地波的传播路径通常包含不同的地形地质，如平原、山脉、沼泽和沙漠等，这些不同的地形地质其地面电参数也不同. 如大地电导率、相对介电常数，地表大气折射率等，对信号接收点的场强与时延的影响也不尽相同，Millington方法只能将传播路径粗略分为均匀光滑若干段，无法考虑地形因素对低频地波的影响. Hufford[5]利用第二格林定理得出了低频地波在不均匀不光滑传播模型中地波场计算的积分方程方法，经过实际测量证明其计算精度较高，目前被广泛应用. 周丽丽等[6]也对其进行了理论研究，并且计算分析了不同形状山脉对低频地波传播性能的影响.
积分方程方法具有很高的计算精度，但是积分方程方法在推导过程中使用了稳定相位积分近似，导致其无法满足极端苛刻的地形. 只适用于地形缓变情况，即传播路径高程数据不含有“噪声”，是一条光滑的曲线. 因此计算过程中首先需要获取高分辨率的传播路径高程数据，一方面可以避免因为积分步长远大于地形采样间隔导致地形数据产生突变；另一方面可以避免计算高程数据的一阶导数和二阶导数误差太大.
地形高程数据通过编写程序处理SRTM3高程数据文件，根据收发点的经纬度和采样点数来获取. 因其固有的离散采样特性与地形本身陡峭等因素，原始路径高程数据不满足积分方程方法缓变地形的要求. 因此我们需要对低频地波传播路径采取一定的“平滑”处理，保留地形的“包络”，即保留地形的几何轮廓，去掉“噪声”. 本文创新性的使用了数学形态法对原始高程数据经行处理，使传播路径高程数据满足积分方程方法的要求.
1. 数学形态学
数学形态学是于1964年由法国数学家Matheron[7]和Serra[8]共同提出的基于集合论、积分几何、拓扑学等数学基础的交叉性学科. 近年来通过数学家不断的研究，数学形态学的理论基础已经逐步完善，相关应用也深入很多领域. 数学形态学的基础算子只有简单的加减法和取极值运算，没有乘除法运算，这样大大提高了运算效率，运算速度很快.
最初数学形态学是应用于图像处理，它是利用结构元素对图像进行变换，应用于去除图像噪声和轮廓边缘检测等方面.
数学形态学在信号处理方面提供了一种可以检测信号几何特征的非线性信号处理方式，被处理信号的几何形状信息可以由作用在信号上的结构元素提取出来，这时的结构元素相当于一个“探针”，在一维或二维信号中不断移动这个探针，就可考察信号波形中的几何关系. 数学形态学已在图像处理、计算机视觉等领域得到广泛应用，目前电力系统、信号处理等领域中也逐渐展开研究与应用.
数学形态学有两个基本运算：腐蚀和膨胀. 腐蚀可以使目标区域变小，造成目标边界收缩，用来消除小且无意义的目标；膨胀会使目标区域变大，造成目标边界扩大，减小目标区域内的谷域以及消除包含在目标区域中的噪声. 将腐蚀运算看作是最小值滤波器，膨胀运算看作是最大值滤波器，它们可分别获得数据的下包络和上包络. 数学形态学基本算子主要包括膨胀、腐蚀及以此为基础构造的开运算、闭运算 4 种运算方法. 定义原始信号
$z(n)$ 为在$Z = (0,1, \cdots, N - $ $ 1)$ 上的离散函数，$g(n)$ 为在$G = (0,1, \cdots ,M - 1)$ 上的离散函数. 腐蚀和膨胀运算定义为：$$\left(z\ominus g\right)\left(n\right)={\rm{min}}[z\left(n+m\right)-g(m\left)\right]{,}$$ (1) $$ \left(z\oplus g\right)\left(n\right)={\rm{max}}[z\left(n+m\right)+g(m\left)\right]{.} $$ (2) $\ominus $ 表示腐蚀运算，$\oplus $ 表示膨胀运算. 开运算“$\circ $ ”和闭运算“$\bullet$ ”定义为：$$z \circ g = (z \ominus g) \oplus g,$$ (3) $$z \bullet g = (z \oplus g) \ominus g.$$ (4) 2. 长波信号传播时延计算
由于低频地波是沿地面传播的，因此接收点的相位延迟不仅与发射站的频率有关，还与发射站与接收站之间的传播距离有关.其中还有一个不可忽略的影响因素即地面. 地面对低频地波传播的影响主要有：一是地面的粗糙度，如山脉等地形起伏的影响；另一个是地质的电磁特性，在地质结构中由于传播介质时空分布不均匀，如低频地波的传播路径在若干年之后由于地表水系或者水土流失等发生改变. 由于地形和大地电导率等因素对低频地波传播路径的影响，实际无线电波信号的传播过程中会发生一系列复杂的变化. 例如：反射、衍射、散射和折射，导致传播速度和相位改变，最终导致导航和授时产生误差，这种误差在授时系统可以达到μm级别，定位系统的误差可能达到数百米到几千米级别.
$$T = {\rm{PF}} + {\rm{SF}} + {\rm{ASPF}},$$ (5) $${\rm{PF}} = \frac{{{n_s}d}}{c} \times {10^6},$$ (6) $${\rm{SF}} + {\rm{ASPF}} = \frac{{\arg ({W_g})}}{\omega } \times {10^6}.$$ (7) 式中：PF为发射站与接收站之间的直线距离和大气折射率对电波传播造成的时延；PF又称为一次相位因子；SF为发射站与接收站之间为纯海水的情况下对电波造成的时延，SF又称为二次相位因子；ASPF为发射站与接收站之间为纯陆地对电波传播造成的时延减去若该路径为纯海水对电波传播造成的时延，即传播路径为陆地相对于纯海水的时延；ASPF又称为附加二次相位因子. PF、SF、ASPF的单位为μs；d是以km为单位的电波传播的距离；arg为求衰减因子幅角；
${n_s}$ 为地表附件大气折射率，约为1.000 338；c为光速；$\omega = 2{\rm{\text{π} }}f$ ，f为电波频率；W为地波衰减函数，是计算接收站地波场强和时延的关键因素，与电波频率，传播距离，大地电导率，传播路径地形等有关，目前低频地波传播信道模型及W的计算主要有以下四种方式：1) 均匀光滑平地面模型
均匀地面指的是电波传播路径上地质结构趋于一致，地面电性参数变化不大，比如大地电导率等可以取一个定值. 光滑地面指的是地形起伏变化程度远小于电波波长，如100 kHz的低频地波波长为3 km，关中平原就可以看作是光滑地面. 当电波传播路径小于60～70 km时，可以把传播路径当做平地面. 继法拉第电磁感应定理的发现到麦克斯韦方程组的建立，垂直电偶极子的辐射场求解问题颇有历史渊源[9]，Sommerfeld[10]通过研究平地面上垂直电偶极子产生的场得出了一个积分表达式，经研究人员不断改进，使其可以在工程中应用.
2) 均匀光滑球面模型
低频地波绕射能力很强，可以沿地面传播很远的距离，当低频地波的传播距离很长时就必须要考虑地球曲率对地波传播的影响. 均匀光滑球地面模型将地球理想化为一个表面光滑，地面电性参数为常数的球体. 地面的发射天线比波长小得多，可以理想化为垂直电偶极子. Watson[11]、Fock[2]等科研人员在上世纪初做了大量研究，得出了垂直电偶极子在球面上的辐射场的积分表达式，并且经过复杂数学运算后将其变成一个收敛相当快的级数表达式[3]，文中取为
${{\rm{e}}^{ - i\omega t}}$ ，均匀光滑球面模型的地波衰减函数计算公式为$${W_g} = {{\rm{e}}^{i\frac{\text{π} }{4}}}\sqrt {\text{π} x} \sum\limits_{s = 1}^\infty {{{\rm{e}}^{ix{t_s}}}} \frac{1}{{{t_s} - {q^2}}}.$$ (8) 式中：
$x = {\left(\displaystyle\frac{{{k_0}a}}{2}\right)^{\frac{1}{3}}}\theta$ ，${k_0}$ 为电磁波在空气中的波数；a为地球半径；$\theta $ 为发射站与接收站之间的大圆角距离；${t_s}$ 为微分方程$\displaystyle\frac{{{\rm{d}}t}}{{{\rm{d}}q}} = \frac{1}{{t - {q^2}}}$ 的根，$q = i{\left(\displaystyle\frac{{{k_0}a}}{2}\right)^{\frac{1}{3}}}{\Delta _g}$ ，${\Delta _g}$ 为归一化地表面阻抗.3) 分段均匀光滑地面模型
通常情况下地面既不均匀又不平坦，因此从发射站到接收站的整个传播路径不能看成是均匀光滑的传播路径. 比如海洋与陆地的电导率差异很大，沙漠、丘陵、平原或沼泽等陆地的大地电导率也不尽相同，这样就不能把整个传播路径看成是均匀的，只能把陆地部分和海洋部分看成两段均匀路径. 陆地部分再根据地形地貌细分为若干段，每一段当做均匀光滑路径，大地电参数在每一段取不同的常数. 分段均匀光滑地面的地波衰减函数计算公式有Wait公式[3]、波模转换法、抛物方程法和Millington公式[4]等. Millington公式简单实用，在工程中广泛应用，传播信道模型如图1所示.
设传播路径分为N段，每段长度为
${d_i}$ ，归一化表面阻抗为${\Delta _i}$ ，地波衰减因子为$$ {W}_{g}=\sqrt{{W}_{\text{正}}{W}_{\text{反}}}.$$ (9) 式中：
$$\begin{aligned}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! {W_{\text{正}}}{\rm{ = }}\frac{{{W_g}({d_1},{\Delta _1}){W_g}({d_1}{\rm{ + }}{d_2},{\Delta _2}) \cdots }}{{{W_g}({d_1},{\Delta _2}){W_g}({d_1}{\rm{ + }}{d_2},{\Delta _3}) \cdots }}{\rm{ }}\frac{{{W_g}({d_1}{\rm{ + }} \cdots {\rm{ + }}{d_{N - 1}},{\Delta _{N - 1}})}}{{{W_g}({d_1}{\rm{ + }} \cdots {\rm{ + }}{d_{N - 1}},{\Delta _N})}}{\rm{ }}{\rm{,}} \end{aligned}$$ (10) $$\begin{split} {W_{\text{反}}}{\rm{ = }}&\frac{{{W_g}({d_N},{\Delta _N}){W_g}({d_N}{\rm{ + }}{d_{N - 1}},{\Delta _{N - 1}}) \cdots }}{{{W_g}({d_N},{\Delta _{N - 1}}){W_g}({d_N}{\rm{ + }}{d_{N - 1}},{\Delta _{N - 2}}) \cdots }}\\ &\frac{{{W_g}({d_N}{\rm{ + }} \cdots {\rm{ + }}{d_1},{\Delta _2})}}{{{W_g}({d_N}{\rm{ + }} \cdots {\rm{ + }}{d_2},{\Delta _1})}}{\rm{ }}{\rm{.}} \end{split}$$ (11) 4) 不均匀不光滑地面模型
实际电波传播路径不仅地质类型不同，同时也存在山脉等地形起伏地区，即不均匀不光滑地面. 积分方程方法适用于均匀光滑平地面模型、均匀光滑球面模型、分段均匀地面模型，与它们具有同样的计算精度，还可以计算出地波传播路径上地形起伏变化产生的影响，是应用非常广泛的一种算法，国内外对复杂路径低频地波传播预测多用此方法. 传播信道模型如图2所示.
假设地面起伏情况如图2中曲线所示，其满足缓变条件. 地波衰减因子计算公式为
$$\begin{split} {{W}_{g}}(P)=&A\Bigg\{1+\displaystyle\frac{i{{k}_{0}}}{2\text{π}}\iint\limits_{{{s}_{0}}}{\frac{{{W}_{g}}(Q){{\rm{e}}^{i{{k}_{0}}({{r}_{1}}+{{r}_{2}}-{{r}_{0}})}}}{{{r}_{1}}{{r}_{2}}}} \\ & \left[\Delta g{{(1+{{(z\text{′})}^{2}})}^{\frac{1}{2}}}+\left(1+\frac{i}{{{k}_{0}}{{r}_{2}}}\right)\frac{\partial {{r}_{2}}}{\partial n}\right]{\rm{d}}s\Bigg\}. \\ \end{split}$$ (12) 信号发射源位于原点O，ldl为垂直电偶极子，接收点为P，地面上的动点为Q，D为从源点到接收点的大圆距离，r0为从源点到接收点P的直线距离，S0为积分区域，∆g为地球归一化表面阻扰，r1为从源点到积分动点Q的直线距离，r2为从积分动点Q到接收点P的直线距离.
3. 计算结果分析
3.1 传播路径地形高程信息提取
航天飞机雷达地形测绘使命（SRTM），是2000年美国国家地理空间情报局（NGA）与美国国家航空航天局（NASA）以及德国和意大利的航天机构的一项联合项目，该项目使用航天飞机进行了为期11天的测量，目的是为地球80%的陆地表面（60°N～56°S的所有陆地）生成数字地形高程数据，在90%置信度下，高程数据的绝对垂直精度为16 m，是目前使用最广泛的高程数据源之一. SRTM数据按照分辨率分为SRTM3和SRTM1两种，SRTM3分辨率为3′′，即90 m，SRTM1分辨率为1′′，即30 m，但是SRTM1数据只包含美国全境，本文所用数据为SRTM3.
3.1.1 SRTM3数据结构及读取
SRTM3数据包含很多文件，每一个文件覆盖地球表面一个纬度乘一个经度块，每个文件里面包含1 201×1 201个采样点的高度数据. 经纬度网格与SRTM3文件对应关系如图3所示.
由图3可知，每个经纬度网格与一个SRTM3文件对应，SRTM3文件名称包含经纬度信息，根据每个经纬度网格左下角的经纬度值构造SRTM3文件名. 比如图3中包含图钉的网格，其左下角坐标为(33°N，107°E)，则这个网格对应的SRTM3文件名称为“N34E107.hgt”. 此文件可以看作为1 201×1 201的矩阵，从左往右的方向为经度增加方向，从下到上的方向为纬度增加的方向，相邻文件最外侧的行与列重叠. 计算机中的存储单位是字节，SRTM3是一种hgt格式的文件，hgt来源于英文单词“hight”，文件中高程数据使用16位有符号整数表示，能取到的最大值为32 767，单位是m. 文件中是以大端模式、顺序结构存储高程数据二维数组，如果采用普通的二进制方式读取文件，则读取后处理比较麻烦，不过python的numpy库只要设置好读取模式就可以轻松读取数据，并且将高程数据转化为二维数组.
3.1.2 传播路径高程提取程序
通常电磁波传播路径几百公里，发射站与接收站之间的高程数据包含多个SRTM3文件，要获取这样传播路径的高程信息比较复杂. 比如图3中的实线表示从陕西省渭南市蒲城发射台到四川省什邡市接收点的电磁波传播路径，路径经过了9个经纬度块，因此需要9个SRTM3文件，按照所需要的SRTM3文件数量对传播路径进行分割，分割后在每个SRTM3文件上提取高程数据，最后再拼接为完整的传播路径. 算法逻辑主要是根据采样点数确定采样步长，然后按照采样步长计算各个采样点的经纬度值，将传播路径离散化，生成采样点经纬度数组. 关键步骤就是构造SRTM3文件名称，各个采样点取其经纬度值整数部分构造其所属的SRTM3文件名称，这一步会产生很多重复名称，只要去除重复文件名称就可以得到传播路径所需的所有SRTM3文件. 由每个SRTM3文件覆盖的经纬度都是1°，经纬度值整数部分相同的采样点自然归属于同一个文件，因此只要取数组中重复的文件名称的下标，就可以得到各个SRTM3文件中的采样开始和采样结束点. 单个文件内提取高程数据的算法逻辑主要为计算偏移量，SRTM3分辨率为3 s，即
$\left(\displaystyle\frac{1}{{1200}}\right)$ °，例如“N34E107.hgt”文件，此文件中左下角为原点，下标从0开始的第i行第j列的点其值代表高程数据，其对应的纬度值为$34 + \displaystyle\frac{i}{{1200}}$ ，经度值为$107 + \displaystyle\frac{j}{{1200}}$ . 程序整体执行过程如下图4.3.2 数学形态学对高程信息处理
选择图3中陕西省渭南市蒲城发射台与四川省什邡市接收点之间的传播路径，使用上节编写的提取高程数据程序，输入发射台与接收点的经纬度和采样点数，得到路径高程信息如图5所示.
由图5可知，路径总长约650 km，从发射站出发，先经过约150 km关中平原，然后翻越秦岭山脉，最后到达约450 km的四川盆地. 整体传播路径的地形变化为先平缓后起伏最后趋于平缓. 关中平原段平缓，满足积分方程计算要求，四川盆地段有少量“噪声”，地形起伏变化主要集中在了秦岭山脉，这两段需要使用数学形态法处理，使其既保留山的轮廓又去除“噪声. 处理结果如图6所示.
由图6可知，经过数学形态学处理后，传播路径在关中平原段没有明显变化，在四川盆地和秦岭段变化较大，四川盆地和秦岭山脉段变得光滑，秦岭山脉轮廓清晰，高度几乎没有损失，证明利用数学形态法处理后的传播路径很好地满足了积分方程方法对地形的要求.
3.3 积分方程计算结果
积分方程方法分别计算上述传播路径ASPF，图7为原始传播路径ASPF计算结果，图8为数学形态学处理过后的传播路径计算结果.
由图7可知，原始路径计算的ASPF在关中平原段比较光滑，四川盆地和秦岭段脉冲噪声较多且噪声幅度较大.
由图8可知，数学形态学处理过后的路径计算的ASPF在四川盆地和秦岭段脉冲噪声明显减少且幅度也显著降低.
4. 总 结
积分方程方法在推导长波信号传播时延过程中采用了稳定相位积分等近似条件，因此其计算过程要求地面是缓变的，地面曲率半径不能太小. 数学形态学处理过后的地形能够保留其原有几何轮廓去除“噪声”，计算速度快，使地面满足缓变的条件. 积分方程方法在处理过后的传播路径上计算的ASPF扰动与跳变减少，趋于平滑，表明数学形态学对于处理低频地波传播路径数据是很有效的.
-
[1] 熊皓. 无线电波传播[M]. 北京: 电子工业出版社, 2000. [2] FOCK V. Diffraction of radio waves around the earth's surface[J]. Journal of physical, 1945(15): 479-496.
[3] WAIT J R, HOUSEHOLDER J. Mixed-path ground wave propagation, II. larger distances[J]. Journal of research of the national bureau of standards, 1957(59): 19-26. DOI: 10.6028/jres.059.003
[4] MILLINGTON G. Ground-wave propagation over an inhomogeneous smooth earth[J]. Journal of the institution of electrical engineers, 1949, 96(39): 53-64. DOI: 10.1049/pi-3.1949.0013
[5] HUFFORD G A. An integral equation approach to the problem of wave propagation over an irregular surface[J]. Quarterly of applied mathematics, 1952, 9(4): 391-404. DOI: 10.1090/qam/44350
[6] 周丽丽, 穆中林, 蒲玉蓉, 等. 不规则地形地波传播衰减因子的改进算法及结果一致性研究[J]. 电子与信息学报, 2015, 37(9): 2254-2259. [7] MATHERON G. Random sets and integral geometry[M]. New York: Wiley, 1975.
[8] SERRA J. Introduction to mathematical morphology[J]. Computer vision, graphics, and image processing, 1986, 35(3): 283-305. DOI: 10.1016/0734-189X(86)90002-2
[9] 潘威炎. 长波超长波及长波传播[M]. 成都: 电子科技大学出版社, 2004. [10] SOMMERFELD A N. The propagation of waves in wireless telegraphy[J]. Annals of physics, 1926, 386（25）: 1135-1153. DOI: 10.1002/andp.19263862516
[11] WATSON G N. The diffraction of radio waves by the earth[J]. Royal sociaty, 1918, 666(95): 83-99. DOI: 10.1098/rspa.1918.0050
-
期刊类型引用(1)
1. 梁堃，胡昀. 一种便携式葡萄剪枝机器设计与实现. 自动化与仪器仪表. 2023(05): 296-300 . 
百度学术
其他类型引用(5)
下载:
计量
- 文章访问数: 464
- HTML全文浏览量: 207
- PDF下载量: 30
- 被引次数: 6
