融入寬范圍物性約束策略的大地電磁和重力的Gramian約束聯(lián)合反演

曾志文, 陳曉, 韓江濤, 郭冬, 鄧居智, 張志勇, 郭一豪

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院, 長春 130026 2 東華理工大學(xué)地球物理與測控技術(shù)學(xué)院, 南昌 330013 3 安徽省勘查技術(shù)院, 合肥 230041

0 引言

地球物理反演的多解性、單一地球物理方法的局限性以及地質(zhì)情況的復(fù)雜性,決定了綜合地球物理研究的必要性,聯(lián)合反演發(fā)揮的優(yōu)勢也越加突顯(楊文采,2002;劉光鼎,2005;底青云等,2020).現(xiàn)階段結(jié)構(gòu)約束聯(lián)合反演發(fā)展比較成熟,最具代表性的如基于交叉梯度方法的聯(lián)合反演(Gallardo and Meju,2003;Wu et al.,2022).巖石物性約束聯(lián)合反演雖然耦合效果強(qiáng),但因?yàn)槲镄约s束方式對(duì)先驗(yàn)信息要求更高,存在難以實(shí)施或產(chǎn)生有偏差的地球物理模型的風(fēng)險(xiǎn)(Tu and Zhdanov,2021),因此研究程度較低.

基于巖石物性約束的聯(lián)合反演,需要建立不同巖石物性之間的理論、經(jīng)驗(yàn)或統(tǒng)計(jì)相關(guān)性的關(guān)聯(lián).雖然有Faust公式(Graul,1987)可以將電阻率和速度聯(lián)系起來,以及Archie公式(Archie,1942)可以將電導(dǎo)率和速度聯(lián)系起來,還有聯(lián)系地震波速和密度之間的伯奇定律(Birch,1960,1961)等,但是在聯(lián)合反演中更多的是利用鉆孔或測井?dāng)?shù)據(jù)來確定巖石物性參數(shù)的關(guān)聯(lián)特征的(Jegen et al., 2009;Dell′Aversan et al.,2016),這也從側(cè)面說明了巖石物性關(guān)聯(lián)特征的復(fù)雜性和模糊性.針對(duì)巖石物性關(guān)聯(lián)信息不易建立、適用性有限的問題,陳曉等(2016,2017)以基于模擬退火算法的 MT 和地震聯(lián)合反演為例,提出了寬范圍物性約束技術(shù),隨后又將寬范圍物性約束技術(shù)融入“多次建模,綜合約束,分步反演”的聯(lián)合反演新框架中.簡而言之,該約束方式將巖石物性關(guān)聯(lián)融入全局優(yōu)化算法的模型生成和更新環(huán)節(jié),直接生成在一定范圍內(nèi)相耦合的地球物理模型參與反演運(yùn)算,進(jìn)而實(shí)現(xiàn)巖石物性參數(shù)的耦合;張磊(2016)針對(duì)巖石物性關(guān)聯(lián)關(guān)系復(fù)雜的情況,提出了隨機(jī)正反比的寬范圍巖石物性約束方案,增加了該技術(shù)的靈活性;郭曼(2018)將寬范圍物性約束技術(shù)引入到基于模擬退火算法的 MT 和重力貝葉斯聯(lián)合反演當(dāng)中,進(jìn)一步擴(kuò)展了該技術(shù)的應(yīng)用范圍.曾志文等(2020)實(shí)現(xiàn)了基于差分進(jìn)化算法的MT和重力的寬范圍物性約束聯(lián)合反演,進(jìn)一步驗(yàn)證了該約束方式在全局優(yōu)化算法中的適用性.

上述分析可以看出,巖石物性統(tǒng)計(jì)特征的復(fù)雜性是客觀存在的.研發(fā)容易實(shí)現(xiàn)的、具有一定容錯(cuò)性的巖石物性約束方式,是該領(lǐng)域的前沿和需求.寬范圍物性約束技術(shù)具有在一定程度上可以降低先驗(yàn)信息的要求、提高巖石物性約束方式的容錯(cuò)率等特點(diǎn),但是該技術(shù)目前只在全局優(yōu)化算法中得以實(shí)現(xiàn).另外,與全局優(yōu)化算法相比,梯度優(yōu)化算法具有計(jì)算效率高、速度快等優(yōu)點(diǎn),實(shí)際中應(yīng)用更為廣泛,然而如何將其推廣至梯度優(yōu)化算法聯(lián)合反演至今沒有明確的策略.在全局優(yōu)化算法中,普遍存在解空間,易于實(shí)現(xiàn)將反演解限定在一定空間內(nèi).而梯度優(yōu)化算法如若強(qiáng)制性給定解空間,則會(huì)直接影響到解的搜索.如何將范圍約束融入到聯(lián)合算法尋優(yōu)過程,以及先驗(yàn)信息如何融入都是亟待解決的科學(xué)問題.

