基于密度模型稀疏表征的重力反演方法

于會(huì)臻， 王金鐸， 王千軍

中國(guó)石化勝利油田分公司勘探開(kāi)發(fā)研究院， 山東東營(yíng) 257000

0 引言

重力反演的目的是由地表觀(guān)測(cè)的重力數(shù)據(jù)恢復(fù)出未知地下空間的密度分布，是礦產(chǎn)資源勘查部署的重要手段(Paterson and Reeves，1985；Oldenburg et al.，1997；管志寧等，1998； Portniaguine and Zhdanov，1999；姚長(zhǎng)利等，2003；孟小紅等，2012).重力反演首先需要對(duì)三維地下半空間進(jìn)行離散化剖分，通常采用的方式是六面體網(wǎng)格方式(Li and Oldenburg，1998).但與地震、電法等地球物理技術(shù)相比，隨著深度的增加，相同大小的密度體所產(chǎn)生的重力異常幅值及頻率衰減更快，導(dǎo)致重力正演核函數(shù)矩陣條件數(shù)較大，同時(shí)受到觀(guān)測(cè)噪聲及異常處理精度的影響，反演結(jié)果分辨率較低、多解性更強(qiáng).

在保證重力異常數(shù)據(jù)擬合的前提下，合理地施加模型約束項(xiàng)是降低重力反演多解性問(wèn)題、獲得可靠密度反演結(jié)果的有效手段.不同的模型約束方法采用了不同的模型假設(shè)，以滿(mǎn)足不同勘探目標(biāo)解釋工作的需求.Li和Oldenburg(1996)利用深度加權(quán)矩陣和L2范數(shù)約束以降低反演結(jié)果的趨膚效應(yīng)；Last和Kubik(1983)引入了基于最小體積約束，Portniaguine和Zhdanov(1999)、Zhdanov(2002)、Zhdanov等(2004)引入最小支撐約束來(lái)獲得具有聚焦特征的反演結(jié)果；Bertete-Aguirre等(2002)引入了最小梯度支撐和總變差(TV)約束來(lái)保證銳化模型反演結(jié)果的邊界梯度；秦朋波和黃大年(2016)、高秀鶴和黃大年(2017)利用聚焦反演方法對(duì)重力及重力梯度數(shù)據(jù)進(jìn)行聯(lián)合反演，通過(guò)融合多維度觀(guān)測(cè)信息提高反演結(jié)果可靠性；Farquharson和Oldenburg(1998)、Farquharson(2008)、Vatankhah等(2017)利用L1范數(shù)來(lái)提高反演分辨率；Sun和Li(2014)、李澤林等(2019)利用Lp范數(shù)來(lái)改善重力反演結(jié)果；彭國(guó)民等(2018)采用柯西分布約束來(lái)保證模型反演結(jié)果的稀疏性.相比來(lái)說(shuō)，具有稀疏特征或稀疏邊界特征的模型約束方法能得到更高的解釋分辨率，在實(shí)際勘探應(yīng)用中更具吸引力.

對(duì)于稀疏或聚焦反演方法來(lái)說(shuō)，獲得可靠的結(jié)果需具備較強(qiáng)的前提條件：一是網(wǎng)格剖分與地質(zhì)體間的匹配性，二是重力剩余異常求取的準(zhǔn)確性.在實(shí)際應(yīng)用中，上述條件很難得到滿(mǎn)足，首先地質(zhì)體發(fā)育模式復(fù)雜、形態(tài)多樣，準(zhǔn)確的網(wǎng)格剖分方案難以確定；其次總的重力異常為地下半空間密度體產(chǎn)生重力異常的疊加，背景異常難以準(zhǔn)確剝離，獲得與聚焦密度體對(duì)應(yīng)的剩余重力異常存在極大的不確定性.這些問(wèn)題都會(huì)導(dǎo)致聚焦或稀疏約束反演密度值和空間分布與真實(shí)情況間存在較大的偏差.為此，通?？梢霚y(cè)井資料或添加光滑、能量最小等多種約束方式(劉展等，2011；朱自強(qiáng)等，2014；Utsugi，2019)來(lái)增強(qiáng)反演結(jié)果的可靠性，獲得介于光滑與聚焦的密度反演結(jié)果.但往往研究區(qū)測(cè)井資料有限，而且光滑與聚焦(稀疏)的聯(lián)合約束也難以描述復(fù)雜的密度空間分布.可以看出，要想滿(mǎn)足高精度重力反演應(yīng)用需求，還需開(kāi)展更深入的模型約束方法研究.

本文提出了一種新的基于密度模型稀疏表征的重力反演方法，引入了可描述典型地質(zhì)體發(fā)育模式及幾何特征的特征模型，并介紹了利用特征模型構(gòu)建模型特征矩陣的流程，在此基礎(chǔ)上，重新推導(dǎo)了重力反演求解方程，將直接求解密度轉(zhuǎn)換為求解與模型特征矩陣對(duì)應(yīng)的分解系數(shù)稀疏求解問(wèn)題，以保證利用最少的、最具有代表性的特征模型組合構(gòu)成期望反演得到的密度模型，從而獲得分辨率更高、可靠性更強(qiáng)的三維密度反演結(jié)果.

