重力和重力梯度基準(zhǔn)圖構(gòu)建的多參量聯(lián)合網(wǎng)格化

袁 園，盧鵬羽，劉 杰，羅偉東

（1.國(guó)家海洋局海底科學(xué)重點(diǎn)實(shí)驗(yàn)室,杭州310012；2.自然資源部第二海洋研究所,杭州310012；3.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,長(zhǎng)春130026；4.吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院,長(zhǎng)春130026；5.中國(guó)地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局,廣州510075）

0 引言

潛艇是現(xiàn)代海軍裝備中隱蔽性最強(qiáng)的裝備體系，憑借搭載的超遠(yuǎn)距離反艦導(dǎo)彈與核導(dǎo)彈，潛艇可以對(duì)敵方構(gòu)成戰(zhàn)術(shù)、戰(zhàn)役、戰(zhàn)略級(jí)別的致命威脅。無(wú)論是作為近海防御還是作為走向遠(yuǎn)洋的安全保障平臺(tái)，潛艇均已成為各國(guó)海上作戰(zhàn)集群中的重要力量之一[1]。

潛艇之所以受到各國(guó)海軍的高度重視，甚至被稱為國(guó)之利器，與潛艇具有的三大特點(diǎn)有重要關(guān)系。隱蔽性強(qiáng)、續(xù)航能力持久、戰(zhàn)斗力強(qiáng)是現(xiàn)代潛艇的主要特征。一般而言，潛艇都是航行在海底幾百米深處，很難被對(duì)方艦艇發(fā)現(xiàn)。只要帶上充足的補(bǔ)給物資，核潛艇就可以在海里潛航幾個(gè)月。大型攻擊型核動(dòng)力潛艇可以攜帶十幾枚甚至幾十枚核彈，戰(zhàn)斗力驚人。在潛艇的三大特征中，最為關(guān)鍵的就是潛艇的隱蔽性能。沒(méi)有良好的隱蔽性，潛艇的作戰(zhàn)效能將大大降低，甚至成為敵方的活靶子。因此，潛艇在深海中航行需要安全可靠、定位精度高、隱蔽性好的導(dǎo)航系統(tǒng)[2-4]。

重力及重力梯度輔助導(dǎo)航利用地球表面隨空間位置不同而連續(xù)變化的重力場(chǎng)來(lái)進(jìn)行導(dǎo)航定位，實(shí)時(shí)校正慣性導(dǎo)航系統(tǒng)因時(shí)間推移而產(chǎn)生的漂移誤差，從而提高慣導(dǎo)系統(tǒng)的導(dǎo)航定位精度[5-6]。重力及重力梯度基準(zhǔn)圖是重力輔助慣性導(dǎo)航系統(tǒng)的重要組成部分，基準(zhǔn)圖的精度和分辨率特征直接決定著定位的精度。從20世紀(jì)70年代末開(kāi)始，海軍、國(guó)家海洋局和中國(guó)地質(zhì)調(diào)查局等測(cè)繪部門(mén)綜合采用航海、衛(wèi)星測(cè)高和航空等測(cè)量手段對(duì)我國(guó)近海和深遠(yuǎn)海部分區(qū)域進(jìn)行了重力測(cè)量作業(yè)，初步建立了海洋重力信息庫(kù)，獲取了多種比例尺的海洋重力測(cè)量數(shù)據(jù)。測(cè)量在局部海區(qū)的分辨率已達(dá)到5km，精度為2mGal。

在實(shí)際生產(chǎn)過(guò)程中，很難獲取大面積的高分辨率原始測(cè)量數(shù)據(jù)，只能通過(guò)網(wǎng)格化插值技術(shù)對(duì)已測(cè)海區(qū)進(jìn)行插值加密。傳統(tǒng)的插值方法包括反距離加權(quán)平均法、樣條函數(shù)法、最小曲率法和Kriging法等[7-12]。雖然并非所有的方法在最開(kāi)始都是針對(duì)地球物理數(shù)據(jù)進(jìn)行設(shè)計(jì)的，但它們?cè)诘厍蛭锢頂?shù)據(jù)的處理與解釋中均已得到了廣泛應(yīng)用。除此之外，Cordell等[13-15]利用等效源法對(duì)位場(chǎng)數(shù)據(jù)進(jìn)行了網(wǎng)格插值研究。郭良輝等[16]將數(shù)據(jù)插值的過(guò)程看作一個(gè)反演過(guò)程，通過(guò)建立插值算子、待插值點(diǎn)及完整數(shù)據(jù)的線性方程，反復(fù)迭代計(jì)算得到了完整的數(shù)據(jù)，進(jìn)而達(dá)到了插值的目的。