此外,Zhdanov等(2012)提出了Gramian約束方式,其本質(zhì)為參與地球物理聯(lián)合反演約束的向量組滿足線性相關(guān).當(dāng)參與地球物理聯(lián)合反演約束的向量組由物性參數(shù)向量組成時(shí),其可視為巖石物性約束.Zhu等(2015)使用Gramian約束求解考慮剩余磁化強(qiáng)度下的總磁場強(qiáng)度(TMI)反演問題.Lin和Zhdanov(2017,2019)將該約束用于鹽丘模型的速度-密度結(jié)構(gòu)反演中,并且應(yīng)用于實(shí)際勘探.Tu和Zhdanov(2020)開展基于Gramian約束的地震和重力聯(lián)合反演方法,在黃石火山口處得到與前人研究成果相同的低密度和低速異常結(jié)果.Malovichko等(2020)通過轉(zhuǎn)換算子將電導(dǎo)率與速度的轉(zhuǎn)換關(guān)系一起寫入Gramian約束中,將電阻率模型作為已知的巖石物理模型指導(dǎo)三維地震全波形反演,但并無涉及模型轉(zhuǎn)換和范圍約束等相關(guān)研究.Gramian約束對(duì)先驗(yàn)信息依賴較低,是比較適用于先驗(yàn)信息儲(chǔ)備低的勘探新區(qū).但若研究區(qū)域勘探程度較高,研究者儲(chǔ)備有一定的先驗(yàn)信息,此時(shí)就需要可以融入先驗(yàn)信息的聯(lián)合反演技術(shù).以Gramian約束為例進(jìn)行分析,Gramian約束是一種“松約束”,不需要明確參與約束的向量之間的關(guān)聯(lián)系數(shù),這一特點(diǎn)降低了Gramian約束應(yīng)用的門檻,但從另一個(gè)角度分析,假如研究區(qū)域有明確的先驗(yàn)物性統(tǒng)計(jì)信息,利用這種“松約束”方式也無法將先驗(yàn)信息確切地融入聯(lián)合反演.但若研究區(qū)域勘探程度較高,此時(shí)就需要可以融入先驗(yàn)信息的聯(lián)合反演技術(shù).

基于此,本文嘗試提出適用于梯度優(yōu)化算法的寬范圍物性約束聯(lián)合反演策略,并將先驗(yàn)信息融入到Gramian約束聯(lián)合反演中,旨在提高先驗(yàn)信息利用率的同時(shí),降低聯(lián)合反演技術(shù)的門檻,進(jìn)而為綜合地球物理聯(lián)合反演提供新的思路.

1 面向梯度優(yōu)化算法的寬范圍物性約束策略

1.1 基本思想

寬范圍物性約束的基本思想是將巖石物性先驗(yàn)信息與所采用的優(yōu)化算法相結(jié)合,不再簡單地將先驗(yàn)物性關(guān)聯(lián)映射獲得的模型直接代入聯(lián)合反演運(yùn)算,而是在一定范圍內(nèi)進(jìn)行再搜索,既能發(fā)揮先驗(yàn)物性信息的導(dǎo)向作用,又可以充分利用優(yōu)化算法的尋優(yōu)能力.

1.2 基于梯度優(yōu)化算法的寬范圍物性約束策略

如何將寬范圍物性約束思維融入梯度優(yōu)化算法中是值得挑戰(zhàn)的.本文提出了適用于梯度優(yōu)化算法的“模型轉(zhuǎn)換、范圍約束、耦合項(xiàng)”的寬范圍物性約束基本策略,可見示意圖1,具體如下:

圖1 寬范圍物性約束策略示意圖Fig.1 Schematic diagram of wide-range petrophysical constraints

第一,先驗(yàn)物性關(guān)聯(lián)為導(dǎo)向.這一點(diǎn)與全局優(yōu)化算法的寬范圍物性約束是一致的,但是引入的方式不同.在全局優(yōu)化算法中,可以直接利用物性關(guān)聯(lián)實(shí)現(xiàn)不同地球物理模型的轉(zhuǎn)換(陳曉等,2016).而在梯度優(yōu)化算法中,可以將物性轉(zhuǎn)換之后模型作為初始模型,或者可以將物性轉(zhuǎn)換之后的模型寫入聯(lián)合反演目標(biāo)函數(shù)中.

