郭一豪，張志勇,2*，陳曉,2，李曼,2，周峰,2，謝尚平,2，曾志文，程三

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

0 引言

由于地球物理場(chǎng)的等效性、觀測(cè)數(shù)據(jù)集的有限性、數(shù)據(jù)觀測(cè)存在誤差等問(wèn)題,地球物理反演不可避免的具有多解性,穩(wěn)定、高精度的反演一直是地球物理反演研究追求的重要目標(biāo).Tikhonov正則化理論（Tikhonov and Arsenin,1977）自提出至今,被公認(rèn)為是可以減少反問(wèn)題不適定性的有效方法之一.地球物理反問(wèn)題,在考慮數(shù)據(jù)誤差的基礎(chǔ)上引入模型穩(wěn)定函數(shù)進(jìn)行正則化可以有效降低反問(wèn)題多解性;同時(shí),通過(guò)穩(wěn)定函數(shù)將先驗(yàn)信息融入反演,也使反演更符合實(shí)際情況、提高了反演的可靠性.

地球物理正則化反演中模型穩(wěn)定函數(shù)的設(shè)計(jì)與選擇直接影響反演的最終效果,其構(gòu)建是反演的關(guān)鍵環(huán)節(jié)之一.根據(jù)反演模型的特征,模型穩(wěn)定函數(shù)通?？煞譃楣饣s束和聚焦約束兩類(lèi).早期地球物理反演研究主要使用最小模型約束函數(shù)（胡正旺等,2021）.Constable等（1987）構(gòu)建了最大化光滑模型穩(wěn)定函數(shù),實(shí)現(xiàn)了大地電磁（Magnetotelluric,MT）數(shù)據(jù)的Occam反演,研究工作表明采用光滑約束可以在一定程度上防止數(shù)據(jù)過(guò)度擬合,獲得與真實(shí)地下物性分布更吻合的結(jié)果.相較光滑反演約束,Rudin等（1992）、Vogel和Oman（1998）在開(kāi)展圖像信息提取研究中提出了一種改進(jìn)的總變化量支撐模型穩(wěn)定函數(shù)（Total Variation,TV）,用于從噪聲數(shù)據(jù)中恢復(fù)圖像.Farquharson和Oldenburg（1998）在一維電磁反演中探討了L-1范數(shù)和L-2范數(shù)模型約束反演的魯棒性,研究表明L-1范數(shù)具有構(gòu)造分段或反演“塊狀”模型能力.Last和Kubik（1983）提出了最小支撐模型函數(shù)（Minimum Support,MS）約束并將其應(yīng)用于重力反演中,結(jié)果表明該約束反演結(jié)果具有致密性好、密度對(duì)比均勻等地質(zhì)特征.Portniaguine和Zhdanov（2002）將MS聚焦函數(shù)引入三維航空磁法反演,實(shí)測(cè)數(shù)據(jù)反演表明,使用MS聚焦函數(shù)可以獲得清晰的、聚焦的反演模型.Zhdanov等（2004）三維重力梯度張量MS聚焦約束反演的研究表明,該方法可以顯著提高反演結(jié)果的分辨率.Zhdanov等（2004）進(jìn)一步研究了基于MS約束的三維MT反演,實(shí)現(xiàn)了非二次最小支撐穩(wěn)定器轉(zhuǎn)化為傳統(tǒng)的二次最小范數(shù)穩(wěn)定器,使反問(wèn)題可以采用統(tǒng)一形式優(yōu)化.Portniaguine和Zhdanov（1999）研究了最小梯度支撐模型函數(shù)（Minimum Gradient Support,MGS）并將其應(yīng)用于重力反演研究中,對(duì)比分析了最光滑約束、總變化量支撐約束、最小支撐和最小梯度支撐的反演結(jié)果,研究表明MGS與加權(quán)函數(shù)相結(jié)合,有利于得到清晰、集中的反演模型.改進(jìn)的TV、MS和MGS約束函數(shù)中均存在一個(gè)小的正實(shí)數(shù),對(duì)于改進(jìn)TV約束函數(shù)來(lái)說(shuō),正實(shí)數(shù)是為了確保TV在零點(diǎn)處可導(dǎo)；對(duì)于MS和MGS約束函數(shù)來(lái)說(shuō),引入的正實(shí)數(shù)被稱(chēng)為聚焦因子,該因子直接影響反演結(jié)果,其選擇問(wèn)題一直是聚焦反演研究的一個(gè)難點(diǎn).Zhao等（2016）提出了一種指數(shù)型模型聚焦函數(shù),重磁反演研究表明該函數(shù)不僅具有聚焦作用,而且更加簡(jiǎn)便、穩(wěn)定,可以避免聚焦因子的選取及函數(shù)梯度計(jì)算困難等問(wèn)題.

