劉 星，鄧居智，陳 曉

(1.東華理工大學(xué)核資源與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室，南昌 330013;2.東華理工大學(xué)地球物理與測(cè)控技術(shù)學(xué)院，南昌 330013)

三維物性反演是重力反演研究的重點(diǎn)，通過(guò)數(shù)據(jù)利用反演方法推算出地下的物性分布情況，從而達(dá)到確定地下地質(zhì)地球物理模型的目的。近年來(lái)，在Tikhonov正則化理論框架下，中外學(xué)者針對(duì)三維反演方法開(kāi)展了大量研究工作，新的反演方法的研究一直是重力反演研究的熱門(mén)。早期反演得到的結(jié)果多是光滑約束的，為了得到更加突出尖銳邊界的結(jié)果，出現(xiàn)了聚焦算法。Last等[1]最早提出最小支撐泛函來(lái)進(jìn)行重力反演的思想，為后來(lái)的聚焦反演奠定了基礎(chǔ)；Portniaguine等[2]通過(guò)引入最小梯度支撐泛函來(lái)實(shí)現(xiàn)聚焦反演，得到清晰且更聚焦的結(jié)果；Zhdanov等[3]在重力梯度數(shù)據(jù)反演中采用三維聚焦反演方法并取得了與理論模型邊界較吻合的反演成果；Rezaie等[4]利用聚焦反演方法實(shí)現(xiàn)了三維重力數(shù)據(jù)的反演，得到一個(gè)能較好吻合真實(shí)異常體形狀和范圍的聚焦解。中國(guó)方面，劉小軍等[5]將聚焦反演方法應(yīng)用于二維大地電磁反演中獲得了比傳統(tǒng)方法更聚焦的反演結(jié)果；粟學(xué)磊等[6]將聚焦反演方法應(yīng)用于重力資料反演中，并引入深度加權(quán)函數(shù)和異常源外擴(kuò)方法來(lái)進(jìn)一步改善反演結(jié)果；陳閆等[7]將三維聚焦反演方法應(yīng)用于重力梯度數(shù)據(jù)中，得到一個(gè)聚焦解，并成功地框定了異常體的分布范圍；楊嬌嬌等[8]實(shí)現(xiàn)了三維重力梯度數(shù)據(jù)的聚焦反演，并與最平滑穩(wěn)定泛函反演結(jié)果進(jìn)行了對(duì)比，結(jié)果表明該方法能獲得更聚焦的結(jié)果；劉?？档萚9]引入粗糙度約束和深度加權(quán)，對(duì)重力數(shù)據(jù)進(jìn)行了三維聚焦反演，較清晰地還原了異常體的邊界；高秀鶴等[10]對(duì)三維重力梯度數(shù)據(jù)聚焦反演方法進(jìn)行了研究，反演結(jié)果顯示異常體密度分布集中且邊界清晰。

