基于深度加權的重力梯度數據聯合相關成像反演

侯振隆，鄭玉君，鞏恩普，程 浩

(1. 東北大學 資源與土木工程學院，遼寧 沈陽 110819;2. 東北大學 深部金屬礦山安全開采教育部重點實驗室，遼寧 沈陽 110819)

與傳統(tǒng)重力異常數據相比，重力梯度數據包含有更多的頻率信息，能夠更好地用于解譯地下礦產的分布情況.利用重力梯度數據可以進行地下目標體的成像與反演，相關成像及其反演是常用的方法之一.該成像方法最初被應用于自然電場數據[1]，是一種利用觀測異常與幾何函數間的相關系數描述地質體空間位置和形態(tài)的算法，能夠應用于不同類型的地球物理場數據，包括重力異常[2]、重力梯度數據[3]、磁異常[4]以及電阻率異常[5]等；在成像基礎上加入物性參數的變化量能夠反演地下目標體的物性分布，國內學者在該方向開展了較深入的研究并取得了良好的效果[6-8].然而，相關成像反演在深度方向上的分辨能力較低；根據其他反演方法的研究，引入深度加權函數雖可能改善結果的縱向分辨率，但當地下存在多個目標體的復雜情況時，結果仍與真實分布有一定差距.

因此，本文將聯合多個重力梯度數據，使用深度加權函數與分區(qū)加權的方法提高相關成像反演的空間分辨率.在數據試驗中，可證實提出的方法大幅提高了縱向分辨能力，且具有一定的可行性.

1 重力梯度數據聯合相關成像反演原理

設研究區(qū)域位于笛卡爾坐標系中，o為坐標原點，xoy為一水平面，z軸的正方向垂直向下；地下空間被劃分為一系列點元，第j個點元在第i個觀測點上引起的重力梯度Vαβ,j為

Vαβ,j(xi,yi,zi)=γρjνjBαβ,j(xi,yi,zi);α,β=x,y,z.

(1)

其中:γ表示萬有引力常數；ρj，νj分別為第j個點元的密度、體積；Bαβ,j表示第j個質點在第i個觀測點上引起的重力梯度數據的幾何函數，則第j個點元的重力梯度相關系數為

(2)

由于相關成像只能計算相關系數，不能反映地下的物性分布.因此，使用相關成像反演方法得到地下物性參數分布情況.首先，建立地下密度分布密度的初始模型，設定密度變化量Δρ，在初始模型的基礎上，利用相關系數與密度變化量的乘積，進行迭代計算，公式為

ρj,iter+1=ρj,iter+Cj×Δρ.

(3)

其中:iter表示迭代的次數;ρj,iter表示第iter次迭代過程中第j個點元的密度值;Cj的意義同公式(2)；Δρ表示密度變化量，根據Hou等的研究[7]，密度變化量可利用觀測數據的最大值除以G×C的最大值來設定:

Δρ=dmax/max[G×C].

(4)

其中，G表示正演矩陣.

通過聯合不同的梯度分量，能夠更加全面地反映地下異常體信息，降低反演問題的多解性，增加結果的可靠性.根據Hou等的研究[7]，多分量重力梯度數據聯合反演中的相關系數為

Cjoint,j=C1,j×C2,j×…×Cn,j.

(5)

當進行聯合反演時，Δρ存在多個不同的值，通常選取最小值參與反演.

如前文提到的，相關成像反演在深度方向上的分辨率有待提高，特別是對目標體底部位置與形態(tài)識別的精度較差.下面引入深度加權函數來提高反演的縱向分辨率，本文采用Commer等提出的空間梯度加權函數[9-10]，該函數能夠引入先驗深度信息改善反演結果，其形式為

(6)

其中:z表示點元中心深度；dz表示研究區(qū)域內地下深度范圍；α為經驗常數，一般為較小的正數；r表示縮放因子；zc1,zc2分別表示模型頂、底部深度,根據先驗地質信息取值.深度加權后的迭代反演公式為

ρj,iter+1=ρj,iter+Cjoint,j×Wj(z)×Δρ.

(7)

上述加權方法對地下空間存在孤立目標體時效果較好，如果地下不同深度上存在多個目標體，不能直接通過公式(7)加權，否則將導致所有異常體處于同一深度范圍內.因此，在處理該情況時，首先根據先驗信息劃分研究區(qū)域，為不同區(qū)域分別設定該區(qū)內目標體的頂、底埋深，再應用公式(7)進行反演.

2 模型試驗

