基于模型降階的貝葉斯方法在三維重力反演中的實踐

劉彥， 呂慶田， 李曉斌， 祁光，趙金花， 嚴加永， 鄧震

1 中國地質科學院礦產(chǎn)資源研究所，國土資源部成礦作用與資源評價重點實驗室， 北京 1000372 中國地質科學院地球深部探測中心， 北京 1000373 中國地質科學院地球物理地球化學勘查研究所， 河北廊坊 0650004 河南理工大學， 河南焦作 454000

?

基于模型降階的貝葉斯方法在三維重力反演中的實踐

劉彥1，2， 呂慶田2，3*， 李曉斌4， 祁光1，2，趙金花1，2， 嚴加永1，2， 鄧震1，2

1 中國地質科學院礦產(chǎn)資源研究所，國土資源部成礦作用與資源評價重點實驗室， 北京 1000372 中國地質科學院地球深部探測中心， 北京 1000373 中國地質科學院地球物理地球化學勘查研究所， 河北廊坊 0650004 河南理工大學， 河南焦作 454000

三維重力反演是地質工作者了解地球深部構造，認知地下結構的重要手段.按照反演單元劃分，三維重力反演有離散多面體（Discrete）反演和網(wǎng)格節(jié)點（Voxels）反演兩種方式.離散多面體反演由于易于吸收先驗地質信息得到的理論場能夠很好地擬合觀測場，因此，在實際重力反演中更受歡迎.目前離散多面體重力反演中初始模型的建立方法繁雜不一，實際應用受到很大的限制.本文本著充分挖掘利用先驗信息和重力觀測數(shù)據(jù)得到豐富可靠的反演結果這一原則，以離散多面體反演技術為基礎，改進建模過程.在初始模型的建立中，吸收貝葉斯算法優(yōu)勢，采用隱馬爾科夫鏈改善樸素貝葉斯方法的分類效果，通過最大似然函數(shù)算法求解，再采取模型降階技術，固定所建模型中幾何體的形態(tài)或密度，達到在幾何體形態(tài)（x,y,z）、密度（σ）和重力值（g）五個參數(shù)中降低維數(shù)目的，從而減小高維不確定性和正演的計算量，由此反演計算的地質體密度和分布范圍相對更準確，更利于重現(xiàn)重力模型結構.通過單位球體和任意形態(tài)幾何體模擬實驗，以及安徽省泥河礦區(qū)三維重力反演實踐，得到非常接近實際的密度或重力值，大幅提高了三維重力反演的精度和效率，說明該方法是有效、實用的.

三維重力反演；離散多面體模擬；貝葉斯方法；最大似然函數(shù)；模型降階

1 引言

礦集區(qū)深部結構、構造關系以及成礦環(huán)境是控制礦床形成的直接因素，因此，無論對于成礦學還是深部資源勘查，厘清深部三維形態(tài)和物性（包括斷裂延伸、巖體的空間展布以及與地層的相互關系）都十分重要.一直以來，地質工作者試圖利用各種已知信息建立地下三維結構模型創(chuàng)造出多種方法，其中，以重力觀測數(shù)據(jù)為基礎以其他先驗地質信息為約束的三維（3 Dimensional，簡稱3D）重力反演，已經(jīng)成為了解地球深部結構的一種重要手段（Oldenburg et al.，1997，2007； Fullagar et al.，2000，2007； Williams et al.，2004； Welford and Hall，2007）.

按照反演單元劃分，三維重力反演通常有網(wǎng)格節(jié)點（Voxels）和離散多面體（Discrete）兩種方式（Li and Oldenburg，1998；姚長利，2002，2007）.離散多面體反演因其既易于吸收先驗地質、構造（包括地層傾向、斷層和礦化體等）信息，又可使反演的理論場與觀測場誤差較小而受到歡迎.對于地質信息約束下的三維重力反演，不少學者做過工作（Roy and Clowes，2000；Jessell，2001；Goleby et al.，2002；Malehmir et al.，2006，2009；Martin et al.，2007；Williams，2008；Wan et al.，2011；呂慶田等，2013）.綜觀過程，雖細節(jié)有所差別，但基本上可歸納為以下三步：①初始模型的建立，②反演的迭代與初始模型的升級，③模型的三維顯示，其中步驟①和②至關重要.對于步驟①初始模型的建立，目前方法繁雜不一，尚未形成統(tǒng)一完美的解決方案，事實上，合理的初始模型能夠減小后續(xù)重力反演的非唯一性，加快反演的收斂速度（McGaughey，2007；McInerney et al.，2007）.步驟②則是對步驟①的進一步優(yōu)化，以獲得符合觀測數(shù)據(jù)的物性參數(shù)空間分布.

貝葉斯方法（Bayes，1763）最初主要用于統(tǒng)計學領域，1987年法國Tarantola教授在其基礎上創(chuàng)立了貝葉斯反演理論（Tarantola，1987），之后，貝葉斯反演逐步應用于地球物理領域并得到廣泛深化和應用（Green，1995；Scales and Snieder，1997；Gao，1997；Gouveia et al.，1998；Grandis et al.，1999；Lasses et al.，2004；Tarantola，2005；Khodja et al.，2010；Stuart，2010；Fernández-Martínez et al.，2013）.貝葉斯反演理論框架建立在概率論的基礎上，它將數(shù)據(jù)信息與模型先驗信息聯(lián)系起來，通過模型的先驗信息來約束模型的后驗參數(shù).由于貝葉斯推理把反演問題與觀測數(shù)據(jù)信息、模型信息以及先驗信息聯(lián)系起來，使得求解約束泛函的極小點不需要計算梯度方向，因此，對約束項沒有光滑性要求；同時，貝葉斯方法不只給出單一的反演值，而是給出解的后驗分布采樣，在此基礎上得到解的條件期望、方差、置信區(qū)間等具有統(tǒng)計意義的結果，這使得貝葉斯方法能夠很好地應用先驗信息，得到豐富可靠的輸出結果.將貝葉斯算法運用到三維重力反演中，應該有助于初始模型的改進.但是，在貝葉斯理論中，測量數(shù)據(jù)和先驗信息包含在后驗概率密度函數(shù)（PPD）中，它可以解釋成模型的單點估計和不確定度等貝葉斯推斷，這些信息的獲取需要對反演問題進行優(yōu)化求最優(yōu)模型和在高維模型空間中對PPD進行采樣積分；而且，從模型的后驗概率密度函數(shù)中提取模型信息，需要計算它的模型估計、參數(shù)不確定度以及參數(shù)間的相關度等屬性.所有這些，無疑會使貝葉斯方法受礙于維數(shù)，并且將為解決正演問題而付出高昂的計算成本.

