基于小波多分辨技術(shù)重力場(chǎng)相關(guān)系數(shù)成像

劉志新，盧 海，2

（1：中國(guó)礦業(yè)大學(xué) 資源與地球科學(xué)學(xué)院，江蘇 徐州221116；2：淮安市水利勘測(cè)設(shè)計(jì)研究院有限公司，江蘇 淮安223001）

重力勘探是利用組成地殼的各種地層、巖體、礦體的密度差異所引起的重力變化而進(jìn)行地質(zhì)勘探的一種方法［1］.常規(guī)重力測(cè)量只觀測(cè)到重力異常，在重力數(shù)據(jù)后期處理時(shí)往往將重力異常值轉(zhuǎn)化成重力梯度異常和延拓重力異常［2-4］.常規(guī)重力異常反演方法是基于反演理論，使目標(biāo)函數(shù)達(dá)到極小的線性或非線性反演［5］，存在多解性.Patclla等在1997年首次提出將概率成像法用于自然電位異常的解釋?zhuān)?］.Abdelrahman等［7］提出了相鄰最小二乘重力異常的歸一化互相關(guān)公式［7］.郭良輝等［8］提出了重力異常三維相關(guān)成像方法和重力梯度數(shù)據(jù)三維相關(guān)成像方法，并提出了基于異常分離的三維相關(guān)成像方法來(lái)提高成像分辨率，通過(guò)合成Y型巖脈模型和合成多個(gè)直立長(zhǎng)方體組合模型的重力異常和重力梯度數(shù)據(jù)試驗(yàn)分析，驗(yàn)證了三維相關(guān)成像方法可顯示出異常地質(zhì)體的空間賦存狀態(tài)和等效剩余質(zhì)量分布，具有良好的縱向和橫向分辨率［8-10］.徐令周等［11］利用概率成像方法對(duì)自然電場(chǎng)和大地電磁數(shù)據(jù)進(jìn)行了研究，取得了一些進(jìn)展和認(rèn)識(shí)［12］.楊文采等［13］利用小波多分辨原理，在對(duì)中國(guó)布格重力異常多尺度分解的基礎(chǔ)上，反演了各種尺度意義下中國(guó)大陸地殼的密度差異，給出了中國(guó)大陸地殼中相對(duì)密度差異的空間分布，利用離散小波變換進(jìn)行重力異常多重分解為研究中國(guó)大陸地殼構(gòu)造，進(jìn)行深部地質(zhì)研究提供了新成果和新方法.

本文提出的重力場(chǎng)參數(shù)相關(guān)系數(shù)成像方法能準(zhǔn)確確定異常體的空間位置和延伸情況，首先給出重力/重力梯度異常相關(guān)系數(shù)的計(jì)算方法，結(jié)合設(shè)置的模型，得出基于小波多分辨技術(shù)的重力參數(shù)相關(guān)系數(shù)成像方法具有很高的橫向、縱向分辨率，且相對(duì)于傳統(tǒng)重力資料成像效果更加直觀和清晰.

1 方法原理

重力異常值是指地下密度異常體在空間中產(chǎn)生的引力的垂直分量，設(shè)地下任意坐標(biāo)為（xq，yq，zq）、體積為Vq的q點(diǎn)的剩余密度為Δσq，則它在空間上任意點(diǎn)（x，y，z）處產(chǎn)生的重力異常Δgq（x，y，z）可表示為

式中，G 為萬(wàn)有引力常數(shù)，G=6.666 7×10-11m3·kg-1·s-2

利用Abdelrahman等［7］提出的重力異常的相關(guān)公式，定義測(cè)區(qū)重力異常與q點(diǎn)剩余質(zhì)量產(chǎn)生的重力異常兩者的相關(guān)公式為

式中：Δg（xi，yi，zi）為測(cè)區(qū)第i個(gè)觀測(cè)點(diǎn)的實(shí)測(cè)重力異常；（xi，yi，zi）為該測(cè)點(diǎn)的坐標(biāo)；Δgq（xi，yi，zi）為測(cè)區(qū)下的q點(diǎn)剩余質(zhì)量在該觀測(cè)點(diǎn)的重力異常；N為總觀測(cè)點(diǎn)數(shù).

