基于q-高斯分布和零階最小熵正則化的三維重力聚焦反演

彭國(guó)民， 劉展

1 齊魯工業(yè)大學(xué)(山東省科學(xué)院)，山東省計(jì)算中心(國(guó)家超級(jí)計(jì)算濟(jì)南中心)，濟(jì)南 250014 2 中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 青島 266580

0 引言

重力勘探已廣泛應(yīng)用于礦產(chǎn)、油氣等資源勘探中，重力反演作為重力勘探的重要定量解釋手段，能夠得到地下密度的空間分布(Nabighian，2005；姚長(zhǎng)利等，2007；劉展等，2011；孟小紅等，2012；王萬(wàn)銀等，2014)，并結(jié)合已有的其他資料進(jìn)行后續(xù)的地質(zhì)和地球物理解釋.眾所周知，重力異常反演具有多解性和不穩(wěn)定性，目前的三維重力異常反演方法主要針對(duì)反演目標(biāo)函數(shù)中的模型項(xiàng)采用有效的正則化技術(shù)來(lái)獲得可靠、穩(wěn)定且具有一定特點(diǎn)的地質(zhì)目標(biāo)體模型(劉銀萍等，2013；高秀鶴和黃大年，2017；李澤林等，2019；周少偉等，2021；李芳等，2021).

從觀測(cè)數(shù)據(jù)擬合角度，三維重力反演問(wèn)題通常構(gòu)建為一最小二乘優(yōu)化問(wèn)題，通過(guò)極小化觀測(cè)數(shù)據(jù)和由當(dāng)前模型正演得到的模擬數(shù)據(jù)之間的誤差，得到能夠擬合觀測(cè)數(shù)據(jù)的地下密度模型.傳統(tǒng)的三維重力反演方法假設(shè)數(shù)據(jù)殘差(觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)之差)服從高斯分布(Li and Oldenburg，1998)，需要可靠地估計(jì)觀測(cè)誤差.此外，野外采集的重力數(shù)據(jù)的觀測(cè)誤差可能較大，甚至含有極端異常值.因此，觀測(cè)誤差的大小和可靠估計(jì)是傳統(tǒng)的三維重力反演方法面臨的挑戰(zhàn)之一.

盡管基于Tikhonov正則化的三維重力異常反演能夠獲得穩(wěn)定解，但是當(dāng)?shù)叵旅芏犬惓ｓw存在尖銳邊界時(shí)，基于L2范數(shù)正則化的三維重力異常反演得到的密度異常體的邊界比較光滑，這給后續(xù)的地質(zhì)和地球物理解釋帶來(lái)一定的困難.針對(duì)此問(wèn)題，眾多學(xué)者提出了基于非L2范數(shù)正則化的非光滑反演方法.大多數(shù)非光滑反演方法采用某一函數(shù)作用于物性模型，通過(guò)極小化目標(biāo)函數(shù)來(lái)獲得具有尖銳邊界的模型解.Last和Kubik(1983)首先基于最小體積約束思想利用最小支撐函數(shù)實(shí)現(xiàn)了重力異常聚焦反演.Portniaguine和Zhdanov(1999)把最小支撐函數(shù)拓展應(yīng)用到模型梯度中，開(kāi)發(fā)了基于最小支撐梯度函數(shù)的三維重力異常聚焦反演.Farquharson和Oldenburg(1998)提出了基于廣義的模型度量的非線(xiàn)性反演方法，能夠利用一定的表征模型復(fù)雜度的函數(shù)產(chǎn)生稀疏解.Sun和Li(2014)、Li 等(2018)基于物性模型的Lp范數(shù)實(shí)現(xiàn)了重磁異常稀疏反演方法.一些學(xué)者分別利用模型的總變差函數(shù)(Bertete-Aguirre et al., 2002)、柯西范數(shù)(Pilkington, 2009；彭國(guó)民等，2018b)、模型指數(shù)變化性質(zhì)(Zhao et al., 2016)、L0擬范數(shù)(Meng, 2018)、最小熵函數(shù)(Rezaie, 2019)實(shí)現(xiàn)了非光滑重磁異常反演問(wèn)題.此外，在重磁異常反演中可以利用已知的物性信息獲得具有尖銳邊界的模型解.Krahenbuhl和Li(2006)利用已知的物性信息開(kāi)發(fā)了二進(jìn)制反演方法來(lái)改善成像結(jié)果.Liu和Jin(2020)實(shí)現(xiàn)了基于改進(jìn)的模糊c均值聚類(lèi)算法的重力異常反演方法.