由于三維反演方法會(huì)出現(xiàn)負(fù)的或者很大的異常值的現(xiàn)象，超出了實(shí)際真實(shí)異常值范圍，使得反演結(jié)果出現(xiàn)多解性問(wèn)題。因此，為了解決這個(gè)問(wèn)題，引入不等式約束讓反演參數(shù)限制在有意義的約束范圍內(nèi)。如何根據(jù)先驗(yàn)物性信息，選擇有效的約束函數(shù)來(lái)將上下限約束融入反演中，是反演領(lǐng)域的研究熱點(diǎn)。目前，各種方法已經(jīng)應(yīng)用于不同的反演方案中以實(shí)現(xiàn)物性不等式約束。絕對(duì)約束法[2,11-12]是比較常見(jiàn)的解決方法，該方法是將迭代過(guò)程中所有超出取值范圍的反演參數(shù)賦值為其接近的上限或下限值，從而使反演得的物性參數(shù)值保持在上下限范圍內(nèi)。Li等[13]引入了對(duì)數(shù)懲罰函數(shù)法(障礙函數(shù)法)將物性約束添加到目標(biāo)函數(shù)中進(jìn)行反演，該方法是沿模型可行域的邊界形成了一個(gè)對(duì)數(shù)屏障以防模型參數(shù)在最小化過(guò)程中越界，從而實(shí)現(xiàn)上下限約束；在Li等[13]的基礎(chǔ)上，Zhang等[14]通過(guò)在模型增廣目標(biāo)函數(shù)中使用拉格朗日乘數(shù)的方法對(duì)三維位場(chǎng)數(shù)據(jù)的反演中不等式約束，將反演參數(shù)搜索區(qū)間限定在有效的取值區(qū)間內(nèi)，從而改善反演解的多解性問(wèn)題；王智等[15]用懲罰函數(shù)法修改目標(biāo)函數(shù)，利用權(quán)重因子將物性約束引入到非線性共軛梯度反演目標(biāo)函數(shù)中，但需要人為確定對(duì)應(yīng)的權(quán)重因子取值；Kim等[16]在最小二乘反演中引入了不等式約束，并利用對(duì)數(shù)變換函數(shù)方法將變化的物性參數(shù)值約束在有效的范圍內(nèi)；Commer[17]在三維重力反演中引入反雙曲線正切變換方法，在上下限范圍內(nèi)約束物性參數(shù)在無(wú)界域中的轉(zhuǎn)換，從而保證反演參數(shù)的合理性；劉斌等[18]通過(guò)對(duì)數(shù)變換方法來(lái)進(jìn)行基于不等式約束的反演，將參數(shù)轉(zhuǎn)換為無(wú)界域的新參數(shù)，利用新的參數(shù)來(lái)完成反演迭代過(guò)程，最終反演的參數(shù)值比未施加不等式約束前更接近真實(shí)異常值，有效減少了一些假異常；劉銀萍等[19]將上下限約束以對(duì)數(shù)變換函數(shù)形式引入到重力三維反演中，將不具有實(shí)際物理意義的反演結(jié)果轉(zhuǎn)化到合理的范圍，有效地限制了反演結(jié)果中剩余密度分布的數(shù)值范圍，還原模型剩余密度；Lin等[20]通過(guò)在三維反演中利用對(duì)數(shù)變換函數(shù)在有效范圍內(nèi)約束密度變化，控制反演解的范圍，從而降低了反演多解性問(wèn)題；Rezaie等[21]對(duì)具有稀疏約束的數(shù)據(jù)空間反演算法進(jìn)行了修改，使得物性邊界約束可以以對(duì)數(shù)變換函數(shù)的形式應(yīng)用于三維重力反演中，虛假構(gòu)造得到減少，反演異常值接近真實(shí)值，異常體更加突出。

可以看出，聚焦反演和不等式約束均是反演領(lǐng)域的研究熱點(diǎn)和前沿。基于此，嘗試將不等式約束引入三維重力聚焦反演，并通過(guò)模型試驗(yàn)來(lái)驗(yàn)證方法的有效性。

1 基于不等式約束的重力三維聚焦反演

1.1 聚焦反演

三維重力聚焦正則化反演目標(biāo)函數(shù)的表達(dá)式可寫(xiě)為

Pα(m)=φ(m)+αS(m)={Wd[A(m)-d]}T{Wd[A(m)-d]}+

α[WeWm(m-mapr)]T[WeWm(m-mapr)]

(1)

式(1)中：Pα(m)表示反演目標(biāo)泛函；φ(m)表示數(shù)據(jù)擬合函數(shù)；S(m)表示模型穩(wěn)定函數(shù)；α為正則化因子；A(m)為正演響應(yīng)函數(shù)；m為模型參數(shù)；d為觀測(cè)數(shù)據(jù)；mapr為模型先驗(yàn)信息；Wd、Wm分別為數(shù)據(jù)和模型加權(quán)矩陣；We為聚焦加權(quán)泛函矩陣。采用文獻(xiàn)[22]提出的自適應(yīng)正則化算法確定正則化因子。為了避免得到傳統(tǒng)反演方法過(guò)于光滑的反演結(jié)果，采用的模型穩(wěn)定泛函為最小支撐泛函[22]，其對(duì)應(yīng)的聚焦加權(quán)泛函矩陣為

(2)

1.2 對(duì)數(shù)不等式約束

為了將反演的參數(shù)控制在合理的物性約束范圍內(nèi)，可以將模型物性參數(shù)取值范圍作為不等式約束條件應(yīng)用到反演過(guò)程中，即

ρmini≤mi≤ρmaxi，i=1,2,…,m

(3)

式(3)中：mi為第i個(gè)網(wǎng)格單元的物性參數(shù)；ρmini和ρmaxi分別為此網(wǎng)格單元的物性下限和上限。