重力異常相關(guān)系數(shù)Aq的物理意義就是用地下q點(diǎn)在地表產(chǎn)生的重力異常值與實(shí)測(cè)重力異常值進(jìn)行互相關(guān)計(jì)算，表示地表重力異常由地下q點(diǎn)產(chǎn)生的可能性的大小，如果Aq為正，則該點(diǎn)質(zhì)量盈余，若為負(fù)，則該點(diǎn)質(zhì)量虧損.

重力梯度異常的相關(guān)系數(shù)為

2 模型設(shè)置及其重力場(chǎng)參數(shù)

計(jì)算模型三維空間范圍為0≤x≤20km、0≤y≤20km、-4km≤z≤0.由淺至深，由小到大共設(shè)置了10個(gè)異常體，分別為深層較大的A1、A2、A3，中層中等范圍的B1、B2、B3、B4，淺層較小的C1、C2、C3.模型中各異常體的具體空間位置及密度參數(shù)見(jiàn)表1.

表1 模型幾何參數(shù)和剩余密度Tab.1 Model geometry parameters and density

圖1為異常體在地面上的投影，圖2為Y=13.0km的斷面切片圖，該斷面上包括C1、C2、B1、B2、B3、A1、A2這7個(gè)異常體.

圖1 模型平面投影Fig.1 Model plane projection map

計(jì)算地質(zhì)體所產(chǎn)生的重力異常，通常是依據(jù)牛頓萬(wàn)有引力定律，首先計(jì)算地質(zhì)體產(chǎn)生的引力位，重力異常是引力位沿重力方向上的導(dǎo)數(shù)，而重力梯度異常是重力異常在空間3個(gè)方向上的導(dǎo)數(shù).如果代入計(jì)算的是地質(zhì)異常體的剩余密度，那么得到的就是重力異常（單位：10-5m·s-2）和重力梯度異常值（單位：10-9s-2）.

圖2 模型Y＝13.0km斷面切片F(xiàn)ig.2 Model Y＝13.0km section slice

圖3a是剩余重力異常圖，可以清楚地反映出A1、A2、A3、B1、B2、B3、B4幾個(gè)較大異常體，C1、C2、C3反映不是很清晰；圖3b是重力梯度異常Vzy，對(duì)各個(gè)級(jí)別的異常反應(yīng)良好，通過(guò)該圖可以區(qū)分出各個(gè)異常體引起的異常.可知，重力梯度異常相比于重力異常更能明顯地突出異常特征，其結(jié)果具有更高的分辨率.

圖3 重力異常及重力梯度異常圖Fig.3 Gravity anomaly and gravity gradient anomaly map

3 重力場(chǎng)參數(shù)相關(guān)系數(shù)

3.1 斷面重力場(chǎng)參數(shù)相關(guān)系數(shù)

圖4 Y＝13.0km處重力異常及重力梯度異常相關(guān)系數(shù)成像圖Fig.4 Correlation coefficients imaging calculated by gravity anomaly and gravity gradient anomaly，Y＝13.0km

首先，選取模型Y=13.0km測(cè)線的重力異常/重力梯度異常曲線，求取其相關(guān)系數(shù)并成像，結(jié)果如圖4所示.分析可知：①重力異常相關(guān)系數(shù)只能反映A1、A2這兩個(gè)深層較大的異常；②Vzx相關(guān)系數(shù)基本上能反映所有異常體的基本情況；③Vzz相關(guān)系數(shù)能反映出A1、B1、B2和B3，對(duì)A2的反映較為模糊；④基于小波多分辨重力異常相關(guān)系數(shù)能清晰地反映 A2、B1、B2、B3、C1和C2，對(duì) A1的反應(yīng)不是很明顯；⑤基于小波多分辨重力梯度異常Vzx相關(guān)系數(shù)準(zhǔn)確地反映 A1、A2，B1、B2、B3、C1和 C2所有異常體的位置、深度以及剩余密度的正負(fù)；⑥基于小波多分辨重力梯度異常Vzz相關(guān)系數(shù)也準(zhǔn)確地反映出A1、A2、B1、B2、B3、C1和C2所有異常體的位置、深度以及剩余密度的正負(fù).