第二,耦合項(xiàng).正如典型的耦合方式,交叉梯度、Gramian約束、余弦相似度(Shi et al.,2018)等,在以往的文獻(xiàn)中學(xué)者們都是使用梯度、共軛梯度或者牛頓法等來實(shí)現(xiàn)關(guān)于耦合項(xiàng)的目標(biāo)函數(shù)的極小化.可以看出,這種帶有耦合項(xiàng)的方式更適合于梯度優(yōu)化算法.因此有必要將典型的耦合方式引入,進(jìn)一步提高梯度優(yōu)化算法聯(lián)合反演結(jié)果的耦合效果.需要指出的是,以Gramian約束為例,傳統(tǒng)的Gramian約束不需要明確參與約束的物性向量之間的相關(guān)系數(shù),這既是其優(yōu)點(diǎn),也同樣是缺點(diǎn).正是因?yàn)樗恍枰_定這些相關(guān)系數(shù),所以即使有明確的先驗(yàn)關(guān)聯(lián)信息,也無法將這些系數(shù)的信息以及實(shí)際中的約束關(guān)系,引入到聯(lián)合反演中來.因此,目前這些典型的耦合方式都無法將明確的先驗(yàn)信息融入聯(lián)合反演.

第三,范圍約束.與全局優(yōu)化算法本身就需要設(shè)置待解參數(shù)的解空間不同,梯度優(yōu)化算法是根據(jù)目標(biāo)函數(shù)的梯度來更新地球物理模型,一般需要確定解更新的方向和步長.梯度優(yōu)化算法在實(shí)現(xiàn)的過程中,如早期采用的最簡單的約束方式,人為地不顧梯度方向,直接限制解的上下限,但這樣的操作會(huì)直接影響解的搜索.故此,基于梯度優(yōu)化算法的寬范圍物性約束需要可以將解限制在一定范圍內(nèi)的方法技術(shù),如:懲罰函數(shù)(Kim et al.,1999),對(duì)數(shù)轉(zhuǎn)換(Commer and Newman,2008)等.

2 寬范圍約束聯(lián)合反演的實(shí)現(xiàn)

2.1 目標(biāo)函數(shù)及實(shí)現(xiàn)

以MT和重力聯(lián)合反演為例,介紹面向梯度優(yōu)化算法的寬范圍物性約束策略的具體實(shí)現(xiàn).基于懲罰函數(shù)的Gramian約束聯(lián)合反演目標(biāo)函數(shù)可寫為:

(1)

其中,Pα(m(1),m(2))表示目標(biāo)函數(shù);m表示模型參數(shù)向量;φ(i)表示數(shù)據(jù)擬合泛函;sMN表示模型穩(wěn)定泛函;α表示正則化因子;β是Gramian約束項(xiàng)的權(quán)重系數(shù);sG是Gramian約束耦合項(xiàng);μ是懲罰函數(shù)項(xiàng)的權(quán)重系數(shù);P(m(i))是懲罰函數(shù)約束項(xiàng).

2.1.1 Gramian約束

(2)

轉(zhuǎn)換算子T可以最簡單地取為單位陣,此時(shí)就是基于巖石物性線性關(guān)聯(lián)的約束.

2.1.2 范圍約束

已有研究表明,雖然懲罰函數(shù)法需要確定懲罰函數(shù)的權(quán)重系數(shù),但更加靈活,適用于物性關(guān)聯(lián)特征復(fù)雜的情況(陳曉等,2023).基于此,本文采用懲罰函數(shù)法來實(shí)現(xiàn)物性參數(shù)的范圍約束.

模型參數(shù)向量滿足約束范圍:

ak≤mk≤bk,

(3)

其中ak和bk分別表示模型參數(shù)mk的最小值和最大值.將其以罰函數(shù)的形式寫入到目標(biāo)函數(shù)當(dāng)中,有:

(4)

其中,P為懲罰函數(shù)約束項(xiàng),q是不等式的總數(shù),對(duì)于上述不等式(3),在這里q=2,hi(mk)的形式為:

h1(mk)=mk-ak,

(5)

h2(mk)=bk-mk.

(6)

顯然,若模型參數(shù)向量滿足式(3),那么不等式約束項(xiàng)等于0,即不起作用.當(dāng)超出界限范圍時(shí),才會(huì)起到約束作用.

2.2 寬范圍約束聯(lián)合反演的流程

MT和重力的寬范圍物性約束聯(lián)合反演的流程圖可參見圖2.本文擬以聯(lián)合反演的重力結(jié)果為例,探討寬范圍物性約束策略的效果.如圖2的右半部分所示,介紹如下:

圖2 聯(lián)合反演流程圖Fig.2 Flow chart of joint inversion

(1) 首先利用MT數(shù)據(jù)進(jìn)行單獨(dú)反演,得到電阻率結(jié)果.

(2) 根據(jù)巖石物性關(guān)聯(lián)的先驗(yàn)信息,由電阻率結(jié)果映射到密度模型,并將此模型作為初始模型引入到密度聯(lián)合反演.