另外，模型參數(shù)物性變換函數(shù)的使用一方面可以結(jié)合先驗(yàn)物性范圍對(duì)反演物性進(jìn)行合理的約束;另一方面,從一定程度上也可起到平衡模型誤差與數(shù)據(jù)誤差的作用,降低正則化因子選擇的難度.Habashy和Abubakar（2004）提出了用于電磁反演的變換函數(shù),確保模型參數(shù)在有實(shí)際物理意義的取值區(qū)間.此外，Loke和Barker（1995）對(duì)直流電阻率、Li和Oldenburg（2003）對(duì)磁法、Commer和Newman（2008）對(duì)三維可控源電磁（Controled Source Electromagnetic Method,CSEM）、Key（2016）對(duì)二維海洋電磁反演問(wèn)題的模型變換函數(shù)進(jìn)行了研究，均取得了理想效果.

隨著勘探程度的提高與物性研究的深入,取得的有效先驗(yàn)信息也越來(lái)越豐富.豐富精確的先驗(yàn)信息可以有效地減少反演多解性,提高反演解釋的正確性.鉆孔和測(cè)井的開(kāi)展提供了局部物性范圍和空間展布的雙重信息,因此融合鉆孔約束的反演被應(yīng)用于地球物理反演中（Bobée et al.,2010;Yan et al.,2017）.Dell′Aversana（2001）將通過(guò)測(cè)井資料確定的電阻率、速度和密度的統(tǒng)計(jì)關(guān)系應(yīng)用于順序反演中,研究了地質(zhì)情況復(fù)雜的逆沖斷層帶.周竹生和周熙襄（1993）將地震資料、測(cè)井信息和先驗(yàn)地質(zhì)知識(shí)相結(jié)合提出了一種寬帶約束反演方法,該方法能提供可與測(cè)井資料進(jìn)行對(duì)比的、具有較高分辨率的絕對(duì)波阻抗參數(shù),為巖性模擬奠定基礎(chǔ).相比測(cè)井及鉆孔信息約束,基于標(biāo)本測(cè)試等方式統(tǒng)計(jì)獲取的信息僅具有物性范圍約束能力,其往往難以直接應(yīng)用于反演約束.

為了將任意多種先驗(yàn)巖石物性約束融入地球物理反演,本文提出一種對(duì)稱(chēng)多項(xiàng)式模型約束函數(shù),該模型約束通過(guò)連乘形式將多種統(tǒng)計(jì)或參考物性引入反演.為實(shí)現(xiàn)反演目標(biāo)函數(shù)優(yōu)化求解,深入研究了該模型約束函數(shù)的梯度與海森矩陣；開(kāi)展了模型和實(shí)測(cè)數(shù)據(jù)反演,并通過(guò)與其他形式正則化模型約束對(duì)比分析,驗(yàn)證了該約束函數(shù)的正確性和有效性.

1 對(duì)稱(chēng)多項(xiàng)式約束

1.1 對(duì)稱(chēng)多項(xiàng)式的定義

定義如下的模型穩(wěn)定函數(shù)：

（1）

1.2 對(duì)稱(chēng)多項(xiàng)式的函數(shù)表示

式（1）中的φi為對(duì)稱(chēng)多項(xiàng)式（為方便討論忽略單元體積項(xiàng)）,

+c1T（mi）k-1+…+ck-1T（mi）1+cn，

（2）

系數(shù)c具有通解（推導(dǎo)見(jiàn)附錄）:

（3）

i=1,2,3,…,k.

（4）

其中,ci為基本對(duì)稱(chēng)多項(xiàng)式,可通過(guò)遞歸算法求解.

