不規(guī)則地形地波傳播衰減因子的改進(jìn)算法及結(jié)果一致性研究

周麗麗 穆中林 蒲玉蓉 席曉莉

1 引言

低頻地波被廣泛用于導(dǎo)航授時(shí)系統(tǒng)。近年有關(guān)低頻地波新的研究應(yīng)用方向主要有：淺水導(dǎo)航、土壤水分探測、新的捕獲降噪方法實(shí)現(xiàn)等方面[13]-。地波信號在傳播過程中受到傳播路徑的大地電參數(shù)、地形起伏及大氣折射指數(shù)等因素的影響，其傳播特性難以精確預(yù)測，從而導(dǎo)致該類系統(tǒng)精度降低[47]-。

目前基于均勻光滑路徑模型的理論預(yù)測方法主要有留數(shù)級數(shù)法，基于復(fù)雜（分段均勻光滑/不均勻不光滑）路徑模型的預(yù)測方法主要有 Millington方法、Wait積分方法、拋物線方程方法以及積分方程方法等，其中能考慮地形影響的算法主要是積分方程方法[8,9]。其最早由文獻(xiàn)[10]1952年基于格林第 2定理提出；1970年文獻(xiàn)[11]根據(jù)拋物線方程得到的另外一種積分形式；1988年文獻(xiàn)[12]則基于補(bǔ)償原理得到更為一般意義的積分形式。研究表明，這些方法在低頻段具有相同的結(jié)果[9]。由于能夠考慮傳播路徑實(shí)際地形起伏變化和電參數(shù)分布，理論上該方法比其它經(jīng)典算法預(yù)測精度更高，且具有普適性。

對于復(fù)雜地形路徑，目前國內(nèi)普遍采用的積分方程公式[1317]-與第1種積分公式相同或相近，衰減因子預(yù)測結(jié)果多是基于此方法。但即使在均勻光滑路徑模型下，文獻(xiàn)[14,15,17]給出的結(jié)果在地面接收距離大于 300 km時(shí)會與留數(shù)級數(shù)法的結(jié)果存在較大偏差，且傳播距離越遠(yuǎn)，誤差越大，在1700 km的距離上二者預(yù)測的傳播相位差別甚至高達(dá) 16 μs[17]。上述文獻(xiàn)中沒有給出造成此項(xiàng)誤差的原因。

本文針對積分算法與其它經(jīng)典算法預(yù)測結(jié)果不一致的問題進(jìn)行討論，在對不同路徑模型地波衰減因子的定義、物理含義及相關(guān)算法公式研究基礎(chǔ)上，給出了積分方程算法的推導(dǎo)過程，分析了誤差產(chǎn)生原因，并對其進(jìn)行改進(jìn)，使之與其它理論算法預(yù)測結(jié)果趨于一致。下文時(shí)間因子統(tǒng)一取 e-jωt。

2 積分方程方法存在問題及算法改進(jìn)

2.1 典型地面模型下地波衰減因子求解算法

在地波傳播理論預(yù)測中，當(dāng)觀察點(diǎn)及場源均離地面不遠(yuǎn)，則觀察點(diǎn)處的電場垂直分量可表示為[8,18]

其中

a為等效地球半徑，dl為電流元長度，θ為觀察點(diǎn)至場源的大圓角距離，Wg稱為地波衰減因子，它是傳播路徑地形及電參數(shù)的復(fù)雜函數(shù)。在實(shí)際工程應(yīng)用中，當(dāng)已知的是發(fā)射臺的輻射功率 Pr，以kW計(jì)。若傳播路徑的角距離 θ ≤π/12，則sinθ≈θ,a≈aθ=d,d為大圓距離，如果d以km計(jì)，則 E0可表示為

E0物理含義是理想導(dǎo)電平面條件下，場源（垂直電偶極子）及接收點(diǎn)均在地面，且傳播距離為d時(shí)接收點(diǎn)處電場垂直于地面的分量。

2.1.1 均勻光滑球地面 若地面為均勻光滑球地面，已知地面相對介電常數(shù)為εr，電導(dǎo)率為σ，留數(shù)級數(shù)算法中地波衰減因子Wg可表示為[4,5]

此處

h1和h2分別為發(fā)射點(diǎn)和接收點(diǎn)離地面高度，w2（ t）是愛里函數(shù)。式中參數(shù)q可表示為

Δg為歸一化地表面阻抗，可表示為

ts是微分方程 w'2（t） - q w2（ t）= 0 的第s個(gè)根。