為此，本文將開展以貝葉斯方法為基礎的離散多面體三維重力反演研究，特別是在初始模型的建立中，吸收貝葉斯反演理論精華，將反演問題同觀測數(shù)據(jù)、模型以及先驗信息聯(lián)系起來，計算中采用模型降階技術克服反演的高維不確定性和復雜的正演成本，快速準確計算出地質體的密度和形態(tài)重現(xiàn)重力模型結構，再通過單位球體和任意形態(tài)幾何體反演模擬實驗以及安徽省泥河礦區(qū)應用實踐，驗證方法的正確性和實用性.

2 方法原理及實現(xiàn)過程

2.1 貝葉斯方法原理

如前所述，貝葉斯反演是綜合利用未知模型參數(shù)的先驗信息、實測數(shù)據(jù)以及由物理規(guī)律得到的正演模型等信息來反演模型參數(shù)，最后得到未知模型參數(shù)的后驗概率分布函數(shù)及其相應的特征量.設m表示未知的M維模型參數(shù)空間，d表示N維數(shù)據(jù)空間，在貝葉斯反演中都看成是隨機變量.實際應用中，數(shù)據(jù)d是已知的，模型參數(shù)化θ也是確定的，為了方便省略θ，得到貝葉斯公式（Carlin and Louis，2000）

（1）

這時P（d）是常量，P（d|m）為已知數(shù)據(jù)d條件下隨模型m變化的函數(shù)，稱為似然函數(shù)，用L（m）表示.上式可以改寫成

圖1 貝葉斯方法原理

（2）

其中P（m|d）表示后驗概率分布函數(shù)（PPD）.其原理可以簡單概括如圖1.

似然函數(shù)L（m）（或P（d|m））表達的是模型參數(shù)m與觀測數(shù)據(jù)d之間的相對概率，在貝葉斯反演計算中，它的大小反映了模型響應與觀測數(shù)據(jù)的適配程度，可見，似然函數(shù)L（m）是貝葉斯反演計算的關鍵參量，模型m的似然函數(shù)可以表示為（Scales and Snieder， 1997）

L（m）=

（3）

其中，g為正演算子（又稱核函數(shù)），d為觀測數(shù)據(jù)矢量，m為模型參數(shù)矢量，T表示轉置，Cd為先驗估計誤差（以下公式中各參數(shù)皆如此）.

引入目標函數(shù)E（m），公式（3）似然函數(shù)L（m）可以改寫為

（4）

其中

（5）

A是一個比例因子.

（6）

對應的協(xié)方差矩陣Cm為

（7）

可見，貝葉斯反演結果最終是由這個不確定度矩陣Cm和最大似然解mML共同定義，當數(shù)據(jù)誤差服從Gaussian分布且是線性的，反演的結果也應滿足Gaussian分布.實際反演通常存在欠參數(shù)化和超參數(shù)化兩種情況：對于欠參數(shù)化情況，一般可采用貝葉斯信息準則（BIC）（Akaike，1980）來判斷與數(shù)據(jù)信息一致的反演層數(shù)；對于超參數(shù)化情況，則采用傾向于簡單結構的帶先驗信息反演，類似于正則化反演方法（TikhonovandArsenin，1977）.本文帶先驗信息的重力反演屬于后者.

通常，重力觀測數(shù)據(jù)或其他先驗地質信息是有噪音的，反演的模型結果也會有誤差，假定這些噪音和誤差滿足Gaussian分布，在重力反演中似然函數(shù)L（m）、模型的先驗概率P（m）、目標函數(shù)E（m）可以類推表達如下：

（8）

（10）

（11）其中（ξ,η,ζ）為地質體體積元的坐標，（x,y,z）為觀測點的坐標，G為常數(shù)，等于6.672×10-11m3/（kg·s2）.

（12）

方程（12）的求解實質是加權二乘擬合，關鍵是在重力觀測數(shù)據(jù)擬合（第一項）和先驗估計（第二項）之間尋求一種折衷，其中權重系數(shù)的選取非常重要.權重系數(shù)越大，反演解越接近先驗估計；反之，反演解越接近于不穩(wěn)定的初始模型.常用方法有L-曲線法（Hansen，1992）和廣義交叉法（Wahba，1990），詳情參考文獻（FarquharsonandOldenburg，2000；HaberandOldenburg，2000）.

2.2 三維重力反演過程

有了以上算法基礎，以離散多面體為反演單元設計三維重力反演，包括：先驗模型的確定、初始模型的建立、重力反演迭代與初始模型的升級、模型的三維顯示，如圖2所示.詳情以下展開.

圖2 三維重力反演流程圖

2.2.1 準備數(shù)據(jù)、確定模型空間

收集已有信息、簡化地質單元、測量巖石物性、電磁或地震剖面以及鉆測井資料的綜合解釋等，定義建模區(qū)域.

2.2.2 根據(jù)已知先驗信息采用貝葉斯算法建立初始模型，開展重力反演迭代與初始模型的升級

樸素貝葉斯分類（Naive Bayes Classifier，或NBC）是貝葉斯反演理論中的一種簡單、穩(wěn)定的分類算法，其思想是：對于給定的待分類項，求解在此項出現(xiàn)的條件下各個類別出現(xiàn)的概率，并按照概率最大劃定待分類項的類別.建立樸素貝葉斯分類模型所需估計的參數(shù)很少，對缺失數(shù)據(jù)不太敏感，算法也比較簡單，運用該方法對先驗信息進行分類.考慮到通過樸素貝葉斯分類建立初始模型，無需了解模型參數(shù)之過程狀態(tài)序列，只需計算模型參數(shù)最終狀態(tài)的概率函數(shù)，因此，采用隱馬爾可夫模型（Hidden Markov Model，或 HMM）更合適.設計如下主要流程：