1.3 對(duì)稱(chēng)多項(xiàng)式函數(shù)曲線(xiàn)

為研究對(duì)稱(chēng)多項(xiàng)式函數(shù)的特點(diǎn),將其與光滑和最小支撐兩種地球物理中常用的模型約束函數(shù)進(jìn)行對(duì)比分析.

光滑約束函數(shù)

（5）

最小支撐約束函數(shù)

（6）

其中,ε為小的正實(shí)數(shù),稱(chēng)為聚焦因子.該因子是聚焦反演的關(guān)鍵參數(shù),Zhdanov和Tolstaya（2004）使用一系列的聚焦因子進(jìn)行了反演,獲得了模型誤差函數(shù)-聚焦因子曲線(xiàn),選取曲線(xiàn)斜率最大點(diǎn)處的聚焦因子作為反演最佳參數(shù),該方法需要在反演前進(jìn)行大量的反演試算,計(jì)算效率較低.

為了定性比較基于對(duì)稱(chēng)多項(xiàng)式和常用的一些模型穩(wěn)定函數(shù),假設(shè)一個(gè)反演單元存在三種可能巖性,物性分別為-1、0、1,繪制模型約束值變化曲線(xiàn)如圖1,圖中紅、綠、藍(lán)、黑線(xiàn)分別為最小模型、聚焦因子為0.01的最小支撐、聚焦因子為0.1的最小支撐以及對(duì)稱(chēng)多項(xiàng)式約束函數(shù)值變化曲線(xiàn).

圖1 三種模型穩(wěn)定函數(shù)曲線(xiàn)Fig.1 The curves of three types model constraint function

分析圖1,基于對(duì)稱(chēng)多項(xiàng)式的約束函數(shù)（黑色曲線(xiàn)）當(dāng)反演物性處在給定參考值的最大和最小值之間時(shí)，函數(shù)取值相對(duì)較小,當(dāng)物性超過(guò)給定約束范圍，則函數(shù)值會(huì)快速增大;當(dāng)物性取值與參考物性一致時(shí)，該函數(shù)取得極小值.對(duì)比最小模型（紅色曲線(xiàn)）、最小支撐（綠、藍(lán)色曲線(xiàn)），多項(xiàng)式約束具有明顯的物性分類(lèi)與取值區(qū)間約束能力.

2 基于對(duì)稱(chēng)多項(xiàng)式的反演

2.1 目標(biāo)函數(shù)

根據(jù)正則化反演理論,構(gòu)造如下目標(biāo)函數(shù)：

（7）

α與τ為用于平衡數(shù)據(jù)誤差與模型誤差的平衡因子,稱(chēng)為正則化因子.可采用經(jīng)驗(yàn)和自適應(yīng)選擇兩種方法確定.其中,經(jīng)驗(yàn)法是利用一定的準(zhǔn)則確定某個(gè)定值作為近似最佳正則化因子,該值一般在整個(gè)反演過(guò)程中保持不變.主要的定值方法包括:無(wú)偏風(fēng)險(xiǎn)估計(jì)法（Wahba,1990）,廣義交叉驗(yàn)證法（Wahba,1977；Golub et al.,1979）和L曲線(xiàn)法（Hansen and O′Leary,1993；Hansen,1994）.自適應(yīng)方法是在反演過(guò)程中根據(jù)一定的準(zhǔn)則自適應(yīng)調(diào)整正則化因子取值的方法.Zhdanov（2002）提出了一種利用初始數(shù)據(jù)擬合泛函和模型穩(wěn)定泛函的比值為正則化因子初值,在數(shù)據(jù)擬合效率下降時(shí)對(duì)正則化因子進(jìn)行衰減的自適應(yīng)方法;陳小斌等（2005）研究了最平緩約束MT反演正則化因子的自適應(yīng)選擇算法.為簡(jiǎn)化研究,本文反演中α與τ的選擇均使用經(jīng)驗(yàn)法.

2.2 對(duì)稱(chēng)多項(xiàng)式函數(shù)的一階導(dǎo)數(shù)向量

由公式（1）可見(jiàn),對(duì)稱(chēng)多項(xiàng)式約束求和中的每一項(xiàng)僅與一個(gè)單元物性相關(guān).因此有