2.1.2 分段均勻光滑球地面 當(dāng)?shù)孛鏋榉侄尉鶆蚬饣虻孛?，設(shè)從發(fā)射點(diǎn)到接收點(diǎn)的傳播路徑分為 N個(gè)有限段，每段長度分別為： d1,d2,…,dN，歸一化表面阻抗分別為Δg1,Δg2,…,ΔgN。對應(yīng)的地波衰減因子求解方法目前主要有 4種：Wait積分法、Millington方法、波模轉(zhuǎn)換法以及拋物線方程方法，下面僅給出拋物線方程方法相關(guān)公式。

設(shè)電波傳播的介質(zhì)環(huán)境沿方位角φ沒有變化，且大氣折射指數(shù)沿高度方向不變，令磁場Hφ=[ε/ （r sin θ）]1/2ejkaθu（ r,θ），進(jìn)行坐標(biāo)變換，則可以推出u（ x, z）滿足式（10）的拋物線方程[19]：

其滿足邊界條件

其中

當(dāng)已知x處地面及空中的u（ x, z），式（10）可以通過分步混合傅里葉變換或有限差分法求解下一段x +Δx 處的任意高度的 u（ x + Δx , z），進(jìn)而不斷迭代得到任意一點(diǎn)的u（ x, z）。再根據(jù)式（13）求解出地波衰減因子

其中 H （x, 0） = - jk Id lejkx/2πx 。迭代初值可以由近0距離的平地面公式求解得到。

2.1.3 不均勻不光滑地面 若電波傳播路徑上地面電導(dǎo)率及地形都有復(fù)雜變化，此時(shí)地波衰減因子可以用積分方程[4,13]的數(shù)值解來求得。記發(fā)射源點(diǎn)的位置為O，觀察點(diǎn)的位置為P，地面上積分動點(diǎn)為Q，傳播路徑幾何示意圖如圖1所示。地波衰減因子滿足的積分方程為[8,18]

圖1 積分方程方法中傳播路徑幾何示意圖

其中 r1表示從源點(diǎn)到地面上積分動點(diǎn)Q之間的直線距離，2r表示從Q點(diǎn)到P的距離， h為所在直角坐標(biāo)系中接收點(diǎn)距離地面的垂直高度。n表示地表面的法線方向。其中x表示觀察點(diǎn)至源點(diǎn)水平距離，L是Q與源點(diǎn)之間的水平距離。地面接收時(shí)A為1.0，空中接收時(shí)A為0.5。

2.2 積分方程方法直接預(yù)測地波衰減因子存在問題

由式（1），式（3）和式（13），光滑路徑和分段均勻光滑路徑模型下相關(guān)算法中地波衰減因子實(shí)質(zhì)上是將實(shí)際接收點(diǎn)電場的垂直分量或磁場的水平分量與良導(dǎo)體平地面上相應(yīng)傳播距離下場分量歸一化的結(jié)果，其考慮了近場感應(yīng)、地球曲率、收發(fā)點(diǎn)高度以及地面有限電導(dǎo)率的影響。為分析采用積分方程方法與其它經(jīng)典算法預(yù)測衰減因子不一致的問題，下面簡述積分方程方法的推導(dǎo)過程[13,17]。

在一個(gè)含有源為垂直電偶子的充滿空氣的封閉區(qū)域V內(nèi)，電磁場的赫茲勢Π滿足標(biāo)量非齊次波動方程

其中J是電流密度矢量，其模用J表示，δ（o）為沖擊函數(shù)，代表源的分布。取直角坐標(biāo)系x,y,z,J的方向?yàn)閦方向。區(qū)域V被大地表面與一個(gè)半徑無窮大的半球面包圍（不包含P點(diǎn)），利用第2格林定理

其中，ψ是一個(gè)格林函數(shù)，定義為

當(dāng)?shù)匦纹鸱淮?，滿足z'2? 1 （z'為地形函數(shù)的一階導(dǎo)數(shù)），根據(jù)邊界阻抗條件

（其中Δg為歸一化地表面阻抗，n為外法線分量），經(jīng)過一系列計(jì)算可以得到

式（19）第1項(xiàng)為源在觀察點(diǎn)貢獻(xiàn)的直達(dá)波，以Π0（P）表示，第2項(xiàng)代表邊界上二次源對觀察點(diǎn)貢獻(xiàn)的場。若地面為一個(gè)無限大良導(dǎo)體平面時(shí)，可以知道此時(shí)P點(diǎn)場為2Π0（P）。此時(shí)衰減因子定義為[13,17]