1 方法原理

1.1 重力反演方法

首先簡(jiǎn)要回顧一下重力反演方法.在笛卡爾坐標(biāo)系下，將地下半空間剖分為三維離散網(wǎng)格，則重力正演可表達(dá)為如下線(xiàn)性方程：

AM=d,

(1)

其中，M∈Rl×1為待求解地下剖分網(wǎng)格內(nèi)的密度模型向量，A∈Rn×l是各個(gè)剖分網(wǎng)格在單位密度下的重力正演核函數(shù)矩陣，d∈Rn×1為重力異常數(shù)據(jù)；n、l分別為重力觀(guān)測(cè)數(shù)據(jù)及反演區(qū)域網(wǎng)格剖分的個(gè)數(shù).

重力反演是正演的逆過(guò)程，為了減少反演多解性，通常會(huì)施加定量的模型約束，可分為線(xiàn)性和非線(xiàn)性模型約束算子兩大類(lèi).前者中最為常用的包括用于保證密度模型能量最小的單位矩陣算子和保證密度模型光滑的拉普拉斯算子等；后者則在迭代反演過(guò)程中因模型的改變而改變，如最小支撐、最小梯度支撐等非線(xiàn)性算子.上述兩類(lèi)模型約束可歸納至統(tǒng)一的重力反演目標(biāo)函數(shù)中，形式如下：

(2)

(3)

第k次迭代的密度模型Mk滿(mǎn)足的目標(biāo)函數(shù)形式如下：

(4)

求解該目標(biāo)函數(shù)的等價(jià)形式為：

(5)

其中，

ο∈Rl×1為一個(gè)零向量.

聚焦反演的目標(biāo)是將密度反演結(jié)果集中在部分網(wǎng)格中，追求的是利用最少的密度網(wǎng)格模型來(lái)擬合重力數(shù)據(jù).但異常分離不準(zhǔn)確時(shí)，會(huì)有一定成分的區(qū)域背景異常干擾，而這部分異常對(duì)應(yīng)的密度體并不符合聚焦算法的應(yīng)用前提.即使異常分離結(jié)果較為準(zhǔn)確，網(wǎng)格剖分過(guò)小或過(guò)大，也難以保證反演結(jié)果的可靠性.

1.2 模型特征矩陣構(gòu)建流程

針對(duì)現(xiàn)有方法存在的問(wèn)題，本文利用模型特征矩陣和分解系數(shù)來(lái)對(duì)密度模型進(jìn)行稀疏表征.下面首先通過(guò)一維模型來(lái)闡述密度模型稀疏表征的含義.如圖1所示，該一維密度模型可分解為三個(gè)特征模型，一是數(shù)值為2 g·cm-3的區(qū)域背景密度模型，二是數(shù)值為0.5 g·cm-3的密度體B1剩余密度模型，三是數(shù)值為0.2 g·cm-3的密度體B2剩余密度模型.

圖1 一維密度模型分解示意圖Fig.1 Diagram of 1D density model decomposition

該一維密度模型M的分解過(guò)程可表達(dá)為模型特征矩陣和分解系數(shù)向量的矩陣相乘：

(6)

其中，將D∈Rl×q稱(chēng)為模型特征矩陣；Γ∈Rq×1稱(chēng)為模型特征矩陣對(duì)應(yīng)的分解系數(shù)向量；q為分解系數(shù)向量的個(gè)數(shù).

矩陣D包含了三個(gè)子矩陣，即

D=[DB1,DB2,Dbg].

(7)

根據(jù)公式(6)的矩陣形式與圖1所示三個(gè)特征模型的對(duì)應(yīng)關(guān)系可看出，矩陣DB1每一列包含了塊體B1的幾何形態(tài)特征，矩陣DB2每一列包含了塊體B2的幾何形態(tài)特征，矩陣Dbg包含了背景信息.為了保證利用最少的特征模型來(lái)組成密度模型，分解向量Γ應(yīng)是稀疏的，其非零點(diǎn)處數(shù)值的大小和位置表示特征模型的密度值大小和空間賦存位置.

模型特征矩陣的構(gòu)建原則來(lái)自于先驗(yàn)地質(zhì)假設(shè)，為更加清晰的描述模型特征矩陣構(gòu)建過(guò)程，下面以二維模型(假設(shè)待反演區(qū)域沿X、Z方向剖分為9×7個(gè)網(wǎng)格)為例進(jìn)行具體說(shuō)明：

(1)根據(jù)已知資料分析，假設(shè)地質(zhì)體發(fā)育模式包含以下兩類(lèi)矩形塊體特征模型，具體參數(shù)為：特征模型Ⅰ(圖2a)：寬W1m、高H1m；特征模型Ⅱ(圖2b)：寬W2m、高H2m；

圖2 構(gòu)建模型特征矩陣的特征模型(a) 特征模型Ⅰ； (b) 特征模型Ⅱ.Fig.2 Feature models for building model feature matrix(a) Feature model Ⅰ; (b) Feature model Ⅱ.