（8）

其中,φi為關(guān)于mi的一元k次多項(xiàng)式函數(shù).進(jìn)一步擴(kuò)展，得到對(duì)稱(chēng)多項(xiàng)式約束的導(dǎo)數(shù)向量為

（9）

2.3 對(duì)稱(chēng)多項(xiàng)式函數(shù)的二階導(dǎo)數(shù)矩陣

結(jié)合一階導(dǎo)數(shù)的推導(dǎo)和公式（1）可知，對(duì)稱(chēng)多項(xiàng)式約束函數(shù)對(duì)物性的二階導(dǎo)數(shù)為

（10）

因此,對(duì)稱(chēng)多項(xiàng)式約束項(xiàng)的海森矩陣為

（11）

2.4 共軛梯度優(yōu)化流程

為簡(jiǎn)化推導(dǎo),參考Zhdanov（2015）,記m=T（m）.對(duì)目標(biāo)函數(shù)（7）的共軛梯度優(yōu)化過(guò)程進(jìn)行推導(dǎo).數(shù)據(jù)項(xiàng)殘差

Rn=A[T-1（mn）]-dobs，

（12）

其中,T-1（·）表示物性變換函數(shù)的反函數(shù),A[T-1（mn）]和dobs分別為計(jì)算數(shù)據(jù)和觀測(cè)數(shù)據(jù),n為當(dāng)前迭代次數(shù).目標(biāo)函數(shù)式（7）的梯度

（13）

（14）

（15）

（16）

（17）

模型更新,

（18）

3 模型試驗(yàn)

分別采用光滑約束、最小支撐約束和對(duì)稱(chēng)多項(xiàng)式約束進(jìn)行重力二、三維模型試驗(yàn),取變換函數(shù)為T(mén)（m）=m.

3.1 模型試驗(yàn)一

設(shè)計(jì)如圖2所示的剩余密度模型,在剩余密度為0 g·cm-3的背景中存在一個(gè)剩余密度為0.1 g·cm-3的塊狀異常體.

圖2 設(shè)計(jì)剩余密度模型Fig.2 Synthetic residual density model

觀測(cè)區(qū)域?yàn)?～10 km,測(cè)點(diǎn)距為500 m,光滑模型約束反演正則化因子選擇0.5,MS聚焦約束反演正則化因子選擇0.05,對(duì)稱(chēng)多項(xiàng)式約束反演取光滑模型約束項(xiàng)權(quán)重因子0.5、多項(xiàng)式約束項(xiàng)權(quán)重因子150,采用矩形網(wǎng)格剖分,模型單元中心位置見(jiàn)示意圖3.反演結(jié)果見(jiàn)圖4,其中圖4a、4b、4c分別為光滑模型、最小支撐（ε=0.03）、對(duì)稱(chēng)多項(xiàng)式約束（參考物性為m*=[0,0.1]）的反演結(jié)果,圖4d、4e、4f為與其對(duì)應(yīng)的誤差曲線(xiàn),圖4g為響應(yīng)曲線(xiàn).

圖3 模型離散單元中心節(jié)點(diǎn)位置Fig.3 Element center position of discrete mesh

對(duì)比圖4a和圖4c可以發(fā)現(xiàn),光滑反演結(jié)果縱向連續(xù)性好,但反演的剩余密度最大僅達(dá)到0.06 g·cm-3,與真實(shí)物性0.1 g·cm-3存在較大差距,而基于對(duì)稱(chēng)多項(xiàng)式約束的反演結(jié)果物性與真實(shí)模型更接近.對(duì)比圖4b和圖4c可以發(fā)現(xiàn),最小支撐函數(shù)約束反演物性可以達(dá)到真實(shí)物性,但存在聚焦作用過(guò)強(qiáng)風(fēng)險(xiǎn),導(dǎo)致反演異常體規(guī)模偏小.由圖4d、4e和4f可見(jiàn)三種方案的反演誤差函數(shù)隨迭代次數(shù)逐漸下降并最終趨于平緩,光滑與對(duì)稱(chēng)多項(xiàng)式約束反演數(shù)據(jù)誤差下降速度比最小支撐模型約束快,最后的擬合差也小于最小支撐模型約束.由圖4g可見(jiàn),三種約束方案均取得了較好的擬合效果,且光滑模型約束與對(duì)稱(chēng)多項(xiàng)式約束擬合效果優(yōu)于最小支撐模型約束結(jié)果.