將式（20）代入式（19），可以得出假定地形沿y方向沒有變化，應(yīng)用穩(wěn)定相位原理及一系列近似，可進(jìn)一步簡化近似為1維積分式（14）。

由式（20）可以看出，積分方程方法中用于求解地波衰減因子的參考場量和以上其它各種方法不同。積分方程方法中參考場量對應(yīng)的電波傳播距離為收發(fā)點(diǎn)之間的直線距離，而其他算法中的參考場量對應(yīng)傳播距離為收發(fā)兩點(diǎn)的地表面上的大圓距離。只有當(dāng)傳播距離很近，且收發(fā)點(diǎn)均在地面上時(shí)，兩種距離才能近似相等。這是積分方程方法與其它算法預(yù)測結(jié)果不一致的首要原因，由此導(dǎo)致即使在均勻光滑球地面模型下，其地波衰減因子預(yù)測結(jié)果與留數(shù)級數(shù)方法存在較大誤差。

再者，采用積分方程方法求解空中接收點(diǎn)時(shí)，離地高度（如圖 1所示）是指直角坐標(biāo)系下離地垂直方向（z方向）高度，而其它算法中求解高空地波衰減因子（如式（6），式（7）和式（10））算法中對應(yīng)的高度均是在球坐標(biāo)系下定義的徑向（r方向）離地高度。該問題會進(jìn)一步造成接收高度越高，與其它理論算法誤差越大的現(xiàn)象，且由于空中預(yù)測算法中離地高度并非海拔高度，算法對于空中傳播預(yù)測結(jié)果應(yīng)用于實(shí)際工程時(shí)難以直接應(yīng)用，需要進(jìn)行高程轉(zhuǎn)換。

2.3 積分方程方法改進(jìn)算法

針對積分方程方法存在的問題，本文提出對積分方法進(jìn)行如下改進(jìn)。

已知傳播路徑隨大圓距離d變化的地形海拔高度 hp（ d），首先將其轉(zhuǎn)換為直角坐標(biāo)系下的地形函數(shù)z（ x） ：其中ha為發(fā)射天線處的海拔高度。將z（ x）代入式（14），求解原積分方程，可得到觀察點(diǎn) Wg（ d, h ）（注意此處h和d分別表示原式（14）中大圓距離和接收高度）。

以 Pg表示改進(jìn)后的衰減因子， hpR表示接收點(diǎn)實(shí)際的海拔高度， dp表示改進(jìn)后對應(yīng)的大圓距離，則

圖2 改進(jìn)算法中地波衰減因子函數(shù)變量幾何示意

從式（24）可以看出本次算法改進(jìn)的方式是對原算法求得的直角坐標(biāo)系中大圓距離為d，離地垂直高度為h處的地波衰減因子，采用球面校正因子和高程轉(zhuǎn)換技術(shù)進(jìn)行校正后得到的實(shí)為大圓距離為dp，海拔高度為 hpR處的地波衰減因子。

由上述推導(dǎo)分析可知直角坐標(biāo)系中相同水平距離下，接收點(diǎn)分別在地面與空中時(shí)，在改進(jìn)算法中對應(yīng)的大圓距離不同。實(shí)際工程應(yīng)用中當(dāng)求解相同經(jīng)緯度，不同海拔高度處的地波衰減因子時(shí)，在原積分算法中對應(yīng)的水平距離或大圓距離也應(yīng)不同。改進(jìn)后的積分算法中，衰減因子包含了近場感應(yīng)、地球曲率、地形起伏、有限電導(dǎo)率以及收發(fā)高度因素等對地波傳播特性的影響。同時(shí)，由于改進(jìn)后的衰減因子中 hpR是觀察點(diǎn)海拔高度，亦更適合工程應(yīng)用。

3 積分方程改進(jìn)算法的數(shù)值結(jié)果驗(yàn)證

該節(jié)對積分方程改進(jìn)算法與其他方法結(jié)果的一致性進(jìn)行了驗(yàn)證。

圖3給出了接收高度 hpR分別為0 km和10 km條件下，采用積分方程、留數(shù)級數(shù)，拋物線方程方法計(jì)算的均勻光滑球地面的結(jié)果，電參數(shù)選擇為：σ= 3 × 1 0-3s/m, εr= 1 5。地球半徑取等效地球半徑 ae= 4 a/3，其中 a = 6 371.12 km。

