重力及其垂直梯度交叉約束聯合密度結構反演方法

馬國慶, 李欣蔚, 王泰涵*, 熊盛青, 高桐

1 吉林大學地球探測科學與技術學院, 長春 130026 2 中國自然資源航空物探遙感中心, 北京 100083

0 引言

重力梯度數據具有較高水平分辨率,相對原始重力場可凸顯較淺地質體.重力梯度測量可通過不同高度測量重力值的方式來獲得,其與重力聯合測量可以探測不同深度的空洞、人防工程(Cordell, 1979; Schlinger, 1990; Bate, 2005; Butler, 2012). 通過重力異常反演地下密度分布可以圈定場源體的分布,主要的反演方法包括光滑約束反演、粗糙約束反演、聚焦反演、協克里金反演、模糊C均值聚類反演等(Li and Oldenburg, 1998; Portniaguine and Zhdanov, 1999; Zhdanov et al., 2004; Pilkington, 2009; Geng et al., 2014; Sun and Li, 2015; 王泰涵等, 2020; 李芳等, 2021).

重力和梯度的聯合反演方法現今主要采用將不同分量數據組合為一個矩陣,通過數據約束來提升反演結果的精度(Tikhonov and Arsenin, 1977; Qin et al., 2016; Wang et al., 2022).Wu等(2013)提出加權最小二乘估計方法,通過自適應權值實現同源數據的聯合.Capriotti和Li(2014)基于靈敏度矩陣平衡核函數衰減速率,并定義了一種重力和張量梯度數據聯合的公式.Qin等(2016)以非線性牛頓法和最小梯度支撐進行重力和重力梯度反演計算,獲得了較清晰的異常邊界分布特征.高秀鶴等(2019)利用密度閾值不斷更新協方差矩陣來實現重力和重力梯度數據協克里金反演.重力與重力梯度聯合反演還可以采用交叉梯度方法,其是從結構上對不同物性參數進行耦合的聯合反演方法,主要用于不同物性源地球物理數據的聯合反演(Gallardo and Meju,2003; Fregoso and Gallardo, 2009; Zhou et al., 2015).重力和重力梯度數據是同源數據,但由于不同參量數據對于深淺場源反演效果的差異,為此采用交叉梯度的結構約束特征可提升反演的分辨率和精度(Zhao et al., 2018; Entezar-Saadat et al., 2020).

為了充分利用重力與重力梯度對地下密度結構的敏感性差異，我們提出重力及其垂直梯度交叉約束聯合密度結構反演方法,以重力與重力梯度單獨反演結果作為權函數來約束交叉梯度聯合反演過程,有效恢復地下不同深度場源的密度結構分布特征.最后通過模型實驗和實際數據證明該方法較常規數據聯合密度反演方法在反演分辨率和精度的優勢，以及在實際應用中的適用性.

1 基本原理

現今重力數據(Vz)和垂直梯度數據(Vzz)密度反演通常將不同數據放置于同一矩陣,采用Tikhonov正則化方法進行反演,其目標函數的表達式為(Tikhonov and Arsenin, 1977; Qin et al., 2016)：

(1)

式中為Φd數據擬合項；Φm為模型約束項；α為正則化參數，用于權衡數據擬合項和模型約束項對正則化方程的影響；Wd為數據加權矩陣；F為重力及其垂直梯度的組合核函數矩陣[Fz,Fzz]T；m為物性(密度)參數；d為重力及其垂直梯度觀測數據的組合矩陣[Vz,Vzz]T；Wz為深度加權矩陣,可以有效抵消核矩陣隨深度的固有衰減,以降低趨膚效應的影響(Li and Oldenburg, 2000).

采用不同大小、不同埋深的異常體A和異常體B來分析重力和重力梯度聯合的Tikhonov正則化反演方法效果,其具體參數如表1所示,分布如圖1a所示.圖1b為重力異常更好地體現較深場源信息,但兩個異常體并未分開.圖1c為垂直梯度異常更好地反映較淺場源的信息,并且兩個異常體的邊界較清晰.圖1d、e分別為利用Tikhonov正則化方法反演得到的重力和重力梯度的密度反演結果在y=1000 m位置的切片圖,可以看出重力結果在深部的分辨率高于重力梯度,重力梯度對于淺部地質目標有更高的反演精度.重力及其垂直梯度異常雖反應地下同一密度結構,但是由于數據敏感度差異，使恢復密度結構易于集中在變化較大區域,模型分辨率產生了差異性.圖1f為重力和重力梯度組合矩陣的Tikhonov正則化方法反演結果,對于淺部地質體的反演效果有一定的提升,但是對于深部地質體依舊較差.