(2)首先計(jì)算特征模型Ⅰ對(duì)應(yīng)的索引向量Index_I.將特征模型I的中心點(diǎn)坐標(biāo)置于二維剖分網(wǎng)格左上角第一個(gè)網(wǎng)格(圖3a)，以此為起始點(diǎn)進(jìn)行第1次搜索，計(jì)算特征模型I所占反演剖分的索引，索引Index_I_1位置的參數(shù)記為整數(shù)1，其余為0(圖3b)；向右平移一個(gè)網(wǎng)格的距離至左上角第2個(gè)網(wǎng)格，計(jì)算特征模型I所占反演網(wǎng)格的索引，并將其記錄為Index_I_2，索引Index_I_2位置的參數(shù)記為整數(shù)1，其余為0；以此類(lèi)推，計(jì)算完第一層之后，從第二層的左側(cè)第一個(gè)網(wǎng)格點(diǎn)為起始點(diǎn)計(jì)算特征模型I的所占網(wǎng)格位置，并進(jìn)行索引標(biāo)記；之后，逐層計(jì)算，設(shè)平移至第22個(gè)網(wǎng)格(X=4、Z=3)(圖3c)進(jìn)行第22次搜索，索引Index_I_22標(biāo)記如圖(圖3d)所示；按此步驟直到遍歷所有網(wǎng)格；

圖3 特征模型Ⅰ對(duì)應(yīng)的模型特征矩陣構(gòu)建過(guò)程(a) 第1次搜索位置; (b) 第1次的索引; (c) 第22次搜索位置; (d) 第22次的索引.Fig.3 Building process of model feature matrix corresponding to feature model Ⅰ(a) The first search location; (b) The first index; (c) The 22nd search location; (d) The 22nd index.

(3)將步驟(2)得到的索引Index_I_k(k=1,…,n)，n為剖分網(wǎng)格個(gè)數(shù)，此次為63，展開(kāi)為列向量并進(jìn)行組合得到特征模型I對(duì)應(yīng)的模型特征矩陣DB1；

(4)特征模型Ⅱ的模型特征矩陣構(gòu)建.重新執(zhí)行第(2)、(3)步，計(jì)算特征模型Ⅱ下的索引向量，并將其進(jìn)行組合得到特征模型Ⅱ?qū)?yīng)的特征矩陣DB2；

(5)背景模型特征矩陣.類(lèi)似特征模型Ⅰ、Ⅱ的方式，選擇具有趨勢(shì)特征的模型作為特征模型，將其生成背景模型特征矩陣Dbg；

(6)完成上述步驟，便可生成二維反演所需的模型特征矩陣D=[DB1,DB2,Dbg].

上述構(gòu)建流程同樣適用于反演不等間隔網(wǎng)格剖分的情況.三維模型特征矩陣的構(gòu)建流程與二維類(lèi)似，在后文三維模型反演實(shí)驗(yàn)中將舉例說(shuō)明，在此不做過(guò)多贅述.

模型特征矩陣所采用的特征模型不限于規(guī)則的長(zhǎng)方體等幾何形態(tài)，也適用于傾斜(巖脈)、球體(孤立火山巖體)等復(fù)雜地質(zhì)模式特征.在實(shí)際應(yīng)用過(guò)程中，可根據(jù)研究區(qū)的先驗(yàn)地質(zhì)認(rèn)識(shí)或其他物探資料來(lái)構(gòu)建模型特征矩陣.在缺少先驗(yàn)地質(zhì)認(rèn)識(shí)的地區(qū)，也可利用重力異常分析來(lái)計(jì)算特征模型的水平和垂向上的尺寸大小及傾角信息.

相比常規(guī)反演，本文將模型特征矩陣D作用于重力正演核函數(shù)A的過(guò)程具有兩大優(yōu)勢(shì)，一是相當(dāng)于對(duì)原始剖分網(wǎng)格按照特征模型的樣式進(jìn)行了重新組合，在組合后的多套網(wǎng)格中分別進(jìn)行分解系數(shù)的稀疏求解，更加符合期望的稀疏假設(shè)條件；二是特征模型可以更高效的將期望滿(mǎn)足的地質(zhì)模式的定性假設(shè)定量化，使得重力反演的過(guò)程更具傾向性，無(wú)論是塊狀、層狀還是傾斜狀等復(fù)雜形態(tài)的地質(zhì)體發(fā)育特征都可以被引入進(jìn)來(lái).

1.3 目標(biāo)函數(shù)及求解算法

將公式(6)代入公式(2)，將原重力反演問(wèn)題的目標(biāo)函數(shù)寫(xiě)為：

(8)

其中L?！蔙q×q為關(guān)于分解系數(shù)Γ線(xiàn)性模型約束算子(取單位矩陣)；φ(Γ)為關(guān)于Γ的非線(xiàn)性模型約束函數(shù).可見(jiàn)，基于重新構(gòu)建的重力反演目標(biāo)函數(shù)(8)，求解目標(biāo)變?yōu)榱朔纸庀禂?shù)Γ.如前文所述，為保證分解系數(shù)Γ的稀疏性，可對(duì)其施加L1范數(shù)約束項(xiàng).則重力反演目標(biāo)函數(shù)可寫(xiě)為：

