重力異常物性正則化反演研究

周宇軒， 雷 宛

(成都理工大學(xué) 地球物理學(xué)院，成都 610059)

?

重力異常物性正則化反演研究

周宇軒， 雷 宛

(成都理工大學(xué) 地球物理學(xué)院，成都 610059)

經(jīng)典位場(chǎng)物性反演是建立在數(shù)與模型之間的線性關(guān)系，利用離散正演公式將每個(gè)剖分單元相對(duì)觀測(cè)點(diǎn)的核函數(shù)進(jìn)行計(jì)算，以尋求觀測(cè)數(shù)據(jù)與核函數(shù)有最好的相關(guān)性。該類(lèi)方法在計(jì)算過(guò)程最大的問(wèn)題就是在解巨大超定方程組時(shí)，導(dǎo)致求解極不穩(wěn)定。在前人的研究基礎(chǔ)上，提出一種基于先驗(yàn)密度模型約束正則化條件下的重力異常視密度求解方法，通過(guò)L曲線法選取最佳的參數(shù)，使求解過(guò)程中有較高的穩(wěn)定性和求解精度。這里以一理論模型作為觀測(cè)數(shù)據(jù)，并引入噪聲作為觀測(cè)誤差，通過(guò)算法的對(duì)比表明，該方法有很高的計(jì)算精度和穩(wěn)定性，對(duì)于需要求解大型超定矩陣的三維模型，在給定先驗(yàn)初始模型的條件下，該方法仍然能夠取得很好的效果。

重力異常； 視密度； 正則化； 先驗(yàn)?zāi)Ｐ停?約束反演

0 引言

重力勘探是地球物理勘探的主要分支之一。目前國(guó)內(nèi)主要通過(guò)地面測(cè)量來(lái)獲取研究區(qū)各點(diǎn)重力場(chǎng)觀測(cè)值，再通過(guò)必要地改正得到重力異常值[1]。重力勘探方法作為重要的勘探手段，已經(jīng)廣泛應(yīng)用于油氣勘探、金屬礦勘查、地質(zhì)填圖及工程與環(huán)境地球物理勘探等領(lǐng)域。隨著儀器、計(jì)算機(jī)、導(dǎo)航定位等技術(shù)的不斷發(fā)展，重力探測(cè)技術(shù)應(yīng)用已經(jīng)擴(kuò)展到衛(wèi)星、航空、海洋等方面。相應(yīng)地，針對(duì)衛(wèi)星、航空、海洋重力實(shí)際資料大面積、大數(shù)據(jù)量以及快速估算解釋的要求[2]，三維正反演也逐漸地發(fā)展起來(lái)。

反演是重力資料定量解釋中的重要環(huán)節(jié)之一[3]。與其他地球物理方法一樣，重力資料解釋的目的在于通過(guò)地面或航空等實(shí)測(cè)數(shù)據(jù)利用某種手段推算出地下的密度分布規(guī)律，從而達(dá)到尋找目標(biāo)地質(zhì)體、推斷地質(zhì)體埋深和地質(zhì)體分布范圍的目的。雖然能很容易找到滿足數(shù)據(jù)擬合空間的解模型，但是由于重力場(chǎng)的體積效應(yīng)、有限觀測(cè)數(shù)據(jù)的不準(zhǔn)確性及反演問(wèn)題的欠定性等幾種因素的存在，反演結(jié)果往往產(chǎn)生多解性，從而很難得到一個(gè)符合實(shí)際情況的解。這就需要在反演過(guò)程中添加一些約束條件及先驗(yàn)信息，最大限度地保證反演結(jié)果的正確性，如王薇等[4]提出的地震波形反演的稀疏約束正則化方法；舒夢(mèng)埕等[5]提出儲(chǔ)層重力密度反演后驗(yàn)約束正則化方法等；Tikhonov[6-7]正則化算法使反演結(jié)果更加可靠。對(duì)于反演模型可用場(chǎng)源目標(biāo)體的形狀或遍歷地下半空間的物性值大小來(lái)表示，這就產(chǎn)生了兩種完全不同的反演方法：①幾何反演；②物性反演。幾何反演是在地下半空間場(chǎng)源體給定物性參數(shù)大小的基礎(chǔ)上，利用地面觀測(cè)異常來(lái)擬合幾何體(如多邊形或多面體)形狀大小，通過(guò)幾何體的形狀大小來(lái)模擬目標(biāo)體的分布規(guī)律。物性反演是將觀測(cè)區(qū)域地下半空間離散化成規(guī)則的網(wǎng)格單元，通過(guò)反演方法確定各離散單元的物性值，由物性的分布確定場(chǎng)源的實(shí)際分布情況。早期以形態(tài)反演為主，主要應(yīng)用于沉積基底界面反演以及形態(tài)簡(jiǎn)單的場(chǎng)源體反演等。隨著計(jì)算機(jī)計(jì)算能力地提高，物性反演后來(lái)居上且從二維發(fā)展為三維，逐漸成為國(guó)內(nèi)、外重力反演的主要方向[8]。

