重力三維密度光滑反演軟件設計與實現

彭衍武 陳文進 譚小龍 洪國慶

1 江西理工大學土木與測繪工程學院,江西省贛州市客家大道1958號,341000

通過重力密度反演分析技術可以更準確地了解地質體結構、形狀及其他特征,為此,國內外學者進行大量工作。Li等[1-2]引入深度加權函數;Portniaguine等[3]引入最小梯度支撐函數約束重力反演;Chasseriau等[4]提出基于隨機方法的反演技術;Zhdanov[5]提出重磁反演和成像新方法;Commer[6]提出深度分辨率增強的加權方案;耿美霞等[7]建立在密度約束下的多變量協同克里金聯合反演方程;張志厚等[8]提出一種基于深度學習的重力異常及重力梯度異常的聯合反演方法。

在三維密度反演領域已經存在諸多算法,但可用于密度反演的公開軟件卻很少?；诖?本文借助 MATLAB語言,開發一個支持三維密度反演、重力場正演、數據分析及可視化等功能的系統軟件,為重力三維密度反演學習和研究提供軟件工具支持。

1 方法原理

1.1 重力場正演

根據牛頓萬有引力定律,由計算點體積質量密度元素產生的萬有引力勢為[9]:

(1)

(2)

式中,r為徑向矢量;ΔVj為第j個棱柱體的體積;ρj為第j個棱柱體的密度;M為長方體總數。本文采用Nagy等[10]推導的位、引力矢量和引力梯度張量表達式用于重力場正演計算。

1.2 三維密度反演

根據待求解密度模型的特點建立合適的反演目標函數:

S=φd+μφm

(3)

式中,φd為數據擬合項。通常認為重力數據觀測噪聲符合高斯分布特點,因而采用L2范數進行構建:

(4)

式中,Wd為每個數據點的權重方陣;dobs為觀測數據點的N×1維向量;dpre為預測數據點的N×1維向量。模型擬合項定義如下:

(5)

(6)

式中,z0為觀測高度;z為單元深度;β為擬合指數項。因此,目標函數可表示為:

φ(ρ)=φ(ρ)d+μφm(ρ)

(7)

最后,通過最小二乘可得:

(8)

2 軟件設計與功能

本文軟件是基于 Windows10系統,開發環境為Intel(R) Core(TM) i5-9300H CPU,8G內存。系統結構和操作流程如圖1所示。該軟件由三維密度反演和數據結果界面兩部分組成(圖2),主要功能包括參數設置、地下模型空間劃分、重力觀測數據選擇、反演結果顯示與預測。

圖1 系統設計框架與操作流程

圖2 三維密度反演主界面與子界面

3 數據測試

3.1 單長方體模型數據測試

測區范圍為1 000 m×1 000 m×500 m,網格劃分為20×20×10,單模型x軸范圍為400～600 m,y軸范圍為400～600 m,z軸范圍為100～300 m,密度為1 g/cm3。模型切面及觀測數據如圖3所示,反演結果如圖4所示。

圖3 模型切面與重力觀測數據

圖4 反演模型在x=500 m、z=200 m處密度切面

從圖4可以看出,反演結果能夠有效地反映地質體的空間分布,驗證了本文軟件的有效性。

3.2 實測數據測試

本文采用文頓鹽丘(Vinton Dome)地區重力梯度實測數據進行反演測試。文頓鹽丘位于美國墨西哥灣沿岸,該地區巖蓋密度約為2.75 g/cm3,鹽丘附近沉積物、圍巖及鹽核密度均約為2.2 g/cm3,即巖蓋剩余密度約為0.55 g/cm3。本文選用的測區范圍為4 000 m×4 000 m×1 000 m,網格劃分為40×40×20。

圖5為實測重力梯度gzz和預測gzz數據,由圖可知,預測數據與實測數據擬合效果較好。圖6為密度反演結果,從圖中可以看出,反演密度集中部分頂部埋深約為200 m,底部埋深約為600 m。該結果與前人研究所得埋深結果200～650 m[11]、200～550 m[12]較為接近,驗證了軟件的可靠性與實用性。

圖5 實測重力梯度數據與密度異常體預測數據

圖6 反演模型在y=2 000 m、z=300 m和z=700 m處密度切面

4 結 語

本文開發出一個重力三維密度光滑反演軟件,具有最優正則化參數選取、數據文件管理、密度反演計算及結果分析和可視化等功能。模擬和實測數據測試結果驗證了該軟件的有效性與實用性,未來可以在計算效率等方面進一步優化。