(3) 利用電阻率結(jié)果和密度模型計(jì)算Gramian約束項(xiàng),并加入范圍約束項(xiàng),將密度的先驗(yàn)信息分布范圍用于約束密度聯(lián)合反演過程.

(4) 密度更新,是否達(dá)到最大迭代次數(shù),否則重復(fù)步驟(3)和(4),是則進(jìn)入下一步.

(5) 輸出密度聯(lián)合反演結(jié)果.

3 模型試驗(yàn)

3.1 “精準(zhǔn)”先驗(yàn)信息模型試驗(yàn)

為了驗(yàn)證寬范圍約束在梯度優(yōu)化算法聯(lián)合反演中的效果,首先設(shè)計(jì)了在先驗(yàn)信息“精準(zhǔn)的”下的模型試驗(yàn).真實(shí)模型見圖3a和3b,在背景值為100 Ωm、0.01 g·cm-3的均勻半空間中,存在兩個(gè)異常體,從左到右剩余密度分別為0.2、0.05 g·cm-3,電阻率值分別為5、20 Ωm,密度約束范圍給定為[0,0.2].圖4e中的黑線是MT和重力巖石物性的先驗(yàn)信息關(guān)聯(lián).可以看到,此時(shí)的先驗(yàn)信息是相對(duì)“精準(zhǔn)的”.首先,對(duì)MT數(shù)據(jù)進(jìn)行40次的單獨(dú)反演,得到電阻率反演結(jié)果圖3d.然后,通過先驗(yàn)物性關(guān)聯(lián)映射得到密度聯(lián)合反演的初始模型,再進(jìn)行聯(lián)合反演.

圖3 二維電阻率和剩余密度簡單塊體模型(a) 剩余密度模型; (b) 電阻率模型; (c) 第一次迭代密度聯(lián)合反演結(jié)果; (d) 電阻率反演結(jié)果.Fig.3 Two-dimensional resistivity and residual density simple block model(a) Residual density model; (b) Resistivity model; (c) Density joint inversion result of the first iteration; (d) Resistivity inversion result.

圖4 精準(zhǔn)先驗(yàn)信息下的4種方案反演結(jié)果(a) 方案1密度聯(lián)合反演結(jié)果; (b) 方案2密度聯(lián)合反演結(jié)果; (c) 方案3密度聯(lián)合反演結(jié)果; (d) 方案4密度聯(lián)合反演結(jié)果; (e) 第一次迭代物性耦合圖; (f) 物性耦合圖; (g) 均方誤差曲線; (h) 重力異常擬合曲線.Fig.4 Inversion results of four schemes with precise prior information(a) Density joint inversion result of Scheme 1; (b) Density joint inversion result of Scheme 2; (c) Density joint inversion result of Scheme 3; (d) Density joint inversion result of Scheme 4; (e) Petrophysical coupling diagram of the first iteration; (f) Petrophysical coupling diagram; (g) Mean square error curves; (h) Gravity anomaly fitting curves.

基于上述基礎(chǔ),設(shè)計(jì)了4種對(duì)比試驗(yàn)方案(見表1).具體而言:方案1,不加Gramian約束項(xiàng)和不加懲罰函數(shù)的密度聯(lián)合反演;方案2,加Gramian約束項(xiàng)和不加懲罰函數(shù)的密度聯(lián)合反演;方案3,不加Gramian約束項(xiàng)和加懲罰函數(shù)的密度聯(lián)合反演;方案4,加Gramian約束項(xiàng)和加懲罰函數(shù)的密度聯(lián)合反演.需要指出的是,為了更單純地對(duì)比懲罰函數(shù)項(xiàng)、Gramian約束項(xiàng)對(duì)聯(lián)合反演效果的影響,本文在電阻率反演中暫不考慮密度結(jié)果的耦合效果,即關(guān)閉了密度對(duì)電阻率的耦合通道,進(jìn)而保證密度聯(lián)合反演在相同的條件下進(jìn)行.

表1 4種聯(lián)合反演方案Table 1 Four joint inversion schemes