1 密度模型正演原理

建立如圖1所示的二維密度模型和觀測(cè)系統(tǒng)，在地表(z=0)處觀測(cè)和接收數(shù)據(jù)，d(i)為在地面第i個(gè)觀測(cè)點(diǎn)所得的數(shù)據(jù)，mj為地下介質(zhì)剖分后第j個(gè) 單元的密度(通常是未知的)。

圖1 密度模型及觀測(cè)系統(tǒng)示意圖Fig.1 Schematic diagram of density model and observation system

根據(jù)圖1的重力觀測(cè)模式，d(i)與mj的定量關(guān)系如式(1)所示。

d(1)=f1,1m1+f1,2m2+f1,3m3…+f1,mmm

(1)

其中：f1,j為重力場(chǎng)觀測(cè)網(wǎng)格化后物理的距離函數(shù)。

剖分矩形截面重力正演公式，如式(2)所示。

(2)

其中：dx、dy分別為地下網(wǎng)格的長(zhǎng)度和寬度；xk,zk為地下觀測(cè)點(diǎn)坐標(biāo)；G表示牛頓萬(wàn)有引力常數(shù)。

為了與實(shí)際重力異常相結(jié)合，此處計(jì)算的是每一個(gè)場(chǎng)源的引力在鉛垂方向上的投影，最后將它們相加。需要計(jì)算一條測(cè)線上面的重力異常值，可以用式(3)的矩陣表示。

(3)

式(3)可以化簡(jiǎn)為式(4)：

F m=d

(4)

其中：F為表示表示觀測(cè)點(diǎn)與場(chǎng)源的核函數(shù)；m為每個(gè)網(wǎng)格單元的密度值；d為測(cè)線上面的重力異常觀測(cè)值。

圖2 三維密度模型和觀測(cè)系統(tǒng)Fig.2 3D density model and observation system

三維模型通常將地下的地質(zhì)體劃分為多個(gè)立方單元體，密度模型和觀測(cè)系統(tǒng)如圖2所示。剖分立方體正演公式,如式(5)所示。

Δg= -G m|||[ξ ln(η+ρ)+

(5)

其中：η、ζ、ξ為網(wǎng)格單元的x、y、z坐標(biāo)值；ρ為觀測(cè)點(diǎn)到單元網(wǎng)格的距離，其表達(dá)式如式(6)所示。

(6)

若F等于式(7)所示的公式：

(7)

則可以將式(5)表示為：Δg=F m，為每個(gè)網(wǎng)格單元體的密度值。

2 正則化反演原理

在反演過(guò)程中，正演擬合參與運(yùn)算的剖分單元數(shù)往往大于觀測(cè)數(shù)據(jù)個(gè)數(shù)，特別是實(shí)際觀測(cè)數(shù)據(jù)總是存在觀測(cè)誤差，導(dǎo)致反算F方程通常是一個(gè)超定矩陣，無(wú)法穩(wěn)定線性的求出精確的解，為了使解更加穩(wěn)定、準(zhǔn)確，一般需要求解一個(gè)最小二乘意義下的近似解式(8)。

‖F(xiàn) m-d‖2→min

(8)

最小二乘本質(zhì)是線性增量修正，由于修正量與解的收斂不確定性，導(dǎo)致求到的解會(huì)出現(xiàn)由穩(wěn)定收斂到逐步發(fā)散的現(xiàn)象，為克服上述問(wèn)題，作者運(yùn)用Tikhonov正則化思想設(shè)計(jì)的正則化反演模型：

(9)

其中：α是正則化參數(shù)；m0為先驗(yàn)的密度模型；通過(guò)L曲線法選取最佳的α參數(shù)。

為降低反演的多解性，需要加入一個(gè)先驗(yàn)?zāi)Ｐ偷募s束。具體求解過(guò)程推導(dǎo)如下：