然而，傳統(tǒng)的船載和機(jī)載重力測(cè)量的測(cè)線間隔通常會(huì)遠(yuǎn)大于測(cè)線上測(cè)點(diǎn)的間距，僅利用單分量的重力數(shù)據(jù)很難恢復(fù)測(cè)線之間被遺漏的一些有效異常信息，且易出現(xiàn)虛假異常，致使出現(xiàn)錯(cuò)誤的成果解釋。隨著全張量重力梯度測(cè)量技術(shù)的發(fā)展，各梯度分量可以同時(shí)反映來(lái)自重力場(chǎng)的橫向和縱向梯度變化。因此，為充分利用水平重力梯度分量在水平方向上的梯度特征，本文將開(kāi)展重力梯度分量多參量的聯(lián)合網(wǎng)格化研究，在避免虛假異常的情況下提高網(wǎng)格數(shù)據(jù)的空間分辨率，以更好地服務(wù)于潛艇水下重力輔助導(dǎo)航系統(tǒng)的重力基準(zhǔn)圖構(gòu)建。

1 重力梯度數(shù)據(jù)聯(lián)合網(wǎng)格化原理

在Descartes坐標(biāo)系中，重力異常為地質(zhì)體的引力位V在z方向上的一階導(dǎo)數(shù)，而重力梯度張量則為V在x、y、z3個(gè)方向上的二階偏導(dǎo)數(shù)。各重力梯度分量組成了重力梯度張量矩陣

其中，

因此，矩陣T為對(duì)稱矩陣。由位場(chǎng)理論可知，引力位V在自由空間滿足Laplace方程

根據(jù)式（2）和式（3），可知T中只含有5個(gè)獨(dú)立的梯度張量成分，即Txx、Txz、Tyx、Tyy、Tyz，而全張量重力梯度儀則提供了Txx、Txz、Tyx、Tyy、Tyz和Tzz共6個(gè)張量成分。

根據(jù)引力位與各重力梯度張量成分在波數(shù)域中的關(guān)系，可得Txx、Txz、Tyx、Tyy、Tyz、Tzz這 6個(gè)分量與引力位在波數(shù)域中的轉(zhuǎn)化公式

式（4）中，kx、ky是x和y方向的波數(shù)，。式（4）中的所有量均代表著其在波數(shù)域中的數(shù)據(jù)。在噪音水平相同的情況下，根據(jù)式（4），可在均衡多個(gè)重力梯度信息的基礎(chǔ)上，得到對(duì)應(yīng)的引力位

將引力位V對(duì)z方向求導(dǎo)，可得到對(duì)應(yīng)網(wǎng)格點(diǎn)在插值后的重力異常。考慮到實(shí)際測(cè)量過(guò)程中噪音的影響，重力張量數(shù)據(jù)在空間域中不能滿足Txx+Tyy+Tzz=Δ2V=0的關(guān)系。因此，在計(jì)算時(shí)也需要考慮Tzz分量的轉(zhuǎn)化。

而在進(jìn)行實(shí)際處理時(shí)，當(dāng)各重力梯度分量的噪音水平不一或需要突出某一個(gè)或某幾個(gè)分量的特性時(shí)，需要按照實(shí)際情況對(duì)不同梯度分量進(jìn)行加權(quán)處理。 設(shè)Txx、Txz、Tyx、Tyy、Tyz和Tzz分量的權(quán)重系數(shù)分別為a、b、c、d、e和f。根據(jù)最小二乘理論，需要使R2值達(dá)到最小

將式（6）對(duì)V進(jìn)行求導(dǎo)，令其導(dǎo)數(shù)為0，可得到加權(quán)后的引力位V的轉(zhuǎn)化公式

