林寶澤　肖　鋒　王明常

(吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林長(zhǎng)春 130000)

1　引言

地球物理的反演問(wèn)題主要包含兩個(gè)方面：一是地質(zhì)體的幾何參數(shù)反演[1,2]； 二是地質(zhì)體的物性參數(shù)反演[3,4]。針對(duì)地質(zhì)體的幾何參數(shù)反演，建模一般分為三個(gè)方案：一是將場(chǎng)源表示為一系列垂直并列僅厚度未知的棱柱體[5-8]； 二是將場(chǎng)源表示為頂點(diǎn)坐標(biāo)未知的多邊形或多面體[9-12]； 三是將包含場(chǎng)源的地下空間剖分成尺寸已知而密度未知的基本矩形單元[13-15]。

徑向反演是一種最早由Silva等[16]提出的場(chǎng)源幾何參數(shù)反演方法； 賈真[17]研究了基于二維模型的重磁梯度分量徑向反演算法； Vanderlei等[18,19]將徑向反演方法擴(kuò)展至重力異常的三維反演以及重力梯度張量的三維聯(lián)合反演。

建模對(duì)反演來(lái)說(shuō)至關(guān)重要，模型過(guò)于復(fù)雜會(huì)導(dǎo)致解的不穩(wěn)定性，可以通過(guò)Tikhonov正則化方法引入約束條件解決這一問(wèn)題。徑向反演的優(yōu)點(diǎn)在于將多邊形頂點(diǎn)位置用極坐標(biāo)表示，反演參數(shù)(即多邊形各頂點(diǎn)對(duì)應(yīng)的半徑)具有同一量綱且反演參數(shù)相對(duì)于直角坐標(biāo)系減少一半。利用Tikhonov正則化引入相關(guān)約束條件(場(chǎng)源幾何形態(tài)等先驗(yàn)信息)時(shí)能取得較好的效果。

徑向反演方法雖然已經(jīng)發(fā)展得較為完善，但在國(guó)內(nèi)關(guān)于該方法的研究和討論幾乎還是空白。此外，該方法的研究多見(jiàn)于模型實(shí)驗(yàn)，對(duì)于影響反演結(jié)果的正則化參數(shù)選取更是缺乏討論。本文介紹了二維徑向反演方法的基本原理，且在此基礎(chǔ)上將自適應(yīng)正則化因子選取引入該方法，并將該方法應(yīng)用到實(shí)際資料處理，取得了較好的效果。

2　方法原理

2.1　正演方法

徑向模型的中心點(diǎn)O(x0,z0)必須設(shè)定在場(chǎng)源內(nèi)部，并給定徑向多邊形的頂點(diǎn)數(shù)M，將圓周平均分成M份，此時(shí)在直角坐標(biāo)系下，多邊形第k個(gè)頂點(diǎn)的坐標(biāo)(xk,zk)可以表示為

(1)

式中：rk為第k個(gè)頂點(diǎn)對(duì)應(yīng)的半徑;θk為第k個(gè)頂點(diǎn)對(duì)應(yīng)的角度(弧度,沿x軸正方向順時(shí)針旋轉(zhuǎn)為正)。模型及參數(shù)如圖1所示。

圖1　徑向模型示意圖

從圖1不難推斷，當(dāng)頂點(diǎn)數(shù)M足夠大時(shí)，徑向模型能取得很好的近似效果。但是從式(1)來(lái)看，應(yīng)用徑向模型時(shí)存在以下限制：一是包含場(chǎng)源的區(qū)域必須是單連通域；二是從模型中心點(diǎn)出發(fā)的射線與場(chǎng)源輪廓線有且僅有一個(gè)交點(diǎn)。

反演過(guò)程中均勻多邊形截面二維體的重力異常正演采用文獻(xiàn)[20]給出的公式

(2)

式中：G為萬(wàn)有引力常數(shù)； Δρ為密度差；Rk、Rk+1分別為測(cè)點(diǎn)到多邊形第k和k+1個(gè)頂點(diǎn)的距離；αk、αk+1分別為測(cè)點(diǎn)與第k和k+1個(gè)頂點(diǎn)連線與水平方向的夾角；φk為多邊形第k條邊與水平方向的夾角，如圖2所示。

圖2　多邊形截面二度體正演模型示意型圖[20]

2.2　約束反演

設(shè)重力異常實(shí)測(cè)數(shù)據(jù)構(gòu)成的N×1維向量為d,徑向多邊形各頂點(diǎn)對(duì)應(yīng)的半徑構(gòu)成的M×1維向量為r,由模型r得到的正演重力異常為N×1維向量g(r)。為了保證反演問(wèn)題的適定性，引入約束條件后，該問(wèn)題轉(zhuǎn)化為求解最優(yōu)化問(wèn)題