相比于高斯分布，基于Tsallis熵最大化原理推導(dǎo)得出的q-高斯分布已應(yīng)用于復(fù)雜系統(tǒng)和交叉學(xué)科領(lǐng)域中，且在擬合實(shí)測(cè)數(shù)據(jù)方面具有更靈活的優(yōu)點(diǎn)(Da Silva et al., 2020).基于最小熵正則化的地球物理反演可以獲得異常體與圍巖之間的尖銳邊界，已被應(yīng)用于大地電磁和重力數(shù)據(jù)反演中(Ramos et al., 1999；Rezaie, 2019).為此，本文提出了一種基于q-高斯分布和零階最小熵正則化的三維重力異常聚焦反演方法，該方法不需要估計(jì)數(shù)據(jù)誤差，且對(duì)較大數(shù)據(jù)誤差的敏感性較弱，同時(shí)能夠獲得聚焦反演結(jié)果，且不需要確定合理的聚焦參數(shù)，并將本文方法應(yīng)用于理論模型和實(shí)際資料以驗(yàn)證本文方法的有效性.

1 重力正演

重力正演采用點(diǎn)元法，首先將地下剖分成若干個(gè)小長(zhǎng)方體單元，每個(gè)剖分單元賦予一剩余密度值，然后計(jì)算每個(gè)剖分單元產(chǎn)生的重力異常，根據(jù)重力場(chǎng)的疊加原理，將每個(gè)剖分單元產(chǎn)生的重力異常進(jìn)行疊加，得到觀測(cè)點(diǎn)位置處的重力異常(Blakely，1995)，即：

Gm=d,

(1)

其中G為重力靈敏度矩陣，m為密度模型向量，d為重力異常數(shù)據(jù)向量.

2 反演方法

2.1 q-高斯分布和零階最小熵正則化

Tsallis熵是統(tǒng)計(jì)物理學(xué)中的一種信息度量方式，能夠?qū)哂胁灰?guī)則形狀的非可加性系統(tǒng)給出較好的物理解釋?zhuān)鋽?shù)學(xué)定義式為(Da Silva et al., 2020)：

(2)

其中k>0，q為一實(shí)參數(shù)，p(x)為概率密度函數(shù).在一定的約束條件下通過(guò)最大化Tsallis熵能夠獲得q-高斯分布，該分布的一個(gè)重要統(tǒng)計(jì)性質(zhì)是更加接近于數(shù)據(jù)殘差滿(mǎn)足的分布特點(diǎn).q-高斯分布函數(shù)可表示為(Da Silva et al., 2020)：

(3)

其中1

基于零階最小熵正則化思想，本文采用如下函數(shù)(Ramos et al., 1999)作為反演目標(biāo)函數(shù)的穩(wěn)定函數(shù)，即：

(4)

(5)

其中：

為一對(duì)角最小熵矩陣，起聚焦反演作用.由于參數(shù)Q是由模型m和參數(shù)δ得到，而參數(shù)δ取值為10-15能夠得到較好的反演結(jié)果，并且具有較好的穩(wěn)定性，因此，該對(duì)角矩陣不需要事先確定合理的聚焦參數(shù).

2.2 目標(biāo)函數(shù)構(gòu)建

傳統(tǒng)的三維重力反演方法假設(shè)數(shù)據(jù)殘差服從高斯分布，基于最小模型范數(shù)正則化的三維重力光滑反演目標(biāo)函數(shù)為：

(6)