為驗證提出方法的有效性，設計立方體組合模型進行數據試驗.模型位于笛卡爾坐標系下，研究區(qū)域設定在坐標軸x，y方向范圍均為0～2 000 m、z方向范圍為0～1 000 m；等分地下空間為20×20×20個立方體，小立方體大小為100 m×100 m×50 m.組合模型由3個立方體A，B和C組成(圖1)，A立方體大小為200 m×300 m×300 m，剩余密度為0.5 g/cm3；B立方體大小為200 m×200 m×200 m，剩余密度為0.7 g/cm3；C立方體大小為500 m×500 m×200 m，剩余密度為0.5 g/cm3；A和B立方體的頂面埋深均為200 m，C立方體的頂面埋深為700 m.為組合模型產生的重力梯度數據加入5%的高斯噪聲驗證算法的抗噪能力；在試驗中使用x=1 500 m，y=1 000 m，z=300 m三個剖面進行結果分析.

根據秦朋波的研究[11]，對比分析不同重力梯度分量的特點，可得出Vxz|Vyz|Vzz組合反演效果最好.因此，本文聯合這三個梯度分量開展試驗，并設定最大迭代次數為100次.為對比深度加權改善的效果，首先使用未改進的相關成像反演(見圖2).經16次迭代計算后，結果的最大值達到0.7 g/cm3，與理論模型一致.在圖2c中，可見反演的橫向分辨率較高，水平方向上的形態(tài)和模型基本一致.但是由圖2a和圖2b中識別出的結果可見，深度方向上的分辨率較低，模型B的下延深度過大，而模型C幾乎無法辨識出.

下面使用基于空間分區(qū)深度加權的方法進行反演.首先對地下空間進行分區(qū)，使分區(qū)后每個立方體獨自占據一個區(qū)域，進而使用深度加權函數為立方體A，B，C分別賦值.立方體A，B的zc1均取200 m，由于下底面深度不同，則模型A的zc2=500 m，模型B的zc2=400 m；對于立方體C，zc1=700 m，zc2=900 m.加權函數中α=0.001，dz=1 000 m，r=20.在16次迭代后，反演結果見圖3，剩余密度最大值達到0.7 g/cm3，與理論模型保持一致.分別對比圖3a和圖2a，圖3b和圖2b，可見本文提出的算法能夠較準確地計算出下底深度，并能同時計算出不同深度的目標體的分布范圍，反演的縱向分辨率和未改進的算法相比大幅提升；結果的數值范圍也與理論模型基本一致.此外，結果還反映出本文方法具有良好的抗噪能力，可用于模擬實測數據的反演研究.

3 文頓鹽丘地區(qū)數據試驗

將本文方法應用于美國文頓鹽丘地區(qū)的航空實測重力梯度數據以驗證其可行性.許多學者都曾對文頓鹽丘地區(qū)開展了大量研究[12-16].該地區(qū)地下存在一剩余密度約為0.55 g/cm3的蓋巖，采集的異常數據便是由它產生的.根據Ennen等[13]、Geng等[14]的研究，蓋巖的深度范圍約在160～650 m.令數據位于WGS84坐標系下，x和y方向范圍分別為440.5～444.5 km,3 332.78～3 336.78 km，研究區(qū)域內均勻分布20×20個測點.設最大反演深度1 km，等分地下空間為20層.在反演前已對數據進行地形校正和高通濾波處理.試驗將選取x=442.56 km，y=3 334.44 km和z=375 m三個剖面進行分析.

根據模型試驗結論，仍使用Vxz|Vyz|Vzz組合進行聯合相關成像反演.在深度加權中，設α=0.001，dz=1 000 m，r=20；結合文頓鹽丘地區(qū)的先驗信息[12],設蓋巖的頂部深度zc1=200 m，底部深度zc2=500 m；根據蓋巖已知范圍，將研究區(qū)域分成“回”字形的兩部分，蓋巖模型占據中間區(qū)域.經過17次迭代得到反演結果(見圖4)，其中最大剩余密度達到0.5 g/cm3，接近蓋巖的已知密度信息；從圖4a和圖4b可見，蓋巖的深度范圍在200～500 m，其形態(tài)也與其他學者研究的結果一致[12-16].因此，認為本文提出的方法能夠用于實測數據的反演解釋.

4 結 語

本文基于深度加權技術，提出了多分量重力梯度數據聯合相關成像反演方法.在數據試驗中，聯合Vxz,Vyz和Vzz三個梯度分量實現反演，結果較為準確地反映出目標體的位置與形態(tài)，密度數值范圍和理論值基本一致；對比未改進的方法，深度分辨率大幅提高；同時，驗證了算法具有抗噪性.在文頓鹽丘地區(qū)的實測數據試驗中，反演出的蓋巖形態(tài)與密度分布均和已知資料一致.上述試驗證明了方法的有效性和可行性.

致謝感謝Bell Geospace公司提供的實測重力梯度數據.