低頻地波傳播時(shí)延計(jì)算的地形優(yōu)化處理方法

朱峙亞，郭偉

( 1.中國(guó)科學(xué)院大學(xué)，北京 100049；2.中國(guó)科學(xué)院國(guó)家授時(shí)中心，西安 710600)

0 引言

近年來(lái)國(guó)際上越來(lái)越重視羅蘭C(Loran-C)陸基遠(yuǎn)程長(zhǎng)波定位導(dǎo)航和授時(shí)系統(tǒng)(PNT)，它與衛(wèi)星導(dǎo)航系統(tǒng)獨(dú)立，同時(shí)也是衛(wèi)星導(dǎo)航系統(tǒng)的備份.Loran-C系統(tǒng)發(fā)射峰值功率達(dá)到2 MW，信號(hào)傳播距離達(dá)到幾千公里，具有很強(qiáng)的抗干擾能力和很廣泛的覆蓋范圍.Loran-C信號(hào)是波長(zhǎng)為3 km 的低頻長(zhǎng)波電磁波，因其沿地球表面?zhèn)鞑ビ址Q為低頻地波.

附加二次相位因子(ASPF)定義為低頻地波傳播路徑，是陸地時(shí)對(duì)傳播時(shí)延產(chǎn)生的影響，ASPF已經(jīng)成為影響Loran-C系統(tǒng)PNT精度的主要因素.計(jì)算ASPF要通過(guò)求解麥克斯韋方程組在滿足邊界條件的解，這是非常復(fù)雜的計(jì)算過(guò)程.對(duì)低頻地波的傳播研究方興未艾，最早由熊皓[1]研究得出均勻光滑平地面模型下地波場(chǎng)的積分表達(dá)式；Fock[2]得出均勻光滑球面模型下地波場(chǎng)的級(jí)數(shù)表達(dá)式；Wait[3]、Millington[4]等得出分段均勻光滑地面模型下的計(jì)算公式，Millington 公式在工程中被廣泛應(yīng)用.《中華人民共和國(guó)電子行業(yè)軍用標(biāo)準(zhǔn)》中長(zhǎng)波地波傳輸信道計(jì)算方法采用的就是以上三種模型與計(jì)算方法.

低頻地波的傳播路徑通常包含不同的地形地質(zhì)，如平原、山脈、沼澤和沙漠等，這些不同的地形地質(zhì)其地面電參數(shù)也不同.如大地電導(dǎo)率、相對(duì)介電常數(shù)，地表大氣折射率等，對(duì)信號(hào)接收點(diǎn)的場(chǎng)強(qiáng)與時(shí)延的影響也不盡相同，Millington 方法只能將傳播路徑粗略分為均勻光滑若干段，無(wú)法考慮地形因素對(duì)低頻地波的影響.Hufford[5]利用第二格林定理得出了低頻地波在不均勻不光滑傳播模型中地波場(chǎng)計(jì)算的積分方程方法，經(jīng)過(guò)實(shí)際測(cè)量證明其計(jì)算精度較高，目前被廣泛應(yīng)用.周麗麗等[6]也對(duì)其進(jìn)行了理論研究，并且計(jì)算分析了不同形狀山脈對(duì)低頻地波傳播性能的影響.

積分方程方法具有很高的計(jì)算精度，但是積分方程方法在推導(dǎo)過(guò)程中使用了穩(wěn)定相位積分近似，導(dǎo)致其無(wú)法滿足極端苛刻的地形.只適用于地形緩變情況，即傳播路徑高程數(shù)據(jù)不含有“噪聲”，是一條光滑的曲線.因此計(jì)算過(guò)程中首先需要獲取高分辨率的傳播路徑高程數(shù)據(jù)，一方面可以避免因?yàn)榉e分步長(zhǎng)遠(yuǎn)大于地形采樣間隔導(dǎo)致地形數(shù)據(jù)產(chǎn)生突變；另一方面可以避免計(jì)算高程數(shù)據(jù)的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)誤差太大.

地形高程數(shù)據(jù)通過(guò)編寫程序處理SRTM3高程數(shù)據(jù)文件，根據(jù)收發(fā)點(diǎn)的經(jīng)緯度和采樣點(diǎn)數(shù)來(lái)獲取.因其固有的離散采樣特性與地形本身陡峭等因素，原始路徑高程數(shù)據(jù)不滿足積分方程方法緩變地形的要求.因此我們需要對(duì)低頻地波傳播路徑采取一定的“平滑”處理，保留地形的“包絡(luò)”，即保留地形的幾何輪廓，去掉“噪聲”.本文創(chuàng)新性的使用了數(shù)學(xué)形態(tài)法對(duì)原始高程數(shù)據(jù)經(jīng)行處理，使傳播路徑高程數(shù)據(jù)滿足積分方程方法的要求.