其中‖·‖2為L(zhǎng)2范數(shù)，dobs為重力觀測(cè)數(shù)據(jù)，Wd為數(shù)據(jù)加權(quán)矩陣，當(dāng)不同觀測(cè)數(shù)據(jù)的誤差互不相關(guān)時(shí)，該矩陣為一對(duì)角矩陣，其對(duì)角線(xiàn)元素為觀測(cè)噪聲的標(biāo)準(zhǔn)差的倒數(shù)，mref為參考密度模型，Wm為一對(duì)角深度加權(quán)矩陣，其對(duì)角線(xiàn)元素為1/zj，zj為第j個(gè)小長(zhǎng)方體單元中心的深度(Peng et al., 2018a)，λ為正則化因子，控制著數(shù)據(jù)擬合項(xiàng)和模型擬合項(xiàng)對(duì)目標(biāo)函數(shù)的相對(duì)貢獻(xiàn).由目標(biāo)函數(shù)(6)可看出，傳統(tǒng)的三維重力光滑反演需要可靠地估計(jì)觀測(cè)誤差以確定數(shù)據(jù)加權(quán)矩陣Wd，此外，基于最小模型范數(shù)正則化的光滑反演得到的異常體邊界比較平滑.

為了避免確定數(shù)據(jù)加權(quán)矩陣以及獲得具有尖銳邊界的反演模型，本文構(gòu)建基于q-高斯分布和零階最小熵正則化的三維重力聚焦反演目標(biāo)函數(shù)，即：

(7)

2.3 優(yōu)化求解

(8)

首先從變換后目標(biāo)函數(shù)(8)的梯度角度分析該方法對(duì)數(shù)據(jù)誤差的敏感性，目標(biāo)函數(shù)(8)的梯度為：

(9)

下面利用高斯牛頓(GN)算法對(duì)變換后的目標(biāo)函數(shù)(8)進(jìn)行最優(yōu)化求解，其法方程為：

(10)

(11)

(12)

其中I為單位矩陣.

(13)

由于構(gòu)建的反演目標(biāo)函數(shù)(7)中涉及到參數(shù)q，本文將在后面的理論模型測(cè)試和實(shí)際資料應(yīng)用中探討q值的選擇對(duì)反演結(jié)果的影響，即反演模型的不確定性分析.

2.4 正則化因子選取

正則化因子的選擇對(duì)于三維重力反演的收斂性及反演模型的可靠性是非常重要的.選擇的正則化因子應(yīng)該能使得反演模型在一定程度上擬合觀測(cè)數(shù)據(jù)，且同時(shí)有最小的結(jié)構(gòu)復(fù)雜度(Farquharson and Oldenburg，2004).本文采用一種呈現(xiàn)遞減規(guī)律變化的正則化因子選取方法，即：

(14)

3 理論模型測(cè)試

3.1 模型1

為了測(cè)試本文反演方法的有效性，設(shè)計(jì)含有兩個(gè)大小相同的規(guī)則異常體的模型(圖1)，模型大小為1000 m×600 m×500 m，兩個(gè)異常體的大小均為200 m×200 m×100 m，其中心位置分別為(300 m,300 m,150 m)、(700 m,300 m,210 m).該模型的背景密度值為0.0 g·cm-3，兩個(gè)異常體的剩余密度均為1.0 g·cm-3.在地表觀測(cè)重力數(shù)據(jù)，其觀測(cè)點(diǎn)數(shù)為39×29=1131，觀測(cè)點(diǎn)間距在x、y方向上分別為25 m、20 m.由于野外采集的重力數(shù)據(jù)不可避免地帶有觀測(cè)誤差，在生成的正演模擬數(shù)據(jù)中加入均值為0、標(biāo)準(zhǔn)差分別為模擬數(shù)據(jù)大小的3%、5%、10%的高斯干擾，得到含有不同強(qiáng)度干擾的重力數(shù)據(jù)(圖2).

圖1 理論模型示意圖Fig.1 Sketch map of synthetic model