圖3c是4種方案的第一次迭代結(jié)果,圖4e中叉號(hào)點(diǎn)是其對(duì)應(yīng)的物性耦合圖.可以看到,由于先驗(yàn)信息是相對(duì)“精準(zhǔn)的”,第一次迭代時(shí)的密度和電阻率已經(jīng)獲得較好的耦合.在方案1(不加Gramian約束和不加懲罰函數(shù)約束)的條件下,密度聯(lián)合反演結(jié)果(圖4a)以及耦合圖(圖4f中的灰圓點(diǎn))顯示,密度和電阻率的耦合呈現(xiàn)線性關(guān)系但偏離了真實(shí)分布,而且存在密度聯(lián)合反演結(jié)果超出了密度先驗(yàn)信息范圍的情況;方案2(只加Gramian約束)的條件下,密度聯(lián)合反演結(jié)果(圖4b)以及耦合圖(圖4f中的矩形點(diǎn))顯示,密度和電阻率分布是相關(guān)的,但也存在密度超過先驗(yàn)信息范圍的現(xiàn)象;在方案3(只加懲罰函數(shù)約束)的條件下,密度聯(lián)合反演結(jié)果(圖4c)以及耦合圖(圖4f中的加號(hào)點(diǎn))顯示,密度值被有效地約束在范圍之內(nèi);在方案4(加Gramian約束和加懲罰函數(shù))的條件下,密度聯(lián)合反演結(jié)果(圖4d)以及耦合圖(4f中的黑圓點(diǎn))顯示,密度和電阻率分布是相關(guān)的,且密度值被有效地約束在范圍之內(nèi).

由上述分析,Gramian約束和懲罰函數(shù)是有效的,在先驗(yàn)信息“精準(zhǔn)的”情況下,可以實(shí)現(xiàn)物性參數(shù)在一定范圍內(nèi)進(jìn)行耦合,也驗(yàn)證了寬范圍約束策略在梯度優(yōu)化算法中的有效性.

3.2 “不精準(zhǔn)”先驗(yàn)信息的簡單模型試驗(yàn)

上述3.1節(jié)驗(yàn)證了寬范圍約束策略在先驗(yàn)信息“精準(zhǔn)的”情況下的有效性,本節(jié)則在先驗(yàn)信息“不精準(zhǔn)的”情況下,驗(yàn)證寬范圍約束策略的效果.設(shè)計(jì)的真實(shí)模型與3.1節(jié)一樣,圖6e中的黑線是MT和重力巖石物性的先驗(yàn)信息關(guān)聯(lián).可以看到,此時(shí)的先驗(yàn)信息是“不精準(zhǔn)的”.同樣,首先對(duì)MT數(shù)據(jù)進(jìn)行40次的單獨(dú)反演,得到電阻率反演結(jié)果圖5b.然后通過先驗(yàn)物性關(guān)聯(lián)映射,得到密度聯(lián)合反演的初始模型,采用和3.1節(jié)相同的4種對(duì)比試驗(yàn)方案.本節(jié)所有圖件色標(biāo)也與3.1節(jié)模型試驗(yàn)一致.

圖5 初步的單獨(dú)反演結(jié)果(a) 第一次迭代密度聯(lián)合反演結(jié)果; (b) 電阻率反演結(jié)果.Fig.5 Preliminary separate inversion results(a) Density joint inversion result of the first iteration; (b) Resistivity inversion result.

圖5a是4種方案的第一次迭代結(jié)果,圖6e是對(duì)應(yīng)的第一次物性耦合圖.此時(shí),由于先驗(yàn)信息是相對(duì)“不精準(zhǔn)的”,第一次迭代時(shí)的密度和電阻率的耦合是偏離真實(shí)物性值分布的.在方案1(不加Gramian約束和不加懲罰函數(shù)約束)的條件下,密度聯(lián)合反演結(jié)果(圖6a)不能較好地反映異常體的位置,耦合圖(圖6f中的灰圓點(diǎn))分布較亂,且偏離了真實(shí)物性的分布;方案2(只加Gramian約束)的條件下,密度聯(lián)合反演結(jié)果(圖6b)較好地還原了兩個(gè)異常體的位置,耦合圖(圖6f中的矩形點(diǎn))顯示密度和電阻率分布是相關(guān)的,但存在密度超過先驗(yàn)信息范圍的現(xiàn)象;

在方案3(只加懲罰函數(shù)約束)的條件下,密度聯(lián)合反演結(jié)果(圖6c)以及耦合圖(圖6f中的加號(hào)點(diǎn))顯示,密度值被有效地約束在范圍之內(nèi),但這個(gè)物性耦合顯然是受“不精準(zhǔn)”的先驗(yàn)關(guān)系所影響;在方案4(加Gramian約束和加懲罰函數(shù))的條件下,密度聯(lián)合反演結(jié)果(圖6d)以及耦合圖(圖6f中的黑圓點(diǎn))顯示,密度和電阻率分布是相關(guān)的,且密度值被有效地約束在范圍之內(nèi),較好地還原了真實(shí)模型.

由上述分析可以看出,寬范圍約束策略在先驗(yàn)信息“不精準(zhǔn)的”情況下也是適用的.它可以更好地促進(jìn)物性參數(shù)在一定范圍內(nèi)進(jìn)行耦合,進(jìn)一步提高了巖石物性關(guān)聯(lián)約束方式在梯度優(yōu)化算法中的適用性,并且具有一定的容錯(cuò)率,可以降低對(duì)先驗(yàn)信息的要求.