1 數(shù)學(xué)形態(tài)學(xué)

數(shù)學(xué)形態(tài)學(xué)是于1964年由法國(guó)數(shù)學(xué)家Matheron[7]和Serra[8]共同提出的基于集合論、積分幾何、拓?fù)鋵W(xué)等數(shù)學(xué)基礎(chǔ)的交叉性學(xué)科.近年來(lái)通過(guò)數(shù)學(xué)家不斷的研究，數(shù)學(xué)形態(tài)學(xué)的理論基礎(chǔ)已經(jīng)逐步完善，相關(guān)應(yīng)用也深入很多領(lǐng)域.數(shù)學(xué)形態(tài)學(xué)的基礎(chǔ)算子只有簡(jiǎn)單的加減法和取極值運(yùn)算，沒(méi)有乘除法運(yùn)算，這樣大大提高了運(yùn)算效率，運(yùn)算速度很快.

最初數(shù)學(xué)形態(tài)學(xué)是應(yīng)用于圖像處理，它是利用結(jié)構(gòu)元素對(duì)圖像進(jìn)行變換，應(yīng)用于去除圖像噪聲和輪廓邊緣檢測(cè)等方面.

數(shù)學(xué)形態(tài)學(xué)在信號(hào)處理方面提供了一種可以檢測(cè)信號(hào)幾何特征的非線性信號(hào)處理方式，被處理信號(hào)的幾何形狀信息可以由作用在信號(hào)上的結(jié)構(gòu)元素提取出來(lái)，這時(shí)的結(jié)構(gòu)元素相當(dāng)于一個(gè)“探針”，在一維或二維信號(hào)中不斷移動(dòng)這個(gè)探針，就可考察信號(hào)波形中的幾何關(guān)系.數(shù)學(xué)形態(tài)學(xué)已在圖像處理、計(jì)算機(jī)視覺(jué)等領(lǐng)域得到廣泛應(yīng)用，目前電力系統(tǒng)、信號(hào)處理等領(lǐng)域中也逐漸展開研究與應(yīng)用.

數(shù)學(xué)形態(tài)學(xué)有兩個(gè)基本運(yùn)算：腐蝕和膨脹.腐蝕可以使目標(biāo)區(qū)域變小，造成目標(biāo)邊界收縮，用來(lái)消除小且無(wú)意義的目標(biāo)；膨脹會(huì)使目標(biāo)區(qū)域變大，造成目標(biāo)邊界擴(kuò)大，減小目標(biāo)區(qū)域內(nèi)的谷域以及消除包含在目標(biāo)區(qū)域中的噪聲.將腐蝕運(yùn)算看作是最小值濾波器，膨脹運(yùn)算看作是最大值濾波器，它們可分別獲得數(shù)據(jù)的下包絡(luò)和上包絡(luò).數(shù)學(xué)形態(tài)學(xué)基本算子主要包括膨脹、腐蝕及以此為基礎(chǔ)構(gòu)造的開運(yùn)算、閉運(yùn)算4種運(yùn)算方法.定義原始信號(hào)z(n) 為 在Z=(0,1,···,N?1) 上的離散函數(shù)，g(n)為在G=(0,1,···,M?1)上的離散函數(shù). 腐蝕和膨脹運(yùn)算定義為：

2 長(zhǎng)波信號(hào)傳播時(shí)延計(jì)算

由于低頻地波是沿地面?zhèn)鞑サ?，因此接收點(diǎn)的相位延遲不僅與發(fā)射站的頻率有關(guān)，還與發(fā)射站與接收站之間的傳播距離有關(guān).其中還有一個(gè)不可忽略的影響因素即地面.地面對(duì)低頻地波傳播的影響主要有：一是地面的粗糙度，如山脈等地形起伏的影響；另一個(gè)是地質(zhì)的電磁特性，在地質(zhì)結(jié)構(gòu)中由于傳播介質(zhì)時(shí)空分布不均勻，如低頻地波的傳播路徑在若干年之后由于地表水系或者水土流失等發(fā)生改變.由于地形和大地電導(dǎo)率等因素對(duì)低頻地波傳播路徑的影響，實(shí)際無(wú)線電波信號(hào)的傳播過(guò)程中會(huì)發(fā)生一系列復(fù)雜的變化.例如：反射、衍射、散射和折射，導(dǎo)致傳播速度和相位改變，最終導(dǎo)致導(dǎo)航和授時(shí)產(chǎn)生誤差，這種誤差在授時(shí)系統(tǒng)可以達(dá)到μm 級(jí)別，定位系統(tǒng)的誤差可能達(dá)到數(shù)百米到幾千米級(jí)別.

