彭國民， 劉展

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

0 引言

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

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

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

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

1 重力正演

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

Gm=d,

(1)

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

2 反演方法

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

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

(2)

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

(3)

其中1

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

(4)

(5)

其中：

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

2.2 目標函數(shù)構建

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

(6)

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

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

(7)

2.3 優(yōu)化求解

(8)

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

(9)

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

(10)

(11)

(12)

其中I為單位矩陣.

(13)

由于構建的反演目標函數(shù)(7)中涉及到參數(shù)q，本文將在后面的理論模型測試和實際資料應用中探討q值的選擇對反演結果的影響，即反演模型的不確定性分析.

2.4 正則化因子選取

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

(14)

3 理論模型測試

3.1 模型1

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

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

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

(15)

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

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

圖3 利用光滑反演方法反演含有不同強度干擾重力數(shù)據(jù)的反演結果 (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ù)的反演結果 (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ù)的反演結果 (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ù)的反演結果 (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

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

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

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

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

表1 使用不同方法反演含有不同強度干擾 重力數(shù)據(jù)的反演結果的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 使用不同方法反演含有不同強度干擾 重力數(shù)據(jù)的反演結果的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 不同強度噪聲的重力數(shù)據(jù) (a) 3%； (b) 5%； (c) 10%.Fig.8 Gravity data with different noise level

圖9 利用光滑反演方法反演含有不同強度干擾重力數(shù)據(jù)的反演結果 (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ù)的反演結果 (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ù)的反演結果 (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ù) 黑色三角為重力觀測點.Fig.13 The residual gravity anomaly of San Nicolas deposit Black triangles indicate gravity observation points.

4 實際資料應用

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

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

5 結論

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

圖14 q=1.1時的反演結果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線的縱向切片； (c) 沿著圖13中N1-N2線的縱向切片； (d) 預測的重力異常數(shù)據(jù).(b)和(c)中的黑色實線框代表礦體真實邊界.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時的反演結果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線的縱向切片； (c) 沿著圖13中N1-N2線的縱向切片； (d) 預測的重力異常數(shù)據(jù).(b)和(c)中的黑色實線框代表礦體真實邊界.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時的反演結果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線的縱向切片； (c) 沿著圖13中N1-N2線的縱向切片； (d) 預測的重力異常數(shù)據(jù).(b)和(c)中的黑色實線框代表礦體真實邊界.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時的反演結果 (a) 深度為300 m處的水平切片； (b) 沿著圖13中E1-E2線的縱向切片； (c) 沿著圖13中N1-N2線的縱向切片； (d) 預測的重力異常數(shù)據(jù).(b)和(c)中的黑色實線框代表礦體真實邊界.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.

致謝感謝加拿大英屬哥倫比亞大學GIF組的公開實測重力數(shù)據(jù)，感謝三位評審專家為本文提出了中肯的修改意見.