從圖3可以看出，經(jīng)改進(jìn)后，積分方程方法結(jié)果與其它方法預(yù)測結(jié)果一致性非常好，均勻光滑地面條件下在1000 km范圍內(nèi)與留數(shù)級數(shù)法幅度誤差小于0.2 dB，相位誤差小于50 ns。

圖4給出了不同接收高度（0 km, 5 km及10 km）條件下，采用積分方程算法、拋物線方程方法、波模轉(zhuǎn)換法以及Millintong算法計(jì)算的陸海兩段均勻光滑球地面的結(jié)果。其中Millintong算法只適用于地面接收的情況，因此只給出了地面接收結(jié)果。第1段路徑長 400 km，電參數(shù)為： σ = 1 0-3s/m,εr=10。第2段路徑長600 km，電參數(shù)為：σ=5 s/m, εr= 8 0。地球半徑仍取等效地球半徑 ae=4a/ 3。

從圖4可以看出，分段均勻光滑條件下改進(jìn)積分方程方法與拋物線方程、波模轉(zhuǎn)換方法預(yù)測結(jié)果一致性也非常好，Millington算法與其它算法誤差略大。

圖 5（a）所示的是一條實(shí)際復(fù)雜傳播路徑上測量點(diǎn)分布區(qū)域圖。測量信號為宣城臺（O 點(diǎn)所示）播發(fā)的導(dǎo)航信號，沿該路徑共測量了8個(gè)點(diǎn)的相位，其中地面上測量了6個(gè)點(diǎn)，海拔為4.26 km的高空測量了1個(gè)點(diǎn)，海拔為6.40 km處測量了1個(gè)點(diǎn)。圖5（b）所示的是沿該路徑地形高程及 8個(gè)測量點(diǎn)空間位置分布，其中1～6為地面測量點(diǎn)，7, 8為兩個(gè)空中測量點(diǎn)。采用改進(jìn)積分方程方法對該路不同高度（地面、4.26 km, 6.40 km）衰減因子相位進(jìn)行預(yù)測，并將預(yù)測結(jié)果與上述8個(gè)點(diǎn)的實(shí)測結(jié)果進(jìn)行比較，如圖 5（c）所示。該路徑中平原地區(qū)電參數(shù)取σ=3 × 1 0-3s/m, εr= 1 3； 山 區(qū)電參 數(shù)取 σ = 4 × 10-4s/m, εr= 7 。

圖3 均勻光滑球地面條件下地波衰減因子算法一致性驗(yàn)證

圖4 分段均勻光滑球地面條件下地波衰減因子算法一致性驗(yàn)證

圖5 復(fù)雜不規(guī)則路徑積分方程方法預(yù)測結(jié)果與實(shí)測結(jié)果比較

從圖5看出對于該復(fù)雜路徑，算法結(jié)果與實(shí)測結(jié)果吻合很好。

4 結(jié)束語

本文主要研究了低頻地波理論衰減因子算法中預(yù)測結(jié)果不一致的問題。給出了地波衰減因子預(yù)測中積分方程方法的改進(jìn)公式，使不同算法在均勻/分段均勻光滑路徑模型下的地波衰減因子預(yù)測結(jié)果趨于一致，提高了電波地面長距離傳播衰減因子的精度，同時(shí)該算法更加適合實(shí)際工程問題的求解。受限于積分方程方法的局限性，諸如模型推導(dǎo)近似、數(shù)值積分近似及奇點(diǎn)處理等，文中改進(jìn)算法對于地形起伏較緩的不規(guī)則路徑預(yù)測精度較高。對于地形起伏劇烈的不規(guī)則路徑，積分方程方法可通過與電磁場數(shù)值計(jì)算方法相結(jié)合等途徑進(jìn)一步提高其預(yù)測精度[20]。