式中：PF為發(fā)射站與接收站之間的直線距離和大氣折射率對(duì)電波傳播造成的時(shí)延；PF又稱為一次相位因子；SF為發(fā)射站與接收站之間為純海水的情況下對(duì)電波造成的時(shí)延，SF又稱為二次相位因子；ASPF為發(fā)射站與接收站之間為純陸地對(duì)電波傳播造成的時(shí)延減去若該路徑為純海水對(duì)電波傳播造成的時(shí)延，即傳播路徑為陸地相對(duì)于純海水的時(shí)延；ASPF又稱為附加二次相位因子.PF、SF、ASPF的單位為μs；d是以km 為單位的電波傳播的距離；arg 為求衰減因子幅角；ns為地表附件大氣折射率，約為1.000 338；c為光速；ω =2πf，f為電波頻率；W為地波衰減函數(shù)，是計(jì)算接收站地波場(chǎng)強(qiáng)和時(shí)延的關(guān)鍵因素，與電波頻率，傳播距離，大地電導(dǎo)率，傳播路徑地形等有關(guān)，目前低頻地波傳播信道模型及W的計(jì)算主要有以下四種方式：

1)均勻光滑平地面模型

均勻地面指的是電波傳播路徑上地質(zhì)結(jié)構(gòu)趨于一致，地面電性參數(shù)變化不大，比如大地電導(dǎo)率等可以取一個(gè)定值.光滑地面指的是地形起伏變化程度遠(yuǎn)小于電波波長(zhǎng)，如100 kHz 的低頻地波波長(zhǎng)為3 km，關(guān)中平原就可以看作是光滑地面.當(dāng)電波傳播路徑小于60～70 km 時(shí)，可以把傳播路徑當(dāng)做平地面.繼法拉第電磁感應(yīng)定理的發(fā)現(xiàn)到麥克斯韋方程組的建立，垂直電偶極子的輻射場(chǎng)求解問(wèn)題頗有歷史淵源[9]，Sommerfeld[10]通過(guò)研究平地面上垂直電偶極子產(chǎn)生的場(chǎng)得出了一個(gè)積分表達(dá)式，經(jīng)研究人員不斷改進(jìn)，使其可以在工程中應(yīng)用.

2)均勻光滑球面模型

低頻地波繞射能力很強(qiáng)，可以沿地面?zhèn)鞑ズ苓h(yuǎn)的距離，當(dāng)?shù)皖l地波的傳播距離很長(zhǎng)時(shí)就必須要考慮地球曲率對(duì)地波傳播的影響.均勻光滑球地面模型將地球理想化為一個(gè)表面光滑，地面電性參數(shù)為常數(shù)的球體.地面的發(fā)射天線比波長(zhǎng)小得多，可以理想化為垂直電偶極子.Watson[11]、Fock[2]等科研人員在上世紀(jì)初做了大量研究，得出了垂直電偶極子在球面上的輻射場(chǎng)的積分表達(dá)式，并且經(jīng)過(guò)復(fù)雜數(shù)學(xué)運(yùn)算后將其變成一個(gè)收斂相當(dāng)快的級(jí)數(shù)表達(dá)式[3]，文中取為e?iωt，均勻光滑球面模型的地波衰減函數(shù)計(jì)算公式為

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

通常情況下地面既不均勻又不平坦，因此從發(fā)射站到接收站的整個(gè)傳播路徑不能看成是均勻光滑的傳播路徑.比如海洋與陸地的電導(dǎo)率差異很大，沙漠、丘陵、平原或沼澤等陸地的大地電導(dǎo)率也不盡相同，這樣就不能把整個(gè)傳播路徑看成是均勻的，只能把陸地部分和海洋部分看成兩段均勻路徑.陸地部分再根據(jù)地形地貌細(xì)分為若干段，每一段當(dāng)做均勻光滑路徑，大地電參數(shù)在每一段取不同的常數(shù).分段均勻光滑地面的地波衰減函數(shù)計(jì)算公式有Wait 公式[3]、波模轉(zhuǎn)換法、拋物方程法和Millington 公式[4]等.Millington 公式簡(jiǎn)單實(shí)用，在工程中廣泛應(yīng)用，傳播信道模型如圖1所示.