minψ=‖d-g(r)‖22+μ(Φ1+Φ2)

(3)

式中

(4)

(5)

(6)

(7)

k=1,2,…,M；j=1,2,…,J

其中：ψ為總目標(biāo)函數(shù)；μ為正則化因子；r0為初始模型多邊形截面各頂點(diǎn)對(duì)應(yīng)的半徑構(gòu)成的M×1維向量；Φ1為相對(duì)鄰近性約束，其作用是使模型r的各個(gè)分量彼此趨于一致，但不能單獨(dú)使用；當(dāng)W為單位矩陣時(shí)，Φ2為絕對(duì)鄰近性約束，作用是穩(wěn)定解；當(dāng)W為非單位矩陣時(shí)，Φ2為會(huì)聚方向約束，使場(chǎng)源物性沿著某一特定方向集中，本文中將場(chǎng)源質(zhì)量會(huì)聚于某一方向，能起到引入與場(chǎng)源幾何形態(tài)有關(guān)的先驗(yàn)信息和穩(wěn)定解的作用；βj為指定的會(huì)聚方向；J表示會(huì)聚方向的總數(shù)；ε為較小的正數(shù)，與場(chǎng)源最小尺寸和最大尺寸之比有關(guān)。顯然，相對(duì)鄰近性約束與會(huì)聚方向約束是相互矛盾的，絕對(duì)鄰近性約束是會(huì)聚方向約束的特殊形式。除上述約束條件外，還有不等性約束和凸性約束。不等性約束的表達(dá)式為

(8)

不等性約束的作用是限制場(chǎng)源的幾何尺寸和保證頂點(diǎn)z坐標(biāo)為非負(fù)。凸性約束就是要在每次迭代中滿足

(9)

式中ck為中心點(diǎn)O到半徑rk與多邊形第k-1、k+1個(gè)頂點(diǎn)連線交點(diǎn)之間的距離，如圖3所示。

圖3　凸性約束示意圖

由式(2)可以看出，g(r)是關(guān)于r的非線性函數(shù)，直接求解式(3)較困難，故采用Gauss-Newton迭代法將該非線性反演問(wèn)題線性化。因此，在r0處對(duì)g(r)進(jìn)行Taylor展開(kāi)并取一級(jí)近似，則有

=g(r0)+JΔr

(10)

式中J為Jacobi矩陣

(11)

令r=r0+Δr,并將式(10)代入式(3)，則

ψ=‖[d-g(r0)]-JΔr‖22+μ‖B(r0+Δr)‖22+

μ‖WΔr‖22

=μ(WΔr)TWΔr+μ[B(r0+Δr)]TB(r0+Δr)+

{[d-g(r0)]-JΔr}T{[d-g(r0)]-JΔr}

(12)

對(duì)式(12)兩端求導(dǎo)，可得

2μ(BTBΔr+BTBr0+WTWΔr)

(13)

令上式等于0，則迭代修正量Δr的表達(dá)式為

Δr=[JTJ+μ(BTB+WTW)]-1×

{JT[d-g(r0)]-μBTBr0}

(14)

2.3　自適應(yīng)正則化因子選取

正則化因子是數(shù)據(jù)擬合泛函和模型約束泛函之間的折衷系數(shù)，故它的選取直接影響反演結(jié)果。在誤差水平未知的情況下，常用的參數(shù)選取方法有GCV準(zhǔn)則和L曲線準(zhǔn)則。但是，對(duì)于非線性問(wèn)題的迭代解法來(lái)說(shuō)，這兩種準(zhǔn)則不太適用。反演時(shí)總是將數(shù)據(jù)擬合置于優(yōu)先的地位，而在迭代過(guò)程中數(shù)據(jù)擬合泛函值越來(lái)越小，故模型約束泛函在總目標(biāo)泛函中的權(quán)重也必須隨之減小，即正則化因子是隨著迭代過(guò)程衰減的。為了適應(yīng)復(fù)雜的反演情況，正則化因子的調(diào)整過(guò)程應(yīng)當(dāng)與反演計(jì)算相聯(lián)系。采用陳小斌[21]提出的自適應(yīng)調(diào)整方案，其表達(dá)式為

(15)