圖1 兩個棱柱體模型的正演異常及不同反演方法結果在y=1000 m處剖面圖(a) 兩個棱柱體組合模型; (b) 重力異常; (c) 重力垂直梯度異常; (d) Tikhonov正則化方法Vz數據反演; (e) Tikhonov正則化方法Vzz數據反演; (f) Tikhonov正則化方法[Vz, Vzz]T數據反演; (g) 交叉梯度Vz數據反演; (h) 交叉梯度Vzz數據反演; (i) 交叉梯度Vzz數據反演結果三維透視圖.Fig.1 Forward response and Section at y=1000 m of inversion results based on different methods of two-prism model(a) The model combined by two prisms; (b) Gravity anomaly; (c) Gravity vertical gradient anomaly; (d) Tikhonov regularization inversion result of Vz data; (e) Tikhonov regularization inversion result of Vzz data; (f) Tikhonov regularization inversion result of [Vz, Vzz]T; (g) Cross-gradient inversion result of Vz data; (h) Cross-gradient inversion result of Vzz data; (i) 3D perspective view of cross-gradient inversion result of Vzz data.

表1 兩個棱柱體模型Table 1 Parameters oftwo-prism model

基于不同重力參量數據對應同源的密度結構,但是不同數據反演結果存在相似性和差異性,采用交叉梯度引入的結構約束可有效提升結果的精度和分辨率.交叉梯度函數(Gallardo and Meju, 2003)的表達式為：

其中，m1和m2分別為不同的物性參數，當兩個梯度方向平行或者不存在時t為0.其在x、y、z三個方向上的分量可以表示為：

(3)

采用交叉梯度的方法進行重力數據和重力梯度聯合反演,其目標函數為：

(5)

式中，λz和λzz為交叉約束項系數.式中參數的下角標z和zz分別代表重力和重力梯度數據相關參數.Wz1和Wz2分別為重力及其梯度數據對應的深度加權矩陣.采用圖1b、c所示的重力及其梯度數據進行交叉梯度反演，得到圖1g、h的兩個反演結果，深部和淺部分辨率有提高,但與實際模型依舊存在較大偏差.

為了提升重力和重力梯度聯合反演的分辨率來獲取不同深度場源的分布,本文提出重力及其垂直梯度交叉約束聯合密度結構反演方法,其首先對重力異常與重力梯度進行單獨Tikhonov正則化反演,然后將各自反演結果形成的歸一化權函數進行交叉互換使用,并采用交叉梯度進行結構約束來實現重力與其梯度的聯合反演,以期達到不同深度場源分辨率的提升,最終取雙方結果的平均值作為最終的密度分布,具體的計算流程如圖2所示.

圖2 交叉約束聯合密度結構反演技術流程Fig.2 The flow of the cross-constrained density joint inversion based on gravity and vertical gradient data

圖2中，Tikhonov正則化方法重力數據反演結果mz1和重力垂直梯度數據反演結果mzz1形成的交叉約束權函數Wρz和Wρzz的表達式為：

(6)

其中(mz1)max、(mzz1)max為矩陣中的最大的元素；(mz1)min、(mzz1)min為矩陣中最小的元素.為了避免產生奇點,令εz和εzz分別等于0.1倍的(mz1)min和0.1倍的(mzz1)min.μ是一個控制交叉約束權在反演中約束影響的常數,通常取值為1.不同參量結果互為權函數的交叉梯度反演表達式為：

最后，取交叉約束權聯合密度反演結果mz2和mzz2的平均值作為最終密度分布mFin.

采用表1中的模型來試驗交叉約束聯合密度結構反演方法的效果,圖3a、b為采用該方法進行反演獲得的重力和重力梯度的反演結果,可以看出異常的位置準確,水平分辨率較高.圖3c為重力和重力梯度的反演結果取平均值融合后得到的最終結果,可以看出相對放置同一矩陣反演結果的精度大,且密度值更接近理論值.圖3d為密度大于0.65 g·cm-3的三維圖,其中藍色邊框為實際模型范圍,本文方法結果形態較前文的其他反演方法更符合實際密度分布,且密度幅值更接近模型的真實密度.為了驗證所提出交叉約束聯合密度結構反演方法的準確性,將不同反演方法結果在y=1000 m處的重力與正演重力異常進行對比,從圖4可以看出在異常變化平緩區域全都擬合較好,而異常變化率較大的位置,本文方法的反演結果的擬合差更小.