第一步，確定模型參數(shù)特征屬性及劃分：根據(jù)先驗信息p確定出模型參數(shù)矢量m.其中，模型參數(shù)矢量包括模型中各地質體的位置、埋深、厚度、密度等物理參數(shù)及其可能變化的范圍.

第二步，獲取所建模型體對應的觀測向量數(shù)據(jù)作為訓練樣本：根據(jù)先驗地質、地球物理信息和建模區(qū)域確定觀測數(shù)據(jù)與地下物理模型參數(shù)之間的函數(shù)關系，即參入反演的核函數(shù)g.再由公式d=gm，通過核函數(shù)g（見公式11）與模型參數(shù)矢量m正演計算出所建模型體的觀測向量數(shù)據(jù)d.

第三步，計算訓練樣本中每個類別的頻率：將模型體生成的觀測向量數(shù)據(jù)劃分為不同的類別.

這里是根據(jù)先驗信息推算模型區(qū)域每個地質體的最大似然解mML.對于重力貝葉斯反演，如公式（12）最大似然模型m的求解，涉及地質體坐標位置（ξ,η,ζ）、剩余質量（σ）和重力異常值（g）五個參數(shù)，計算維數(shù)較高，計算成本也高，而且影響反演解的收斂性增加不確定度.因此，借助模型降階技術采用固定部分參數(shù)減小計算維數(shù)的方式，固定模型中地質體的坐標位置使得類別參數(shù)降到最少.這樣，當體積確定時，根據(jù)“質量=密度×體積”，通過分段劃分地質體與區(qū)域圍巖的相對密度變化區(qū)間，等同于劃分剩余質量的變化區(qū)間，再根據(jù)萬有引力公式就可計算出模型地質體的重力異常數(shù)值范圍.

第四步，計算這一位置固定的地質體的相對密度可能概率.

第五步，使用分類器鑒別所述模型地質體的最大似然解并將其分類.按照上述最大似然解算法，依據(jù)公式（12）計算模型，要求模型參數(shù)最大限度滿足重力觀測值與先驗信息這一分類原則，確定模型中各個地質體（巖體）的坐標位置、埋深、厚度、密度等物理參數(shù).

圖3 貝葉斯算法程序

計算過程如圖3所示，概括為：根據(jù)已知的先驗信息選擇初始模型；根據(jù)所選擇的初始模型確定重力異常訓練樣本；計算訓練樣本中每個類別的頻率（歸結為計算模型中地質體剩余質量的變化區(qū)間）；通過貝葉斯分類器對觀測向量數(shù)據(jù)進行分類鑒別（運用模型地質體的最大似然解），確定初始模型參數(shù)的逼近區(qū)間.若逼近區(qū)間符合模型重力與觀測重力的擬合方差要求，則輸出反演結果；若該逼近區(qū)間不符合擬合方差要求，則根據(jù)先驗信息重新選擇初始模型.通過上述六步，就可針對性地給出模型中各地質體的密度值和分布范圍，得到較為準確的初始地質模型.具體程序見“基于貝葉斯反演理論的重力建模軟件”（劉彥等，2015）.

2.2.3 采用離散體模擬方法將初始地質模型構建成2.5D地質模型

對每一個賦予初始密度的模型體，采用人機交互法修改建立2D剖面模型，完成建模區(qū)域內所有2D剖面的重力模擬后，將每條2D剖面的模型走向方向長度延長同剖面長度，初始模型變?yōu)?.5D地質模型.再對2.5D模型進行重力曲線擬合、模型修改處理，直到達到滿意的數(shù)據(jù)擬合和獲得合理的地質模型.可采用如澳大利亞EncomTechnologyPtyLtd.公司開發(fā)的專業(yè)軟件EncomModelVisionProTM等重力數(shù)據(jù)處理、反演軟件完成（祁光等，2012）.

2.2.4 將2.5D地質模型拼接擬合成3D地質模型

得到系列2.5D地質模型后，按照剖面的空間順序依次將2.5D模型拼合成3D模型，計算3D模型的理論異常，并與實際異常對比，擬合誤差較大的區(qū)域，需返回到2.5D剖面重新修改，如此調整得到3D地質模型.

2.2.5 對3D地質模型進行可視化和結構解譯處理

采用三維可視化模擬軟件（如：GOCAD或EncomPA）對3D地質模型進行可視化和結構解譯處理，經(jīng)過可視化及結構解譯處理后的模型能夠形象直觀地反映地下三維密度結構.

3 模擬實驗

為驗證上述貝葉斯算法開展重力反演的有效性，選擇單位球體和任意形態(tài)幾何體兩種模型體開展模擬實驗.重力值可依據(jù)牛頓萬有引力定律計算，球體的計算公式為（曾華霖，2005）

（14）

其中Δg為球體的重力異常值，其單位取μGal，1μGal=10-8m·s-2，G為萬有引力常量，取值6.672×10-11m3/（kg·s2），M為場源質量，h為場源中心深度或埋深，R為球體半徑，X為觀測平面上測點到場源中心的水平距離，這里指測線距離，h、R、X的單位均為m，σ為模型體相對圍巖的密度差或相對密度，單位為kg·m-3.對于任意形態(tài)幾何體，可通過三角剖分法將其劃分成若干個極小單元的四面體，根據(jù)這些最小單元四面體底面中心點的坐標（x，y，z）與三角形面積計算每個四面體的體積，再與密度乘積得到四面體的質量.為簡化計算，任意形態(tài)幾何體剖分至百個數(shù)量級單元體，剖分得到的小四面體近似于單位質量球體，重力異常值可采用公式（14）來計算，最后對所有剖分體的重力異常值求和，得到任意形態(tài)幾何體的重力異常總值.

3.1 單位球體數(shù)值模擬

以單位球體為例，測點到場源中心的水平距離X和埋深h為固定值，采用上述貝葉斯算法，計算球體的密度范圍.假設單位球體處于測點坐標X0=1000 m處，埋深h0為800 m（如圖4a所示），球體與圍巖的相對密度（即密度差）σ為30 kg·m-3，測線總長2000 m，測點間距為10 m，加入正態(tài)分布的觀測噪聲，依據(jù)公式（14）正演計算得到重力觀測值，如圖4b所示.