3.3 “不精準(zhǔn)”先驗(yàn)信息的復(fù)雜模型試驗(yàn)

為了進(jìn)一步驗(yàn)證寬范圍約束策略在復(fù)雜模型中的效果,本節(jié)設(shè)計(jì)了復(fù)雜的4個(gè)塊體試驗(yàn).在背景值為200 Ωm、0.01 g·cm-3的均勻半空間中,存在4個(gè)異常體,從左到右剩余密度分別為0.4、0.1、0.05、0.04 g·cm-3,電阻率值分別為5、20、40、50 Ωm,密度約束范圍給定為[0,0.4].圖8e中黑線是MT和重力巖石物性的先驗(yàn)信息關(guān)聯(lián),模擬的是實(shí)際中可能獲得的“不精準(zhǔn)”的物性關(guān)聯(lián)的情況.與上文一樣,對(duì)MT數(shù)據(jù)進(jìn)行40次的單獨(dú)反演,得到電阻率反演結(jié)果圖7d.然后,通過黑線的先驗(yàn)物性關(guān)聯(lián)映射,得到密度聯(lián)合反演的初始模型.和3.2節(jié)一樣,設(shè)計(jì)相同的4種對(duì)比試驗(yàn)方案.

圖8 復(fù)雜模型的4種方案反演結(jié)果(a) 方案1密度聯(lián)合反演結(jié)果; (b) 方案2密度聯(lián)合反演結(jié)果; (c) 方案3密度聯(lián)合反演結(jié)果; (d) 方案4密度聯(lián)合反演結(jié)果; (e) 第一次迭代物性耦合圖; (f) 物性耦合圖; (g) 均方誤差曲線; (h) 重力異常擬合曲線.Fig.8 Inversion results for four schemes of complex models(a) Density joint inversion result of Scheme 1; (b) Density joint inversion result of Scheme 2; (c) Density joint inversion result of Scheme 3; (d) Density joint inversion result of Scheme 4; (e) Petrophysical coupling diagram of the first iteration; (f) Petrophysical coupling diagram; (g) Mean square error curves; (h) Gravity anomaly fitting curves.

圖7c是4種方案的第一次迭代密度聯(lián)合反演結(jié)果,圖8e是對(duì)應(yīng)的第一次迭代的物性耦合圖.由于先驗(yàn)信息是相對(duì)“不精準(zhǔn)的”,第一次迭代時(shí)的密度和電阻率的耦合也是偏離真實(shí)物性值分布的.與3.2節(jié)的模型試驗(yàn)規(guī)律類似,在沒有Gramian約束的情況下,方案1(圖8f灰圓點(diǎn))和方案3(圖8f加號(hào)點(diǎn))的物性耦合較差,偏離了真實(shí)值分布;在有Gramian約束的情況下,方案2(圖8f矩形點(diǎn))和方案4(圖8f黑圓點(diǎn))密度和電阻率的耦合關(guān)系得到了加強(qiáng),且分布在真實(shí)值附近,而再相比于方案2,加上懲罰函數(shù)的方案4密度值被有效地約束在范圍之內(nèi).但圖中部分?jǐn)?shù)據(jù)點(diǎn)也顯示出嚴(yán)重偏離,這是因?yàn)楸疚牟捎玫氖菃蜗虻拿芏软樞蚵?lián)合反演,電阻率結(jié)果在密度聯(lián)合反演時(shí)不會(huì)更新.輸入不準(zhǔn)的電阻率信息反演結(jié)果,導(dǎo)致了部分?jǐn)?shù)據(jù)的偏離.

可以看出,寬范圍約束策略適用于復(fù)雜的地球物理模型.總而言之,本文嘗試?yán)脦r石物性關(guān)聯(lián)獲得初始模型,并利用懲罰函數(shù)開展范圍約束,成功實(shí)現(xiàn)了基于梯度優(yōu)化算法的寬范圍物性約束聯(lián)合反演.

4 實(shí)際資料處理

茶亭礦區(qū)位于南陵—宣城區(qū)域盆地內(nèi),是該區(qū)域近些年來找礦的重大發(fā)現(xiàn).燕山期強(qiáng)烈的構(gòu)造-巖漿活動(dòng)使該區(qū)受到巖漿-熱液作用的疊加改造,為區(qū)內(nèi)礦床的形成創(chuàng)造了有利條件(肖慶玲等,2018;陶龍等,2019).礦區(qū)大部為中分村組火山巖覆蓋,主要圍巖成分為花崗閃長巖,還有各類侵入巖體等.根據(jù)參考文獻(xiàn)(洪大軍等,2019),茶亭礦區(qū)的部分礦石平均密度和平均電阻率見表2.圖9是通過巖石物性統(tǒng)計(jì)表2建立的電導(dǎo)率和剩余密度巖石物性關(guān)聯(lián)圖,鑒于此區(qū)域的巖石物性關(guān)系比較雜,本文利用了兩種不同函數(shù)關(guān)系進(jìn)行擬合.