[1] Feng Y and Astin I. Remote sensing of soil moisture using the propagation of Loran-C navigation signals[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12（1）:195-198.

[2] Lei W and Chen B. Low complexity sparse noise reduction method for Loran-C skywave delay estimation[J]. Electronics Letters, 2013, 49（24）: 1572-1574.

[3] 李實(shí)鋒, 王玉林, 華宇, 等. 羅蘭-C信號抗干擾的快速檢測方法[J]. 西安交通大學(xué)學(xué)報(bào), 2013, 47（10）: 91-96.Li Shi-feng, Wang Yu-lin, Hua Yu, et al.. A fast anti-interference detection method for Loran-C signal[J].Journal of Xi’an Jiaotong University, 2013, 47（10）: 91-96.

[4] Zhou Li-li, Xi Xiao-li, Zhang Jin-sheng, et al.. A new method for Loran-C ASF calculation over irregular terrain[J]. IEEE Transactions on Aerospace and Electronic Systems, 2013,49（3）: 1738-1744.

[5] ?afa? J, Williams P, and Vejra?ka F. Accuracy performance of e-Loran receivers under cross-rate interference conditions[J]. Annual of Navigation, 2013, 19（1）: 133-148.

[6] 李實(shí)鋒, 王玉林, 華宇, 等. 羅蘭-C信號快速捕獲方法及其性能分析[J]. 電子與信息學(xué)報(bào), 2013, 35（9）: 2175-2179.Li Shi-feng, Wang Yu-lin, Hua Yu, et al.. Loran-C signal fast acquisition method and its performance analysis[J]. Journal of Electronics & Information Technology, 2013, 35（9）:2175-2179.

[7] Seo J and Kim M. eLoran in Korea - current status and future plans[C]. Proceedings of the European Navigation Conference, Vienna, 2013: 23-27.

[8] 潘威炎. 長波超長波極長波傳播[M]. 成都: 電子科技大學(xué)出版社, 2004: 106-173.

[9] Gesny H and Ravard O. Propagation over irregular terrain in the VHF band a review of integral equation models[C].Proceedings of the National Conference on Antennas and Propagation, York, UK, 1999: 61-64.

[10] Hufford G A. An integral equation approach to the problem of wave propagation over an irregular surface[J]. Quarterly Journal of Applied Mathematics, 1952, 9（4）: 391-404.

[11] Ott R H. An alternative integral equation for propagation over irregular terrain, 2[J]. Radio Science, 1971, 6（4）:429-435.

[12] Wu Z, Maclean T, Bagwell D, et al.. Propagation over an inhomogeneous irregular surface[J]. Radio Science, 1988,23（1）: 33-40.

[13] 梁仲寰, 李濱. 山區(qū)地形對100 kHz地波相位影響的實(shí)驗(yàn)研究[J]. 陜西天文臺臺刊, 1996, 19（1）: 88-97.Liang Zhong-huan and Li Bin. The experiments on the effects of mountain terrain upon 100 kHz ground wave phase[J].CSAO Publications, 1996, 19（1）: 88-97.

[14] 王東文. 復(fù)雜地面條件下低頻地波傳播的三維預(yù)測[J]. 電波科學(xué)學(xué)報(bào), 1988, 3（1）: 50-63.Wang Dong-wen. Three-dimensional prediction for the propagation of low-frequency groundwave over inhomogeneous irregular terrain[J]. Chinese Journal of Radio Science, 1988, 3（1）: 50-53.

[15] 溫芳茹. 地形對地波衰減因子的影響[J]. 電波科學(xué)學(xué)報(bào), 1987,2（1）: 27-38.Wen Fang-ru. The effect of the terrain on the groundwave attenuation factor[J]. Chinese Journal of Radio Science, 1987,2（1）: 27-38.

[16] 彭懷云, 樊文生, 潘威炎, 等. 沿不規(guī)則不均勻地面?zhèn)鞑サ膫?cè)面波[J]. 電波科學(xué)學(xué)報(bào), 2006, 21（4）: 497-502.Peng Huai-yun, Fan Wen-sheng, Pan Wei-yan, et al.. Lateral wave propagation along irregular and inhomogeneous ground[J]. Chinese Journal of Radio Science, 2006, 21（4）: 497-502.

[17] 陳秀明. 羅蘭C組合導(dǎo)航中地波ASF修正研究[D]. [碩士論文],中國科學(xué)院, 2007.Chen Xiu-ming. Study on ASF corrections in integrated navigation systems of Loran-C[D]. [Master dissertation],Chinese Academy of Sciences, 2007.

[18] 熊皓. 無線電波傳播[M]. 北京: 電子工業(yè)出版社, 2000:331-363.

[19] 潘威炎, 彭懷云, 張紅旗. 非均勻光滑球面地波衰減因子的拋物方程算法[J]. 電波科學(xué)學(xué)報(bào), 2006, 21（1）: 37-42.Pan Wei-yan, Peng Huai-yun, and Zhang Hong-qi. Parabolic equation algorithm of wave attenuation along inhomogeneous smooth ground[J]. Chinese Journal of Radio Science, 2006,21（1）: 37-42.

[20] Xi X, Zhou L, Zhang J, et al.. Combined IE-FDTD algorithm for long-range Loran-C ground-wave propagation[J]. IEEE Transactions on Antennas and Propagation, 2012, 60（8）:3802-3808.