圖1 分段均勻光滑地面模型圖

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

實(shí)際電波傳播路徑不僅地質(zhì)類型不同，同時(shí)也存在山脈等地形起伏地區(qū)，即不均勻不光滑地面.積分方程方法適用于均勻光滑平地面模型、均勻光滑球面模型、分段均勻地面模型，與它們具有同樣的計(jì)算精度，還可以計(jì)算出地波傳播路徑上地形起伏變化產(chǎn)生的影響，是應(yīng)用非常廣泛的一種算法，國(guó)內(nèi)外對(duì)復(fù)雜路徑低頻地波傳播預(yù)測(cè)多用此方法.傳播信道模型如圖2所示.

圖2 不光滑不均勻地面模型圖

假設(shè)地面起伏情況如圖2中曲線所示，其滿足緩變條件.地波衰減因子計(jì)算公式為

信號(hào)發(fā)射源位于原點(diǎn)O，ldl 為垂直電偶極子，接收點(diǎn)為P，地面上的動(dòng)點(diǎn)為Q，D為從源點(diǎn)到接收點(diǎn)的大圓距離，r0為從源點(diǎn)到接收點(diǎn)P的直線距離，S0為積分區(qū)域，Δg為地球歸一化表面阻擾，r1為從源點(diǎn)到積分動(dòng)點(diǎn)Q的直線距離，r2為從積分動(dòng)點(diǎn)Q到接收點(diǎn)P的直線距離.

3 計(jì)算結(jié)果分析

3.1 傳播路徑地形高程信息提取

航天飛機(jī)雷達(dá)地形測(cè)繪使命(SRTM)，是2000年美國(guó)國(guó)家地理空間情報(bào)局(NGA)與美國(guó)國(guó)家航空航天局(NASA)以及德國(guó)和意大利的航天機(jī)構(gòu)的一項(xiàng)聯(lián)合項(xiàng)目，該項(xiàng)目使用航天飛機(jī)進(jìn)行了為期11天的測(cè)量，目的是為地球80%的陸地表面(60°N～56°S的所有陸地)生成數(shù)字地形高程數(shù)據(jù)，在90%置信度下，高程數(shù)據(jù)的絕對(duì)垂直精度為16 m，是目前使用最廣泛的高程數(shù)據(jù)源之一.SRTM 數(shù)據(jù)按照分辨率分為SRTM3和SRTM1兩種，SRTM3分辨率為3′′，即90 m，SRTM1分辨率為1′′，即30 m，但是SRTM1數(shù)據(jù)只包含美國(guó)全境，本文所用數(shù)據(jù)為SRTM3.

3.1.1 SRTM3數(shù)據(jù)結(jié)構(gòu)及讀取

SRTM3數(shù)據(jù)包含很多文件，每一個(gè)文件覆蓋地球表面一個(gè)緯度乘一個(gè)經(jīng)度塊，每個(gè)文件里面包含1 201×1 201個(gè)采樣點(diǎn)的高度數(shù)據(jù).經(jīng)緯度網(wǎng)格與SRTM3文件對(duì)應(yīng)關(guān)系如圖3所示.

圖3 經(jīng)緯度網(wǎng)格與SRTM 3文件對(duì)應(yīng)圖