3.2 模型試驗(yàn)二

設(shè)計(jì)如圖5所示的剩余密度模型,在剩余密度為0 g·cm-3的背景中存在剩余密度分別為0.1 g·cm-3和-0.1 g·cm-3的兩個(gè)塊狀異常體.

觀測(cè)區(qū)域?yàn)?～10 km,測(cè)點(diǎn)距為500 m,光滑反演正則化因子選擇0.3,MS聚焦約束反演正則化因子選擇0.05,對(duì)稱(chēng)多項(xiàng)式模型約束反演光滑項(xiàng)權(quán)重0.3、多項(xiàng)式項(xiàng)權(quán)重3000,單元剖分同圖3,反演結(jié)果見(jiàn)圖6,其中,圖6a,6b,6c分別為光滑、最小支撐（ε=0.1）、對(duì)稱(chēng)多項(xiàng)式約束（參考物性為m*=[0,0.1,-0.1]）的反演結(jié)果,圖6d,6e，6f為對(duì)應(yīng)的反演誤差曲線(xiàn),圖6g為響應(yīng)曲線(xiàn).

對(duì)比圖6a和圖6c可以發(fā)現(xiàn),光滑反演結(jié)果連續(xù)性好,但反演邊界模糊,物性仍然小于真實(shí)物性,而基于對(duì)稱(chēng)多項(xiàng)式約束的反演結(jié)果物性與真實(shí)模型更接近,且邊界更清晰.對(duì)比圖6b和圖6c可以發(fā)現(xiàn),最小支撐函數(shù)約束反演物性可以達(dá)到真實(shí)物性,但反演聚焦效果強(qiáng),異常體規(guī)模偏小.分析圖6d、6e和6f可以發(fā)現(xiàn)三種方案的反演誤差函數(shù)隨迭代次數(shù)逐漸下降并最終趨于平緩.分析圖6g可以發(fā)現(xiàn),三種方案的反演響應(yīng)與實(shí)測(cè)響應(yīng)基本吻合.

圖4 合成模型反演（a） 光滑模型約束反演密度斷面;（b） 最小支撐模型約束反演密度斷面;（c） 對(duì)稱(chēng)多項(xiàng)式約束反演密度斷面;（d） 光滑模型約束反演誤差曲線(xiàn);（e） 最小支撐模型約束反演誤差曲線(xiàn);（f） 對(duì)稱(chēng)多項(xiàng)式約束反演誤差曲線(xiàn);（g） 響應(yīng)擬合曲線(xiàn).Fig.4 Synthetic model inversion（a） The density section of smooth model constraint inversion;（b） The density section of minimum support model constraint inversion;（c） The density section of symmetric polynomial model constrain inversion;（d） The smooth model constraint inversion error curves;（e） The minimum support model constraint inversion error curves;（f） The symmetric polynomial model constraint inversion error curves;（g） Data fitting curves.

圖5 含兩個(gè)異常體的剩余密度模型Fig.5 Synthetic residual density model with two abnormal blocks

3.3 模型試驗(yàn)三

設(shè)計(jì)如圖7所示的三維剩余密度模型,在剩余密度為0 g·cm-3的背景中存在兩個(gè)剩余密度分別為0.1 g·cm-3和-0.1 g·cm-3的塊狀異常體.

觀測(cè)區(qū)域?yàn)殚L(zhǎng)寬各10 km的正方形,測(cè)點(diǎn)距為500 m,測(cè)線(xiàn)距500 m,使用均勻長(zhǎng)方體網(wǎng)格剖分,網(wǎng)格長(zhǎng)寬高分別為500 m、500 m和200 m,光滑約束反演正則化因子選擇5.0,MS聚焦約束反演正則化因子選擇0.001,對(duì)稱(chēng)多項(xiàng)式約束反演光滑項(xiàng)權(quán)重為5.0、多項(xiàng)式項(xiàng)權(quán)重為1000,反演結(jié)果見(jiàn)圖8、9和10.其中,圖8為光滑反演結(jié)果,圖9為最小支撐（ε=0.1）反演結(jié)果,圖10為對(duì)稱(chēng)多項(xiàng)式約束（參考物性為m*=[0,0.1,-0.1]）的反演結(jié)果.各圖中（a）為反演模型,（b）為反演物性統(tǒng)計(jì)圖,（c）為誤差曲線(xiàn)圖.

