二維重力數據徑向反演及應用

林寶澤　肖　鋒　王明常

(吉林大學地球探測科學與技術學院，吉林長春 130000)

1　引言

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

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

建模對反演來說至關重要，模型過于復雜會導致解的不穩定性，可以通過Tikhonov正則化方法引入約束條件解決這一問題。徑向反演的優點在于將多邊形頂點位置用極坐標表示，反演參數(即多邊形各頂點對應的半徑)具有同一量綱且反演參數相對于直角坐標系減少一半。利用Tikhonov正則化引入相關約束條件(場源幾何形態等先驗信息)時能取得較好的效果。

徑向反演方法雖然已經發展得較為完善，但在國內關于該方法的研究和討論幾乎還是空白。此外，該方法的研究多見于模型實驗，對于影響反演結果的正則化參數選取更是缺乏討論。本文介紹了二維徑向反演方法的基本原理，且在此基礎上將自適應正則化因子選取引入該方法，并將該方法應用到實際資料處理，取得了較好的效果。

2　方法原理

2.1　正演方法

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

(1)

式中：rk為第k個頂點對應的半徑;θk為第k個頂點對應的角度(弧度,沿x軸正方向順時針旋轉為正)。模型及參數如圖1所示。

圖1　徑向模型示意圖

從圖1不難推斷，當頂點數M足夠大時，徑向模型能取得很好的近似效果。但是從式(1)來看，應用徑向模型時存在以下限制：一是包含場源的區域必須是單連通域；二是從模型中心點出發的射線與場源輪廓線有且僅有一個交點。

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

(2)

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

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

2.2　約束反演

設重力異常實測數據構成的N×1維向量為d,徑向多邊形各頂點對應的半徑構成的M×1維向量為r,由模型r得到的正演重力異常為N×1維向量g(r)。為了保證反演問題的適定性，引入約束條件后，該問題轉化為求解最優化問題

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

(3)

式中

(4)

(5)

(6)

(7)

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

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

(8)

不等性約束的作用是限制場源的幾何尺寸和保證頂點z坐標為非負。凸性約束就是要在每次迭代中滿足

(9)

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

圖3　凸性約束示意圖

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

=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)

對式(12)兩端求導，可得

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

(13)

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

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

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

(14)

2.3　自適應正則化因子選取

正則化因子是數據擬合泛函和模型約束泛函之間的折衷系數，故它的選取直接影響反演結果。在誤差水平未知的情況下，常用的參數選取方法有GCV準則和L曲線準則。但是，對于非線性問題的迭代解法來說，這兩種準則不太適用。反演時總是將數據擬合置于優先的地位，而在迭代過程中數據擬合泛函值越來越小，故模型約束泛函在總目標泛函中的權重也必須隨之減小，即正則化因子是隨著迭代過程衰減的。為了適應復雜的反演情況，正則化因子的調整過程應當與反演計算相聯系。采用陳小斌[21]提出的自適應調整方案，其表達式為

(15)

式中：μi為正則化因子經驗初值；k為迭代次數；φdf、φmc分別為數據擬合泛函與模型約束泛函。上述正則化因子的選取方式雖然顯著加快了收斂速度，但在將數據擬合置于優先地位的情況下，迭代次數過多會導致過擬合，故本文將擬合標準差不再明顯改變時對應的迭代次數作為最佳迭代次數，所以當迭代次數達到最佳迭代次數或擬合標準差達到閾值時，迭代中止。

3　巖蓋模型實驗

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

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

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

圖4　絕對鄰近約束下無噪反演結果

圖5　絕對鄰近約束下有噪反演結果

圖6　絕對鄰近約束和相對鄰近約束下無噪反演結果

圖7　絕對鄰近約束和相對鄰近約束下有噪反演結果

圖8　會聚約束下無噪反演結果

圖9　會聚約束下有噪反演結果

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

上文已經提到自適應正則化方法雖然在一定程度上提高了數據擬合度，但是可能會造成過擬合。由于正則化因子伴隨著迭代修正不斷改變，但不是一味衰減而是伴隨著振蕩，這也就導致了擬合誤差在迭代到一定次數后也會發生振蕩。這種振蕩可能達到多個數量級，這也就導致了反演結果出現很大的變化，這種變化的好壞是未知的，所以需要確定最優迭代次數。下面的例子僅僅是反演實驗中得到的一個特例。

圖10　重力異常擬合標準差隨迭代次數變化

圖11　巖蓋模型迭代100次反演結果

圖10和圖11是絕對和相對鄰近聯合約束下得到的誤差收斂曲線及無噪數據誤差收斂曲線和反演結果。在25次迭代之后，擬合誤差出現了劇烈振蕩。對于迭代100次得到的結果，已經接近于會聚約束下才能得到的結果。但這僅僅是一個特例，筆者認為迭代中止原則還是遵循上文給出的最佳迭代次數原則為宜。

4　應用實例

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

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

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

從圖13a中的異常擬合圖可見，反演得到的模型正演異常很好地擬合了實測數據。從反演的場源形態看，反演結果較好地反映了場源的深度范圍，頂部和左側較好地反映了場源的邊界；右側底部與實際差異較大， 可能是由于實測異常并不是嚴格對稱的，從而導致反演結果在底部整體向左偏移。借鑒界面反演的相關做法，計算反演多邊形的半徑各頂點標準差為17.1186m，說明反演結果較可靠。

圖13　會聚及凸性約束下反演結果

5　結論

徑向反演只能應用于單連通域，且初始模型中心必須位于場源內部。對于不太復雜的單個場源，在利用歐拉反褶積及相關成像等方法得到場源中心近似位置后，利用徑向反演可以圈定場源的大致形態；如果能進一步獲得關于場源分布的半定量信息，利用會聚約束能得到較好的反演結果。

正則化因子的自適應調整能提高數據擬合精度和收斂速度，但是由于實際數據的復雜性和反演的多解性，片面追求數據擬合效果可能會造成過擬合，利用最佳迭代次數選擇原則雖然能解決這一問題，但不可避免地引入人為因素。

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

[1]谷文彬,陳清禮,王余泉等.饒陽凹陷虎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]李婧,崔永謙,宋景明等.重力—地震聯合建模正演技術在SH地區潛山勘探中的應用.石油地球物理勘探,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]耿美霞,楊慶節.應用RBF神經網絡反演二維重力密度分布.石油地球物理勘探,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]王浩然,陳超,杜勁松.重力梯度張量數據的三維反演方法與應用.石油地球物理勘探,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]賈真.均勻二度體位場正演與徑向反演方法[學位論文].吉林長春:吉林大學,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]賈真,孟令順.均勻多邊形截面二度體重力異常計算公式的改進.地球物理學進展,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]陳小斌.大地電磁正反演新算法研究及資料處理與解釋的可視化集成系統開發[學位論文].北京:中國地

震局地質研究所,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.