高偉，魯國瑞

( 天津城建大學(xué)地質(zhì)與測繪學(xué)院, 天津 300384 )

0 引言

地球物理場和慣性導(dǎo)航系統(tǒng)(inertial navigation system，INS)聯(lián)合的輔助導(dǎo)航技術(shù)可有效遏制INS誤差累積問題，始終是國內(nèi)外水下導(dǎo)航領(lǐng)域的研究難點和熱點. 水下輔助導(dǎo)航技術(shù)主要包括地磁匹配導(dǎo)航、地形匹配導(dǎo)航和重力匹配導(dǎo)航等. 而水下重力匹配導(dǎo)航通過測量重力信息與地球重力場進(jìn)行匹配實現(xiàn)導(dǎo)航定位，在重力輔助導(dǎo)航系統(tǒng)運行時，采集重力特征顯著區(qū)域的重力信息與重力基準(zhǔn)圖對比得到位置信息，實現(xiàn)天空海一體化水下潛器慣性/重力組合導(dǎo)航系統(tǒng)重調(diào)[1-4].

地球形狀的不規(guī)則性和密度的不均勻性導(dǎo)致地球各點的重力場并不相同，表現(xiàn)為空間位置(經(jīng)度、緯度、高度)的函數(shù)[5]. 因此潛航器可在航行過程中，通過重力測量儀器采集航線途經(jīng)的重力數(shù)據(jù)，并與預(yù)先存儲的重力數(shù)據(jù)相匹配，獲得潛航器當(dāng)前定位信息，進(jìn)而對INS 累積的位置誤差進(jìn)行校正[6-7]. 重力輔助導(dǎo)航系統(tǒng)在測量重力場數(shù)據(jù)過程中，潛航器無需露出或接近水面，測量儀器也無需向外部發(fā)射信號或接收外部信號，系統(tǒng)可進(jìn)行無源、隱蔽的導(dǎo)航定位，潛航器在衛(wèi)星、無線電定位系統(tǒng)失效或遭受破壞的特殊情況下仍可達(dá)到自主隱蔽導(dǎo)航的目的[8-10].

隨著重力測量儀器及空間測量技術(shù)的發(fā)展[11]，在全球范圍內(nèi)快速、準(zhǔn)確地獲取重力數(shù)據(jù)成為現(xiàn)實，使得重力輔助導(dǎo)航系統(tǒng)具備校正INS 積累位置誤差的能力[12-13]. 重力輔助導(dǎo)航系統(tǒng)主要有兩種代表性系統(tǒng)，均由Bell Aerospace 研發(fā)，分別為于1990 年研制出的重力梯度導(dǎo)航系統(tǒng)[14]，以及于1991 年研制出的重力輔助INS[15]. 重力輔助導(dǎo)航系統(tǒng)由INS、重力測量儀器、數(shù)字重力基準(zhǔn)圖[16]和匹配定位算法等四個主要部分組成，各組成部分的性能差異對重力輔助定位系統(tǒng)的性能有重要影響. 數(shù)字重力基準(zhǔn)圖(重力模型)是重力輔助定位的基礎(chǔ)[17-18]，重力基準(zhǔn)圖描述是否準(zhǔn)確，包含的重力特征是否豐富，分辨率是否滿足需求，都會影響重力輔助導(dǎo)航系統(tǒng)的性能[19].

為減小線性化誤差，對重力異常參考地圖進(jìn)行預(yù)處理，選擇合適的計算模型對地圖進(jìn)行插值逼近，得到高分辨率、小格網(wǎng)地圖，以減小線性化誤差對位置觀測精度的影響. 在重力異常數(shù)據(jù)匱乏的地區(qū)，隨著相鄰點間的插值數(shù)的增加，擬合函數(shù)和Kriging 方法無法滿足逼近精度要求. 其中，擬合函數(shù)法跟插值的規(guī)模、控制點數(shù)目、控制點統(tǒng)計特征和地圖統(tǒng)計特征差異、地圖統(tǒng)計特征相關(guān)，不適宜重力異常數(shù)據(jù)匱乏的區(qū)域，而Kriging 方法不適合大尺度條件下的重力插值，因此采用分形重構(gòu)模擬真實的重力環(huán)境，為系統(tǒng)提供良好的仿真環(huán)境. 分形通常被用在模擬真實地表和其他自然形態(tài)的描述中，它用分維數(shù)表達(dá)數(shù)據(jù)場的統(tǒng)計特征. 與Kriging 方法的空間變異函數(shù)相似，低分辨率地圖不利于長距離的統(tǒng)計特征描述[20-22]. 本文考慮在這種條件下，提出利用分形插值的方法對重力異常隨機數(shù)據(jù)場進(jìn)行重構(gòu)，重構(gòu)模擬將為水下重力輔助INS 提供更接近真實的重力環(huán)境[23].