從以上分析結(jié)果可知，經(jīng)過(guò)基于小波多分辨技術(shù)的重力異常分離后，成像準(zhǔn)確度有明顯提高.尤其是重力異常Δg，本來(lái)只能辨別A1和A2這兩個(gè)大的異常，經(jīng)過(guò)小波多分辨分解后，已經(jīng)能辨別出所有異常體的分布情況，另外對(duì)Vzz的相關(guān)系數(shù)成像分辨率的提高也很明顯.可見(jiàn)基于小波多分辨分析的相關(guān)系數(shù)在提高解釋縱向和橫向分辨率方面可行且有效.

3.2 小波函數(shù)及分解層數(shù)的選定

利用小波多分辨技術(shù)進(jìn)行重力異常分離，重點(diǎn)是小波函數(shù)的選取.本文討論使用較多且效果較好的幾種小波函數(shù)，主要包括Haar小波、Daubechies小波系、Biorthogonal小波系、Coiflet小波系、Symlets小波系.判斷一個(gè)小波的優(yōu)劣主要根據(jù)處理結(jié)果和理論結(jié)果的誤差ε大小.

式中：ΔgF為重力異常Δg的分解值；ΔgL為Δg的理論計(jì)算值.

ε作為選定小波的依據(jù)，分別用重力異常Δg曲線和重力梯度異常Vzx曲線的分解結(jié)果與理論結(jié)果的分析比較，分別給出適合重力異常和重力梯度異常的小波函數(shù).

對(duì)于重力異常Δg曲線，模型Y=13km剖面上的重力異常曲線如圖5a所示，單獨(dú)計(jì)算A1和A2在地面引起的異常值即理論異常值，如圖5b所示，選取不同小波函數(shù)在剖面上的重力異常中分解出的A1和A2在地面引起的異常值即分解異常值，求取兩者的誤差值，可以得到各小波函數(shù)分解誤差的大小，認(rèn)為誤差小的小波函數(shù)即為適合重力異常曲線分解的小波，表2給出其誤差值.

表2給出了Haar小波、Daubechies小波系、Biorthogonal小波系、Coiflet小波系、Symlets小波系中的部分小波函數(shù)分解結(jié)果的誤差值，從中可以看出，db8小波函數(shù)分解結(jié)果的誤差值為3.729，明顯小于其他小波函數(shù)分解結(jié)果的誤差值.因此，認(rèn)為db8小波函數(shù)最適應(yīng)重力異常曲線，選定db8小波函數(shù)對(duì)重力異常數(shù)據(jù)進(jìn)行多分辨分析.

圖5 重力異常理論曲線Fig.5 Theoretical curves of gravity anomaly

表2 不同小波函數(shù)重力異常分解結(jié)果誤差值Tab.2 Error values of different wavelet function decomposition values of gravity anomaly

表3給出了Haar小波、Daubechies小波系、Biorthogonal小波系、Coiflet小波系、Symlets小波系中的部分小波函數(shù)Vzx分解結(jié)果的誤差值，從中可以看出，db9、bior6.8和sym8小波函數(shù)分解結(jié)果的誤差值相近，且遠(yuǎn)小于其他小波函數(shù)分解結(jié)果的誤差值，但db9小波函數(shù)分解誤差更小.認(rèn)為db9小波函數(shù)最適應(yīng)重力梯度異常曲線，選定db9小波函數(shù)對(duì)重力梯度異常數(shù)據(jù)進(jìn)行多分辨分析.

表3 不同小波函數(shù)Vzx分解結(jié)果誤差值Tab.3 Error values of different wavelet function decomposition values of Vzx

小波分解階數(shù)的選取依據(jù)下述規(guī)則：對(duì)于模型而言，如果某一階分解結(jié)果最大值出現(xiàn)比上一階分解結(jié)果最大值小兩個(gè)數(shù)量級(jí)，就不再分解，認(rèn)為分解結(jié)果是小波波形和模型數(shù)據(jù)曲線不一致造成的誤差；對(duì)于實(shí)測(cè)數(shù)據(jù)而言，如果某一階分解結(jié)果最大值小于實(shí)測(cè)時(shí)的觀測(cè)均方誤差值，就不再分解，認(rèn)為分解結(jié)果是隨機(jī)誤差.