(10)

利用最小二乘原理將上訴問(wèn)題轉(zhuǎn)化為正則方程，得到的結(jié)果如式(11)所示。

(GTG+α2I)m=GTd+α2Im0

(GTG+α2I)(m-m0)=GT(d-Gm0)

(11)

令m*=m-m0,d*=d-Gm0，將其代入式(11)中，則可得到式(12)。

(GTG+α2I)m*=GTd*

(12)

求解上述方程，對(duì)矩陣G進(jìn)行奇異值(SVD)分解：

(VSTUTUSVT+α2I)m*=VSTUTd*

(VSTSVT+α2I)m*=VSTUTd*

假設(shè)x*=VTm*，m=Vx,VVT=I

則有：

(STS+α2I)x*=STUTd*

(13)

從而可求到：

(14)

令m-m0=m*=Vx*,d=d*+Gm0，代入式(14)，可以得到：

(15)

其中：S、U、V由矩陣G的SVD分解得到；Gm0、m0都是已知量，因此，作者采用L-curve曲線方法[9](圖3)求解。

圖3 L-curve曲線Fig.3 L-curve curve

正則化參數(shù)α的取值規(guī)則是從小取到大，依次用不同α算出的‖G m-d‖、‖m-m0‖取對(duì)數(shù)顯示成L曲線。通過(guò)L曲線可以看出，當(dāng)α很小時(shí)‖G m-d‖的值較小，而‖m-m0‖較大，這說(shuō)明取較小的α?xí)r數(shù)據(jù)擬合效果較好，但反演結(jié)果與初始后驗(yàn)?zāi)Ｐ拖啾葦_動(dòng)較大。當(dāng)α取較大時(shí)擬合誤差‖G m-d‖變大，而反演結(jié)果與后驗(yàn)?zāi)Ｐ偷木嚯x‖m-m0‖變小。因此選擇L曲線曲率最大處作為最優(yōu)有α參數(shù)，這樣既能平衡擬合誤差，又能兼顧后驗(yàn)?zāi)Ｐ蚚10-12]。

求出α以后，代入到公式(15)中，即可求出m。

3 數(shù)值試驗(yàn)

設(shè)計(jì)一個(gè)二維模型(圖4)，其參數(shù)設(shè)置為：中間方形的密度為ρ=2 g/cm3，旁側(cè)的密度設(shè)為 g/cm3；長(zhǎng)寬分別為300 m、300 m，中心埋深為500 m；數(shù)據(jù)網(wǎng)格為10 m×20 m，測(cè)線上的數(shù)據(jù)點(diǎn)采樣為2 000個(gè)，計(jì)算重力正演異常(圖5)。

為模擬觀測(cè)數(shù)據(jù)的隨機(jī)觀測(cè)誤差，對(duì)該重力異常進(jìn)行隨機(jī)噪聲添加(圖6)。

對(duì)圖6所得重力異常數(shù)據(jù)進(jìn)行最小二乘優(yōu)化反演，得到結(jié)果如圖7所示，反演的視密度比較接近真實(shí)值2，但從視密度成像幾何分布特征來(lái)看，與設(shè)計(jì)理論模型圖4存在很明顯的差異。在給定先驗(yàn)密度模型約束條件下，采用正則化方法，對(duì)含噪聲重力異常數(shù)據(jù)進(jìn)行反演，其視密度成像結(jié)果如圖8所示。從圖8中可以看出,無(wú)論密度反演值還是密度分布成像特征，均較最小二乘法有很好的改進(jìn)，與設(shè)計(jì)模型(圖4)有很好的對(duì)比性，表明基于正則化條件下的視密度反演算法，在噪聲條件下，反演解未陷于局部收斂，其解有很高的穩(wěn)定性和求解精度。