(9)

其中，‖·‖1代表L1范數(shù)稀疏約束.

稀疏求解問(wèn)題在目前信號(hào)處理領(lǐng)域得到了廣泛的應(yīng)用，求解算法包括基追蹤(Chen et al.，2001)、交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)(Boyd et al., 2010)、正交匹配追蹤(orthogonal matching pursuit)(Sahoo and Makur，2015)等.為了獲得更加高的稀疏度、更加準(zhǔn)確的分解系數(shù)，本文采用了Majorization-Minimization(優(yōu)化最小化，簡(jiǎn)稱(chēng)MM)優(yōu)化框架來(lái)加速稀疏分解系數(shù)Γ的求解過(guò)程.MM框架(Figueiredo et al.，2007；Selesnick and Bayram，2014；Nguyen，2017)是用來(lái)求解非凸函數(shù)的一種性能良好的優(yōu)化算法框架.MM算法通過(guò)將原來(lái)的復(fù)雜優(yōu)化問(wèn)題分解為一系列簡(jiǎn)單優(yōu)化問(wèn)題，對(duì)于分解系數(shù)進(jìn)行稀疏求解(具體推導(dǎo)過(guò)程見(jiàn)附錄A)，其第k次迭代的Γk公式如下：

(10)

WΓk-1為第k-1次迭代求解結(jié)果Γk-1對(duì)應(yīng)的對(duì)角加權(quán)矩陣；φ′(Γk-1)為模型約束函數(shù)φ(Γ)關(guān)于Γk-1的一階偏導(dǎo)，其計(jì)算公式為：

(11)

為避免反演出現(xiàn)嚴(yán)重的趨膚現(xiàn)象，對(duì)公式(10)中的重力異常正演核函數(shù)進(jìn)行深度加權(quán)，則：

(12)

其中，

Wdepth∈Rl×l為深度加權(quán)矩陣，N為深度加權(quán)參數(shù)，υ∈Rq×1為與Γ大小相等的零向量.

由于期望獲得Γ是稀疏的，即Γ中必然會(huì)出現(xiàn)一定數(shù)量的0值，而WΓk-1中的Γk-1位于分母位置，易出現(xiàn)奇異值，可對(duì)公式(10)中的求逆項(xiàng)進(jìn)行重新推導(dǎo)得到如下公式.

(13)

其中，

綜合前文所述的模型特征矩陣構(gòu)建及分解系數(shù)求解過(guò)程，基于密度模型稀疏表征的重力反演方法主要包含以下步驟：

(1)根據(jù)實(shí)際需求剖分原始反演網(wǎng)格，計(jì)算重力正演核函數(shù)矩陣A；

(2)根據(jù)工區(qū)實(shí)際情況、重力異常形態(tài)等信息設(shè)置特征模型(可包括塊體、傾斜狀巖脈、球體等各類(lèi)典型地質(zhì)體)，根據(jù)1.2節(jié)介紹的過(guò)程構(gòu)建模型特征矩陣D；

(3)設(shè)置正則化參數(shù)λ1、λ2及深度加權(quán)參數(shù)N(一般N=4)，設(shè)置最大迭代次數(shù)IterMax；

(4)根據(jù)模型特征矩陣列向量個(gè)數(shù)設(shè)置稀疏分解系數(shù)初始模型向量Γk-1=ΓInit(k=1)，注意ΓInit不要含有0值；

(6)利用公式(10)和(13)計(jì)算第k次的分解系數(shù)Γk；

(7)利用公式(6)得到第k次迭代對(duì)應(yīng)的密度模型Mk=WdepthDΓk；

(8)獲得公式(9)的計(jì)算結(jié)果，若滿(mǎn)足目標(biāo)函數(shù)收斂條件則執(zhí)行步驟(9)；若不滿(mǎn)足，重復(fù)步驟(5)—(7)，修改反演參數(shù)；

(9)評(píng)價(jià)(8)獲得的反演結(jié)果，若符合地質(zhì)認(rèn)識(shí)，則執(zhí)行步驟(10)；若不滿(mǎn)足，重復(fù)(2)—(8)，修改特征模型，重新構(gòu)建模型特征矩陣并求取反演結(jié)果；

(10)最終的密度反演結(jié)果輸出.

2 理論模型實(shí)驗(yàn)分析

為了驗(yàn)證本文所提重力反演方法的有效性，分別開(kāi)展了二維及三維模型實(shí)驗(yàn)，并與聚焦反演算法進(jìn)行了對(duì)比分析.在實(shí)驗(yàn)過(guò)程中均未加入密度閾值約束，即不在反演過(guò)程中限定密度的最大值和最小值.

2.1 二維模型