將該模型剖分成40×30×25=30000個(gè)剖分單元，每個(gè)剖分單元的大小為25 m×20 m×20 m.本文將密度模型的最小值、最大值分別設(shè)置為0.0 g·cm-3、1.0 g·cm-3.將本文反演方法應(yīng)用于這三種帶有干擾的重力數(shù)據(jù)，以探討本文反演方法對(duì)含有不同強(qiáng)度干擾的重力數(shù)據(jù)的適用性，并與傳統(tǒng)的三維重力光滑反演方法進(jìn)行對(duì)比.為了定量評(píng)價(jià)反演結(jié)果，本文采用如下兩種統(tǒng)計(jì)度量，即模型擬合度mRMS和Pearson相關(guān)系數(shù)R，其定義式分別為：

(15)

圖3為利用光滑反演方法反演含有不同強(qiáng)度干擾重力數(shù)據(jù)的反演結(jié)果，圖4—圖6分別為使用不同q值反演含有3%、5%、10%干擾重力數(shù)據(jù)的反演結(jié)果，圖中水平切片圖的深度均為150 m，縱向剖面圖均沿著x方向且y=300 m，黑色實(shí)線(xiàn)為異常體的真實(shí)邊界.表1為使用不同方法反演含有不同強(qiáng)度干擾重力數(shù)據(jù)的反演結(jié)果的mRMS和R值.對(duì)于傳統(tǒng)的重力光滑反演方法，從圖3和表1可看出：①盡管異常體的位置能夠反演出來(lái)，但是異常體的邊界不是很清晰，且反演的密度值也不能較好地反映異常體的真實(shí)值；②隨著干擾強(qiáng)度變大，反演模型的mRMS值變大，R值變小，反演結(jié)果變差.對(duì)于本文聚焦反演方法，從圖4—圖6和表1可看出：①當(dāng)q取值為1.1、1.5、2.0時(shí)，隨著干擾強(qiáng)度變大，反演模型仍具有比較清晰的異常體邊界，且反演密度值也比較接近于真實(shí)值；②當(dāng)q取值為2.5時(shí)，盡管本文方法得到的反演結(jié)果相對(duì)變差，但是也優(yōu)于傳統(tǒng)的重力光滑反演結(jié)果；③在該理論模型測(cè)試中，當(dāng)干擾強(qiáng)度為3%時(shí)，q取值為1.5反演結(jié)果具有最小的mRMS值和最大的R值，q取值為2.0，次之；當(dāng)干擾強(qiáng)度為5%時(shí)，q取值為1.5反演結(jié)果具有最小的mRMS值和最大的R值，q取值為1.1，次之；當(dāng)干擾強(qiáng)度為10%時(shí)，q取值為1.1反演結(jié)果具有最小的mRMS值和最大的R值，q取值為1.5，次之.

圖2 含有不同強(qiáng)度干擾的重力數(shù)據(jù) (a) 3%； (b) 5%； (c) 10%.Fig.2 Gravity data with noise of different level

圖3 利用光滑反演方法反演含有不同強(qiáng)度干擾重力數(shù)據(jù)的反演結(jié)果 (a) 3%，平面圖； (b) 3%，剖面圖； (c) 5%，平面圖； (d) 5%，剖面圖； (e) 10%，平面圖； (f) 10%，剖面圖.Fig.3 Inversion results from inversion of gravity data with noise of different level using smooth inversion method (a) 3%, horizontal slice; (b) 3%, vertical section; (c) 5%, horizontal slice; (d) 5%, vertical section; (e) 10%, horizontal slice; (f) 10%, vertical section.

圖4 使用不同q值反演含有3%干擾重力數(shù)據(jù)的反演結(jié)果 (a) q=1.1，平面圖； (b) q=1.1，剖面圖； (c) q=1.5，平面圖； (d) q=1.5，剖面圖； (e) q=2.0，平面圖； (f) q=2.0，剖面圖； (g) q=2.5，平面圖； (h) q=2.5，剖面圖.Fig.4 Inversion results from inversion of gravity data with 3% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

圖5 使用不同q值反演含有5%干擾重力數(shù)據(jù)的反演結(jié)果 (a) q=1.1，平面圖； (b) q=1.1，剖面圖； (c) q=1.5，平面圖； (d) q=1.5，剖面圖； (e) q=2.0，平面圖； (f) q=2.0，剖面圖； (g) q=2.5，平面圖； (h) q=2.5，剖面圖.Fig.5 Inversion results from inversion of gravity data with 5% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