圖5 正演重力異常Fig.5 Positive gravity anomaly

圖6 加噪音正演重力異常Fig.6 Positive gravity anomaly with noise

將上述算法推廣到三維模型，設(shè)計(jì)如圖9的模型，其中模型密度為ρ=2 g/cm3，圍巖密度為ρ=0.4 g/cm3；長(zhǎng)、寬、高分別為300 m、300 m、300 m，中心埋深為500 m；數(shù)據(jù)網(wǎng)格為10 m×10 m×10 m，數(shù)據(jù)節(jié)點(diǎn)為1 000個(gè)；觀測(cè)面數(shù)據(jù)采樣為50×50=2 500個(gè)；其正演重力異常見(jiàn)圖10。 設(shè)置如圖11所示的先驗(yàn)初始模型，采用該先驗(yàn)?zāi)Ｐ瓦M(jìn)行正則化條件下的重力異常數(shù)據(jù)三維物性反演，得到的結(jié)果如圖12所示。從圖12中可以看出，反演的視密度大小及幾何形狀分布與圖9有很高的吻合度。對(duì)比二維和三維反演的計(jì)算效率(表1)。

由表1可以看出，三維正則化物性反演雖然在迭代時(shí)間上消耗較多，但迭代次數(shù)相對(duì)二維反演較小，且能達(dá)到同等的效果。表明該方法在三維大型數(shù)據(jù)運(yùn)算中仍具有很高的求解精度和穩(wěn)定的求算解。

圖7 最小二乘算法反演結(jié)果Fig.7 Least square algorithm inversion results

圖8 加噪音正則化算法反演結(jié)果Fig.8 Inversion results with noise

圖9 三維正演模型Fig.9 3D forward modeling

圖10 三維正演平面圖Fig.10 3D positive results plan

圖11 先驗(yàn)初始模型Fig.11 A priori initial model

圖12 反演結(jié)果Fig.12 Inversion results

網(wǎng)格點(diǎn)數(shù)迭代次數(shù)時(shí)間消耗二維反演10*20150.34三維正則化反演10*10*10100.53

4 結(jié)論

1)重力異常物性反演在解巨大超定方程組時(shí)，導(dǎo)致求解極不穩(wěn)定，作者提出基于正則化條件下的約束反演，對(duì)于解決這類(lèi)問(wèn)題有很好的效果。

2)二維噪聲模型和三維密度模型數(shù)值模擬結(jié)果表明，基于正則化條件下的反演算法，特別是在給定先驗(yàn)初始密度模型約束條件下，在求解的精度和求解的穩(wěn)定性有其相對(duì)優(yōu)勢(shì)。

[1] GROETSCH C W. Inverse Problems[M]. Washington: The Mathematical Association of America, 1999.

[2] 陳召曦,孟小紅,郭良輝.重磁數(shù)據(jù)三維物性反演方法進(jìn)展[J].地球物理學(xué)進(jìn)展,2012,27(12)：503-511.

CHEN Z X, MENG X H, GUO L H. Review of 3D property inversion of gravity and magneticdata[J]. 2012,27(12)：503-511.(In Chinese)

[3] 陳桂廷,羅強(qiáng),周宇軒,等.三維位場(chǎng)正則化向下延拓成像技術(shù)[J].科技創(chuàng)新導(dǎo)報(bào),2014(26):46-47.

CHEN G T, LUO Q, ZHOU Y X ,et al.Down-ward extension imaging techn ique of threedimensional potential field[J]. Scienc Science and Technology Innovation Herald,2014(26):46-47. (In Chinese)

[4] 肖庭延. 反問(wèn)題的數(shù)值解法[M]. 北京：科學(xué)出版社, 2003.

