劉彥， 呂慶田， 李曉斌， 祁光，趙金花， 嚴(yán)加永， 鄧震

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

?

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

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

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

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

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

1 引言

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

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

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

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

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

2.1 貝葉斯方法原理

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

（1）

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

圖1 貝葉斯方法原理

（2）

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

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

L（m）=

（3）

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

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

（4）

其中

（5）

A是一個(gè)比例因子.

（6）

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

（7）

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

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

（8）

（10）

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

（12）

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

2.2 三維重力反演過程

有了以上算法基礎(chǔ)，以離散多面體為反演單元設(shè)計(jì)三維重力反演，包括：先驗(yàn)?zāi)Ｐ偷拇_定、初始模型的建立、重力反演迭代與初始模型的升級、模型的三維顯示，如圖2所示.詳情以下展開.

圖2 三維重力反演流程圖

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

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

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

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

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

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

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

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

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

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

圖3 貝葉斯算法程序

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

2.2.3 采用離散體模擬方法將初始地質(zhì)模型構(gòu)建成2.5D地質(zhì)模型

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

2.2.4 將2.5D地質(zhì)模型拼接擬合成3D地質(zhì)模型

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

2.2.5 對3D地質(zhì)模型進(jìn)行可視化和結(jié)構(gòu)解譯處理

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

3 模擬實(shí)驗(yàn)

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

（14）

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

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

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

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

表1 單位球體模擬貝葉斯計(jì)算結(jié)果Table 1 Bayesian calculation results of unit sphere simulation

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

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

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

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

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

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

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

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

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

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

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

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

運(yùn)用貝葉斯反演程序建立初始模型.采用離散多面體建模需要依次對每條剖面開展重力反演工作.以礦區(qū)9線剖面為例，依據(jù)已有地質(zhì)信息、鉆孔資料以及電磁地震等認(rèn)識，勾畫出先驗(yàn)?zāi)Ｐ腿鐖D7a所示，該先驗(yàn)?zāi)Ｐ突窘缍ǖ刭|(zhì)體的坐標(biāo)位置和埋深，這等同于固定地質(zhì)體的體積.再計(jì)算各地質(zhì)體的密度區(qū)間，采用劉彥等（2015）設(shè)計(jì)的貝葉斯反演程序.結(jié)合測試的882塊物性數(shù)據(jù)，給出各地層和巖體的密度區(qū)間范圍，利用程序快速算出初始模型中各地質(zhì)體的密度（見表2），省卻常規(guī)離散體反演建模過程中反復(fù)調(diào)試密度的困頓，可大大提高建模效率.從地表起由上至下，第四系砂礫巖（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 采用貝葉斯計(jì)算的泥河礦區(qū)初始模型中各地質(zhì)體的密度

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

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

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

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

5 結(jié)論與認(rèn)識

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

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

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

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.

附中文參考文獻(xiàn)

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

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

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

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

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

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

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

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

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

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

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

張昆,嚴(yán)加永,呂慶田等.2014.安徽泥河玢巖鐵礦電磁法探測試驗(yàn).地質(zhì)學(xué)報(bào),88（4）:498-506.

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

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

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

周濤發(fā),范裕,袁峰等.2014.安徽廬樅盆地泥河鐵礦床與膏鹽層的成因聯(lián)系及礦床成礦模式.地質(zhì)學(xué)報(bào),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

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

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

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

10.6038/cjg20151233.

10.6038/cjg20151233

P631

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

劉彥， 呂慶田， 李曉斌等. 2015. 基于模型降階的貝葉斯方法在三維重力反演中的實(shí)踐.地球物理學(xué)報(bào)，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.