式中：μi為正則化因子經(jīng)驗(yàn)初值；k為迭代次數(shù)；φdf、φmc分別為數(shù)據(jù)擬合泛函與模型約束泛函。上述正則化因子的選取方式雖然顯著加快了收斂速度，但在將數(shù)據(jù)擬合置于優(yōu)先地位的情況下，迭代次數(shù)過(guò)多會(huì)導(dǎo)致過(guò)擬合，故本文將擬合標(biāo)準(zhǔn)差不再明顯改變時(shí)對(duì)應(yīng)的迭代次數(shù)作為最佳迭代次數(shù)，所以當(dāng)?shù)螖?shù)達(dá)到最佳迭代次數(shù)或擬合標(biāo)準(zhǔn)差達(dá)到閾值時(shí)，迭代中止。

3　巖蓋模型實(shí)驗(yàn)

采用巖蓋模型[17]分析不同約束條件對(duì)重力異常反演的影響。在有噪反演中，所有模型數(shù)據(jù)采用的噪聲變化范圍是原始異常幅值的±2%。

巖蓋模型大致呈菌狀穹窿狀，傘蓋部分沿剖面方向?qū)捈s1km，傘柄部分與剖面方向大致成45°向下延伸。觀測(cè)剖面垂直于模型走向(走向垂直于紙面)，其中心位于傘蓋中心正上方，沿剖面方向向兩側(cè)各延伸1000m，總長(zhǎng)度為2000m，點(diǎn)距為20m，密度差為0.5g/cm3，徑向模型參數(shù)個(gè)數(shù)為100，初始正則化因子為0.1，會(huì)聚方向由虛線表示，與剖面方向分別成0°、45°、180°，會(huì)聚因子ε=0.167。圖4～圖9為巖蓋模型的反演結(jié)果。

對(duì)比圖4與圖5、圖6與圖7、圖8與圖9可見(jiàn)，隨機(jī)噪聲對(duì)反演結(jié)果的影響不大。從圖4～圖7可以看出，絕對(duì)鄰近約束條件下的反演結(jié)果能反映模型上邊界的大致范圍，但是整體上反映的上邊界偏深，對(duì)下邊界特征則完全沒(méi)有反映。對(duì)比圖4與圖6、圖5與7可見(jiàn)，絕對(duì)和相對(duì)鄰近聯(lián)合約束下得到的反演結(jié)果比絕對(duì)鄰近約束下得到的結(jié)果更光滑，但是整體形態(tài)并沒(méi)有明顯變化。從圖8與圖9可以看出，會(huì)聚約束條件下的反演結(jié)果較為準(zhǔn)確地反映了模型的真實(shí)形態(tài)，準(zhǔn)確的先驗(yàn)信息有助于得到更好的反演結(jié)果。

圖4　絕對(duì)鄰近約束下無(wú)噪反演結(jié)果

圖5　絕對(duì)鄰近約束下有噪反演結(jié)果

圖6　絕對(duì)鄰近約束和相對(duì)鄰近約束下無(wú)噪反演結(jié)果

圖7　絕對(duì)鄰近約束和相對(duì)鄰近約束下有噪反演結(jié)果

圖8　會(huì)聚約束下無(wú)噪反演結(jié)果

圖9　會(huì)聚約束下有噪反演結(jié)果

以上模型實(shí)驗(yàn)表明，在會(huì)聚方向約束條件下，徑向反演方法效果最好；絕對(duì)鄰近約束條件下，反演結(jié)果能大致反映模型的上邊界，而下邊界總是趨于弧狀；相對(duì)鄰近約束能讓絕對(duì)鄰近約束條件下得到的反演結(jié)果更加光滑；凸性約束在控制反演結(jié)果的凹凸性上能起到一定的效果。

上文已經(jīng)提到自適應(yīng)正則化方法雖然在一定程度上提高了數(shù)據(jù)擬合度，但是可能會(huì)造成過(guò)擬合。由于正則化因子伴隨著迭代修正不斷改變，但不是一味衰減而是伴隨著振蕩，這也就導(dǎo)致了擬合誤差在迭代到一定次數(shù)后也會(huì)發(fā)生振蕩。這種振蕩可能達(dá)到多個(gè)數(shù)量級(jí)，這也就導(dǎo)致了反演結(jié)果出現(xiàn)很大的變化，這種變化的好壞是未知的，所以需要確定最優(yōu)迭代次數(shù)。下面的例子僅僅是反演實(shí)驗(yàn)中得到的一個(gè)特例。

圖10　重力異常擬合標(biāo)準(zhǔn)差隨迭代次數(shù)變化

圖11　巖蓋模型迭代100次反演結(jié)果