圖6 兩個(gè)異常體合成模型反演（a） 光滑模型約束反演密度斷面;（b） 最小支撐模型約束反演密度斷面;（c） 對(duì)稱(chēng)多項(xiàng)式約束反演密度斷面;（d） 光滑模型約束反演誤差曲線(xiàn);（e） 最小支撐模型約束反演誤差曲線(xiàn);（f） 對(duì)稱(chēng)多項(xiàng)式約束反演誤差曲線(xiàn);（g） 響應(yīng)擬合曲線(xiàn).Fig.6 Synthetic model inversion with two abnormal blocks（a） The density section of smooth model constraint inversion;（b） The density section of minimum support model constraint inversion;（c） The density section of symmetric polynomial model constrain inversion;（d）The smooth model constraint inversion error curves;（e） The minimum support model constraint inversion error curves;（f） The symmetric polynomial model constraint inversion error curves;（g） Data fitting curves.

對(duì)比分析三種方案的反演結(jié)果,可以發(fā)現(xiàn):光滑約束反演,在異常體兩側(cè)存在假異常,反演邊界模糊,物性小于真實(shí)物性；最小支撐反演結(jié)果異常體更緊湊,異常體規(guī)模與真實(shí)模型更接近,但物性范圍僅為[-0.07,0.07] g·cm-3，與真實(shí)物性依舊存在差異；基于對(duì)稱(chēng)多項(xiàng)式約束的反演結(jié)果物性值與真實(shí)模型更接近,且邊界更清晰,但縱向分辨率略低,分析可能與光滑約束權(quán)重選擇有關(guān).由物性統(tǒng)計(jì)直方圖可見(jiàn),對(duì)稱(chēng)多項(xiàng)式約束形成三個(gè)物性統(tǒng)計(jì)中心,反映對(duì)稱(chēng)多項(xiàng)式約束有較好的聚焦效果.分析誤差曲線(xiàn)圖,可以發(fā)現(xiàn)三種方案的反演誤差函數(shù)隨迭代次數(shù)逐漸下降并最終趨于平緩.

圖7 兩個(gè)異常體三維剩余密度模型Fig.7 Synthetic 3D model with two residual density blocks

圖8 光滑約束反演結(jié)果（a） 反演剩余密度模型;（b） 反演物性統(tǒng)計(jì)圖;（c） 反演誤差曲線(xiàn).Fig.8 Smooth constraint inversion results（a） Inversion residual density model;（b） Statistical chart of inversion physical property;（c） Inversion error curves.

圖9 最小支撐約束反演結(jié)果（a） 反演剩余密度模型;（b） 反演物性統(tǒng)計(jì)圖;（c） 反演誤差曲線(xiàn).Fig.9 Minimum support constraint inversion results（a） Inversion residual density model;（b） Statistical chart of inversion physical property;（c） Inversion error curves.

圖10 對(duì)稱(chēng)多項(xiàng)式約束反演結(jié)果（a） 反演剩余密度模型;（b） 反演物性統(tǒng)計(jì)圖;（c） 反演誤差曲線(xiàn).Fig.10 Symmetric polynomial constrain inversion results（a） Inversion residual density model;（b） Statistical chart of inversion physical property;（c） Inversion error curves.

綜合分析圖4、圖6、圖8、圖9和圖10,可以發(fā)現(xiàn),基于對(duì)稱(chēng)多項(xiàng)式約束的反演可以將任意數(shù)量的物性統(tǒng)計(jì)信息融入反演,并能獲得類(lèi)似聚焦反演的效果,反演物性處于給定參考值范圍內(nèi).

4 實(shí)測(cè)數(shù)據(jù)反演