圖3 交叉約束聯合密度結構反演結果(a) Vz數據反演結果在y=1000 m處剖面圖; (b) Vzz數據反演結果在y=1000 m處剖面圖; (c) 反演最終結果; (d) 密度大于0.65 g·cm-3反演結果三維透視圖.Fig.3 Inversion results of the cross-constrained joint inversion method(a) Section at y=1000 m of inversion result of Vz data; (b) Section at y=1000 m of inversion result of Vzz data; (c) Section at y=1000 m of final inversion result; (d) 3D perspective view of inversion result with density larger than 0.65 g·cm-3.

圖4 不同方法反演結果在y=1000m處的與觀測數據擬合曲線Fig.4 The fitting curve map of inversion results of different methodsand observed data at section of y=1000 m

2 模型試驗

為了檢驗本文方法的抗噪能力,對圖5a所示的模型的重力和梯度數據加入最大值為10%的高斯噪聲,含噪聲的數據正演異常如圖5b、c所示.表2統計了重力數據進行正則化反演、重力及其梯度數據正則化聯合反演、重力及梯度數據交叉梯度聯合反演以及本文方法的預測數據與觀測數據的擬合差進行對比，可以看出本文方法較其他反演方法擬合差更小,證明了本文方法抗噪性良好.

圖5 含噪數據正演異常和不同方法反演結果在y=1000 m處剖面圖(a) 兩個棱柱體組合模型; (b) 重力異常; (c) 重力垂直梯度異常; (d) Tikhonov正則化方法Vz數據反演結果; (e) Tikhonov正則化方法Vzz數據反演結果; (f) Tikhonov正則化方法[Vz, Vzz]T反演結果; (g) 交叉梯度Vz數據反演結果; (h) 交叉梯度Vzz數據反演結果; (i) 交叉約束Vz數據反演結果; (j) 交叉約束Vzz數據反演結果; (k) 交叉約束方法反演最終結果; (l) 交叉約束方法結果三維透視圖.Fig.5 Forward anomaly and sections at y=1000 m of inversion results of different methods of data with noise(a) The model combined by two prisms; (b) Gravity anomaly; (c) Gravity vertical gradient anomaly; (d) Tikhonov regularization inversion result of Vz data; (e) Tikhonov regularization inversion result of Vzz data; (f) Tikhonov regularization inversion result of [Vz, Vzz]T; (g) Cross-gradient inversion result of Vz data; (h) Cross-gradient inversion result of Vzz data; (i) Cross-constrained inversion result of Vz data; (j) Cross-constrained inversion result of Vzz data; (k) Final inversion result of cross-constrained method; (l) 3D perspective view of cross-constrained inversion result.

表2 含噪數據不同反演方法擬合差Table 2 The misfit of the different inversion methods for data with noise

對比重力異常圖5b和重力梯度異常圖5c,重力梯度數據水平分辨率高于重力數據.重力數據較梯度數據對于深部異常更為敏感,重力梯度數據對于淺部異常更敏感.含噪聲模型實驗的Tikhonov正則化反演方法梯度數據結果圖5e淺部異常的分辨率和水平分辨率優于重力數據反演結果圖5d,圖5d深部異常的分辨率優于圖5e，驗證了重力和梯度數據的數據特征.并且梯度數據對于高頻信號較為敏感,噪聲屬于高頻信號,可以看到有梯度數據參與反演的Tikhonov正則化方法Vzz數據反演結果圖5e、Tikhonov正則化方法[Vz,Vzz]T反演結果圖5f和交叉梯度Vzz數據反演結果圖5h在淺部都一定程度上受噪聲干擾.采用交叉約束聯合密度結構反演方法,可以能夠降低噪聲的影響,使反演結果更集中于密度高值分布的區域.交叉約束聯合密度結構反演方法的結果圖5i、j在兩個異常體定位和分辨率方面明顯優于圖5d—f所示的Tikhonov正則化方法反演結果以及圖5g、h交叉梯度的反演結果.圖5l為密度值大于0.65 g·cm-3的部分的三維密度結果,比較符合真實密度分布.

采用不同密度的異常體組合模型來驗證交叉約束密度結構反演方法對密度差異區間較大的復雜地質體的適用性,模型參數見表3,圖6a為異常體的密度真實分布.

圖6 三個棱柱體模型及其正演異常(a) 三個棱柱體組成的模型; (b) 重力異常; (c) 重力垂直梯度.Fig.6 Three-prism model and forward anomaly(a) The model combined by three prisms; (b) Gravity anomaly; (c) Gravity vertical gradient anomaly.

表3 三個棱柱體模型參數Table 3 Parameters of three-prism model