利用式（5）和式（7）對(duì)重力張量數(shù)據(jù)進(jìn)行插值。首先，在離測(cè)點(diǎn)較近的位置使用通用的插值方法進(jìn)行賦值，而將遠(yuǎn)離測(cè)點(diǎn)位置的值預(yù)設(shè)為0；隨后，根據(jù)式（5）和式（7）得到引力位，再利用式（4）計(jì)算各重力張量值。此時(shí)，若將近測(cè)點(diǎn)的位置賦替換為零值，而遠(yuǎn)離測(cè)點(diǎn)的值保持不變，重復(fù)上述計(jì)算過(guò)程，便可以減少由于數(shù)據(jù)的初始任意賦值而引發(fā)的轉(zhuǎn)化結(jié)果不連續(xù)的問(wèn)題。反復(fù)迭代上述過(guò)程，當(dāng)在測(cè)點(diǎn)位置得到的結(jié)果收斂于觀測(cè)結(jié)果時(shí)，終止迭代過(guò)程，可得到最終的插值結(jié)果。

2 模型試驗(yàn)

為驗(yàn)證本文方法的正確性，需建立一個(gè)簡(jiǎn)單的球體模型，球體埋深為25m，球體直徑為5m，剩余密度為1g/cm3，觀測(cè)網(wǎng)格的間距為30m。球體位于觀測(cè)點(diǎn)之間，球體位置與測(cè)點(diǎn)分布如圖1所示。首先，通過(guò)正演計(jì)算得到256×256的網(wǎng)格數(shù)據(jù)，并將其作為理論數(shù)據(jù)。然后，對(duì)數(shù)據(jù)抽稀得到如測(cè)點(diǎn)分布的模擬測(cè)量數(shù)據(jù)。球體的理論正演值如圖2所示。

圖1 球體模型與測(cè)點(diǎn)分布平面圖Fig.1 Sphere model and measuring point distribution

圖2 球體模型正演重力值Fig.2 Forward gravity anomaly of sphere model

采用兩種不同的方式進(jìn)行插值試驗(yàn)。第1種方式為，直接利用3次樣條spline函數(shù)對(duì)重力數(shù)據(jù)進(jìn)行插值。第2種方式為，先利用spline方法對(duì)重力張量數(shù)據(jù)進(jìn)行插值，然后再將得到的重力張量梯度值利用式（5）進(jìn)行計(jì)算，最后將結(jié)果對(duì)z方向進(jìn)行求導(dǎo)。得到的兩組重力插值結(jié)果如圖3所示，綠色圓形代表著球體的位置。從圖3可以看出，直接利用spline進(jìn)行插值得到的重力異常呈不對(duì)稱分布，且異常中心點(diǎn)偏離球體正確位置。本文方法的插值結(jié)果不僅能夠準(zhǔn)確定位異常峰值的位置，而且插值的異常結(jié)果更加接近理論值。雖然由于Fourier變換處理的相關(guān)影響，圖3（b）中的球體周?chē)嬖诓痪鶆虻姆植记闆r，但即便如此，該方法仍能體現(xiàn)較好的插值結(jié)果。經(jīng)過(guò)兩種方法插值得到的重力值與理論重力值的均方根誤差分別為2.3×10-3gu和 1.9×10-3gu（gu 為國(guó)際重力單位，1gu=10-6m/s2）。針對(duì)以上的球體模型，由于它遠(yuǎn)離各測(cè)點(diǎn)，因而從理論上而言，使用不涉及z方向的重力分量會(huì)有更好的反映。通過(guò)增加水平梯度分量的權(quán)重，會(huì)取得更好的結(jié)果。當(dāng)Txy權(quán)重為56%、Txx為12%、Tyy為12%、Tzx和Tzy分別為6%時(shí)，根據(jù)式（7）計(jì)算取得的均方根誤差為1.75×10-3gu。

為驗(yàn)證方法的有效性，需建立一個(gè)復(fù)雜的地質(zhì)體模型，該模型由多個(gè)棱柱體和一個(gè)球體組成。棱柱體走向與測(cè)網(wǎng)之間存在夾角，為了體現(xiàn)重力張量梯度數(shù)據(jù)的優(yōu)勢(shì)，在模型深部加入了一個(gè)較大球體。其中，棱柱體的剩余密度為1g/cm3，長(zhǎng)和寬均為4m，高為10m，上底面深度為20m；球體的剩余密度為-2g/cm3，中心埋深為100m，球體半徑為10m。復(fù)雜模型的整體立體三維圖如圖4所示，模型的正演重力異常圖如圖5所示。使用與簡(jiǎn)單球體模型同樣的處理流程，但對(duì)2組流程采取不同的采樣間隔，得到的插值結(jié)果如圖6所示。