利用db8小波和db9小波分別對(duì)重力異常Δgq和重力梯度異常Vzx進(jìn)行四階分解，四階低頻逼近即為深層大范圍的密度體引起的異常，三階、四階高頻細(xì)節(jié)疊加就為中層中等范圍的異常，一階、二階高頻細(xì)節(jié)疊加就為淺層小范圍的異常.將分解的重力異常和重力梯度異常求取相關(guān)系數(shù)并進(jìn)行疊加，得到基于小波多分辨技術(shù)的相關(guān)系數(shù)結(jié)果如圖4d～4f所示，對(duì)比圖4a～4c和圖4d～4f發(fā)現(xiàn)Vzx和Vzz在縱向和橫向上對(duì)模型有一個(gè)比較好的分辨率，梯度異常的相關(guān)系數(shù)比之重力異常的相關(guān)系數(shù)分辨率更高.可以明顯看出經(jīng)過(guò)基于小波多分辨技術(shù)的重力異常分離后，成像準(zhǔn)確度有明顯提高.尤其是剩余重力異常Δgq，本來(lái)只能辨別兩個(gè)大的異常，經(jīng)過(guò)小波多分辨分解后，已經(jīng)能辨別出所有異常體的分布情況，另外對(duì)Vzz的相關(guān)系數(shù)成像分辨率的提高也很明顯.

4 模型結(jié)果

經(jīng)過(guò)3.1節(jié)對(duì)選取測(cè)線重力數(shù)據(jù)相關(guān)系數(shù)的檢驗(yàn)，利用基于小波多分辨分解的重力參數(shù)相關(guān)系數(shù)分析更加準(zhǔn)確.用db8二維小波對(duì)平面重力異常進(jìn)行分解以及db9二維小波對(duì)重力梯度異常Vzx進(jìn)行分解，可以直觀地發(fā)現(xiàn)已經(jīng)能將異常分解成不同范圍和量級(jí)的異常，將三維空間按100×100網(wǎng)格化，分解出的不同范圍的局部異常，最后將相關(guān)系數(shù)疊加，最終重力異常和重力梯度異常Vzx相關(guān)系數(shù)成像結(jié)果如圖6所示.

圖6三維效果圖和表1中模型的幾何參數(shù)相比較，發(fā)現(xiàn)利用相關(guān)系數(shù)值圈定的異常曲面，能反映出模型各異常體的位置和范圍，而且空間各點(diǎn)相關(guān)系數(shù)值的正負(fù)，表征了空間點(diǎn)剩余密度的正負(fù)，所以認(rèn)為該種基于小波多分辨分析的重力/重力梯度異常相關(guān)系數(shù)成像方法取得了非常直觀和清晰的結(jié)果.

圖6 基于小波多分辨的三維重力異常和重力梯度異常相關(guān)系數(shù)圖Fig.6 Correlation coefficient map of gravity anomaly and gravity gradient anomaly based on wavelet multi-resolution

5 結(jié)論

本文是利用基于小波多分辨分析的相關(guān)系數(shù)方法對(duì)模型重力資料進(jìn)行再處理，得出結(jié)論如下：

（1）重力梯度異常相關(guān)系數(shù)成像相比于重力異常相關(guān)系數(shù)成像具有更高的分辨率，這是因?yàn)橹亓μ荻犬惓１旧砭褪菍?duì)重力異常的求導(dǎo)，相比重力異常更能突出異常位置.

（2）選取模型Y=13.0km剖面的重力參數(shù)，分別求取了未經(jīng)小波多分辨分解的重力參數(shù)的相關(guān)系數(shù)值以及經(jīng)過(guò)小波多分辨分解的重力參數(shù)的相關(guān)系數(shù)，并將它們的結(jié)果進(jìn)行了對(duì)比分析.發(fā)現(xiàn)基于小波多分辨重力參數(shù)的相關(guān)系數(shù)在反映異常體的分辨率上有顯著的提高，所以利用基于小波多分辨分解的重力參數(shù)計(jì)算相關(guān)系數(shù)是非常必要的.

（3）基于小波多分辨技術(shù)的模型重力異常三維相關(guān)系數(shù)成像結(jié)果，基本能反映出所有異常體的空間位置及分布情況，與模型參數(shù)對(duì)應(yīng)良好，重力場(chǎng)參數(shù)相關(guān)系數(shù)成像效果明顯.

