劉繁明 姚劍奇 荊 心

1 哈爾濱工程大學自動化學院，哈爾濱市南通大街145號，150001

隨著重力測量精度的不斷提高，出現了以地球重力場特征隨地理位置變化的特點進行定位的重力輔助導航。相對衛(wèi)星／無線電導航定位方法，重力輔助導航依靠重力測量的無源性，可保證潛器航行的隱蔽性與導航定位的可靠性［1］。重力輔助導航是基于重力基準圖的有圖匹配定位方法，基準圖的精確性將影響匹配定位精度。

由于經濟成本及其他客觀原因，目前難以獲得大范圍、高分辨率的實測重力數據。而低分辨率數據又難以滿足導航定位的要求，所以常用空間插值的方法加密實測重力數據［2－3］。因此，插值精度對重力基準圖的構建較為重要。

本文首先對比研究了幾種重力場插值效果較好的插值算法；而后提出以基準圖小波細節(jié)信息為轉移概率，由隨機采樣方案在基準圖高頻區(qū)域補充測點，從而提高插值算法精度的策略；最后進行了實驗驗證。

1 重力基準圖插值算法

經過大量的對比實驗得知，重力場插值效果較為優(yōu)秀的空間插值算法有Kriging算法、改進的Shepard算法、徑向基函數算法、格林函數樣條算法。

1.1 Kriging算法

Kriging算法［2］是以空間區(qū)域變量理論為基礎的無偏最優(yōu)估計算法。該算法將觀測量視為依賴空間具體位置的隨機場，其中隨位置變化的變量為區(qū)域化變量。并結合變異函數理論，由一定范圍內空間變量的相關性對未知位置變量進行估計。

1.2 改進的Shepard算法

改進的Shepard算法［3］通過局部區(qū)域最佳二次曲面擬合數據代替已知數據，由未知點一定鄰域內的擬合數據，按其與未知點的距離加權進行插值計算。具體的計算形式為：

式中，p0為未知點，pi為p0鄰域內的已知測點，di為pi與p0的 距 離，Qi為 擬 合 數 據，Rw為p0的鄰域半徑，其選取方法可參見文獻［4］。

1.3 徑向基函數算法

徑向基函數插值［5］為利用已知基準數據集｛pi，zi，i＝1，…，L｝，選取徑向基函數φ并通過平移構造基函數系｛φ（‖p－pi‖）｝進行插值的算法。具體計算形式為：

其中wi為系數，‖·‖為歐氏距離的模值。常用的基函數φ的類型有逆多重二次曲面函數、多重二次曲面函數、高斯函數及薄板樣條函數，基函數中含有平滑因子c。以多重二次曲面為例，其形式為：

基函數及平滑因子的選取對插值精度影響較大，該影響與目標數據類型、已知基準數據分布、計算精度密切相關。經實驗得知，高斯基函數的重力場插值精度較低，而文獻［6］提出的優(yōu)化方法可確定較優(yōu)的平滑因子取值。

1.4 格林函數樣條插值算法

格林函數樣條插值［7］由以已知基準點為中心位置的格林函數線性疊加，形成插值曲面?；鶞蕯祿蛱荻刃畔⒖刂聘骰鶞庶c處的格林函數系數。此外，引入張力以限制樣條函數在基準數據間的過度擺動。具體計算形式為：

式中，wi為系數，K0為0階的第二類修正Bessel函數，t為樣條中的張力。

1.5 重力基準圖插值實驗

以Topex衛(wèi)星重力數據（范圍2°×2°，分辨率1′×1′）作為原始基準圖進行插值實驗，如圖1所示。圖1（a）、（b）的重力變化范圍分別為－294～260mGal及－142～139mGal。

降采樣基準圖至3′×3′，由插值算法恢復至原分辨率并驗證插值精度。各算法誤差如表1所示（“AShepard”為改進Shepard算法，“RBF”為徑向基函數算法）。Kriging采用球狀變異函數模型，徑向基函數算法采用多重二次曲面基函數，平滑因子為1.8。

圖1 重力基準圖Fig.1 Gravity reference map