采用設計的貝葉斯算法程序經(jīng)過10次迭代運算誤差小于0.02，滿足精度要求，得到模型的相對密度值為31.7578125±5.2734375 kg·m-3，運算過程如表1所示.再采用隱馬爾科夫鏈分類器鑒別的反演計算，得到模型與圍巖的相對密度值為30.3107438614098 kg·m-3，方差為1.87447932360261×10-5，更接近實際密度差30 kg·m-3.可見，隱馬爾科夫鏈分類器鑒別的貝葉斯反演結果相對更準確.

表1 單位球體模擬貝葉斯計算結果Table 1 Bayesian calculation results of unit sphere simulation

圖4 單位球體數(shù)值模擬示意圖

3.2 任意形態(tài)幾何體數(shù)值模擬

選擇與圍巖的相對密度σ為30 kg·m-3的任意形態(tài)幾何體正演計算重力數(shù)據(jù).為貼近實際，幾何體定義為長軸a等于300 m、短軸b等于30 m的透鏡體，三角剖分后如圖5a所示.透鏡體位于測點坐標X0=1000 m處，埋深h0為800 m，測線總長2000 m，測點間距為10 m.加入正態(tài)分布的噪聲，正演計算得到的重力值如圖5b所示.

同樣，采用上述貝葉斯算法程序計算，經(jīng)過10次迭代運算誤差小于0.02，滿足精度要求，得到相對密度值為31.7578125±5.2734375 kg·m-3，計算結果與單個球體相同.可見多個球體疊加運算對于貝葉斯算法影響并不大，同時也證明貝葉斯算法具有一定的穩(wěn)定性，對梯度方向沒有要求.采用隱馬爾科夫鏈分類器鑒別反演計算，得到剩余密度值為30.0626173207599 kg·m-3，方差為1.72748273715972×10-5，結果與實際給定的相對密度30 kg·m-3非常接近，說明本文設計的貝葉斯反演方法對復雜形體同樣有效.

圖5 任意形態(tài)幾何體數(shù)值模擬示意圖

4 安徽泥河礦區(qū)三維重力反演實踐

泥河礦床位于安徽省廬樅火山巖盆地西北邊緣，是長江中下游地區(qū)實施第二空間找礦以來在廬樅礦集區(qū)500～1500 m深度范圍內最先發(fā)現(xiàn)的大型“玢巖型”鐵礦床.礦區(qū)處于羅河—黃屯北東向成礦帶上，南西距羅河鐵礦床3 km，北東距龍橋鐵礦床13 km（吳明安等，2011）.泥河深部礦體的發(fā)現(xiàn)引起眾多的關注，由此產(chǎn)生一系列關于“玢巖型”鐵礦床成礦與控礦模式以及如何區(qū)分礦體與非礦體異常的思考.欲解決這些關乎找礦成敗的關鍵問題，需要精細解剖地球物理異常性質，深入認識礦區(qū)深部結構，明晰區(qū)域地質狀況.

早在2007年初，高銳等在廬樅礦集區(qū)及外圍布設兩個剖面總長達153.28 km的十字形深地震反射測線（高銳等，2010），2008年呂慶田等部署兩條穿過泥河、羅河、大包莊等主要礦體的10 km長的地震剖面測線（呂慶田等，2010），2009—2010年呂慶田等又在廬樅地區(qū)完成5條總長達326 km的深地震反射剖面（呂慶田等，2014，2015），這些反射地震成像很好地顯現(xiàn)出礦區(qū)所在區(qū)域的深部構造.不僅如此，針對泥河礦床更是實施完成不少具體工作：吳明安（2011）總結了泥河鐵礦的發(fā)現(xiàn)過程及區(qū)域找礦方向，趙文廣（2011）查明礦床的地質特征并分析了成因，認為輝石閃長玢巖是泥河鐵礦的成礦母巖，正長斑巖穿切火山巖地層和礦體，是成礦期后形成的脈巖；劉彥（2012）、祁光（2012）以泥河礦床為例開展了基于重力異常分離方法尋找深部隱伏鐵礦的實踐以及地質信息約束下的三維重磁建模工作；張昆（2014）開展了音頻大地電磁、可控源音頻大地電磁、瞬變電磁和復電阻率四種電磁法探測試驗，獲得泥河鐵礦區(qū)三維電性結構模型；范裕（2014）開展了礦床成巖成礦年代學研究，認為泥河鐵礦床的成巖成礦作用幾乎同時發(fā)生，盆地內130 Ma左右的輝石閃長玢巖侵入體是尋找泥河式玢巖型鐵礦床的勘探靶區(qū)；張舒（2015）分析測試了泥河鐵礦主成礦期礦石礦物的稀土元素、硫同位素及鉛同位素，認為泥河鐵礦是由次火山巖體演化產(chǎn)生的含礦高溫熱液在閃長玢巖穹窿頂部，通過交代充填作用形成的玢巖型鐵硫礦床；張明明（2014）以鉆孔資料為依據(jù)，采用離散光滑插值法建立成礦巖體的精細三維模型，重現(xiàn)閃長玢巖體的形態(tài)特征；周濤發(fā)（2014）則系統(tǒng)總結了泥河鐵礦床的地質地球化學特征，提出泥河礦床的形成除了與輝石閃長玢巖有關外，還與三疊系膏鹽層有著重要的成因關系，并建立起泥河鐵礦床兩期成礦模式.有了以上深入細致的工作，泥河礦區(qū)結構與成礦控礦特征漸趨明朗，積累了大量寶貴的資料，這些信息資料是三維重力反演試驗成功的重要保證.

利用收集到的資料和數(shù)據(jù)如上述設計程序開展礦區(qū)三維重力反演，通過試驗效果檢驗方法的效力.資料包括：1∶10000地形地質圖， 1∶50000高精度重力觀測數(shù)據(jù)，74個鉆孔累計8.33萬米的鉆井資料（含69個鉆孔測井數(shù)據(jù)），882塊鉆孔巖芯物性測試數(shù)據(jù)，以及廬樅礦集區(qū)5條電磁剖面和9條地震剖面.