圖3 簡(jiǎn)單球體模型的兩種重力插值結(jié)果的對(duì)比圖Fig.3 Comparison the gravity anomalies of two different interpolation methods for simple spherical model

圖4 復(fù)雜模型整體立體三維圖Fig.4 3D view of complex model

圖5 復(fù)雜模型正演重力結(jié)果Fig.5 Forward gravity anomaly of complex model

從圖6可以看出，在測(cè)點(diǎn)間距為32m時(shí)，直接利用spline插值得到的結(jié)果不能明顯體現(xiàn)模型中棱柱體所產(chǎn)生的異常，而考慮了重力張量梯度數(shù)據(jù)的插值會(huì)取得更好的結(jié)果。進(jìn)一步增大采樣間隔，當(dāng)測(cè)點(diǎn)間距為64m時(shí)，雖然兩種方法均已無(wú)法取得較準(zhǔn)確的結(jié)果，但重力張量梯度聯(lián)合插值的結(jié)果更加接近理論正演結(jié)果。其中，spline插值結(jié)果的均方根誤差為1.5×10-2gu，重力梯度聯(lián)合插值的均方根誤差為1.2×10-2gu，優(yōu)于直接進(jìn)行spline插值的結(jié)果。提高Txy的權(quán)重到20%，插值結(jié)果的均方根誤差為1.0×10-2gu。

圖6 復(fù)雜模型的插值結(jié)果對(duì)比Fig.6 Interpolated gravity anomalies results of complex model with different methods

為驗(yàn)證噪聲對(duì)重力梯度聯(lián)合插值的影響，在測(cè)點(diǎn)間隔為32m的數(shù)據(jù)中加入白噪聲進(jìn)行插值試驗(yàn)。對(duì)Tzx分量加入異常峰值為10%的白噪聲，對(duì)其他梯度分量加入對(duì)應(yīng)峰值為5%的白噪聲。圖7給出的是直接spline插值與權(quán)重相同的重力張量插值的結(jié)果。從圖7可以看出，直接插值得到的重力值分布已經(jīng)出現(xiàn)了變形，而基于重力梯度聯(lián)合的插值結(jié)果基本保持了不變，驗(yàn)證了方法的抗噪性。若將Tzx分量的加權(quán)適當(dāng)減小，會(huì)取得更好的結(jié)果。選擇利用200m～350m范圍內(nèi)的數(shù)據(jù)進(jìn)行插值分析，spline插值結(jié)果的均方根誤差為7.0×10-2gu，重力梯度聯(lián)合插值結(jié)果的均方根誤差為6.6×10-2gu。若將Tzx分量的權(quán)重降低至10%，插值結(jié)果與理論重力值的均方根誤差為6.4×10-2gu。

圖7 含噪情況下的插值結(jié)果對(duì)比Fig.7 Interpolated gravity anomalies results of complex model with noisy data

3 結(jié)論

潛艇水下重力輔助匹配導(dǎo)航技術(shù)的發(fā)展，對(duì)重力及重力梯度基準(zhǔn)圖的精度和分辨率提出了更高的需求。然而，目前我國(guó)海上重力測(cè)量數(shù)據(jù)的比例尺還屬于中小比例尺，不能完全滿足重力輔助匹配導(dǎo)航的分辨率要求，需要對(duì)其進(jìn)行網(wǎng)格化插值。傳統(tǒng)的插值方法在稀疏測(cè)線區(qū)域經(jīng)常會(huì)導(dǎo)致虛假異常和異常位置偏離等情況。因此，本文避開(kāi)了傳統(tǒng)數(shù)據(jù)網(wǎng)格化的思想，提出了重力梯度多分量的聯(lián)合網(wǎng)格化方法，并通過(guò)多個(gè)模型數(shù)據(jù)驗(yàn)證了方法的準(zhǔn)確性，證實(shí)了該方法相對(duì)傳統(tǒng)網(wǎng)格化技術(shù)存在一定優(yōu)勢(shì)。