重力異常圖6b和重力垂直梯度異常6c中異常體③的分辨率較其他兩個異常體高,梯度數據原始異常的異常體②和③的邊界較重力異常的邊界更為清晰.取x=1400 m,y=1000 m和z=400 m處的反演結果圖對比不同方法的反演結果.由于Tikhonov和交叉梯度法的反演結果幅值較低,為了方便對比不同反演方法的結果,將兩種反演方法的結果統一色標柱以便對比.

圖7d—f梯度數據的Tikhonov正則化反演結果較重力Tikhonov正則化反演結果圖7a—c整體分辨率較高.Tikhonov正則化方法[Vz,Vzz]T反演結果圖7h以及交叉梯度方法的圖7k和圖7n中異常體②水平方向密度分布較Tikhonov正則化法的反演結果圖7b、e更為收斂.Tikhonov正則化方法和交叉梯度方法的在z=400 m 反演結果圖7c、f、i、l、o在水平方向上的分辨率都表現良好,三個異常體的邊界都比較明晰.由本文方法的水平方向剖面圖7r可以看出密度高值集中于異常體中心,其垂向剖面圖圖7p、q較Tikhonov正則化方法和交叉梯度方法的垂向剖面圖7a、b、d、e、g、h、j、k、m、n密度分布更符合實際情況,說明本文方法在提升縱向分辨率上優于Tikhonov正則化方法和交叉梯度方法.由圖7r本文方法的反演結果得到的三個異常體中心的密度值分別為0.11 g·cm-3、0.46 g·cm-3、0.73g·cm-3，幅值較接近理論模型的值,證明了本文方法在反演情況復雜的異常體時,對異常體位置定位準確,密度分布更加接近實際情況,三個異常體的邊界都很清晰,并且有較高的縱向分辨率.

圖7 不同反演方法在x=1400 m, y=1000 m, z=400 m反演結果圖(a),(b)及(c)為Tikhonov正則化方法Vz數據反演; (d), (e)及(f) Tikhonov正則化方法Vzz 數據反演; (g), (h)及(i) Tikhonov正則化方法[Vz, Vzz]T反演; (j), (k)及(l) 交叉梯度Vz數據反演; (m),(n)及(o) 交叉梯度Vzz數據反演; (p), (q)及(r) 交叉約束方法反演.Fig.7 Slices at x=1400 m, y=1000 m and z=400 m of inversion results of diffrent methods(a), (b) and (c) Tikhonov regularization inversion of Vz data; (d), (e) and (f) Tikhonov regularization inversion of Vzz data; (g), (h) and (i) Tikhonov regularization inversion of [Vz, Vzz]T; (j), (k) and (l) Cross-gradient inversion of Vz data; (m), (n) and (o) Cross-gradient inversion of Vzz data; (p), (q) and (r) Cross-constrained method.

3 實際應用

3.1 遼源礦區重力及垂直梯度聯合反演

遼源煤礦位于長白山余脈與松嫩平原交界,礦產資源豐富.由于長期的無序開采,地面塌陷給人民生活和城市建設帶來安全隱患 (武欣等,2022).圖8為遼源礦區圖,圖中標注了開采的礦點、已廢棄礦點以及地表塌陷區域,黑色線框內區域為研究區域,面積為3 km×2 km.我們將目標區域劃分為30×20×10個棱柱體單元,每個單元的大小為100 m×100 m×30 m.通過重力及其梯度交叉約束聯合密度結構反演方法,確定采空區的位置,以圈定測區內采空區范圍.

圖8 遼源礦區Fig.8 Liaoyuan mining area

采空區引起的重力異常相對較小,為了尋找采空區采用CG-5重力儀對測區進行高精度的微重力測量(圖9a、c).儀器讀數分辨率1 μGal,測量的均方誤差為0.087 mGal.對實測重力異常去除背景場,并對重力數據和垂直梯度數據做了地形改正(陳善, 1986; 葉周潤等, 2011),得到的剩余異常和重力梯度異常如圖9b、d所示.圍巖密度為2.1 g·cm-3,采空區密度為1.0 g·cm-3,采空區引起的異常密度差接近-1.1 g·cm-3.

圖9 野外重力測量及重力異常(a) 重力梯度測量; (b) 重力異常; (c) CG-5重力儀; (d) 重力垂直梯度.Fig.9 Gravity measurement of field exploration and gravity anomaly(a) Gravity gradiometry; (b) Gravity anomaly; (c) CG-5 Gravimeter; (d) Gravity vertical gradient anomaly.