為了驗(yàn)證對(duì)稱(chēng)多項(xiàng)式約束函數(shù)的實(shí)用性,選擇某測(cè)線(xiàn)實(shí)測(cè)資料進(jìn)行反演驗(yàn)證.測(cè)線(xiàn)長(zhǎng)9 km,點(diǎn)距200 m,已知研究深度內(nèi)存在侏羅系、志留系、奧陶系地層,其中侏羅系以凝灰?guī)r為主,標(biāo)本測(cè)試密度為2.5～2.55 g·cm-3,志留系以砂巖、砂質(zhì)泥巖及頁(yè)巖為主，密度為2.65～2.67 g·cm-3,奧陶系以鈣質(zhì)泥巖、灰?guī)r為主，密度為2.68～2.7 g·cm-3.

分別使用光滑約束、最小支撐約束（ε=0.05）和對(duì)稱(chēng)多項(xiàng)式約束（參考物性為m*=[2.53,2.66,2.7]）進(jìn)行反演.計(jì)算采用矩形網(wǎng)格剖分,單元中心節(jié)點(diǎn)位置示意圖見(jiàn)圖11.為便于與先驗(yàn)物性對(duì)比,采用絕對(duì)密度表示反演結(jié)果,見(jiàn)圖12.其中,圖12a、12b、12c分別為光滑、最小支撐、對(duì)稱(chēng)多項(xiàng)式約束反演結(jié)果,其中黑線(xiàn)是根據(jù)反演模型給出的推測(cè)地層界面；圖12d、12e、12f為對(duì)應(yīng)的反演誤差曲線(xiàn),圖12g為響應(yīng)曲線(xiàn).

對(duì)比圖12a、12b和12c可以發(fā)現(xiàn)，光滑反演可以獲得平滑過(guò)渡的反演模型,但反演物性相對(duì)已知物性普遍偏小,邊界分辨能力很差.最小支撐反演可以獲得淺層物性分布信息,但對(duì)深部信息反演存在不足.對(duì)稱(chēng)多項(xiàng)式約束反演結(jié)果具有清晰邊界且反演物性接近給定的先驗(yàn)物性信息.分析圖12d、12e、12f和12g可以發(fā)現(xiàn),三種反演方案數(shù)據(jù)誤差均穩(wěn)定下降并趨于平緩,觀測(cè)數(shù)據(jù)與計(jì)算數(shù)據(jù)擬合程度很高.

為了直觀比較三種方案的反演物性分布,繪制反演物性的頻數(shù)分布直方圖,圖13a、13b和13c分別為光滑、最小支撐和對(duì)稱(chēng)多項(xiàng)式約束反演的物性統(tǒng)計(jì)結(jié)果,圖中的藍(lán)色虛線(xiàn)為先驗(yàn)的密度統(tǒng)計(jì)值.

從反演的頻數(shù)分布直方圖中可以發(fā)現(xiàn),光滑反演物性存在明顯的過(guò)渡,物性連續(xù)性好.MS反演結(jié)果物性分布相對(duì)孤立,反演結(jié)果存在明顯的“分塊”.但光滑反演和MS反演結(jié)果物性分布峰值與先驗(yàn)的物性統(tǒng)計(jì)信息存在明顯偏差.基于對(duì)稱(chēng)多項(xiàng)式約束的反演結(jié)果物性分布峰值與先驗(yàn)統(tǒng)計(jì)物性更加接近,一定程度證明了對(duì)稱(chēng)多項(xiàng)式的約束能力.

實(shí)測(cè)數(shù)據(jù)反演表明對(duì)稱(chēng)多項(xiàng)式約束函數(shù)可以充分利用先驗(yàn)物性信息,通過(guò)將任意多個(gè)先驗(yàn)物性信息融入目標(biāo)函數(shù),一定程度上克服了復(fù)雜區(qū)域先驗(yàn)統(tǒng)計(jì)信息應(yīng)用困難的問(wèn)題.

5 結(jié)論與展望

本文提出了一種基于對(duì)稱(chēng)多項(xiàng)式形式的模型約束函數(shù),并將其應(yīng)用于重力數(shù)據(jù)正則化反演,通過(guò)模型試驗(yàn)和實(shí)測(cè)數(shù)據(jù)處理，形成以下結(jié)論.