分形就是對不規(guī)則和支離破碎的自然形狀的描述，分形理論的發(fā)展促進(jìn)了利用它對各種自然形態(tài)描述的研究. Goodchild 在1980 年考慮用分形理論分析地文表面，后來Clark 和Schweitzer 定義了具有更好魯棒性的對地形表面分形估計. 從1991 年開始，Kenichi Arakawa 和Eric Krotkov 對“利用分形重構(gòu)自然地形”進(jìn)行了研究，利用分形布朗運動(fractal Brownian motion, FBM)在不規(guī)則采樣點集上估計自然地形的分形結(jié)構(gòu)，并以此為基礎(chǔ)構(gòu)造出自然地形，Naokazu 等給出基于FBM 模型的插值方法并不斷進(jìn)行完善[24].

本文提出采用FBM 分形插值方法對重力異常隨機場進(jìn)行重構(gòu)，構(gòu)造可用于水下重力輔助導(dǎo)航系統(tǒng)的重力異?；鶞?zhǔn)圖，通過實驗驗證分形插值方法的可行性，將為水下重力輔助導(dǎo)航系統(tǒng)的實現(xiàn)提供理論參考.

1 FBM 函數(shù)及分維數(shù)估計

布朗運動(Brownian Motion，BM)是1982 年英國植物學(xué)家R. Brown 發(fā)現(xiàn)的，1923 年德國數(shù)學(xué)家N.Wiener 建立了BM 的數(shù)學(xué)模型. BM 是一種隨機過程B(t)，滿足以下兩個條件：

式中：C為常數(shù)；隨機過程的任意兩點之間(B(b)-B(a)和B(c)-B(b),a≤b≤c) 的增量是相互獨立的.

B(t)被寫成

顯然，它具有自緊密性.

Mandelbrot 提出FBM 函數(shù)，由它產(chǎn)生BM 過程.

Penland 定義了FBM 函數(shù)統(tǒng)計自相似性特征

式中：f(t) 為FBM 函數(shù)；g(x) 為累計分布函數(shù)(g(x)～N(0,σ2))；H為自仿射參數(shù)，它與f(t) 的分維數(shù)D之間的關(guān)系為

n是幾何拓樸維數(shù)，在重力異常地圖中，n=2 .

由式(6)可得

又

可以得到

兩邊取對數(shù)

在對二維隨機數(shù)據(jù)場進(jìn)行分維數(shù)估計時， Δt用地圖上兩點距離代替，則

2 重力異常地圖分維數(shù)估計

重力異常地圖的分維數(shù)的估計也是以式(11)為基礎(chǔ)進(jìn)行的.

在計算時需要得到在不同 ΔX下的E[|f(X+ΔX)-f(X)|]，顯然在重力異常地圖陣列的基礎(chǔ)上，可以得到的 ΔX是離散的，它的取值與地圖網(wǎng)格單位距離和地圖陣列的規(guī)模有關(guān)，如圖1 所示.

圖1 地圖網(wǎng)格示意圖

能取的距離為ΔX

式中，i=1,2,···,mx；j=1,2,···,my；這里僅討論ux=uy=u的情形.

顯然，ΔX增加，在計算時，可參與統(tǒng)計點數(shù)目會急劇減少. 基于該思路，設(shè)定距離容差，可以更好地反映地圖的統(tǒng)計特征.

即，設(shè) εx為距離容差，當(dāng)?shù)貓D上兩點距離滿足

那么，兩點重力異常增量可以計入與ΔX對應(yīng)的計算中.其中，ΔXP是兩點之間的實際網(wǎng)格距離.

按此方法分別對2 個地塊區(qū)域的重力異常地圖進(jìn)行了分維數(shù)計算，結(jié)果如圖2 所示.

圖2 重力異常地圖分形特征曲線

用最小二乘法擬合直線，得到的結(jié)果如表1 所示.

表1 分形特征計算結(jié)果

3 FBM 分形插值方法

用FBM 分形插值，生成新的小尺度重力異常地圖的基本步驟如圖3 所示.

圖3 分形插值方法基本流程

分形特征一個有限的尺度范圍內(nèi)才能充分體現(xiàn)，這個范圍稱為無標(biāo)度區(qū)間 [ΔXmin,ΔXmax] .

由地圖得到的分形特征無標(biāo)度區(qū)間的下限由地圖間隔決定

其上限為分形特征曲線線性擬合區(qū)域的上限.

分形插值后得到地圖的最小尺度uc小于u，因此分形插值所能得到的最小尺度不能由原始地圖確定.本文所討論的重力異常地圖分形插值只能用于重力異常地圖的粗插值.

