重力及其垂直梯度交叉約束聯(lián)合密度結(jié)構(gòu)反演方法

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

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

0 引言

重力梯度數(shù)據(jù)具有較高水平分辨率,相對原始重力場可凸顯較淺地質(zhì)體.重力梯度測量可通過不同高度測量重力值的方式來獲得,其與重力聯(lián)合測量可以探測不同深度的空洞、人防工程(Cordell, 1979; Schlinger, 1990; Bate, 2005; Butler, 2012). 通過重力異常反演地下密度分布可以圈定場源體的分布,主要的反演方法包括光滑約束反演、粗糙約束反演、聚焦反演、協(xié)克里金反演、模糊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).

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

為了充分利用重力與重力梯度對地下密度結(jié)構(gòu)的敏感性差異，我們提出重力及其垂直梯度交叉約束聯(lián)合密度結(jié)構(gòu)反演方法,以重力與重力梯度單獨(dú)反演結(jié)果作為權(quán)函數(shù)來約束交叉梯度聯(lián)合反演過程,有效恢復(fù)地下不同深度場源的密度結(jié)構(gòu)分布特征.最后通過模型實(shí)驗(yàn)和實(shí)際數(shù)據(jù)證明該方法較常規(guī)數(shù)據(jù)聯(lián)合密度反演方法在反演分辨率和精度的優(yōu)勢，以及在實(shí)際應(yīng)用中的適用性.

1 基本原理

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

(1)

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

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

圖1 兩個(gè)棱柱體模型的正演異常及不同反演方法結(jié)果在y=1000 m處剖面圖(a) 兩個(gè)棱柱體組合模型; (b) 重力異常; (c) 重力垂直梯度異常; (d) Tikhonov正則化方法Vz數(shù)據(jù)反演; (e) Tikhonov正則化方法Vzz數(shù)據(jù)反演; (f) Tikhonov正則化方法[Vz, Vzz]T數(shù)據(jù)反演; (g) 交叉梯度Vz數(shù)據(jù)反演; (h) 交叉梯度Vzz數(shù)據(jù)反演; (i) 交叉梯度Vzz數(shù)據(jù)反演結(jié)果三維透視圖.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 兩個(gè)棱柱體模型Table 1 Parameters oftwo-prism model

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

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

(3)

采用交叉梯度的方法進(jìn)行重力數(shù)據(jù)和重力梯度聯(lián)合反演,其目標(biāo)函數(shù)為：

(5)

式中，λz和λzz為交叉約束項(xiàng)系數(shù).式中參數(shù)的下角標(biāo)z和zz分別代表重力和重力梯度數(shù)據(jù)相關(guān)參數(shù).Wz1和Wz2分別為重力及其梯度數(shù)據(jù)對應(yīng)的深度加權(quán)矩陣.采用圖1b、c所示的重力及其梯度數(shù)據(jù)進(jìn)行交叉梯度反演，得到圖1g、h的兩個(gè)反演結(jié)果，深部和淺部分辨率有提高,但與實(shí)際模型依舊存在較大偏差.

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

圖2 交叉約束聯(lián)合密度結(jié)構(gòu)反演技術(shù)流程Fig.2 The flow of the cross-constrained density joint inversion based on gravity and vertical gradient data

圖2中，Tikhonov正則化方法重力數(shù)據(jù)反演結(jié)果mz1和重力垂直梯度數(shù)據(jù)反演結(jié)果mzz1形成的交叉約束權(quán)函數(shù)Wρz和Wρzz的表達(dá)式為：

(6)

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

最后，取交叉約束權(quán)聯(lián)合密度反演結(jié)果mz2和mzz2的平均值作為最終密度分布mFin.

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

圖3 交叉約束聯(lián)合密度結(jié)構(gòu)反演結(jié)果(a) Vz數(shù)據(jù)反演結(jié)果在y=1000 m處剖面圖; (b) Vzz數(shù)據(jù)反演結(jié)果在y=1000 m處剖面圖; (c) 反演最終結(jié)果; (d) 密度大于0.65 g·cm-3反演結(jié)果三維透視圖.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 不同方法反演結(jié)果在y=1000m處的與觀測數(shù)據(jù)擬合曲線Fig.4 The fitting curve map of inversion results of different methodsand observed data at section of y=1000 m

2 模型試驗(yàn)

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

圖5 含噪數(shù)據(jù)正演異常和不同方法反演結(jié)果在y=1000 m處剖面圖(a) 兩個(gè)棱柱體組合模型; (b) 重力異常; (c) 重力垂直梯度異常; (d) Tikhonov正則化方法Vz數(shù)據(jù)反演結(jié)果; (e) Tikhonov正則化方法Vzz數(shù)據(jù)反演結(jié)果; (f) Tikhonov正則化方法[Vz, Vzz]T反演結(jié)果; (g) 交叉梯度Vz數(shù)據(jù)反演結(jié)果; (h) 交叉梯度Vzz數(shù)據(jù)反演結(jié)果; (i) 交叉約束Vz數(shù)據(jù)反演結(jié)果; (j) 交叉約束Vzz數(shù)據(jù)反演結(jié)果; (k) 交叉約束方法反演最終結(jié)果; (l) 交叉約束方法結(jié)果三維透視圖.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 含噪數(shù)據(jù)不同反演方法擬合差Table 2 The misfit of the different inversion methods for data with noise

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

采用不同密度的異常體組合模型來驗(yàn)證交叉約束密度結(jié)構(gòu)反演方法對密度差異區(qū)間較大的復(fù)雜地質(zhì)體的適用性,模型參數(shù)見表3,圖6a為異常體的密度真實(shí)分布.

圖6 三個(gè)棱柱體模型及其正演異常(a) 三個(gè)棱柱體組成的模型; (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 三個(gè)棱柱體模型參數(shù)Table 3 Parameters of three-prism model

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

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

圖7 不同反演方法在x=1400 m, y=1000 m, z=400 m反演結(jié)果圖(a),(b)及(c)為Tikhonov正則化方法Vz數(shù)據(jù)反演; (d), (e)及(f) Tikhonov正則化方法Vzz 數(shù)據(jù)反演; (g), (h)及(i) Tikhonov正則化方法[Vz, Vzz]T反演; (j), (k)及(l) 交叉梯度Vz數(shù)據(jù)反演; (m),(n)及(o) 交叉梯度Vzz數(shù)據(jù)反演; (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 實(shí)際應(yīng)用

3.1 遼源礦區(qū)重力及垂直梯度聯(lián)合反演

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

圖8 遼源礦區(qū)Fig.8 Liaoyuan mining area

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

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

圖10 反演密度分布結(jié)果與高密度電阻率法結(jié)果對比(a) 本文交叉約束方法反演結(jié)果; (b) 交叉約束反演結(jié)果在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 長白山火山區(qū)重力及垂直梯度聯(lián)合反演

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

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

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

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

圖12 長白山區(qū)域反演結(jié)果(a) 交叉約束方法結(jié)果三維透視圖; (b) 反演結(jié)果在L1處的剖面圖; (c) 反演結(jié)果在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 結(jié)論

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

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