二維反演包括三個(gè)模型實(shí)驗(yàn)，分別是無(wú)背景模型干擾的兩個(gè)地質(zhì)體、存在背景模型干擾的兩個(gè)地質(zhì)體和存在背景模型干擾的傾斜狀地質(zhì)體.三個(gè)模型實(shí)驗(yàn)均采用相同的觀(guān)測(cè)系統(tǒng)，具體參數(shù)為：觀(guān)測(cè)點(diǎn)個(gè)數(shù)為30個(gè)，間距為200 m，觀(guān)測(cè)點(diǎn)均位于海拔0 m.反演網(wǎng)格采用長(zhǎng)方形剖分，X、Z方向剖分的個(gè)數(shù)為30×28個(gè)，網(wǎng)格大小為200 m×50 m，從海拔0 m向下剖分，觀(guān)測(cè)點(diǎn)位于第一層網(wǎng)格的中心位置.

(1)模型一：無(wú)背景模型干擾的兩個(gè)地質(zhì)體

該模型包括兩個(gè)孤立的、不同大小的長(zhǎng)方體(圖4b)，參數(shù)見(jiàn)表1.假設(shè)構(gòu)造背景異常被徹底分離，即反演區(qū)域的背景密度為0 g·cm-3，模型參數(shù)見(jiàn)表1.觀(guān)測(cè)數(shù)據(jù)加入了3%的高斯噪聲(圖4a).

圖4 反演結(jié)果對(duì)比(a) 重力異常(黑線(xiàn))及擬合數(shù)據(jù)(紅線(xiàn)); (b) 密度模型; (c) 常規(guī)聚焦反演結(jié)果; (d) 本文反演結(jié)果.Fig.4 Comparison of inversion results(a) Gravity anomaly (black line) and fitting data (red line); (b) Density model; (c) The results of conventional focus inversion; (d) The results of this paper.

表1 二維密度模型參數(shù)Table 1 Parameters of 2D density model

構(gòu)建模型特征矩陣采用了兩類(lèi)模型特征矩陣構(gòu)建方式，一類(lèi)是選擇一個(gè)特征模型，X、Z方向尺寸為600 m×200 m的矩形密度體；另一類(lèi)是選擇三個(gè)幾何形態(tài)差異較大的特征模型，X、Z方向尺寸分別為200 m×200 m，600 m×200 m，1600 m×100 m時(shí).設(shè)置反演參數(shù)(λ1=0.1、λ2=2，最大迭代次數(shù)IterMax=150).

從該模型實(shí)驗(yàn)可看出，在沒(méi)有密度閾值約束前提條件下，常規(guī)聚焦反演結(jié)果(圖4c)雖然可以較好地反映地質(zhì)體質(zhì)心的埋深，但是聚焦程度的控制比較難以施加，易導(dǎo)致反演結(jié)果過(guò)于聚焦而出現(xiàn)較大的密度值.利用本文所提方法進(jìn)行反演，采用兩類(lèi)模型特征矩陣構(gòu)建方式都會(huì)得到相同的結(jié)果，如圖4d所示.可見(jiàn)，即使采用第二類(lèi)模型特征矩陣構(gòu)建方式，在分解系數(shù)稀疏約束的控制下，仍在三個(gè)特征模型中選擇了(600 m×200 m)特征模型作為待反演密度異常體的主要構(gòu)成單元，與稀疏分解系數(shù)進(jìn)行組合，共同恢復(fù)出了形態(tài)與密度值更接近于真實(shí)密度模型的反演結(jié)果.

(2)模型二：存在背景模型干擾的兩個(gè)地質(zhì)體

該模型包括三個(gè)密度體(圖5b)，其中兩個(gè)為孤立的、不同大小、不同埋深的長(zhǎng)方形密度體，另外一個(gè)為埋深較大的地質(zhì)體(密度為0.2 g·cm-3)作為背景干擾.與模型一不同的是，兩個(gè)長(zhǎng)方形的剩余密度體值有正、有負(fù)，具體參數(shù)見(jiàn)表1.在模型正演數(shù)據(jù)中加入了3%的高斯噪聲作為觀(guān)測(cè)數(shù)據(jù)(圖5a).

圖5 反演結(jié)果對(duì)比(a) 重力異常(黑線(xiàn))及擬合數(shù)據(jù)(紅線(xiàn)); (b) 密度模型; (c) 常規(guī)聚焦反演結(jié)果; (d) 本文反演結(jié)果.Fig.5 Comparison of inversion results(a) Gravity anomaly (black line) and fitting data (red line); (b) Density model; (c) The results of conventional focus inversion; (d) The results of this paper.

構(gòu)建模型特征矩陣采用了四個(gè)幾何形態(tài)差異較大的特征模型，X、Z方向尺寸分別為200 m×200 m，600 m×150 m，1000 m×250 m，15000 m×100 m.設(shè)置反演參數(shù)(λ1=0.2、λ2=2，最大迭代次數(shù)IterMax=150).