圖6 使用不同q值反演含有10%干擾重力數(shù)據(jù)的反演結(jié)果 (a) q=1.1，平面圖； (b) q=1.1，剖面圖； (c) q=1.5，平面圖； (d) q=1.5，剖面圖； (e) q=2.0，平面圖； (f) q=2.0，剖面圖； (g) q=2.5，平面圖； (h) q=2.5，剖面圖.Fig.6 Inversion results from inversion of gravity data with 10% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

3.2 模型2

設(shè)計(jì)含有兩個(gè)方向相反、大小不同的傾斜脈狀體的模型(圖7)，模型大小為600 m×600 m×300 m，背景密度值為0.0 g·cm-3.兩個(gè)傾斜脈狀體的頂部埋深均為地表以下40 m位置處，左邊傾斜脈狀體的寬度為80 m，且向下延伸100 m，而右邊傾斜脈狀體的寬度為100 m，且向下延伸120 m，每個(gè)傾斜脈狀體的剩余密度均為1.0 g·cm-3.在地表觀測(cè)重力數(shù)據(jù)，其觀測(cè)點(diǎn)數(shù)為29×29=841，觀測(cè)點(diǎn)間距在x、y方向上均為20 m.由于野外采集的重力數(shù)據(jù)不可避免地帶有觀測(cè)誤差，在生成的正演模擬數(shù)據(jù)中加入均值為0、標(biāo)準(zhǔn)差分別為模擬數(shù)據(jù)大小的3%、5%、10%的高斯干擾，得到含有不同強(qiáng)度干擾的重力數(shù)據(jù)(圖8).

圖7 理論模型示意圖Fig.7 Sketch map of synthetic model

將該模型剖分成30×30×15=13500個(gè)剖分單元，每個(gè)剖分單元的大小為20 m×20 m×20 m.將本文反演方法應(yīng)用于這三種帶有干擾的重力數(shù)據(jù)，以探討本文反演方法對(duì)含有不同強(qiáng)度干擾的重力數(shù)據(jù)的適用性，并與傳統(tǒng)的三維重力光滑反演方法進(jìn)行對(duì)比.

圖9為利用光滑反演方法反演含有不同強(qiáng)度干擾重力數(shù)據(jù)的反演結(jié)果，圖10—圖12分別為使用不同q值反演含有3%、5%、10%干擾重力數(shù)據(jù)的反演結(jié)果，圖中水平切片圖的深度均為60 m，縱向剖面圖均沿著x方向且y=300 m，黑色實(shí)線(xiàn)為傾斜脈狀體的真實(shí)邊界.表2為使用不同方法反演含有不同強(qiáng)度干擾重力數(shù)據(jù)的反演結(jié)果的mRMS和R值.對(duì)于傳統(tǒng)的重力光滑反演方法，從圖9和表2可看出：①盡管兩個(gè)傾斜脈狀體的位置和傾斜方向均能夠反演出來(lái)，但是反演出的傾斜脈狀體的邊界不是很清晰，且反演的密度值不能較好地反映傾斜脈狀體的真實(shí)值；②隨著干擾強(qiáng)度變大，反演模型的mRMS值變大，R值變小，反演結(jié)果變差.對(duì)于本文聚焦反演方法，從圖10—圖12和表2可看出：①隨著干擾強(qiáng)度變大，反演模型也具有比較清晰的異常體邊界，且反演密度值也比較接近于真實(shí)值；②對(duì)于干擾強(qiáng)度為10%的重力數(shù)據(jù)，本文方法仍能較好地得到傾斜脈狀體的形態(tài)和密度值，反演結(jié)果優(yōu)于傳統(tǒng)的重力光滑反演結(jié)果；③在該理論模型測(cè)試中，q取值為2.0反演結(jié)果具有最小的mRMS值和最大的R值，q取值為1.5，次之.

表1 使用不同方法反演含有不同強(qiáng)度干擾 重力數(shù)據(jù)的反演結(jié)果的mRMS和R值Table 1 The values of mRMS and R of inversion results from inversion of gravity data with noise of different level using different inversion methods