重力及其梯度交叉約束聯合密度結構反演方法結果如圖10a所示,采空區密度為負值,通過本文方法計算得到的采空區密度最低值為-0.88 g·cm-3與推測的剩余密度值較為符合.在反演結果150 m的切片結果上根據異常邊界圈定了Ⅰ-Ⅴ號采空區,Ⅵ號區域附近為塌陷區,所以也將其認定為一個采空區(圖10a).根據圖8遼源礦區圖,將研究區域內的已知礦點畫在150 m的反演結果上,10個礦點的位置,都位于圈定的采空區內,表明本文方法的反演結果與以往已知礦點位置對應良好.采空區的深度范圍為0～250 m,平均深度200 m.圖10c為測線LL′測線的高密度電阻率剖面,推測在該剖面上采空區深度大約為170 m.圖10b為反演結果在測線同一位置的剖面圖,該剖面三處區域的形狀及位置與高密度電阻率剖面上的范圍有較好的位置對應,驗證了反演結果的準確性.

圖10 反演密度分布結果與高密度電阻率法結果對比(a) 本文交叉約束方法反演結果; (b) 交叉約束反演結果在LL′測線的剖面圖; (c) LL′測線的高密度電阻率法剖面圖.Fig.10 Comparison of inversion results of density distribution and high density resistivity method(a) Inversion result of cross-constrained method; (b) Slice of Line LL′ of the cross-constrained method; (c) Slice of Line LL′ of the high-density resistivity method.

3.2 長白山火山區重力及垂直梯度聯合反演

長白山火山區是新生代多成因復合火山,地球物理勘探發現長白山天池火山下方有巖漿囊的存在,火山有再次噴發的風險,因此研究長白山火山巖漿囊的位置和范圍對火山的監測具有重大價值 (郭履燦等, 1996; 呂政等, 2007).

本文選擇430 km×430 km 的長白山火山區衛星重力數據,去除沉積層和莫霍面起伏產生的影響,得到圖11a剩余重力異常(Blakely,1995).對重力異常求導得到圖11b垂直梯度異常來進行重力及其垂直梯度聯合反演.

圖11 長白山研究區域重力異常(a) 重力異常; (b) 重力垂直梯度異常.Fig.11 Gravity anomaly of Changbai Mountain survey area(a) Gravity anomaly; (b) Gravity vertical gradient anomaly.

將地下網格剖分為43×43×10個單元,每個單元大小為10 km×10 km×3 km.巖漿囊在高溫熔融狀態較圍巖呈低密度,重力異常顯示為低值.重力及其梯度交叉約束聯合密度結構反演方法的三維密度結果如圖12a.圖12b、c為穿過巖漿囊的兩條剖面.根據通過巖漿囊中心的L1和L2剖面圖,巖漿囊呈拱形,位于北向250～370 km,東向190～260 km,呈東北-西南方向分布.本文方法估計的巖漿囊的中心埋深約15 km,深度范圍約為7～20 km.湯吉等(2001)根據大地電磁測深反演剖面結果顯示10～30 km的深度范圍內存在低阻體,認為其可能為活動的巖漿囊; 管彥武等(2020)利用重力剖面數據及地質信息進行人機交互建立的密度模型推斷深度7～15 km處存在巖漿囊,本文推斷的巖漿囊深度整體范圍及中心埋深與這些前人研究比較吻合,進一步驗證了本文方法的適用性.

圖12 長白山區域反演結果(a) 交叉約束方法結果三維透視圖; (b) 反演結果在L1處的剖面圖; (c) 反演結果在L2處的剖面圖.Fig.12 Inversion result of Changbai Moutain area(a) 3D perspective view of the cross-constrained inversion result; (b) Slice of Line L1 of the inversion result; (c) Slice of Line L2 of the inversion result.

4 結論

本文提出了重力及其垂直梯度聯合的交叉約束聯合密度結構反演方法,通過建立交叉權函數的進行交叉梯度聯合反演,有效提升了結果的分辨率和精度.通過模型試驗證明本文方法能高分辨率地獲取不同深度場源的密度變化特征,尤其是相對以往方法對于較深場源的反演精度大大提升.此外,本文方法的抗噪性較強,對復雜地質體適用性良好.將本文方法應用于遼源礦區的實際數據反演,圈定了6個采空區,其位置與電法資料以及已知礦點完全吻合,且描繪出了大致范圍.通過本文方法對長白山火山區的衛星數據的應用,我們估計了巖漿囊的位置以及范圍,其結果與前人研究結果基本一致.實際數據應用證明了本文方法的實用性,為區域規劃建設及火山災害監測提供了準確的依據.

致謝感謝審稿專家提出的修改意見.