表2 茶亭礦區(qū)部分巖石物性參數(shù)Table 2 Some petrophysical parameters in the Chating mine area

圖9 電導(dǎo)率和密度的巖石物性關(guān)聯(lián)圖Fig.9 Petrophysical parameters correlation diagram of the conductivity and density

4.1 先驗(yàn)地質(zhì)地球物理模型

圖10a是由安徽省勘查技術(shù)院提供的茶亭區(qū)域的地質(zhì)解譯圖,該結(jié)果是結(jié)合重磁反演、地質(zhì)模型、鉆孔數(shù)據(jù)的解譯而來.謝巧勤等(2020)指出茶亭銅金礦床的礦體主要賦存于石英閃長玢巖侵入體(殼?；煸?且以幔源為主)包孕的角礫巖筒中,并根據(jù)地球化學(xué)的手段分析其為角礫巖型礦床.田自強(qiáng)(2020)指出茶亭銅金礦床的石英閃長玢巖是在侵位之后短時(shí)間內(nèi)被隆升至近地表,再覆蓋中分村組火山巖.此外,鉆孔資料顯示稻山村逆斷層上盤的花崗斑巖為無根巖體,指示該無根巖體不是原地侵位的產(chǎn)物,應(yīng)是逆斷層活動(dòng)過程中從深部運(yùn)移而來的某個(gè)巖體的一部分.

圖10 先驗(yàn)?zāi)Ｐ?a) 地質(zhì)解譯圖; (b) 密度參考模型.Fig.10 Priori model(a) Geological interpretation map; (b) Density reference model.

可以看出,圖10a基本符合實(shí)際地質(zhì)地球物理情況,具有較高的合理性且有明顯的地質(zhì)含義,可以作為本文開展地球物理聯(lián)合反演的先驗(yàn)地質(zhì)地球物理模型.圖10b是依據(jù)地質(zhì)解譯圖及表2的巖石物性統(tǒng)計(jì)表建立的密度參考模型.但此模型只有淺部5 km以上的密度構(gòu)造,并沒有5 km以下的信息.

圖11a是對(duì)該測線的MT數(shù)據(jù)進(jìn)行二維單獨(dú)反演得到的電阻率結(jié)果.結(jié)合參考模型5 km以上的信息分析(圖中的黑線代表構(gòu)造線),可以看出,在19～24 km處的淺部兩個(gè)侵入巖體處大致呈現(xiàn)中高阻反映,但茶亭礦床處侵入巖體的電性連通性不強(qiáng),并且30～36 km處的稻山村的淺部和沙溪深部的巖體對(duì)應(yīng)性有待進(jìn)一步驗(yàn)證.

圖11 聯(lián)合反演結(jié)果(a) 電阻率單獨(dú)反演結(jié)果; (b) 電阻率聯(lián)合反演結(jié)果; (c) 密度初始模型; (d) 密度聯(lián)合反演結(jié)果; (e) 均方誤差曲線; (f) 重力異常擬合曲線.Fig.11 Joint inversion results(a) Separate inversion result of resistivity; (b) Joint inversion result of resistivity; (c) Initial density model; (d) Density joint inversion result; (e) Mean square error curve; (f) Gravity anomaly fitting curve.

綜上分析,目前該測線還存在許多待解決的問題:茶亭礦床附近鉆孔不足2 km, 茶亭銅金礦床的“根”在哪,在深部該礦床是如何展布的還是未知;另外無“根”的花崗斑巖的源在哪;此外,先驗(yàn)地質(zhì)地球物理模型深度僅有5 km,在探索深部構(gòu)造的同時(shí),能否把先驗(yàn)?zāi)Ｐ鸵踩谌脒M(jìn)去.下一步,本文將針對(duì)以上問題,運(yùn)用第2節(jié)所提出的新技術(shù)來處理該測線,嘗試獲得新的認(rèn)識(shí).

4.2 聯(lián)合反演

為了將先驗(yàn)地質(zhì)地球物理模型融入聯(lián)合反演,結(jié)合上文所提出的寬范圍約束聯(lián)合反演新策略,本文設(shè)計(jì)了如下的實(shí)際資料反演方案.首先,將圖11a的MT單獨(dú)反演結(jié)果作為聯(lián)合反演的電阻率初始模型,開展基于分區(qū)域的Gramian約束MT和重力聯(lián)合反演(郭一豪,2020).所謂分區(qū)域,是指電阻率模型和密度參考模型在5 km以上做Gramian約束,而5 km以下則不做要求.在迭代40次之后,得到此次電阻率聯(lián)合反演結(jié)果圖11b.

