朱峙亞，郭偉

( 1.中國科學院大學，北京 100049；2.中國科學院國家授時中心，西安 710600)

0 引言

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

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

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

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

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

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

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

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

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

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

2 長波信號傳播時延計算

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

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

1)均勻光滑平地面模型

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

2)均勻光滑球面模型

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

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

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

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

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

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

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

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

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

3 計算結(jié)果分析

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

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

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

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

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

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

3.1.2傳播路徑高程提取程序

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

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

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

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

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

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

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

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

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

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

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

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

4 總結(jié)

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