［1］ 曾華霖.重力場(chǎng)與重力勘探［M］.北京：地質(zhì)出版社，2005.ZENG Hualin.Gravity field and gravity exploration ［M］.Beijing：Geological Press，2005.

［2］ 霍志周.重力資料預(yù)處理系統(tǒng)研制［D］.西安：長(zhǎng)安大學(xué)，2006.HUO Zhizhou.Study of pretreatment system of gravity data［D］.Xi’an：Chang’an University，2006

［3］ 許德樹(shù)，曾華霖.優(yōu)選延拓技術(shù)及其在中國(guó)布格重力異常圖處理上的應(yīng)用［J］.現(xiàn)代地質(zhì)，2000，14（2）：215.XU Deshu， ZENG Hualin.Optimization continuation technology and its application of processing china bouguer gravity anomaly map［J］.Geosciences，2000，14（2）：215.

［4］ 曾華霖，許德樹(shù).最佳向上延拓高度的估計(jì)［J］.地學(xué)前緣.2002，9（2）：499.ZENG Hualin，XU Deshu.Estimate of the best upward continuation height［J］.Geoscience Frontiers，2002，9（2）：499.

［5］ 陳石.利用重力異常進(jìn)行三維斷層參數(shù)反演方法研究［D］.西安：長(zhǎng)安大學(xué)，2005.CHEN Shi.Research on the inversion method of 3D fault parameters by gravity anomaly ［D］.Xi’an：Chang’an University，2005

［6］ Mauriello P，Patclla D.Principles of probability tomography for natural-sources electromagnetic induction fields ［J］.Goephysics，1999，64：1404.

［7］ Abdelrahman E M，Bayoumi A I，Abdelhady Y E.Gravity interpretation using correlation factors between successive least squares residual anomalies［J］.Geophysics，1989，54：1614.

［8］ 郭良輝，孟小紅，石磊.重力異常分離的相關(guān)法［J］.地球物理學(xué)進(jìn)展，2008，23（5）：1425.GUO Lianghui，MENG Xiaohong，SHI Lei.The correlation method for gravity anomaly separation ［J］.Progress in Geophysics，2008，23（5）：1425.

［9］ 陳召曦，孟小紅，郭良輝.重磁數(shù)據(jù)三維物性反演方法進(jìn)展［J］.地球物理學(xué)進(jìn)展，2012，27（2）：503.CHEN Zhaoxi，MENG Xiaohong，GUO Lianghui.Review of 3D property inversion of gravity and magnetic data［J］.Progress in Geophysics，2012，27（2）：503.

［10］ 陳召曦，孟小紅，郭良輝.基于GPU并行的重力、重力梯度三維正演快速計(jì)算及反演策略［J］.地球物理學(xué)報(bào)，2012，55（12）：4069.CHEN Zhaoxi，MENG Xiaohong，GUO Lianghui.Threedimensional fast forward modeling and the inversion strategy for large scale gravity and gravimetry data based on GPU［J］.Chinese Journal of Geophysics，2012，55（12）：4069.

［11］ 徐令周，關(guān)繼騰，房文靜.高次導(dǎo)數(shù)的概率成像原理［J］.青島大學(xué)學(xué)報(bào)，2003，16（4）：32.XU Lingzhou，GUAN Jiteng，F(xiàn)ANG Wenjing.Probability tomography theory of higher derivative［J］.Journal of Qingdao University，2003，16（4）：32.

［12］ 毛立峰，王緒本，高永才.大地電磁概率成像的效果評(píng)價(jià)［J］.地球物理學(xué)報(bào)，2005，48（6）：679.MAO Lifeng， WANG Xuben， GAO Yongcai.Effect evaluation of electromagnetic probability imaging［J］.Chinese Journal of Geophysics，2005，48（6）：679.

［13］ 楊文采，施志群，侯遵澤.離散小波變換與重力異常多重分解［J］.地球物理學(xué)報(bào)，2001，44（4）：534.YANG Wencai，SHI Zhiqun，HOU Zhunze，et al.Discrete wavelet transform and gravity anomalies multiple decomposition［J］.Chinese Journal of Geophysics，2001，44（4）：534.