XIAO T Y.Numerical solution of the inverse problem[M].Beijing:Science Press, 2003. (In Chinese)

[5] 舒夢(mèng)埕,王彥飛.儲(chǔ)層重力密度反演后驗(yàn)約束正則化方法[J].地球物理學(xué),2015,58(6):2079-2086.SHU M C,WANG Y F.The posterior constrained regularization method for reservoirdensity inversion[J].ActaGeophysicaSinica,2015,58(6):2079-2086.(In Chinese)

[6] 徐新禹，李建成，王正濤,等.Tiknohov正則化方法在GOCE重力場(chǎng)中的模擬研究[J].測(cè)繪學(xué)報(bào)，2010,39(5)：466-470.

XU X Y, LI J C, WANG Z T, et al.The simulation research on the Tikhonov regularization applied in Gravity Field[J].Acta Geodaeticaet Cartograph-icaSinica, 2010,39(5)：466-470. (In Chinese)

[7] 傅初黎，洪芳，向團(tuán).不適定問(wèn)題的Tikhonov正則化方法[J].計(jì)算數(shù)學(xué),2006,28(3):237-246.

FU C L,HONG F,XIANG T.Iterated Tikhonov regularization for ill-poseproblems[J].Mathematica Numerica Sinica,2006,28(3)：237-246.(In Chinese)

[8] 王邦華,王理.重磁位場(chǎng)的正則化向下延拓[J].物探化探計(jì)算技術(shù),1998,20(1):30-35.

WANG B H. WANG L.A new normalized thod of down ward extrapolation for potential field[J]. Computing Techniques for Geophysical and Geochemical Exploration,1998,20(1):30-35. (In Chinese)

[9] TVEITO A, LANGTANGEN H P, NIELSEN B F,et al. Parameter estimation and Inverse problems[M].Pittsburgh:Academix Press,2013.

[10]馮寶賓. 不適定問(wèn)題的正則化解算方法設(shè)計(jì)及應(yīng)用[D]. 成都：成都理工大學(xué), 2010.

FENG B B.Design and application of theregularization method for solving ill posed problems[D].Chengdu：Chengdu University of Technology,2010.(In Chinese)

[11]ARSENIN, IA. V. Solutions of ill-posed problems [M]. Washington:Winston,1977.

[12]劉銀萍,王祝文,杜曉娟,等. 基于Tikhonov正則化算法的重力數(shù)據(jù)三維約束反演[J]. 地球物理學(xué)報(bào), 2013, 56(5):1650-1659.

LIU Y P, WANG Z W, DU X J, et al. 3D constrained inversion of gravity data based on Tikhonov regularization algorithm[J].Chinese Journal Geophysics,2013,56(5):1650-1659. (In Chinese)

Study on physical properties regularization inversion of gravity anomaly

ZHOU Yu-xuan， LEI Wan

(College of Geophysics,Chengdu University of Technology,Chengdu 610059，China)

Classical potential field physical properties inversion is based on the linear relationship between the number and the model.We calculate each cross-sectional unit relative to observation point kernel function by using discrete forward formula to seek the best correlation in the observational data and the kernel function. The biggest problem of this kind of method calculation is unstable system to solving when we solve the large set of equations. On the basis of previous research,a kind of density algorithm is put forward based on the prior density model regularization constraint conditions of gravity anomaly in this study.The solution process has high stability and accuracy by using the L-curve method to select optimal parameters. In this paper,we establish a theoretical model as the observation data and introduce a noise as the observation error.Calculation results show that the method has high accuracy and stability through the comparison of different algorithms.Given a priori model of the initial conditions,the method can still achieve good results for the need for solving large over determined matrix of 3D model.

gravity anomaly; apparent density; regularization; a priori model; constrained inversion

2016-04-13 改回日期：2016-06-14

中國(guó)地質(zhì)調(diào)查項(xiàng)目(12130113095100)

周宇軒(1992-),男，碩士，從事地球物理勘探方法研究，E-mail:zhouyuxuan@163.com。

1001-1749(2016)05-0579-05

P 631.1

A

10.3969/j.issn.1001-1749.2016.05.01