首先簡化地質單元，提取重要信息.將1∶10000地形地質圖簡化（見圖6），泥河礦區(qū)存在第四系（Q）砂礫巖、白堊統(tǒng)楊灣組（K1y）砂巖、白堊統(tǒng)雙廟組（K1sh）和白堊統(tǒng)磚橋組（K1zh）火山巖4套地層.地層產(chǎn)狀平緩，褶皺不發(fā)育，往深部地層產(chǎn)狀略有起伏變化.構造簡單，僅有火山巖地層的北西傾斜單斜構造和成礦期后的淺層斷裂.需要重點關注的是，泥河鐵礦體主要賦存在閃長玢巖體的頂部與磚橋組下段的火山碎屑巖中，閃長玢巖體的穹狀隆起部位是賦礦的有利構造（吳明安等，2011；趙文廣等，2011）.

接著備好數(shù)據(jù)，確定模型空間.重力資料采用1∶50000高精度觀測數(shù)據(jù)（源自安徽省地質調查院），需開展位場分離工作保留剩余重力值（劉彥，2012）.反演區(qū)域是在綜合礦區(qū)地層、巖體及整體結構并充分考慮到鉆孔資料可明確的1200 m巖性深度這些已知信息基礎上劃定的.建模2D剖面部署如圖6方框內線條所示，平面面積為5.6 km2（2.8 km×2.0 km），剖面垂直構造走向取北偏西40°方向，間距為100 m，共劃分28條平行剖面，區(qū)域重力場的背景密度設為2.62 g·cm-3.

圖6 泥河礦區(qū)地質簡化圖

運用貝葉斯反演程序建立初始模型.采用離散多面體建模需要依次對每條剖面開展重力反演工作.以礦區(qū)9線剖面為例，依據(jù)已有地質信息、鉆孔資料以及電磁地震等認識，勾畫出先驗模型如圖7a所示，該先驗模型基本界定地質體的坐標位置和埋深，這等同于固定地質體的體積.再計算各地質體的密度區(qū)間，采用劉彥等（2015）設計的貝葉斯反演程序.結合測試的882塊物性數(shù)據(jù)，給出各地層和巖體的密度區(qū)間范圍，利用程序快速算出初始模型中各地質體的密度（見表2），省卻常規(guī)離散體反演建模過程中反復調試密度的困頓，可大大提高建模效率.從地表起由上至下，第四系砂礫巖（Q）密度為2.10±0.03 g·cm-3，白堊統(tǒng)楊灣組（K1y）地層密度為2.30±0.05 g·cm-3，白堊統(tǒng)雙廟組上段（K1sh2）地層密度為2.55±0.04 g·cm-3，白堊統(tǒng)雙廟組下段（K1sh1）地層密度為2.57±0.03 g·cm-3，白堊統(tǒng)磚橋組（K1zh）地層密度為2.65±0.03 g·cm-3；其它巖體及礦體的密度范圍：閃長玢巖（δμ）巖體密度為2.90±0.04 g·cm-3，磁鐵礦（Mt）、黃鐵礦（Py）以及石膏礦（Ah）密度分別是3.48±0.09 g·cm-3、3.35±0.09 g·cm-3、2.93±0.07 g·cm-3.

圖7 泥河礦區(qū)9線剖面重力反演示意圖

表2 采用貝葉斯計算的泥河礦區(qū)初始模型中各地質體的密度

2.5D模型的構建.有了確定的密度區(qū)間再調整修改地質體的形態(tài)，即進行重力的反演迭代與初始模型的升級，直到模型曲線和實際重力曲線擬合完好，此時模型生成的重力值與實際重力觀測值方差滿足允許范圍，由此建立剖面的2.5D初始模型（如圖7b所示）.同樣方式，依次做好全部28條剖面的重力反演迭代與初始模型的升級工作.

圖8 安徽省泥河礦區(qū)三維重力反演模型

3D模型的生成與顯示.選用離散多面體反演建模專業(yè)軟件Encom ModelVision ProTM拼接擬合，完成整個泥河礦區(qū)的三維重力反演.反演最后得到的模型加載到三維可視化專業(yè)軟件Encom PA，實現(xiàn)多方位成像和結構解譯，最終模型如圖8所示.

圖8所示的3D模型顯示：泥河礦區(qū)內的主要礦體為磁鐵礦、黃鐵礦以及石膏礦；礦體整體呈北東向展布，延伸至東北部時礦體稍有抬升；礦區(qū)磁鐵礦、黃鐵礦較多，石膏礦相對較少；磁鐵礦主要位于礦區(qū)的西南部，呈透鏡狀；黃鐵礦主要集中在礦區(qū)的東北部，中部少量的石膏礦；黃鐵礦和石膏礦埋深相對于磁鐵礦較淺，礦體深度大致在地下600～1100 m范圍內.礦區(qū)東北部礦體主要是垂直方向上的兩層黃鐵礦，上層礦體體積較小，呈層狀分布，礦體在西南部較小，東北部較大，平均寬度約為245 m，埋深約為600 m，厚度約40 m；下層礦體體積較大，呈透鏡狀，埋深在800～1050 m，最大寬度約為680 m.三維重力反演結果還揭示出泥河礦床的控礦特征：礦區(qū)內地層由淺到深主要為第四系蓋層（Q）、白堊統(tǒng)楊灣組（K1y）、白堊統(tǒng)雙廟組（K1sh）和白堊統(tǒng)磚橋組（K1zh）；較淺層的黃鐵礦和石膏礦多位于磚橋組上段地層，體積較小，呈層狀似層狀分布；礦體主要集中在磚橋組下段侵入巖中，侵入巖主要為輝石閃長玢巖和脈巖；侵入體頂面形成隆起，礦區(qū)西南部隆起陡立，礦區(qū)東北部隆起寬緩.該模型刻畫展示的礦體和巖體的空間關系，與吳明安（2011）、趙文廣（2011）等運用地質學，范裕（2014）等運用成礦年代學，周濤發(fā)（2014）、張舒（2015）等運用地球化學，張明明（2014）等運用離散光滑插值法，張昆（2014）等運用電磁法等推論結果吻合一致，說明本方法開展三維重力反演是實用有效的.

5 結論與認識