誤差分布如圖2所示（只顯示誤差的絕對值，以改善視覺效果）?？梢姴煌逯邓惴ň哂邢嗨频恼`差分布情況。經觀察發(fā)現，插值誤差較大的區(qū)域往往存在幅度較大的重力場高頻信息。

圖2 插值算法誤差分布Fig.2 The error distribution of interpolation algorithm

2 插值精度提升策略

2.1 插值誤差分析

插值算法實質是利用已知數據估計未知位置信息的過程，所以除算法本身的特點外，已知數據的分布特征必然影響插值精度。根據Nyquist采樣定理，已知數據的采樣頻率限制了其能反映的空間頻率上限，而因實測基準圖的空間分辨率有限，無法充分反映重力場細節(jié)變化，所以在高頻區(qū)域產生相對較大的插值誤差是難免的。

含大幅度重力場高頻信息的局部區(qū)域具有弱相關性，是重力輔助導航的重點區(qū)域。由以上分析知，高頻區(qū)域插值誤差較大的原因為區(qū)域內的已知數據不足以刻畫高頻細節(jié)信息。若可在該類區(qū)域補充基準數據量（即補充測點），則可以較小的測量代價，達到有效提高插值精度的目的。

2.2 基于小波的基準圖分析

小波分析具有時／頻域局部分析能力，是重／磁位場數據分析的有力工具。根據相關文獻［8］的研究可知，重力場的多分辨率小波細節(jié)信息可感知不同深度地質異常對應的重力異常的空間頻率信息。而通過觀察重力場空間頻率分布發(fā)現，淺層異常產生的高頻重力信息通常對應一定頻段內的連續(xù)分量，若已知基準重力數據的小波細節(jié)可感知淺層異常對應頻段的部分頻率信息，則可確定可能含有高頻重力信息的區(qū)域。

考察不同分辨率基準圖小波細節(jié)信息分布，將原基準圖分別降采樣為2′×2′、3′×3′及4′×4′，采用重力場分析效果較好的bior3.5小波母函數［9］。因目標為確定高頻區(qū)域，所以僅進行1階小波分解。圖3所示流程可將不同分辨率基準圖小波細節(jié)轉化為相同大小，以便比較小波細節(jié)分布情況。圖3中的DHL、DLH、DHH分別為水平方向、垂直方向及對角線方向的小波細節(jié)信息。

圖3 小波細節(jié)轉換流程Fig.3 The flow chart of wavelet detail conversion

3種分辨率小波細節(jié)如圖4。2′×2′基準圖小波細節(jié)清晰，與插值誤差分布最相符；3′×3′的小波細節(jié)較模糊，但與誤差分布基本相符；4′×4′的細節(jié)非常模糊，與誤差分布相符度最差。因此，原基準數據須具有一定的密度才能探測高頻區(qū)域。

2.3 基于隨機采樣方案的測點分布

由以上分析知，小波細節(jié)可反映基準圖局部區(qū)域具有大插值誤差的可能程度，則可由隨機采樣策略設計新增測點的分布。Metropolis采樣可生成以任意分布為不變分布的馬爾科夫鏈，則可設計如下新增測點的隨機采樣方案：

1）以小波細節(jié)均值為閾值，在細節(jié)圖中篩選新增測點的候選位置。

2）以候選位置的小波細節(jié)為轉移概率，由Metropolis采樣產生測點位置。設新增測點數量為已知測點數量的20%。

圖4 基準圖小波細節(jié)分布Fig.4 The wavelet detail distribution of reference map

3）為防止測點過于集中，若新測點與已知測點重合或其8鄰域內已存在測點，放棄該點并繼續(xù)抽樣。

4）因孤立測點會降低測量工作的效率，在抽樣完畢后，刪除一定鄰域內無其余測點的孤立點并繼續(xù)抽樣，直至所有測點滿足以上條件。鄰域半徑設為15′。

基于以上策略，以3′×3′基準圖的小波細節(jié)生成的測點分布如圖5所示。

增加測點后的插值算法誤差分布如圖6，誤差統(tǒng)計如表2。因增加測點后，各插值算法仍為無偏估計，即誤差均接近0，表2中省去均值而只給出標準差。此外，設各算法原插值誤差標準差為a，新增測點后的誤差標準差為b，則表2中“精度提升”為（a－b）／a×100%。

圖5 新增測點分布Fig.5 The distribution of incremental measurement points

圖6 增加測點的插值誤差分布Fig.6 The error distribution of interpolation algorithm for measurement points increased

表2 增加測點后的插值誤差統(tǒng)計Tab.2 Statistical data of interpolation error for measurement points increased

可見，各算法在高頻區(qū)域的誤差有所減弱。由統(tǒng)計數據知，誤差的范圍、標準差均有一定程度的改善，而徑向基函數算法的改善程度最為顯著。

2.4 對比方案

在整個基準圖范圍內，按均勻分布生成測點（新測點不與已知測點重疊）。設新增測點數為1 000（相當于已知測點數量的61%），該方案的插值誤差統(tǒng)計數據如表3所示。

通過對比兩種方案可知，在插值精度提升程度基本相同的情況下，在高頻區(qū)域增加測點的方案所需的測點數量較少，而且誤差的極值較小。

表3 新測點均勻分布的插值算法誤差統(tǒng)計Tab.3 Statistical data of interpolation error for measurement points distributing uniformity

3 結 語

本文通過插值實驗指出，插值算法在重力場高頻區(qū)域往往含有較大的誤差。而后利用小波分析可感知淺層地質異常產生的重力異?？臻g頻率信息的特點，提出以基準圖小波分解的細節(jié)信息為轉移概率，由隨機采樣方法在重力場高頻區(qū)域增加測點，提高對高頻信息的刻畫能力，從而提高插值精度的策略。插值實驗表明，相對均勻分布新增測點的方案，在插值精度提升程度相同的情況下，本文策略所需的測點數較少，誤差的極值較小，以較小的代價獲得了較優(yōu)的插值效果。

［1］Lowreys J A.Passive Navigation Using Inertial Navigation Sensors and Maps［J］.Naval Engineers Journal，1997（5）：245－251

［2］成怡，劉繁明，郝燕玲.海洋重力圖加密仿真實驗研究［J］.中國慣性技術學報，2007，15（2）：206－208（Cheng Yi，Liu Fanming，Hao Yanling.Simulation Experiment Research on Marine Gravity Map Interpolation［J］.Journal of Chinese Inertial Technology，2007，15（2）：206－208）

［3］吳太旗，黃謨濤，陸秀平，等.重力場匹配導航的重力圖生成技術［J］.中國慣性技術學報，2007，15（4）：438－441（Wu Taiqi，Huang Motao，Lu Xiuping，et al.Gravity Map Creating Technology in Gravity Matching Navigation［J］.Journal of Chinese Inertial Technology，2007，15（4）：438－441）

［4］Renka R J.Multivariate Interpolation of Large Sets of Scattered Data［J］.ACM Transaction on Mathematical Software，1988，14（2）：139－148

［5］Wright G B.Radial Basis Function Interpolation：Numerical and Analytical Developments［D］.University of Colorado，2003

［6］Rippa S.An Algorithm for Selecting a Good Value for the Parametercin Radial Basis Function Interpolation［J］.Advances in Computational Mathematics，1999，11：193－210

［7］Wessel P.A General－Purpose Green’s Function－Based Interpolator［J］.Computer ＆ Geosciences，2009，35（6）：1 247－1 254.

［8］楊文采，施志群，侯遵澤，等.離散小波變換與重力異常多重分解［J］.地球物理學報，2001，44（4）：534－541（Yang Wencai，Shi Zhiqun，Hou Zunze，et al.Discrete Wavelet Transform for Multiple Decomposition of Gravity Anomalies［J］.Chinese Journal of Geophysics，2001，44（4）：534－541）

［9］李健，周云軒，許惠平.重力場數據處理中小波母函數的選擇［J］.物探與化探，2001，25（6）：410－417（Li Jian，Zhou Yunxuan，Xu Huiping.The Selection of Wavelet Generating Functions in Data－Processing of Gravity Field［J］.Geophysical ＆Geochemical Exploration，2001，25（6）：410－416）