從該模型實(shí)驗(yàn)可看出，利用公式(5)求解得到的常規(guī)聚焦反演結(jié)果(圖5c)與原始模型偏差較大，原因主要是由于深部地質(zhì)體產(chǎn)生的背景異常起到了干擾作用，聚焦后左側(cè)的正密度體埋深加大，結(jié)果偏向與深層背景地質(zhì)體與淺層地質(zhì)體之間，而右側(cè)的負(fù)密度體反演結(jié)果較淺；而從本文方法獲得反演結(jié)果(圖5d)可看出，由于模型特征矩陣構(gòu)建過(guò)程中包含了水平方向長(zhǎng)條形特征模型(15000 m×100 m)的引入，在對(duì)應(yīng)的分解系數(shù)求解過(guò)程中一定程度上抵消了背景異常的干擾，相當(dāng)于一個(gè)異常自動(dòng)分離的過(guò)程，而其余的部分又由于包含了更接近于原始模型中塊體特征的特征模型(1000 m×250 m)，雖然結(jié)果的密度值較原有模型略大，但反演結(jié)果的空間分布更加逼近于原始模型.

(3)模型三：存在背景模型干擾的傾斜狀地質(zhì)體

該模型包括一個(gè)傾斜狀的地質(zhì)體，具體參數(shù)見(jiàn)表1.同時(shí)，與模型二類(lèi)似，模擬剩余異常仍包含一定的深層信息的情形，在深部(埋深在1150～1400 m范圍內(nèi)的網(wǎng)格)加入一個(gè)密度沿著X方向變化的橫向非均勻地質(zhì)體(密度值介于0.2～0.35 g·cm-3).在模型正演數(shù)據(jù)中加入了3%的高斯噪聲作為觀(guān)測(cè)數(shù)據(jù).

構(gòu)建模型特征矩陣采用了兩種特征模型，一是傾斜地質(zhì)體(圖6)，另外一類(lèi)是長(zhǎng)方形(X、Z方向尺寸為5000 m×100 m).設(shè)置反演參數(shù)(λ1=0.2、λ2=2，最大迭代次數(shù)IterMax=150).

圖6 構(gòu)建模型特征矩陣的二維特征模型Fig.6 2D feature model for building model feature matrix

從該模型實(shí)驗(yàn)可看出，利用公式(5)獲得的常規(guī)反演結(jié)果(圖7c)雖然在一定程度上體現(xiàn)了傾斜狀地質(zhì)體的特征，但由于深部地質(zhì)體產(chǎn)生的背景異常的干擾，為同時(shí)保證異常的擬合以及反演模型的連續(xù)聚焦目標(biāo)，使得結(jié)果埋深較大；而反觀(guān)利用本文方法得到的反演結(jié)果(圖7d)，在層狀特征模型的模型特征矩陣的聯(lián)合控制下，背景密度異常體的位置得到較好的恢復(fù).同時(shí)，具有傾斜形狀的特征模型(與模型大小并不相同)的引入又大大增加了逼近原始模型的可能性，使得即使沒(méi)有在反演過(guò)程中控制密度閾值范圍，恢復(fù)的淺層傾斜狀地質(zhì)體密度值及空間分布也與原始模型十分接近.

圖7 反演結(jié)果對(duì)比(a) 重力異常(黑線(xiàn))及擬合數(shù)據(jù)(紅線(xiàn)); (b) 密度模型; (c) 常規(guī)聚焦反演結(jié)果; (d) 本文反演結(jié)果.Fig.7 Comparison of inversion results(a) Gravity anomaly (black line) and fitting data (red line); (b) Density model; (c) The results of conventional focus inversion; (d) The results of this paper.

2.2 三維模型

三維模型反演實(shí)驗(yàn)中，重力觀(guān)測(cè)系統(tǒng)在X、Y方向的間距為200 m，觀(guān)測(cè)網(wǎng)格X、Y方向的個(gè)數(shù)為20×20=400個(gè).反演網(wǎng)格采用長(zhǎng)方體剖分，剖分網(wǎng)格在X、Y、Z方向的大小為200 m×200 m×200 m，反演網(wǎng)格在X、Y、Z方向剖分為20×20×10=4000個(gè)網(wǎng)格單元.三維模型形態(tài)呈“Y”字型，由2個(gè)傾斜狀地質(zhì)體構(gòu)成，如圖8a所示，參數(shù)見(jiàn)表2.三維密度模型正演模擬獲得的數(shù)據(jù)添加了3%的噪聲作為觀(guān)測(cè)數(shù)據(jù)，如圖8b所示.

表2 “Y”字型地質(zhì)體密度模型參數(shù)Table 2 Parameters of Y-type geological body density model

圖8 三維密度模型及重力異常(a) 三維密度模型； (b) 重力正演異常(含3%高斯噪聲).Fig.8 3D density model and gravity anomaly(a) 3D density model; (b) Gravity forward anomaly (including 3% Gaussian noise).

構(gòu)建模型特征矩陣采用了兩個(gè)三維特征模型，一是沿X軸正方向的向上右傾的三維特征模型a，如圖9a所示，另一個(gè)是沿X軸反方向的向上左傾的三維特征模型b，如圖9b所示.設(shè)置反演參數(shù)(λ1=1、λ2=1，最大迭代次數(shù)IterMax=200).