（1） 本文發(fā)揮離散多面體反演技術和貝葉斯方法優(yōu)勢，通過理論分析、綜合研究，設計程序完成模型試驗和礦區(qū)實踐，改進了傳統(tǒng)的三維重力反演方式，得到如下認識：將貝葉斯方法引入重力反演初始模型的建立中，可以定量地利用已知先驗信息、觀測數(shù)據(jù)以及物理規(guī)律來確定模型參數(shù)具體數(shù)值，節(jié)省模型選擇時間，提高反演效率；最大似然函數(shù)算法應用于分類器的鑒別中改善了傳統(tǒng)樸素貝葉斯分類方法，采取模型降階技術固定所建模型中幾何體的形態(tài)或密度，達到在幾何體形態(tài)、密度和重力值五個參數(shù)中降低維數(shù)，可以減小高維不確定性和正演的計算量，有利于降低多解性，提高反演的準確度；任意形態(tài)幾何體模擬試驗以及泥河礦區(qū)反演實踐證明，這種基于貝葉斯算法的三維重力反演方式，能夠充分吸收先驗信息，改善反演精度，提高建模效率.

（2） 重點主要集中于改善反演初始模型的建立方式，三維重力反演的實現(xiàn)是通過本次設計的貝葉斯方法模塊與其它成熟反演建模軟件（如：Encom ModelVisionProTM）整合完成的.在貝葉斯反演計算中，可以通過固定密度、變化體積這種方式準確圈定出異常體的展布形態(tài)，這對于實際勘探生產(chǎn)將更具有應用價值.

致謝 本文得到中國地質大學（北京）薛濤老師及學生張永強的大力支持和幫助，審稿專家富有建設性的意見對于本文的完善起了很好的作用，稿件修改期間還得到中國地質科學院礦產(chǎn)資源研究所史大年研究員的指導以及與四川大學周永道博士的有益探討，所有這些一并致謝！

Akaike H. 1980. Likelihood and the Bayes procedure. Bernardo J M, DeGroot M H, Lindley D V , et al eds.Bayesian Statistics. Valencia:University Press, 141-166.

Carlin B P, Louis T A. 2000.Bayes and Empirical Bayes Methods for Data Analysis.2nd ed. New York:Chapman & Hall.

Dosso S E, Wilmut MJ. 2006. Data uncertainty estimation in matched-fieldgeoacousticinversion.IEEEJournalofOceanicEngineering, 31（2）: 470-479. Farquharson CG, Oldenburg DW. 2000. Automatic estimation of the trade-off parameter in nonlinear inverse problems using the GCVand L-curve criteria. 70th SEG Program Expanded Meeting.Expanded Abstracts, 265-268.

Fernández-Martínez J L, Fernández-Muiz Z, Pallero J L G, et al. 2013.From Bayes to Tarantola: New insights to understand uncertainty in inverse problems.JournalofAppliedGeophysics, 98（11）: 62-72.

Fullagar P K, Hughes N A, Paine J. 2000. Drilling-constrained 3D gravity interpretation.ExplorationGeophysics, 31（2）: 17-23.

Fullagar P K, Pears G A. 2007. Towards geologically realistic inversion. Milkereit Bed. Proceedings of the Fifth Decennial Conference on Mineral Exploration.ExpandedAbstracts, 445-460.

Gao R, Lu Z W, Liu J K, et al. 2010. A result of interpreting from deep seismic reflection profile: Revealing fine structure of the crust and tracing deep process of the mineralization in Luzong deposit area.ActaPetrologicaSinica（in Chinese）, 26（9）: 2543-2552.

Gao S. 1997. A Bayesian nonlinear inversion of seismic body-wave attenuation factors.Bull.Seism.Soc.Am., 87（4）: 961-970.

Goleby B R, Korsch R J, FominT, et al. 2002. Preliminary 3-D geological model of the Kalgoorlie region, YilgarnCraton, Western Australia, based on deep seismic-reflection and potential-field data.AustralianJournalofEarthSciences, 49（6）: 917-933.

Gouveia W P, Scales J A. 1998. Bayesian seismic waveform inversion: Parameter estimation and uncertainty analysis.J.Geophys.Res., 103（2）: 2759-2779.

Grandis H, Menvielle M,Roussignol M. 1999. Bayesian inversion with Markov chains：The magnetotelluricone-dimensional case.GeophysicalJournalInternational, 138（3）: 757-768.

Green P J. 1995. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.Biometrika, 82（4）: 711-732.

Haber E, Oldenburg D. 2000. A GCV based method for nonlinear ill-posed problems.ComputationalGeosciences, 4（1）: 41-63.

Hansen PC. 1992.Analysis of discrete ill-posed problems by means of the L-curve.SIAMReview, 34（4）: 561-580.

Jessell M. 2001. Three-dimensional geological modelling of potential-field data.Computers&Geosciences, 27（4）: 455-465.

Khodja M R, PrangeM D, DjikpesseH A. 2010. Guided Bayesian optimal experimental design.InverseProblems, 26（5）: 055008.

Lasses M, Siltanen S. 2004. Can one use total variation prior for edge-preserving Bayesian inversion?.InverseProblems, 20（5）: 1537-1563.

Li Y G, Oldenburg D W. 1998. 3-D inversion of gravity data.Geophysics, 63（1）: 109-119.

Liu Y, Yan J Y, Wu M A, et al. 2012. Exploring deep concealed ore bodies based on gravity anomaly separation methods: A case study of the Nihe iron deposit.ChineseJ.Geophys.（in Chinese）, 55（12）: 4181-4193, doi: 10.6038/j.issn.0001-5733.2012.12.030.

Lü Q T, Han L G, Yan J Y, et al. 2010. Seismic imaging of volcanic hydrothermal iron-sulfur deposits and its hosting structure in Luzong ore district.ActaPetrologicaSinica（in Chinese）, 26（9）: 2598-2612.

Lü Q T, Qi G, Yan J Y. 2013. 3D geologic model of Shizishan ore field constrained by gravity and magnetic interactive modeling: A case history.Geophysics, 78（1）: B25-B35.

Lü Q T, Liu Z D, Tang J T, et al. 2014. Upper crustal structure and deformation of Lu-Zong ore district: Constraints from integrated geophysical data.ActaGeologicalSinica（in Chinese）, 88（4）: 447-465.

Lü Q T, Liu Z D, Yan J Y, et al. 2015.Crustal-scale structure and deformation of Lu-Zong ore district: Joint interpretation from integrated geophysical data.Interpretation, 3（2）: SL39-SL61.

Malehmir A, Tryggvason A, Juhlin C, et al. 2006. Seismic imaging and potential field modelling to delineate structures hosting VHMS deposits in the Skellefte Ore District, Northern Sweden.Tectonophysics, 426（3-4）: 319-334.

