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

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

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

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

1 重力梯度數(shù)據(jù)聯(lián)合相關(guān)成像反演原理

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

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

(1)

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

(2)

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

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

(3)

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

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

(4)

其中，G表示正演矩陣.

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

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

(5)

當(dāng)進(jìn)行聯(lián)合反演時，Δρ存在多個不同的值，通常選取最小值參與反演.

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

(6)

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

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

(7)

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

2 模型試驗(yàn)

為驗(yàn)證提出方法的有效性，設(shè)計(jì)立方體組合模型進(jìn)行數(shù)據(jù)試驗(yàn).模型位于笛卡爾坐標(biāo)系下，研究區(qū)域設(shè)定在坐標(biāo)軸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.為組合模型產(chǎn)生的重力梯度數(shù)據(jù)加入5%的高斯噪聲驗(yàn)證算法的抗噪能力；在試驗(yàn)中使用x=1 500 m，y=1 000 m，z=300 m三個剖面進(jìn)行結(jié)果分析.

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

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

3 文頓鹽丘地區(qū)數(shù)據(jù)試驗(yàn)

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

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

4 結(jié) 語

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

致謝感謝Bell Geospace公司提供的實(shí)測重力梯度數(shù)據(jù).