圖9 構(gòu)建模型特征矩陣的三維特征模型(a) 三維特征模型a； (b) 三維特征模型b.Fig.9 3D feature model for building model feature matrix(a) 3D feature model a; (b) 3D feature model b.

從該模型實(shí)驗(yàn)可看出，利用公式(5)獲得的常規(guī)反演結(jié)果(圖10c、圖10d)雖然在中心埋深上與模型比較接近，且一定程度上體現(xiàn)了“Y”字型特征，但反演結(jié)果分辨率較低，且傾斜體II的形態(tài)與埋深均得不到較好的恢復(fù)；而觀(guān)察本文方法得到的反演結(jié)果(圖10e、圖10f)，雖然構(gòu)建模型特征矩陣的兩類(lèi)傾斜形狀特征模型并非與真實(shí)模型的傾斜異常體大小相同，但包含接近的傾斜狀發(fā)育特征，該約束信息的引入使得反演過(guò)程中更加注重傾斜狀密度異常體的生成，使得反演結(jié)果逼近原始模型的可能性大大增加，即便沒(méi)有施加密度閾值約束，模型的分辨率與形態(tài)特征都與原始模型保持更高的一致性.

圖10 反演結(jié)果對(duì)比(a) 理論模型A-A′垂直剖面； (b) 理論模型B-B′垂直剖面； (c) 聚焦反演A-A′垂直剖面； (d) 聚焦反演B-B′垂直剖面； (e) 本文方法A-A′垂直剖面； (f) 本文方法B-B′垂直剖面.Fig.10 Comparison of inversion results(a) A-A′ profile of theoretical model; (b) B-B′ profile of theoretical model; (c) A-A′ profile of focus inversion; (d) B-B′ profile of focus inversion; (e) A-A′ profile of proposed method; (f) B-B′ profile of proposed method.

通過(guò)二維及三維模型實(shí)驗(yàn)可看出，相比常規(guī)反演方法，利用本文所提方法得到的反演結(jié)果在分辨率和可靠性方面都得到了有效提升.

雖然在進(jìn)行分解系數(shù)求解時(shí)，未知數(shù)的個(gè)數(shù)會(huì)有所增加，但是特征模型的引入反而在一定程度減少了密度模型的求解空間，這是由于期望獲得的密度模型不僅需要考慮正演結(jié)果與實(shí)測(cè)數(shù)據(jù)的擬合，還需要滿(mǎn)足由盡可能少的特征模型來(lái)組合得到；而且從重力異常解釋效率的角度來(lái)看，本文所提方法可更加靈活、快速地施加密度模型約束信息，無(wú)需反復(fù)調(diào)整網(wǎng)格剖分方案及過(guò)多地調(diào)整重力異常分離參數(shù)，在一定程度上克服了求解分解系數(shù)帶來(lái)的計(jì)算量增加的問(wèn)題.

3 實(shí)際資料應(yīng)用測(cè)試

現(xiàn)將本文方法應(yīng)用于墨西哥薩卡特卡斯州圣尼古拉斯硫化物銅鋅礦區(qū)的重力實(shí)際資料.研究區(qū)的礦床賦存于鎂鐵質(zhì)和長(zhǎng)英質(zhì)火山巖中，根據(jù)鉆井及巖心資料可知該硫化物礦床具有高密度、高磁化率、高極化率和低電阻率的特點(diǎn)(Phillips et al.，2001)，其中密度可達(dá)3.5 g·cm-3，較圍巖大概有1.1～1.4 g·cm-3的密度差，且埋深較淺，便于利用重力勘探技術(shù)預(yù)測(cè)礦床的空間分布.工區(qū)重力測(cè)點(diǎn)數(shù)為198個(gè)，將其利用克里金插值為均勻網(wǎng)格異常(圖11a)，X(東西，-2700～-600 m)、Y(南北-1100～600 m)方向間距分別為100 m、100 m，觀(guān)測(cè)網(wǎng)格X、Y方向的個(gè)數(shù)為22×18=396個(gè).

利用本文方法開(kāi)展重力三維反演.首先將地下半空間沿X、Y、Z三個(gè)方向剖分為22×18×15=4000個(gè)網(wǎng)格單元，尺寸為100 m×100 m×100 m.第一個(gè)網(wǎng)格的位置X、Y、Z坐標(biāo)位于(-2750 m，-1150 m，0 m).構(gòu)建模型特征矩陣采用了三個(gè)長(zhǎng)方體特征模型，在X、Y、Z方向的尺寸分別為300 m×300 m×300 m、800 m×200 m×400 m和1800 m×1200 m×200 m，前兩個(gè)特征模型用于構(gòu)建局部密度異常體，第三個(gè)特征模型用于擬合均勻的背景剩余密度.給定的分解系數(shù)初始模型為0.1，設(shè)置反演參數(shù)(λ1=1、λ2=1，最大迭代次數(shù)IterMax=200).