圖10和圖11是絕對(duì)和相對(duì)鄰近聯(lián)合約束下得到的誤差收斂曲線及無(wú)噪數(shù)據(jù)誤差收斂曲線和反演結(jié)果。在25次迭代之后，擬合誤差出現(xiàn)了劇烈振蕩。對(duì)于迭代100次得到的結(jié)果，已經(jīng)接近于會(huì)聚約束下才能得到的結(jié)果。但這僅僅是一個(gè)特例，筆者認(rèn)為迭代中止原則還是遵循上文給出的最佳迭代次數(shù)原則為宜。

4　應(yīng)用實(shí)例

實(shí)際資料選取的是在新墨西哥州里奧格蘭德峽谷大橋(圖12)上測(cè)量得到的重力數(shù)據(jù)[22]。里奧格蘭德峽谷大橋長(zhǎng)約400m，高約175m。峽谷周邊地層大致可分為三種，分別為玄武巖地層，沖積沉積層，夾雜極薄湖成沉積層，其中玄武巖密度約為2.70g/cm3、沉積巖密度約為2.17g/cm3，具體情況如圖12所示。

圖12　奧格蘭德峽谷大橋截面圖

假定峽谷的周邊地層的85%為玄武巖、 15%為沖積沉積物，則峽谷與圍巖及空氣的密度差為-2.62g/cm3。觀測(cè)剖面長(zhǎng)度為800m，點(diǎn)距為20m，徑向模型參數(shù)個(gè)數(shù)為40，初始正則化因子為0.1，初始模型半徑為25m，會(huì)聚因子ε=0.32。反演結(jié)果如圖13所示。

從圖13a中的異常擬合圖可見(jiàn)，反演得到的模型正演異常很好地?cái)M合了實(shí)測(cè)數(shù)據(jù)。從反演的場(chǎng)源形態(tài)看，反演結(jié)果較好地反映了場(chǎng)源的深度范圍，頂部和左側(cè)較好地反映了場(chǎng)源的邊界；右側(cè)底部與實(shí)際差異較大， 可能是由于實(shí)測(cè)異常并不是嚴(yán)格對(duì)稱的，從而導(dǎo)致反演結(jié)果在底部整體向左偏移。借鑒界面反演的相關(guān)做法，計(jì)算反演多邊形的半徑各頂點(diǎn)標(biāo)準(zhǔn)差為17.1186m，說(shuō)明反演結(jié)果較可靠。

圖13　會(huì)聚及凸性約束下反演結(jié)果

5　結(jié)論

徑向反演只能應(yīng)用于單連通域，且初始模型中心必須位于場(chǎng)源內(nèi)部。對(duì)于不太復(fù)雜的單個(gè)場(chǎng)源，在利用歐拉反褶積及相關(guān)成像等方法得到場(chǎng)源中心近似位置后，利用徑向反演可以圈定場(chǎng)源的大致形態(tài)；如果能進(jìn)一步獲得關(guān)于場(chǎng)源分布的半定量信息，利用會(huì)聚約束能得到較好的反演結(jié)果。

正則化因子的自適應(yīng)調(diào)整能提高數(shù)據(jù)擬合精度和收斂速度，但是由于實(shí)際數(shù)據(jù)的復(fù)雜性和反演的多解性，片面追求數(shù)據(jù)擬合效果可能會(huì)造成過(guò)擬合，利用最佳迭代次數(shù)選擇原則雖然能解決這一問(wèn)題，但不可避免地引入人為因素。

此外，由于本文采用的是Gauss-Newton迭代法，故模型初值的選擇對(duì)反演結(jié)果有一定的影響。本方法的初始模型一般采用給定中心埋深的無(wú)限長(zhǎng)水平圓柱體。在實(shí)際數(shù)據(jù)反演中，地質(zhì)體中心埋深往往是利用其他方法(如歐拉反褶積、特征點(diǎn)法等)得到的估計(jì)值，估計(jì)值越接近真實(shí)值，最終反演結(jié)果精度也越高。采用更加合理、精度更高的計(jì)算中心埋深方法與本方法結(jié)合將是后續(xù)研究的方向。當(dāng)然，如果有更多的地球物理資料并建立較無(wú)限長(zhǎng)水平圓柱體更加復(fù)雜的初始模型，反演的精度會(huì)進(jìn)一步得到提高。

[1]谷文彬,陳清禮,王余泉等.饒陽(yáng)凹陷虎8北潛山重力三維松約束反演.石油地球物理勘探,2016,51(6):1219-1226.

Gu Wenbin,Chen Qingli,Wang Yuquan et al.Part-constrained 3D gravity inversion for the Hubabei buried hill in Raoyang Sag.OGP,2016,51(6):1219-1226.