表2 使用不同方法反演含有不同強(qiáng)度干擾 重力數(shù)據(jù)的反演結(jié)果的mRMS和R值Table 2 The values of mRMS and R of inversion results from inversion of gravity data with noise of different level using different inversion methods

圖8 不同強(qiáng)度噪聲的重力數(shù)據(jù) (a) 3%； (b) 5%； (c) 10%.Fig.8 Gravity data with different noise level

圖9 利用光滑反演方法反演含有不同強(qiáng)度干擾重力數(shù)據(jù)的反演結(jié)果 (a) 3%，平面圖； (b) 3%，剖面圖； (c) 5%，平面圖； (d) 5%，剖面圖； (e) 10%，平面圖； (f) 10%，剖面圖.Fig.9 Inversion results from inversion of gravity data with noise of different level using smooth inversion method (a) 3%, horizontal slice; (b) 3%, vertical section; (c) 5%, horizontal slice; (d) 5%, vertical section; (e) 10%, horizontal slice; (f) 10%, vertical section.

圖10 使用不同q值反演含有3%干擾重力數(shù)據(jù)的反演結(jié)果 (a) q=1.1，平面圖； (b) q=1.1，剖面圖； (c) q=1.5，平面圖； (d) q=1.5，剖面圖； (e) q=2.0，平面圖； (f) q=2.0，剖面圖； (g) q=2.5，平面圖； (h) q=2.5，剖面圖.Fig.10 Inversion results from inversion of gravity data with 3% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

圖11 使用本文方法基于不同q值反演含有5%干擾的重力數(shù)據(jù)的反演結(jié)果 (a) q=1.1，平面圖； (b) q=1.1，剖面圖； (c) q=1.5，平面圖； (d) q=1.5，剖面圖； (e) q=2.0，平面圖； (f) q=2.0，剖面圖； (g) q=2.5，平面圖； (h) q=2.5，剖面圖.Fig.11 Inversion results from inversion of gravity data with 5% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

圖13 San Nicolas礦區(qū)的剩余重力異常數(shù)據(jù) 黑色三角為重力觀測(cè)點(diǎn).Fig.13 The residual gravity anomaly of San Nicolas deposit Black triangles indicate gravity observation points.

4 實(shí)際資料應(yīng)用

為了驗(yàn)證本文提出方法的實(shí)用性，現(xiàn)將本文方法應(yīng)用于墨西哥薩卡特卡斯州圣尼古拉斯(San Nicolas)礦區(qū)的重力數(shù)據(jù).San Nicolas礦區(qū)上侏羅系到下白堊系的長(zhǎng)英質(zhì)和鐵鎂質(zhì)火山巖中發(fā)育有火山成因塊狀硫化物，火山基巖被第三系的火山碎屑角礫巖覆蓋，其厚度為50～150 m.San Nicolas礦區(qū)塊狀硫化物礦體的最高密度為3.5 g·cm-3，圍巖的密度為2.3～2.7 g·cm-3，因此該礦區(qū)的密度差最大可達(dá)1.2 g·cm-3(Phillips et al.，2001)，這為利用重力勘探方法直接找塊狀硫化物礦體提供了物性基礎(chǔ).圖13為San Nicolas礦區(qū)的剩余重力異常圖，由鉆井信息知大約中間位置處的高重力異常是由地下的塊狀硫化物礦體引起的.