Malehmir A, ThunehedH,Tryggvason A. 2009. The paleoproterozoic Kristineberg mining area, northern Sweden: Results from integrated 3D geophysical and geologic modeling, and implications for targeting ore deposits.Geophysics, 74（1）: B9-B22.

Martin L, PerronG, MassonM. 2007. Discovery from 3D visualization and quantitative modelling. Milkereit Bed.Proceedings of the Fifth International Conference on Mineral Exploration.ExpandedAbstracts, 543-550.

McGaughey J. 2007. Geological models, rock properties, and the 3D inversion of geophysical data. MilkereitBed.Proceedings of Exploration: Fifth Decennial International Conference on Mineral Exploration. Expanded Abstracts, 473-483.

McInerney P, Goldberg A, Calcagno P, et al. 2007. Improved 3D geology modelling using an implicit function interpolator and forward modelling of potential field data. MilkereitBed.Proceedings of Exploration 2007: Fifth Decennial International Conference on Mineral Exploration. Expanded Abstracts, 919-922.

Mitsuhata Y. 2004. Adjustment of regularization in ill-posed linear inverse problems by the empirical Bayes approach.GeophysicalProspecting, 52（3）: 213-239.

Oldenburg D W, Li Y G, Ellis R G. 1997.Inversion of geophysical data over a copper gold porphyry deposit: A case history for Mt. Milligan.Geophysics, 62（5）: 1419-1431.

Oldenburg D W, Pratt D A. 2007.Geophysical inversion for mineral exploration: A decade of progress in theory and practice. MilkereitBed.Proceedings of Exploration Fifth Decennial International Conference on Mineral Exploration. Expanded Abstracts, 61-95.

Qi G, LüQT, Yan J Y, et al. 2012.Geologic constrained 3D gravity and magnetic modeling of Nihe deposit-Acasestudy.ChineseJ.Geophys. （in Chinese）, 55（12）: 4194-4206, doi: 10.6038/j.issn.0001-5733.2012.12.031.

Roy B, Clowes R M. 2000. Seismic and potential-field imaging of the Guichon Creek batholith, Brithsh Columia, Canada, to delineate structures hosting porphyry copper deposits.Geophysics, 65（5）: 1418-1434. Scales J A, Snieder R. 1997.To Bayes or not to Bayes?Geophysics, 62（4）: 1045-1046.

Stuart A M. 2010. Inverse problems: A Bayesian perspective.ActaNumerica, 19: 451-559.

Tarantola A. 2005. Inverse Problem Theoryand Methods for Model Parameter Estimation.Philadelphia, PA:SIAM. Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-posed Problems.Washington, DC: John Wiley & Sons, Inc.Wahba G. 1990. Spline Models for Observational Data.Philadelphia, PA:Society of Industrial and Applied Mathematics.Wang G W, ZhangS T, Yan C H, et al. 2011. Mineral potential targeting and resource assessment based on 3D geological modeling in Luanchuan region,China.Computers&Geosciences, 37（12）1976-1988.

Welford K J, HallJ. 2007. Crustal structure of the Newfoundland rifted continental margin from constrained 3-D gravity inversion.GeophysicalJournalInternational, 171（2）: 890-908.

Williams N C, Lane R, Lyons P. 2004.Regional constrained 3D inversion of potential field data from the Olympic Cu-Au province, South Australia.ASEGExtendedAbstracts,（1）: 1-4.

Williams N C. 2008.Geologically-constrained UBC-GIF gravity and magnetic inversions with examples from the Agnew-Wiluna greenstone belt, Western Australia[Ph. D. thesis]. Vancouver,Canada: The University of British Columbia.Wu M A, Wang Q S, Zheng G W, et al. 2011. Discovery of the Nihe iron deposit in Lujiang, Anhui, and its exploration significance.ActaGeologicaSinica（in Chinese）, 85（5）: 802-809.

Yang W C. 1997.Theory and Methods in Geophysical Inversion （in Chinese）.Beijing: Geological Publish House.

Yao C L, Hao T Y,Guan Z N. 2002.Restrictions in gravity and magnetic inversions and technical strategy of 3D properties inversion.Geophysical&GeochemicalExploration（in Chinese）, 26（4）: 253-257.

Yao C L, Zheng Y M,Zhang Y W. 2007. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces.ChineseJournalofGeophysics（in Chinese）, 50（5）: 1576-1583, doi: 10.3321/j.issn:0001-5733.2007.05.035.

Zeng H L. 2005. Gravity Field and Gravity Exploration （in Chinese）.Beijing: Geological Publish House.

Zhang K, Yan J Y, Lü QT, et al.2014.The electromagnetic exploration experimentation of Nihe porphyry iron ore in anhui.ActaGeologicaSinica（in Chinese）, 88（4）:498-506.

Zhang M M, Zhou T F, Yuan F, et al. 2014.Three-dimensional morphological analysis method for mineralization related intrusion and prospecting indicators of nihe iron deposit in luzong basin.ActaGeologicaSinica（in Chinese）, 88（4）：574-583.

Zhang S, Wu M A, Zhao W G, et al. 2014.Geochemistry characteristics of Nihe iron deposit in Lujiang, Anhui Province and their constrains to ore genesis.ActaPetrologicaSinica（in Chinese）, 30（5）:1382-1396.

Zhao W G, Wu M A, Zhang Y Y, et al. 2011. Geological characteristics and genesis of the Nihe Fe-S deposit, Lujiangcountry, Anhui province.ActaGeologicaSinica（in Chinese）, 85（5）: 789-801.

Zhou T F, Fan Y, Yuan F, et al.2014.The metallogenic model of nihe iron deposit in Lu-Zong basin and genetic relationship between gypsum-salt layer and deposit.ActaGeologicaSinica（in Chinese）, 88（4）:562-573.

附中文參考文獻

范裕,劉一男,周濤發(fā)等.2014.安徽廬樅盆地泥河鐵礦床年代學研究及其意義.巖石學報，30（5）:1369-1381.

高銳,盧占武,劉金凱等．2010.廬樅金屬礦集區(qū)深地震反射剖面解釋結果-揭露地殼精細結構，追蹤成礦深部過程. 巖石學報,26（9）:2543-2552.