（1）對(duì)稱(chēng)多項(xiàng)式可以將任意多個(gè)先驗(yàn)物性信息融入目標(biāo)函數(shù),一定程度上克服了復(fù)雜區(qū)域先驗(yàn)統(tǒng)計(jì)信息應(yīng)用困難的問(wèn)題,提高了反演的實(shí)際應(yīng)用效果.

圖11 模型離散單元中心節(jié)點(diǎn)位置Fig.11 Element center position of discrete mesh

（2）對(duì)稱(chēng)多項(xiàng)式模型約束函數(shù)極小值取值明確,反演具有較強(qiáng)的聚焦能力.

（3）在給定參考物性時(shí),目標(biāo)函數(shù)在最小與最大參考物性范圍外，對(duì)稱(chēng)多項(xiàng)式約束函數(shù)取值快速增加,表明該約束隱含有物性范圍約束能力.

（4）對(duì)稱(chēng)多項(xiàng)式模型約束函數(shù)的梯度與海森矩陣具有解析形式,反演目標(biāo)函數(shù)優(yōu)化求解簡(jiǎn)單.

此外,智能算法在地球物理中的應(yīng)用受到了越來(lái)越多的關(guān)注,以聚類(lèi)算法為代表的反演方法已成功應(yīng)用于地球物理反演.大多數(shù)聚類(lèi)算法可以有效地獲得聚類(lèi)中心,但沒(méi)有明確的代價(jià)函數(shù),因此當(dāng)前只有具有明確代價(jià)函數(shù)的模糊C均值聚類(lèi)在地球物理反演中獲得了廣泛的應(yīng)用.基于對(duì)稱(chēng)多項(xiàng)式約束可以有效地結(jié)合各種聚類(lèi)算法獲得的聚類(lèi)中心,實(shí)現(xiàn)聚類(lèi)反演,如何將對(duì)稱(chēng)多項(xiàng)式與聚類(lèi)算法相結(jié)合將是今后的研究方向.

致謝感謝審稿專(zhuān)家對(duì)本文詳細(xì)專(zhuān)業(yè)的修改及編輯部的支持.

附錄

設(shè)有一函數(shù)為

（A1）

可以發(fā)現(xiàn):

（A2）

因?yàn)?/p>

1≤i

（A3）

所以Φ為對(duì)稱(chēng)多項(xiàng)式,有

（A4）

其中,σi為基本對(duì)稱(chēng)多項(xiàng)式,其表達(dá)式如下:

（A5）

圖12 實(shí)測(cè)數(shù)據(jù)反演（a） 光滑模型約束反演密度斷面;（b） 最小支撐模型約束反演密度斷面;（c） 對(duì)稱(chēng)多項(xiàng)式約束反演密度斷面;（d） 光滑模型約束反演誤差曲線(xiàn);（e） 最小支撐模型約束反演誤差曲線(xiàn);（f） 對(duì)稱(chēng)多項(xiàng)式約束反演誤差曲線(xiàn);（g） 響應(yīng)擬合曲線(xiàn).Fig.12 Field data inversion（a） The density section of smooth model constraint inversion;（b） The density section of minimum support model constraint inversion;（c） The density section of symmetric polynomial model constrain inversion;（d） The smooth model constraint inversion error curves;（e） The minimum support model constraint inversion error curves;（f） The symmetric polynomial model constraint inversion error curves;（g） Response fitting curves.

圖13 重力反演密度統(tǒng)計(jì)圖（a） 光滑反演;（b） 最小支撐反演;（c） 對(duì)稱(chēng)多項(xiàng)式約束反演.Fig.13 Gravity inversion density statistics（a） Smooth inversion;（b） Minimum support inversion;（c） Symmetric polynomial constrained inversion.

代入φ（A+Bx）有:

φ（A+Bx）=c0xk+c1xk-1+…+ck-1x1+ck.

（A6）

根據(jù)函數(shù)定義結(jié)合基本對(duì)稱(chēng)多項(xiàng)式定義可知,c的通解為

（A7）

（A8）

最終,函數(shù)f可以表示為如下的多項(xiàng)式函數(shù)乘積形式:

f=（c0xk+c1xk-1+…+ck-1x1+ck）

×（c0xk+c1xk-1+…+ck-1x1+ck）.

（A9）