由圖3可知，每個(gè)經(jīng)緯度網(wǎng)格與一個(gè)SRTM3文件對(duì)應(yīng)，SRTM3文件名稱包含經(jīng)緯度信息，根據(jù)每個(gè)經(jīng)緯度網(wǎng)格左下角的經(jīng)緯度值構(gòu)造SRTM3文件名.比如圖3中包含圖釘?shù)木W(wǎng)格，其左下角坐標(biāo)為(33°N，107°E)，則這個(gè)網(wǎng)格對(duì)應(yīng)的SRTM3文件名稱為“N34E107.hgt”.此文件可以看作為1 201×1 201的矩陣，從左往右的方向?yàn)榻?jīng)度增加方向，從下到上的方向?yàn)榫暥仍黾拥姆较?，相鄰文件最外?cè)的行與列重疊.計(jì)算機(jī)中的存儲(chǔ)單位是字節(jié)，SRTM3是一種hgt 格式的文件，hgt 來(lái)源于英文單詞“hight”，文件中高程數(shù)據(jù)使用16位有符號(hào)整數(shù)表示，能取到的最大值為32 767，單位是m.文件中是以大端模式、順序結(jié)構(gòu)存儲(chǔ)高程數(shù)據(jù)二維數(shù)組，如果采用普通的二進(jìn)制方式讀取文件，則讀取后處理比較麻煩，不過(guò)python 的numpy 庫(kù)只要設(shè)置好讀取模式就可以輕松讀取數(shù)據(jù)，并且將高程數(shù)據(jù)轉(zhuǎn)化為二維數(shù)組.

3.1.2傳播路徑高程提取程序

圖4 高程信息提取程序執(zhí)行框圖

3.2 數(shù)學(xué)形態(tài)學(xué)對(duì)高程信息處理

選擇圖3中陜西省渭南市蒲城發(fā)射臺(tái)與四川省什邡市接收點(diǎn)之間的傳播路徑，使用上節(jié)編寫的提取高程數(shù)據(jù)程序，輸入發(fā)射臺(tái)與接收點(diǎn)的經(jīng)緯度和采樣點(diǎn)數(shù)，得到路徑高程信息如圖5所示.由圖5可知，路徑總長(zhǎng)約650 km，從發(fā)射站出發(fā)，先經(jīng)過(guò)約150 km 關(guān)中平原，然后翻越秦嶺山脈，最后到達(dá)約450 km 的四川盆地.整體傳播路徑的地形變化為先平緩后起伏最后趨于平緩.關(guān)中平原段平緩，滿足積分方程計(jì)算要求，四川盆地段有少量“噪聲”，地形起伏變化主要集中在了秦嶺山脈，這兩段需要使用數(shù)學(xué)形態(tài)法處理，使其既保留山的輪廓又去除“噪聲.處理結(jié)果如圖6所示.

圖5 原始傳播路徑高程信息圖

圖6 數(shù)學(xué)形態(tài)法優(yōu)化后路徑高程信息圖

由圖6可知，經(jīng)過(guò)數(shù)學(xué)形態(tài)學(xué)處理后，傳播路徑在關(guān)中平原段沒(méi)有明顯變化，在四川盆地和秦嶺段變化較大，四川盆地和秦嶺山脈段變得光滑，秦嶺山脈輪廓清晰，高度幾乎沒(méi)有損失，證明利用數(shù)學(xué)形態(tài)法處理后的傳播路徑很好地滿足了積分方程方法對(duì)地形的要求.

3.3 積分方程計(jì)算結(jié)果

積分方程方法分別計(jì)算上述傳播路徑ASPF，圖7為原始傳播路徑ASPF計(jì)算結(jié)果，圖8為數(shù)學(xué)形態(tài)學(xué)處理過(guò)后的傳播路徑計(jì)算結(jié)果.

圖7 原始路徑計(jì)算結(jié)果圖

圖8 優(yōu)化后路徑計(jì)算結(jié)果圖

由圖7可知，原始路徑計(jì)算的ASPF在關(guān)中平原段比較光滑，四川盆地和秦嶺段脈沖噪聲較多且噪聲幅度較大.

由圖8可知，數(shù)學(xué)形態(tài)學(xué)處理過(guò)后的路徑計(jì)算的ASPF在四川盆地和秦嶺段脈沖噪聲明顯減少且幅度也顯著降低.

4 總結(jié)

積分方程方法在推導(dǎo)長(zhǎng)波信號(hào)傳播時(shí)延過(guò)程中采用了穩(wěn)定相位積分等近似條件，因此其計(jì)算過(guò)程要求地面是緩變的，地面曲率半徑不能太小.數(shù)學(xué)形態(tài)學(xué)處理過(guò)后的地形能夠保留其原有幾何輪廓去除“噪聲”，計(jì)算速度快，使地面滿足緩變的條件.積分方程方法在處理過(guò)后的傳播路徑上計(jì)算的ASPF擾動(dòng)與跳變減少，趨于平滑，表明數(shù)學(xué)形態(tài)學(xué)對(duì)于處理低頻地波傳播路徑數(shù)據(jù)是很有效的.