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

1 引言

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

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

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

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

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

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

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

其中

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

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

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

此處

h1和h2分別為發(fā)射點(diǎn)和接收點(diǎn)離地面高度，w2（ t）是愛(ài)里函數(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è)有限段，每段長(zhǎng)度分別為： d1,d2,…,dN，歸一化表面阻抗分別為Δg1,Δg2,…,ΔgN。對(duì)應(yīng)的地波衰減因子求解方法目前主要有 4種：Wait積分法、Millington方法、波模轉(zhuǎn)換法以及拋物線(xiàn)方程方法，下面僅給出拋物線(xiàn)方程方法相關(guān)公式。

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

其滿(mǎn)足邊界條件

其中

當(dāng)已知x處地面及空中的u（ x, z），式（10）可以通過(guò)分步混合傅里葉變換或有限差分法求解下一段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ù)值解來(lái)求得。記發(fā)射源點(diǎn)的位置為O，觀(guān)察點(diǎn)的位置為P，地面上積分動(dòng)點(diǎn)為Q，傳播路徑幾何示意圖如圖1所示。地波衰減因子滿(mǎn)足的積分方程為[8,18]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3給出了接收高度 hpR分別為0 km和10 km條件下，采用積分方程、留數(shù)級(jí)數(shù)，拋物線(xiàn)方程方法計(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ù)測(cè)結(jié)果一致性非常好，均勻光滑地面條件下在1000 km范圍內(nèi)與留數(shù)級(jí)數(shù)法幅度誤差小于0.2 dB，相位誤差小于50 ns。

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

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

圖 5（a）所示的是一條實(shí)際復(fù)雜傳播路徑上測(cè)量點(diǎn)分布區(qū)域圖。測(cè)量信號(hào)為宣城臺(tái)（O 點(diǎn)所示）播發(fā)的導(dǎo)航信號(hào)，沿該路徑共測(cè)量了8個(gè)點(diǎn)的相位，其中地面上測(cè)量了6個(gè)點(diǎn)，海拔為4.26 km的高空測(cè)量了1個(gè)點(diǎn)，海拔為6.40 km處測(cè)量了1個(gè)點(diǎn)。圖5（b）所示的是沿該路徑地形高程及 8個(gè)測(cè)量點(diǎn)空間位置分布，其中1～6為地面測(cè)量點(diǎn)，7, 8為兩個(gè)空中測(cè)量點(diǎn)。采用改進(jìn)積分方程方法對(duì)該路不同高度（地面、4.26 km, 6.40 km）衰減因子相位進(jìn)行預(yù)測(cè)，并將預(yù)測(cè)結(jié)果與上述8個(gè)點(diǎn)的實(shí)測(cè)結(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ù)測(cè)結(jié)果與實(shí)測(cè)結(jié)果比較

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

4 結(jié)束語(yǔ)

本文主要研究了低頻地波理論衰減因子算法中預(yù)測(cè)結(jié)果不一致的問(wèn)題。給出了地波衰減因子預(yù)測(cè)中積分方程方法的改進(jìn)公式，使不同算法在均勻/分段均勻光滑路徑模型下的地波衰減因子預(yù)測(cè)結(jié)果趨于一致，提高了電波地面長(zhǎng)距離傳播衰減因子的精度，同時(shí)該算法更加適合實(shí)際工程問(wèn)題的求解。受限于積分方程方法的局限性，諸如模型推導(dǎo)近似、數(shù)值積分近似及奇點(diǎn)處理等，文中改進(jìn)算法對(duì)于地形起伏較緩的不規(guī)則路徑預(yù)測(cè)精度較高。對(duì)于地形起伏劇烈的不規(guī)則路徑，積分方程方法可通過(guò)與電磁場(chǎng)數(shù)值計(jì)算方法相結(jié)合等途徑進(jìn)一步提高其預(yù)測(cè)精度[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信號(hào)抗干擾的快速檢測(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信號(hào)快速捕獲方法及其性能分析[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] 潘威炎. 長(zhǎng)波超長(zhǎng)波極長(zhǎng)波傳播[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ū)地形對(duì)100 kHz地波相位影響的實(shí)驗(yàn)研究[J]. 陜西天文臺(tái)臺(tái)刊, 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ù)測(cè)[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] 溫芳茹. 地形對(duì)地波衰減因子的影響[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]. [碩士論文],中國(guó)科學(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] 熊皓. 無(wú)線(xiàn)電波傳播[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.