(4)

(5)

(6)

(7)

因此，新的模型參數(shù)的靈敏度矩陣表示為

(8)

1.3 共軛梯度法

共軛梯度(conjugate gradient，CG)算法解決目標(biāo)函數(shù)最小化的問(wèn)題，CG算法的計(jì)算流程為

(9)

βn=‖ln‖2/‖ln-1‖2

(10)

(11)

(12)

(13)

(14)

2 模型試驗(yàn)

為了驗(yàn)證基于不等式約束的聚焦反演的有效性，以重力反演為例，設(shè)計(jì)了兩個(gè)密度模型，采用規(guī)則六面體單元剖分網(wǎng)格，通過(guò)利用文獻(xiàn)[24]中的直立長(zhǎng)方體重力異常積分解公式進(jìn)行正演計(jì)算，再對(duì)模型進(jìn)行三維反演。觀測(cè)區(qū)剖分為21×21 個(gè)觀測(cè)點(diǎn)，將地下均勻半空間剖分為21×21×10 個(gè)網(wǎng)格，網(wǎng)格間距均為50 m，x方向?yàn)檎龞|方向，y方向?yàn)檎狈较颉?/p>

2.1 單塊體模型

單塊體模型分布如圖1所示，在密度設(shè)置為0的地下均勻半空間中，存在一個(gè)中心坐標(biāo)為(475、450、200 m)的長(zhǎng)方體異常體，正東向邊長(zhǎng)為350 m，正北向邊長(zhǎng)為300 m，異常體頂部埋深為100 m，底部埋深為300 m，其剩余密度為0.5 g/cm3。

圖1 單塊體模型切片和立體圖Fig.1 Single cube model slices and 3D structure diagram

設(shè)計(jì)了如下的對(duì)比試驗(yàn)：方案一是無(wú)正則化項(xiàng)反演；方案二是聚焦反演；方案三是基于不等式約束的聚焦反演，該方案給所有的模型均勻半空間網(wǎng)格單元體賦予密度約束范圍為0～0.6 g/cm3；方案四是基于分區(qū)域不等式約束的聚焦反演，在深度范圍100 m≤z≤300 m的區(qū)域?qū)⒉坏仁郊s束范圍設(shè)為0～0.6 g/m3，而其余區(qū)域網(wǎng)格設(shè)置為0～0.2 g/cm3。為了能更好地進(jìn)行對(duì)比，4種反演方案的反演參數(shù)均是一致的，其中模型先驗(yàn)信息mapr=0，反演初始模型選取為0，反演迭代收斂終止條件為均方誤差(root mean square，RMS)=0.01。

圖2～圖5分別為無(wú)正則化項(xiàng)反演、未施加不等式約束聚焦反演、不等式約束條件下的聚焦反演和分區(qū)域不等式約束聚焦反演結(jié)果。對(duì)比3種反演方案的反演結(jié)果可看出，無(wú)正則化項(xiàng)反演所得的密度分布立體圖與理論模型偏差較大[圖2(a)]，反演密度集中分布于地表附近0～100 m處，圖2(b)～圖2(d)顯示異常值最大值為0.3 g/cm3，相對(duì)原模型較小，無(wú)法還原出原模型的位置分布及密度。未施加不等式約束的聚焦反演結(jié)果如圖3所示，圖3 (a)中顯示異常體與理論模型在位置上相近，但異常體的分布范圍比理論模型大，圖3(b)所示的正東方向切片圖表明在100～300 m深處較好地吻合理論模型，在深度300～500 m處，反演結(jié)果存在0.2～0.5 g/cm3的假異常。圖3(c)表明反演的模型分布范圍在淺層100～200 m基本與原模型吻合，水平分辨率也相對(duì)較高，聚焦效果明顯，深部邊界擬合相對(duì)較差。從圖3(b)～圖3(d)可以看出，反演結(jié)果出現(xiàn)了一些負(fù)密度，且異常體中心部位出現(xiàn)大于真實(shí)異常值的密度，異常體剩余密度為0.3～0.53 g/cm3。

圖2 無(wú)正則化項(xiàng)反演結(jié)果切片和立體結(jié)構(gòu)Fig.2 The slices and 3D structure of the inversion results without regularization