劉彥, 嚴加永, 吳明安等. 2012. 基于重力異常分離方法尋找深部隱伏鐵礦—以安徽泥河鐵礦為例. 地球物理學報, 55（12）: 4181-4193, doi: 10.6038/j.issn.0001-5733.2012.12.030.

劉彥,張永強等．2015. 基于貝葉斯反演理論的重力建模程序軟件．中華人民共和國國家版權局：軟著登字第1044862號．

呂慶田,韓立國,嚴加永等. 2010. 廬樅礦集區(qū)火山氣液型鐵、硫礦床及控礦構造的反射地震成像．巖石學報，26（9）:2598-2612.

呂慶田,劉振東,湯井田等. 2014. 廬樅礦集區(qū)上地殼結構與變形：綜合地球物理探測結果. 地質學報, 88（4）: 447-465.

祁光, 呂慶田, 嚴加永等. 2012. 先驗地質信息約束下的三維重磁反演建模研究—以安徽泥河鐵礦為例. 地球物理學報, 55（12）: 4194-4206, doi: 10.6038/j.issn.0001-5733.2012.12.031. 吳明安, 汪青松, 鄭光文等. 2011. 安徽廬江泥河鐵礦的發(fā)現(xiàn)及意義. 地質學報, 85（5）: 802-809.

楊文采. 1997. 地球物理反演的理論與方法. 北京: 地質出版社.

姚長利, 郝天珧, 管志寧. 2002. 重磁反演約束條件及三維物性反演技術策略. 物探與化探, 26（4）: 253-257.

姚長利, 鄭元滿, 張聿文. 2007. 重磁異常三維物性反演隨機子域法方法技術. 地球物理學報, 50（5）: 1576-1583, doi: 10.3321/j.issn:0001-5733.2007.05.035.

曾華霖. 2005. 重力場與重力勘探. 北京: 地質出版社.

張昆,嚴加永,呂慶田等.2014.安徽泥河玢巖鐵礦電磁法探測試驗.地質學報,88（4）:498-506.

張明明,周濤發(fā),袁鋒等.2014.廬樅盆地泥河鐵礦床成礦巖體三維形態(tài)學分析及找礦指標研究.地質學報,88（4）:574-583.

張舒,吳明安,趙文廣等.2014.安徽廬江泥河鐵礦礦床地球化學特征及其對成因的制約.巖石學報,30（5）:1382-1396.

趙文廣, 吳明安, 張宜勇等. 2011. 安徽省廬江縣泥河鐵硫礦床地質特征及成因初步分析. 地質學報,85（5）:789-801.

周濤發(fā),范裕,袁峰等.2014.安徽廬樅盆地泥河鐵礦床與膏鹽層的成因聯(lián)系及礦床成礦模式.地質學報,88（4）:562-573.

（本文編輯 劉少華）

3D gravity inversion based on Bayesian method with model order reduction

LIU Yan1，2, Lü Qing-Tian2，3*, LI Xiao-Bin4, QI Guang1，2, ZHAO Jin-Hua1，2, YAN Jia-Yong1，2, DENG Zhen1，2

1InstituteofMineralResourcesChineseAcademyofGeologicalSciences,MLRKeyLaboratoryofMetallogenyandMineralAssessment,Beijing100037，China2ChinaDeepExplorationCenter—SinoProbeCenter,ChineseAcademyofGeologicalSciences,Beijing100037,China3InstituteofGeophysicalandGeochemicalExplorationChineseAcademyofGeologicalSciences,HebeiLangfang065000，China4HenanPolytechnicUniversity,HenanJiaozuo454000，China

3D gravity inversion is an important way for geologist to understand the deep structure of earth. As the key part of 3D gravity inversion, inversion unit includes discrete simulation and voxels. In fact discrete simulation inversion method is welcome, which is easy to absorb prior geological information, and can make the theory gravity field fitting observation gravity field. This paper summarizes the problems of discrete simulation inversion method in 3D gravity inversion and improves its modeling process. In order to build an initial accurate model, we develop the feature of Bayesian inversion theory. The algorithm of maximum-likelihood function is used and hidden Markov chain is added to build a naive Bayes classifier. Based on probability theory, Bayesian method can relate measure data into the priori information, and constrain posterior parameters through the prior information model. While the priori information is fully utilized, output results is abundant and reliable. In the process of Bayesian method, model order reduction techniques are taken to reduce parameter dimension with fixing the model body geometry shape or density in the five parameters of geometry shape （x,y,z）, density （σ） and gravity （g）.Thereby high-dimensional uncertainty and the forward calculation are also decreased. As density and distribution of geological body can be accurately obtained, gravity model structure is fine also. To test the method of 3D gravity inversion based on those Bayes algorithm, simulation experiments of unit sphere and any geometry with normal distribution noise carry on, and 3D gravity inversion experiment of Nihe mining area in Anhui province have performed. The density and gravity value of inversion calculation is very close to the reality. It is proved that the Bayesian method is practical and effective.

3D gravity inversion； Discrete body simulation； Bayesian method； Maximum-likelihood function； Model order reduction

“十二五”國家科技支撐計劃課題（2011BAB04B01）,國家深部探測專項（SinoProbe-03）,國家自然科學基金項目（41304100）,中國地質調查項目（12120114019401）,中央級公益性科研院所基本科研業(yè)務費專項資金項目（K1317）聯(lián)合資助.

劉彥，女，1975年生，高級工程師，理學博士，主要從事地球物理勘探技術和深部探測工作.E-mail：liuy@cags.ac.cn

*通訊作者 呂慶田，男，1964年生，研究員，博士生導師.主要從事深部探測和金屬礦勘查技術方法研究.E-mail：lqt@cags.ac.cn

10.6038/cjg20151233.

10.6038/cjg20151233

P631

2015-02-12，2015-10-10收修定稿

劉彥， 呂慶田， 李曉斌等. 2015. 基于模型降階的貝葉斯方法在三維重力反演中的實踐.地球物理學報，58（12）:4727-4739,

Liu Y, Lü Q T, Li X B, et al. 2015. 3D gravity inversion based on Bayesian method with model order reduction.ChineseJ.Geophys. （in Chinese）,58（12）:4727-4739,doi:10.6038/cjg20151233.