[2]李婧,崔永謙,宋景明等.重力—地震聯(lián)合建模正演技術(shù)在SH地區(qū)潛山勘探中的應(yīng)用.石油地球物理勘探,2016,51(增刊1):152-160.

Li Jing,Cui Yongqian,Song Jingming et al.Gravity forward modeling with model built based on drilling and seismic data:an example of buried hill exploration in Area SH.OGP,2016,51(S1):152-160.

[3]耿美霞,楊慶節(jié).應(yīng)用RBF神經(jīng)網(wǎng)絡(luò)反演二維重力密度分布.石油地球物理勘探,2013,48(4):651-657.

Geng Meixia,Yang Qingjie.2-D density inversion with the RBF neural network method.OGP,2013,48(4):651-657.

[4]王浩然,陳超,杜勁松.重力梯度張量數(shù)據(jù)的三維反演方法與應(yīng)用.石油地球物理勘探,2013,48(3):474-481.

Wang Haoran,Chen Chao,Du Jinsong.3-D inversion of gravity gradient tensor data and its applications.OGP,2013,48(3):474-481.

[5]Cordell L and Henderson R G.Iterative three-dimensional solution of gravity anomaly data using a digital computer.Geophysics,1968,33(4):596-601.

[6]Dyrelius D and Vogel A.Improvement of convergence in iterative gravity interpretation.Geophysical Journal of Royal Astronomical Society,1972,27(2):195-205.

[7]Pedersen L B.Interpretation of potential field data - A generalized inverse approach.Geophysical Prospecting,1977,25(2):199-230.

[8]Barbosa V C F,Silva J B C and Medeiros W E.Stable inversion of gravity anomalies of sedimentary basins with non-smooth basement reliefs and arbitrary density contrast variations.Geophysics,1999,64(3):754-764.

[9]Corbato CE.A least-squares procedure for gravity interpretation.Geophysics,1965,30(2):228-233.

[10]Al-Chalabi M.Interpretation of gravity anomalies by nonlinear optimization.Geophysical Prospecting,1972,20(1):1-16.

[11]Moraes R A V and Hansen R O.Constrained inver-sion of gravity fields for complex 3D structures.Geophysics,2001,66(2):501-510.

[12]Wildman R A and Gazonas G A.Gravitational and magnetic anomaly inversion using a tree-based geometry representation.Geophysics,2009,74(3):I23-I35.

[13]Li Y and Oldenburg D W.3-D inversion of gravity data.Geophysics,1998,63(1):109-119.

[14]Portniaguine O and Zhdanov M S.Focusing geophysical inversion images.Geophysics,1999,64(3):874-887.

[15]Silva J B C and Barbosa V C F.Interactive gravity inversion.Geophysics,2006,71(1):1-9.

[16]Silva J B C and Barbosa V C F.Generalized radial inversion of 2D potential field data.Geophysics,2004,69(6):1405-1413.

[17]賈真.均勻二度體位場(chǎng)正演與徑向反演方法[學(xué)位論文].吉林長(zhǎng)春:吉林大學(xué),2009.

Jia Zhen.Forward and Radial Inverse Modeling of Potential Fields Produced by a 2D Homogeneous Source[D].Jilin University,Changchun,Jilin,2009.

[18]Vanderlei C,Oliveira J,Valeria C F et al.Source geo-metry estimation using the mass excess criterion to constrain 3-D radial inversion of gravity data.Geophysical Journal International,2011,187(2):754-772.

[19]Vanderlei C,Oliveira J,Valeria C F et al.3-D radial gravity gradient inversion.Geophysical Journal International,2013,195(2):883-902.

[20]賈真,孟令順.均勻多邊形截面二度體重力異常計(jì)算公式的改進(jìn).地球物理學(xué)進(jìn)展,2009,24(2):462-467.

Jia Zhen,Meng Lingshun.Some improvements on the formula for calculating the gravity anomaly due to a

2D homogeneous polygonal source.Progress in Geophysics,2009,24(2):462-467.

[21]陳小斌.大地電磁正反演新算法研究及資料處理與解釋的可視化集成系統(tǒng)開(kāi)發(fā)[學(xué)位論文].北京:中國(guó)地

震局地質(zhì)研究所,2003.

Chen Xiaobin.New Forward and Inversion Algorithms and a Visual Integrated System for MT data [D].Institute of Geology,China Earthquake Administration,Beijing,2003.

[22]Titus W J,Titus S J and Davis J R.A Bayesian approach to modeling 2D gravity data using polygons.Geophysics,2017,82(1):G1-G21.