圖3 未施加不等式約束的聚焦反演結(jié)果切片和立體結(jié)構(gòu)Fig.3 The slices and 3D structure of thefocusing inversion results without inequality constraints

圖4 不等式約束條件下的聚焦反演結(jié)果切片和立體結(jié)構(gòu)Fig.4 The slices and 3D structure of thefocusing inversion results with inequality constraints

圖5 分區(qū)域不等式約束條件下的聚焦反演結(jié)果切片和立體圖Fig.5 The slices and 3D structure of thefocusing inversion results with subregional inequality constraints

相較未施加不等式約束的反演結(jié)果而言，圖4的施加不等式約束聚焦反演的密度模型與理論模型位置、分布形狀基本一致[圖4(a)]，能夠較好地還原原模型分布形態(tài)，并且由圖4(b)～圖4(d)的切片圖像可以看出，沒(méi)有出現(xiàn)密度范圍之外的值，異常值不超過(guò)0.5 g/cm3，假異常得到了很大的減少，深度方向的分辨率得到了提高，驗(yàn)證了不等式約束能夠減少解的多解性問(wèn)題，提高聚焦反演解的質(zhì)量。雖然還是存在異常體的底部邊界較模糊的現(xiàn)象，在深度300～350 m處仍然存在0.3～0.35 g/cm3的虛假構(gòu)造，但總體而言，基于不等式約束聚焦反演方法是可行的。

由圖5可以看出，與不等式約束聚焦反演結(jié)果相比有進(jìn)一步的改善，圖5(a)的反演的密度模型與原模型位置、規(guī)模范圍、分布形狀基本吻合，圖5(b)～圖5(d)的切片圖顯示異常體剩余密度值范圍為0.3～0.5 g/cm3，地下異常體深部邊界信息清晰，進(jìn)一步驗(yàn)證了分區(qū)域不等式約束可以提高聚焦反演結(jié)果精度從而還原真實(shí)異常體。圖6為4種反演方案的RMS迭代誤差擬合曲線，從中可看出，RMS在反演迭代初期迅速減小，隨著迭代次數(shù)的增加而緩慢減小，直到達(dá)到設(shè)定的擬合值迭代結(jié)束。

圖6 迭代誤差曲線Fig.6 Curves of RMS Iterative error

2.2 臺(tái)階狀密度模型

設(shè)有一個(gè)臺(tái)階狀地下異常體，其模型分布如圖7所示，是5個(gè)250 m×300 m×50 m的長(zhǎng)方體疊加交錯(cuò)50 m拼接成的階梯狀模型，埋深范圍為50～300 m，異常體剩余密度為Δρ=1.0 g/cm3，圍巖密度均為0。

圖7 臺(tái)階狀模型切片和立體圖Fig.7 Inclined plate model slices and 3D structure diagram

與模型一類似，設(shè)計(jì)了如下的對(duì)比試驗(yàn)：方案一是無(wú)正則化項(xiàng)反演；方案二是聚焦反演；方案三是基于不等式約束的聚焦反演，該方案給所有的模型均勻半空間網(wǎng)格單元體賦予密度約束范圍為0～1.1 g/cm3；方案四是基于分區(qū)域不等式約束的聚焦反演，在深度范圍50 m≤z≤300 m的區(qū)域?qū)⒉坏仁郊s束范圍設(shè)為0～1.1 g/cm3，而其余區(qū)域網(wǎng)格設(shè)置為0～0.5 g/cm3。為了更好地進(jìn)行對(duì)比，4種反演方案的反演參數(shù)均是一致的，其中模型先驗(yàn)信息mapr=0，反演初始模型選取為0，反演迭代收斂終止條件為RMS=0.01。