利用本文方法獲得的密度反演三維結(jié)果如圖11b所示，經(jīng)過(guò)礦床的南北A-A′的垂向剖面如圖11c所示，東西向垂直剖面B-B′如圖11d所示，密度反演結(jié)果與根據(jù)鉆井資料獲得的硫化物礦床地質(zhì)認(rèn)識(shí)(Phillips et al.，2001)十分吻合.對(duì)比國(guó)內(nèi)外研究學(xué)者在該礦區(qū)開(kāi)展的工作(Lelivere and Oldenburg，2009；李澤林等，2019)，反演結(jié)果在其他區(qū)域也恢復(fù)出了具有較高分辨率且更易解釋的密度分布，有效驗(yàn)證了本文方法的適用性.分析原因，模型特征矩陣中引入的背景特征模型使得在三維反演結(jié)果中產(chǎn)生了較為合理的背景密度模型，在一定程度上消除了區(qū)域異常分離不準(zhǔn)確的問(wèn)題.與此同時(shí)，在分解系數(shù)的稀疏求解過(guò)程中，不同尺寸的特征模型得到優(yōu)選，并用來(lái)逼近目標(biāo)密度體，使得反演的可靠性得到了保證.進(jìn)一步的，在實(shí)際資料應(yīng)用中，針對(duì)不斷細(xì)化的勘探工作需求，可嘗試構(gòu)建更精細(xì)的特征模型及模型特征矩陣，從而不斷提升重力反演分辨率，逐步恢復(fù)出更符合真實(shí)情況的地質(zhì)體密度空間分布.

圖11 實(shí)際資料反演結(jié)果(a) 剩余重力異常; (b) 密度三維反演結(jié)果; (c) A-A′垂直剖面反演結(jié)果； (d) B-B′垂直剖面反演結(jié)果. (c)—(d)中黑色多邊形代表硫化物礦床的位置.Fig.11 Inversion results of actual data(a) Residual gravity anomaly; (b) 3D inversion result of density; (c) The result of A-A′ vertical profile; (d) The result of B-B′ vertical profile.The sulphide location is indicated by the black polygon in (c)—(d).

4 結(jié)論

本文提出了一種基于密度模型稀疏表征的重力反演方法.首先，將密度模型表征為模型特征矩陣和分解系數(shù)的乘積，其次，給出了由特征模型構(gòu)建模型特征矩陣的技術(shù)流程，再次，重新構(gòu)建重力反演目標(biāo)函數(shù)，將直接求解密度模型問(wèn)題轉(zhuǎn)換為求解分解系數(shù)，最后，給出了分解系數(shù)的稀疏求解方法，實(shí)現(xiàn)了可有效融合地質(zhì)模式信息、提高重力反演分辨率和可靠性的目的.理論模型實(shí)驗(yàn)和實(shí)際資料應(yīng)用測(cè)試表明，相比常規(guī)聚焦反演算法，在反演目標(biāo)包含多個(gè)地質(zhì)體或剩余異常無(wú)法準(zhǔn)確求取時(shí)，本文方法都可取得較好的應(yīng)用效果.此外，本文所提出的模型約束方法并不僅限于重力反演，也為磁力、電法、地震等其他地球物理反演方法提供了一種靈活添加地質(zhì)約束信息的有效途徑.

致謝感謝三位匿名評(píng)審專(zhuān)家和期刊編輯對(duì)本文提出的寶貴修改意見(jiàn).

附錄A 基于Majorization-Minimization優(yōu)化框架的分解系數(shù)稀疏求解

根據(jù)正文1.3節(jié)中公式(8)，原重力反演問(wèn)題的目標(biāo)函數(shù)寫(xiě)為關(guān)于分解系數(shù)Γ的最優(yōu)化目標(biāo)函數(shù)：

(A1)

由于在稀疏求解時(shí)，φ(Γ)通常為一不可微的凸函數(shù)，現(xiàn)定義函數(shù)G(Γ)形式為：

G(Γ)=a·Γ2+b,

(A2)

其中，a、b為常數(shù).

可知函數(shù)G(Γ)是關(guān)于Γ的嚴(yán)格凸函數(shù).

根據(jù)MM框架理論，可知在待求的Γ=Γv處，G(Γv)應(yīng)滿(mǎn)足如下假設(shè)：

(A3)

其中，Γv表示第v次迭代所對(duì)應(yīng)的Γ值.

將公式(A3)代入公式(A2)，則

(A4)

其中，a和b為常數(shù)項(xiàng).

求解方程組(A4)，可得a和b.

(A5)

將公式(A5)代入公式(A2)，G(Γ)可寫(xiě)為：

(A6)

若Γ、Γv都為大小為q×1的向量，則根據(jù)MM思想：

(A7)

設(shè)

(A8)

其中，WΓv為與第v次迭代求解結(jié)果對(duì)應(yīng)的對(duì)角加權(quán)矩陣，c為與Γv有關(guān)的數(shù)值.

(A9)

因此，原目標(biāo)函數(shù)(A1)可表示為：

(A10)

對(duì)其求最小二乘解，可得：

(A11)