低频地波传播时延计算的地形优化处理方法

朱峙亚，郭伟

( 1.中国科学院大学，北京 100049；2.中国科学院国家授时中心，西安 710600)

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,···,N−1) 上的离散函数，g(n)为在G=(0,1,···,M−1)上的离散函数. 腐蚀和膨胀运算定义为：

2 长波信号传播时延计算

由于低频地波是沿地面传播的，因此接收点的相位延迟不仅与发射站的频率有关，还与发射站与接收站之间的传播距离有关.其中还有一个不可忽略的影响因素即地面.地面对低频地波传播的影响主要有：一是地面的粗糙度，如山脉等地形起伏的影响；另一个是地质的电磁特性，在地质结构中由于传播介质时空分布不均匀，如低频地波的传播路径在若干年之后由于地表水系或者水土流失等发生改变.由于地形和大地电导率等因素对低频地波传播路径的影响，实际无线电波信号的传播过程中会发生一系列复杂的变化.例如：反射、衍射、散射和折射，导致传播速度和相位改变，最终导致导航和授时产生误差，这种误差在授时系统可以达到µm 级别，定位系统的误差可能达到数百米到几千米级别.

式中：PF为发射站与接收站之间的直线距离和大气折射率对电波传播造成的时延；PF又称为一次相位因子；SF为发射站与接收站之间为纯海水的情况下对电波造成的时延，SF又称为二次相位因子；ASPF为发射站与接收站之间为纯陆地对电波传播造成的时延减去若该路径为纯海水对电波传播造成的时延，即传播路径为陆地相对于纯海水的时延；ASPF又称为附加二次相位因子.PF、SF、ASPF的单位为µs；d是以km 为单位的电波传播的距离；arg 为求衰减因子幅角；ns为地表附件大气折射率，约为1.000 338；c为光速；ω =2πf，f为电波频率；W为地波衰减函数，是计算接收站地波场强和时延的关键因素，与电波频率，传播距离，大地电导率，传播路径地形等有关，目前低频地波传播信道模型及W的计算主要有以下四种方式：

1)均匀光滑平地面模型

均匀地面指的是电波传播路径上地质结构趋于一致，地面电性参数变化不大，比如大地电导率等可以取一个定值.光滑地面指的是地形起伏变化程度远小于电波波长，如100 kHz 的低频地波波长为3 km，关中平原就可以看作是光滑地面.当电波传播路径小于60～70 km 时，可以把传播路径当做平地面.继法拉第电磁感应定理的发现到麦克斯韦方程组的建立，垂直电偶极子的辐射场求解问题颇有历史渊源[9]，Sommerfeld[10]通过研究平地面上垂直电偶极子产生的场得出了一个积分表达式，经研究人员不断改进，使其可以在工程中应用.

2)均匀光滑球面模型

低频地波绕射能力很强，可以沿地面传播很远的距离，当低频地波的传播距离很长时就必须要考虑地球曲率对地波传播的影响.均匀光滑球地面模型将地球理想化为一个表面光滑，地面电性参数为常数的球体.地面的发射天线比波长小得多，可以理想化为垂直电偶极子.Watson[11]、Fock[2]等科研人员在上世纪初做了大量研究，得出了垂直电偶极子在球面上的辐射场的积分表达式，并且经过复杂数学运算后将其变成一个收敛相当快的级数表达式[3]，文中取为e−iωt，均匀光滑球面模型的地波衰减函数计算公式为

3)分段均匀光滑地面模型

通常情况下地面既不均匀又不平坦，因此从发射站到接收站的整个传播路径不能看成是均匀光滑的传播路径.比如海洋与陆地的电导率差异很大，沙漠、丘陵、平原或沼泽等陆地的大地电导率也不尽相同，这样就不能把整个传播路径看成是均匀的，只能把陆地部分和海洋部分看成两段均匀路径.陆地部分再根据地形地貌细分为若干段，每一段当做均匀光滑路径，大地电参数在每一段取不同的常数.分段均匀光滑地面的地波衰减函数计算公式有Wait 公式[3]、波模转换法、抛物方程法和Millington 公式[4]等.Millington 公式简单实用，在工程中广泛应用，传播信道模型如图1所示.

图1 分段均匀光滑地面模型图

4)不均匀不光滑地面模型

实际电波传播路径不仅地质类型不同，同时也存在山脉等地形起伏地区，即不均匀不光滑地面.积分方程方法适用于均匀光滑平地面模型、均匀光滑球面模型、分段均匀地面模型，与它们具有同样的计算精度，还可以计算出地波传播路径上地形起伏变化产生的影响，是应用非常广泛的一种算法，国内外对复杂路径低频地波传播预测多用此方法.传播信道模型如图2所示.

图2 不光滑不均匀地面模型图

假设地面起伏情况如图2中曲线所示，其满足缓变条件.地波衰减因子计算公式为

信号发射源位于原点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 经纬度网格与SRTM 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传播路径高程提取程序

图4 高程信息提取程序执行框图

3.2 数学形态学对高程信息处理

选择图3中陕西省渭南市蒲城发射台与四川省什邡市接收点之间的传播路径，使用上节编写的提取高程数据程序，输入发射台与接收点的经纬度和采样点数，得到路径高程信息如图5所示.由图5可知，路径总长约650 km，从发射站出发，先经过约150 km 关中平原，然后翻越秦岭山脉，最后到达约450 km 的四川盆地.整体传播路径的地形变化为先平缓后起伏最后趋于平缓.关中平原段平缓，满足积分方程计算要求，四川盆地段有少量“噪声”，地形起伏变化主要集中在了秦岭山脉，这两段需要使用数学形态法处理，使其既保留山的轮廓又去除“噪声.处理结果如图6所示.

图5 原始传播路径高程信息图

图6 数学形态法优化后路径高程信息图

由图6可知，经过数学形态学处理后，传播路径在关中平原段没有明显变化，在四川盆地和秦岭段变化较大，四川盆地和秦岭山脉段变得光滑，秦岭山脉轮廓清晰，高度几乎没有损失，证明利用数学形态法处理后的传播路径很好地满足了积分方程方法对地形的要求.

3.3 积分方程计算结果

积分方程方法分别计算上述传播路径ASPF，图7为原始传播路径ASPF计算结果，图8为数学形态学处理过后的传播路径计算结果.

图7 原始路径计算结果图

图8 优化后路径计算结果图

由图7可知，原始路径计算的ASPF在关中平原段比较光滑，四川盆地和秦岭段脉冲噪声较多且噪声幅度较大.

由图8可知，数学形态学处理过后的路径计算的ASPF在四川盆地和秦岭段脉冲噪声明显减少且幅度也显著降低.

4 总结

积分方程方法在推导长波信号传播时延过程中采用了稳定相位积分等近似条件，因此其计算过程要求地面是缓变的，地面曲率半径不能太小.数学形态学处理过后的地形能够保留其原有几何轮廓去除“噪声”，计算速度快，使地面满足缓变的条件.积分方程方法在处理过后的传播路径上计算的ASPF扰动与跳变减少，趋于平滑，表明数学形态学对于处理低频地波传播路径数据是很有效的.