接下來,按照新技術(shù)反演流程圖2的步驟,將上述的電阻率反演結(jié)果作為聯(lián)合反演的電阻率初始模型,通過圖9的先驗(yàn)關(guān)系映射得到密度的初始模型圖11c.密度模型的約束范圍如下:5 km以上的給予緊約束,將初始模型密度值上下擾動(dòng)±0.02 g·cm-3;5 km以下則給定比較寬松的約束范圍±0.2 g·cm-3.最后,結(jié)合圖11b的電阻率反演結(jié)果,開展基于寬范圍物性約束的重力和MT聯(lián)合反演,迭代80次.

對(duì)比圖11a和11b可以看到,圖11b的19～20 km處,侵入巖體形態(tài)更加連續(xù)和清晰,在23～24 km處,銅金礦床在淺部5 km以上連通起來,這與先驗(yàn)信息相符.另外,在30～36 km、深度0～5 km處存在兩個(gè)電阻率相對(duì)應(yīng)的區(qū)域,結(jié)合密度聯(lián)合反演結(jié)果圖11d,推覆體的密度與沙溪深部的密度分布也是對(duì)應(yīng)的.據(jù)以上分析,此處構(gòu)造巖體應(yīng)該是推覆構(gòu)造的產(chǎn)物,推測與沙溪深部巖體是同源的.此外,侵入巖體的5 km以下的形態(tài)也得到了揭示,呈現(xiàn)出倒立的長鐘狀.

綜上分析,本文提出的聯(lián)合反演新策略,成功地將先驗(yàn)?zāi)Ｐ腿谌刖C合地球物理解釋.聯(lián)合反演結(jié)果所揭示的礫巖筒型銅金礦床在淺部的分布符合先驗(yàn)信息,驗(yàn)證了稻山村逆斷層上盤無根的花崗斑巖為區(qū)域的推覆運(yùn)動(dòng)的產(chǎn)物,勾勒了角礫巖筒礦床在深部的空間分布形態(tài),為茶亭銅金礦床深部找礦提供了深部模型.

5 結(jié)論

本文提出了適用于在梯度優(yōu)化算法的寬范圍物性約束策略,通過巖石物性關(guān)聯(lián)生成初始模型的方式引入先驗(yàn)信息,在范圍約束項(xiàng)和耦合項(xiàng)的共同作用下進(jìn)行聯(lián)合反演,使得聯(lián)合反演既依賴于先驗(yàn)?zāi)Ｐ托畔?又可以在一定程度上“擺脫”先驗(yàn)信息的嚴(yán)格控制,從而使巖石物性參數(shù)在真實(shí)值附近實(shí)現(xiàn)耦合.基于此,開展了模型試驗(yàn)和實(shí)測資料反饋,驗(yàn)證了新策略的適用性和實(shí)用性,具體而言:

(1) 寬范圍物性約束不僅是一種技術(shù),更是一種思維.寬范圍約束策略可以提高先驗(yàn)巖石物性信息利用率,并降低不精確先驗(yàn)信息帶來的風(fēng)險(xiǎn).可以在一定程度上將先驗(yàn)巖石物性關(guān)聯(lián)融入聯(lián)合反演的迭代尋優(yōu)過程,且使巖石物性耦合在一定范圍內(nèi).

(2) 本文提出的“巖石物性關(guān)聯(lián)+范圍約束+耦合項(xiàng)”的思維具有較強(qiáng)的適用性,并且具有結(jié)合其他優(yōu)化算法和約束方式進(jìn)一步擴(kuò)展的潛質(zhì).

(3) 本文提出的寬范圍物性約束策略具有一定的實(shí)用性,新方法驗(yàn)證了茶亭銅金礦床淺部的分布形態(tài),無根的花崗斑巖為推覆構(gòu)造的結(jié)果,勾勒了茶亭銅金礦床的深部結(jié)構(gòu)的空間展布.

需要指出的是,第一,如何將寬范圍物性約束思維拓展到其他的地球物理聯(lián)合反演耦合方式中去,還需要進(jìn)一步試驗(yàn),以及如何將新的地球物理聯(lián)合反演技術(shù)推廣至多地球物理方法的聯(lián)合反演也值得進(jìn)一步研究;第二,實(shí)測資料處理顯示,Gramian約束可以實(shí)現(xiàn)先驗(yàn)信息的融入,也為該地區(qū)的地質(zhì)解譯提供了新視角和新思路.但是,如何評(píng)判先驗(yàn)信息可靠性以及如何更合理地實(shí)現(xiàn)先驗(yàn)信息的融入是值得進(jìn)一步探討的.