圖8～圖11分別為無(wú)正則化項(xiàng)反演、未施加不等式約束聚焦反演、不等式約束條件下聚焦反演和分區(qū)域不等式約束條件下聚焦反演的結(jié)果。對(duì)比四種反演方案的反演結(jié)果可看出，圖8(a)所示的無(wú)正則化項(xiàng)反演得到三維模型的密度集中分布于臺(tái)階狀模型垂直上方的地表附近150 m內(nèi)，無(wú)法還原出真實(shí)模型的真實(shí)形態(tài)，且圖9(b)～圖9(d)的切片圖顯示反演所得剩余密度異常范圍小，除臺(tái)階狀模型頂層階梯對(duì)應(yīng)的橫向分布范圍的近地表0～50 m深處出現(xiàn)0.9 g/cm3左右的異常值，其余范圍內(nèi)的剩余密度值均小于0.6 g/cm3。未施加不等式約束的聚焦反演由圖9(a)表明在異常體位置出現(xiàn)形狀與真實(shí)模型相差較大的異常體，在淺部50～100 m深度處能夠較好地反映異常體的階梯形狀和邊界，深部反映效果差，圖9(b)～圖9(d)的切片圖顯示在淺部50～100 m深度處能夠較準(zhǔn)確地反演出異常體形狀和邊界，水平分辨率也相對(duì)較高，但在深部分辨率很差，在理論模型邊界外部出現(xiàn)了較多的假密度分布，異常體底部邊界延伸到400 m，且反演的背景場(chǎng)密度有出現(xiàn)較大的負(fù)的數(shù)值-0.3 g/cm3，異常體剩余密度值在0.5～1.1 g/cm3。

圖8 無(wú)正則化項(xiàng)反演結(jié)果切片和立體圖Fig.8 The slices and 3D structure of the inversion results without regularization

圖9 未施加不等式約束的聚焦反演結(jié)果切片和立體結(jié)構(gòu)Fig.9 The slices and 3D structure of thefocusing inversion results without inequality constraints

相較未施加不等式約束的聚焦反演結(jié)果而言，圖10(a)為施加不等式約束聚焦反演的三維密度模型，該模型與圖9 (a)對(duì)比，能較好地?cái)M合于理論模型位置、分布范圍和階梯形狀，圖10(b)～圖10(d)所示的切片圖顯示，該方法的反演異常體剩余密度范圍為0.5～1.0 g/cm3，有效抑制了密度范圍之外的值，水平和深度方向分辨率也得到很大改善，減少了虛假異常，也證明了基于不等式約束聚焦反演的可行性。

圖10 不等式約束條件下的聚焦反演結(jié)果切片和立體結(jié)構(gòu)Fig.10 The slices and 3D structure of thefocusing inversion results with inequality constraints

圖11所示的分區(qū)域不等式約束條件下的聚焦反演結(jié)果與圖10相比有進(jìn)一步的改善，圖11(a)顯示的三維密度模型與理論模型位置、規(guī)模范圍、分布形狀基本吻合，與圖10(a)相比，深度方向300～350 m處的缺陷得到補(bǔ)充，圖11(b)～圖11(d)的切片圖顯示異常體剩余密度值范圍為0.5～1.0 g/cm3，異常體深部邊界信息更加清晰，較好地還原了模型的真實(shí)分布特征,進(jìn)一步驗(yàn)證了分區(qū)域不等式約束可以提高聚焦反演結(jié)果精度從而提高異常體還原度。圖12為4種反演方案的RMS迭代誤差擬合曲線，從圖2可以看出，4種反演方案都能穩(wěn)定收斂。

圖11 分區(qū)域不等式約束條件下的聚焦反演結(jié)果切片和立體圖Fig.11 The slices and 3D structure of thefocusing inversion results with subregional inequality constraints

圖12 迭代誤差曲線Fig.12 Curves of RMS iterative error

3 結(jié)論

將先驗(yàn)物性信息的上限和下限作為約束條件，以對(duì)數(shù)轉(zhuǎn)換函數(shù)形式的將不等式約束融入到三維重力聚焦反演，模型試驗(yàn)表明以下結(jié)論。

(1)聚焦反演在淺部能夠較準(zhǔn)確地框定異常體形狀和范圍，但在深部存在一些虛假異常。

(2)基于不等式約束的三維聚焦反演可以將解控制在合理的物性范圍之內(nèi)，減少虛假異常構(gòu)造的產(chǎn)生，此外，分區(qū)域不等式約束的加入可以提高先驗(yàn)物性信息應(yīng)用的靈活性，反演結(jié)果更加逼近真實(shí)模型。

需要指出的，理論模型反演中給定的不等式約束的范圍是較為理想的，但是在實(shí)際情況中異常體分布信息更為復(fù)雜，在將一些豐富的地質(zhì)信息轉(zhuǎn)化成合理的約束條件的問(wèn)題上還存在困難，因此如何獲得比較合理準(zhǔn)確的邊界范圍以及物性的上下限約束條件需要進(jìn)一步進(jìn)行研究。