將該礦區(qū)地下空間離散成43×36×20=30960個(gè)剖分單元，每個(gè)剖分單元的大小為50 m×50 m×50 m，在反演過(guò)程中設(shè)置密度差的上、下邊界約束分別為1.2 g·cm-3、-0.1 g·cm-3(Phillips et al.，2001).圖14—圖17分別為q取值為1.1、1.5、2.0、2.5時(shí)的反演結(jié)果，每個(gè)圖中的子圖(a)為反演結(jié)果在深度為300 m位置處的水平切片，子圖(b)為反演結(jié)果沿著圖13中的E1-E2線(xiàn)的縱向切片，子圖(c)為反演結(jié)果沿著圖13中的N1-N2線(xiàn)的縱向切片，子圖(d)為由反演結(jié)果進(jìn)行正演得到的預(yù)測(cè)重力異常數(shù)據(jù)，圖中的黑色實(shí)線(xiàn)框代表礦體的真實(shí)邊界.從這些圖中可看出，q取值為1.1(對(duì)應(yīng)反演結(jié)果圖14)或1.5(對(duì)應(yīng)反演結(jié)果圖15)時(shí)，本文方法能夠得到與已知信息比較一致的反演模型.在圖14a和15a中，中間的密度異常體為San Nicolas礦體，位于東北方向的異常體是由鐵鎂質(zhì)火山巖引起的(Phillips et al.，2001)；在圖14b、c和圖15b、c中，已知的鉆井信息表明，塊狀硫化物礦體的頂部埋深在地表以下150～220 m，礦體向下延伸到地表以下400～450 m，反演的塊狀硫化物礦體的位置與已知信息比較一致，且與礦體的真實(shí)邊界也比較吻合；從圖14d和圖15d可看出，反演結(jié)果也能夠較好地?cái)M合觀測(cè)數(shù)據(jù).

5 結(jié)論

本文提出了一種基于q-高斯分布和零階最小熵正則化的三維重力聚焦反演方法，該方法不需要估計(jì)數(shù)據(jù)誤差，且對(duì)較大數(shù)據(jù)誤差的敏感性較弱，同時(shí)能夠獲得聚焦反演結(jié)果，也不需要事先確定合理的聚焦參數(shù).理論模型測(cè)試表明，相比于傳統(tǒng)的光滑反演方法，本文方法能夠較好地刻畫(huà)地質(zhì)體邊界和密度值，且對(duì)含有較強(qiáng)干擾的重力數(shù)據(jù)具有一定的適用性.墨西哥圣尼古拉斯礦區(qū)的實(shí)測(cè)重力數(shù)據(jù)應(yīng)用表明，本文方法反演的硫化物礦體與已知資料具有較好的一致性.理論模型和實(shí)際資料測(cè)試中的反演模型不確定性分析表明，q取值為1.5能夠得到令人滿(mǎn)意的反演結(jié)果.

圖14 q=1.1時(shí)的反演結(jié)果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線(xiàn)的縱向切片； (c) 沿著圖13中N1-N2線(xiàn)的縱向切片； (d) 預(yù)測(cè)的重力異常數(shù)據(jù).(b)和(c)中的黑色實(shí)線(xiàn)框代表礦體真實(shí)邊界.Fig.14 Inversion results for q=1.1 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

圖15 q=1.5時(shí)的反演結(jié)果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線(xiàn)的縱向切片； (c) 沿著圖13中N1-N2線(xiàn)的縱向切片； (d) 預(yù)測(cè)的重力異常數(shù)據(jù).(b)和(c)中的黑色實(shí)線(xiàn)框代表礦體真實(shí)邊界.Fig.15 Inversion results for q=1.5 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

圖16 q=2.0時(shí)的反演結(jié)果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線(xiàn)的縱向切片； (c) 沿著圖13中N1-N2線(xiàn)的縱向切片； (d) 預(yù)測(cè)的重力異常數(shù)據(jù).(b)和(c)中的黑色實(shí)線(xiàn)框代表礦體真實(shí)邊界.Fig.16 Inversion results for q=2.0 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

圖17 q=2.5時(shí)的反演結(jié)果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線(xiàn)的縱向切片； (c) 沿著圖13中N1-N2線(xiàn)的縱向切片； (d) 預(yù)測(cè)的重力異常數(shù)據(jù).(b)和(c)中的黑色實(shí)線(xiàn)框代表礦體真實(shí)邊界.Fig.17 Inversion results for q=2.5 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

致謝感謝加拿大英屬哥倫比亞大學(xué)GIF組的公開(kāi)實(shí)測(cè)重力數(shù)據(jù)，感謝三位評(píng)審專(zhuān)家為本文提出了中肯的修改意見(jiàn).