分形插值采用隨機中點位移法，即插入點附近區(qū)域的點理解為一個符合 (0,σ2) 的隨機數(shù)據(jù)與區(qū)域的平均值的和，它的實現(xiàn)如圖4 所示.

圖4 中點位移插值法

在圖4 中，填充點為原始地圖的控制點，虛線是插值點；插值點分為兩類，標(biāo)號為1 的插值點 (x,y) ，由其周圍的4 個原地圖控制點為參照進(jìn)行插值

標(biāo)號為2 的插值點，所用到的參照點由兩個原地圖控制點和標(biāo)號為1 的兩個插值點構(gòu)成

其中， Gauss 是符合均值為0，標(biāo)準(zhǔn)差為1 的隨機信號.

4 實驗結(jié)果與分析

4.1 與原始地圖比較

針對某塊 60×60 的重力異常地圖進(jìn)行插值，在原始地圖上進(jìn)行等間隔抽樣得到一塊 30×30 的地圖，然后用上述分形插值方法進(jìn)行插值，并對插值地圖和原始地圖進(jìn)行比較. 抽樣后的地圖分形特征如表2所示.

表2 分形特征參數(shù)表

插值后的重力異常地圖與原始地圖的比較如圖5所示.

圖5 在原始數(shù)據(jù)上抽樣后分形插值結(jié)果

插值后的精度分析如表3 所示，采用誤差平均值(E(ge))、絕對誤差平均值、絕對誤差最大值(max(|ge|))、誤差均方值和絕對誤差標(biāo)準(zhǔn)差(D(ge) )等參數(shù)來評定.

表3 分形插值精度分析

4.2 與原重力異常地圖比較

在原重力異常地圖上進(jìn)行插值，分析插值后地圖的分形特征的變化. 圖6 是插值后的等值線圖比較.圖7 是插值后地圖分形特征曲線. 其特征參數(shù)見表2第2 行數(shù)據(jù).

圖6 原始數(shù)據(jù)直接分形插值結(jié)果

圖7 插值后地圖分形特征曲線

顯然，插值后的地圖與插值前地圖的分形特征變化很小. 分形插值的精度與標(biāo)準(zhǔn)差相關(guān)，標(biāo)準(zhǔn)差越大，精度越低；分形插值引入隨機信號，因此插值點的重力異常值具有不確定性；采用中點位移法的插值，插值點位置固定，即uc=u/2 .

分形插值作為地形自然表面的重構(gòu)方法，具有很強的真實感. 重力異常與地形之間的關(guān)聯(lián)性，產(chǎn)生了利用分形對重力異常進(jìn)行插值的聯(lián)想，為得到可靠的重力異常地圖的分維數(shù)，在分維數(shù)定義中加入距離容差，改善了分維數(shù)在長空間距離下的統(tǒng)計特征描述.分形插值的精度雖然低，但插值后的地圖仍保持原重力異常參考地圖的統(tǒng)計特征，該方法為仿真提供了很好的環(huán)境模擬. 由于條件的限制，雖然沒有對插值結(jié)果進(jìn)行實際測量的驗證，但從計算機仿真實驗結(jié)果看，分形插值作為一種粗插值是可行的.

本文是針對水下重力數(shù)據(jù)匱乏區(qū)域的重力場仿真研究，由于數(shù)據(jù)異常和條件差，分形插值的實驗結(jié)果僅與原始數(shù)據(jù)的計算結(jié)果進(jìn)行了對比研究. 由于重力數(shù)據(jù)的匱乏，采用其他插值方法的效果不太理想，限于篇幅和數(shù)據(jù)條件，未進(jìn)行詳細(xì)分析，但后續(xù)仍需要通過獲取大量實測數(shù)據(jù)來進(jìn)行反復(fù)驗證與探討.

5 結(jié)束語

本文針對重力異常數(shù)據(jù)匱乏的地區(qū)，提出了采用FBM 分形插值方法對重力異常隨機場進(jìn)行重構(gòu)，實現(xiàn)可用于水下重力輔助導(dǎo)航的重力異常基準(zhǔn)圖生成. 實驗結(jié)果表明，分形插值后的地圖與插值前地圖的分形特征變化很小. 分形插值的精度與標(biāo)準(zhǔn)差相關(guān)，標(biāo)準(zhǔn)差越大，精度越低；分形插值引入隨機信號，因此插值點的重力異常值具有不確定性；采用中點位移法的插值，插值點位置可固定. 分形插值作為地形自然表面的重構(gòu)方法，具有很強的真實感. 從數(shù)值模擬結(jié)果來看，分形插值方法作為一種粗插值是可行的和有效的，可為水下重力輔助INS 提供更接近真實的重力環(huán)境，為水下重力輔助導(dǎo)航系統(tǒng)的實現(xiàn)提供理參考.