連晟,程正璞,羅旋,李敬杰,田蒲源

(1.中國地質(zhì)調(diào)查局 水文地質(zhì)環(huán)境地質(zhì)調(diào)查中心,河北 保定 071051; 2.天津市地熱資源勘查開發(fā)工程研究中心,天津 300300)

0 引言

作為深地資源勘探中重要的定量解釋工作,地球物理反演通過地球物理觀測數(shù)據(jù)(重力數(shù)據(jù)、磁法數(shù)據(jù)、地震數(shù)據(jù)和電磁數(shù)據(jù))反推地下場源的空間分布情況,提供目標地質(zhì)構(gòu)造的物性和幾何形態(tài)等特征信息。然而由于反演的多解性和單一地球物理方法的局限性等問題,導致了反演獲得的結(jié)果與真實地下結(jié)構(gòu)存在不同程度的差異。地球物理反演多解性問題主要是由于有限的觀測數(shù)據(jù)量、觀測數(shù)據(jù)的離散化、觀測數(shù)據(jù)包含的隨機噪聲以及地球物理場自身的等效性。為改善反演多解性和單一方法的局限性,獲得可靠的反演結(jié)果,融合多種地球物理場信息,并充分考慮巖石地球物理測試結(jié)果,進行多種地球物理場聯(lián)合反演被認為是符合實際規(guī)律的一種解決方案。

基于空間分布結(jié)構(gòu)性耦合的聯(lián)合反演方法,由于強調(diào)結(jié)構(gòu)上的相似性,不依賴于巖石物性間關(guān)系,應用更加廣泛。Gallardo等[1]提出了交叉梯度函數(shù),該方法通過對不同物性梯度之間叉乘來識別結(jié)構(gòu)邊界,實現(xiàn)了重、磁、電、震4種方法的交叉梯度聯(lián)合反演;Moorkamp等[2]實現(xiàn)了重力、大地電磁和地震交叉梯度聯(lián)合反演和巖石物性約束反演,并對二者進行了對比研究;李桐林等[3]實現(xiàn)了部分區(qū)域約束下的交叉梯度多種地球物理數(shù)據(jù)聯(lián)合反演。

通過巖石物性實驗室測試得到的地球物理參數(shù)是最為精確的方式之一,在實驗基礎(chǔ)上通過理論推導,部分巖石物性參數(shù)之間存在明確的函數(shù)關(guān)系,例如密度和地震波速度之間存在明確的正比關(guān)系,但大量的地球物理參數(shù)之間并不存在明確的函數(shù)關(guān)系,僅巖石物性測試尺度中存在統(tǒng)計學相關(guān)性。如何將該種來自巖石物理測試獲得的地球物理參數(shù)之間的統(tǒng)計學關(guān)系應用到地球物理聯(lián)合反演中,學者們做了大量研究工作。張磊[4]基于正則化思想利用寬范圍巖石物性約束方案實現(xiàn)了MT與重力聯(lián)合反演。Sun等[5]將巖石地球物理先驗信息在參數(shù)域利用模糊C聚類方案以參考聚類中心的方式融入反演過程,有效提高反演效果。此外,模糊C聚類算法在磁法、重力等數(shù)據(jù)反演中的應用均取得了良好效果[5-6]。易柯等[7]在經(jīng)典最小結(jié)構(gòu)模擬正則化約束的基礎(chǔ)上開展了多地球物理參數(shù)的聯(lián)合反演研究工作。隨著大數(shù)據(jù)的發(fā)展考慮巖性約束的聯(lián)合反演方案應用越來越廣泛,其中利用有限混合模型,將不同巖石樣本物性測試得到的統(tǒng)計學分布特征用簡單的概率密度函數(shù)模擬,使得分析結(jié)果具有較好的空間連續(xù)性和穩(wěn)定性[8]。學術(shù)及實踐中發(fā)展了大量類型的有限混合模型,其中最具代表性的是高斯混合模型。高斯混合模型憑借形式簡單、計算方便等優(yōu)勢,成為廣泛應用于科學研究有限混合模型[9],Astic等[9]使用高斯混合模型實現(xiàn)了用巖石先驗信息引導重力、磁法聯(lián)合一維反演。Di-Giuseppe等[10-11]基于共網(wǎng)格技術(shù)提出了一種后反演(post-inverstion)方案,該方案將k均值聚類分析技術(shù)應用于先前獲得的單變量反演的二維地震折射層析和可控源音頻大地電磁反演剖面,通過定量分析不同參數(shù)間的相關(guān)程度,識別了聚類后的數(shù)據(jù)集與巖性的對應關(guān)系。Li等[6]將聚類技術(shù)應用于大地電磁三維反演結(jié)果,識別和歸類了反演結(jié)果中的地質(zhì)構(gòu)造,達到了地質(zhì)分異(geology differentiation)的目的。

本研究實現(xiàn)了高斯混合模型進行巖石物理和地質(zhì)引導的多種地球物理場聯(lián)合反演以及基于共網(wǎng)格的單變量反演模型的多元物理信息k均值聚類的地質(zhì)分異,達到了多種地球物理場聯(lián)合反演及數(shù)據(jù)融合分析的目的。本文僅就聯(lián)合反演部分進行展開,將高斯混合模型作為巖石地球物理的約束方案,在地震勘查的結(jié)構(gòu)基礎(chǔ)上,聯(lián)合重磁、大地電磁法數(shù)據(jù)進行聯(lián)合反演。

1 重力、磁法和大地電磁正反演理論

本文重力正演計算采用曾華霖[12]給出的解析算法,磁法正演使用管志寧[13]給出的算法。大地電磁正演使用基于六面體網(wǎng)格剖分的矢量有限元法[14]。

φ=φd+λφm,

(1)

其中:φd為數(shù)據(jù)擬合項,表達式如下:

(2)

φm為模型平滑約束項:

(3)

正則化反演的實現(xiàn)可以認為是一個尋優(yōu)過程,其目標是找到一個使目標函數(shù)最小的地球物理模型m,該地球物理優(yōu)化問題的形式為[13]

(4)

(5)

由圖6可知，降低購房者購買被動房的增量成本可有效降低購房者購買普通房的概率，即提高購房者購買被動房的意愿，間接說明可提高對購房者激勵的有效性，反之，若購買被動房的增量成本增加會降低被動房的購買意愿。

2 高斯混合模型及巖性約束

為實現(xiàn)巖石物性對地球物理反演的約束,首先需要將巖石物理測量的物性參數(shù)如密度、磁導率和電導率等測量結(jié)果進行數(shù)字化表達,本次研究使用了高斯分布來描述某種巖石物性參數(shù)的概率分布規(guī)律,同時利用不同巖性的高斯分布相疊加后得到多巖性單地球物理參數(shù)的高斯混合模型。若M個巖石樣本均測量了4個地球物理參數(shù)(密度、磁化率、電阻率、速度),這些巖石樣本可以分為泥巖、砂巖、花崗巖3種類型,對應每個地球物理參數(shù)會有4個高斯混合模型,每個高斯混合模型有K=3個子高斯模型,其概率密度函數(shù)[9]為:

(6)

式中:K代表混合模型中子高斯模型的數(shù)量,即所測量巖石類別的個數(shù);πk是第k個高斯分量的先驗權(quán)值,該權(quán)值滿足條件:

(7)

φ(x∣θk)是第k個高斯分量的高斯分布密度函數(shù),表達式如下:

(8)

3 多地球物理巖性約束反演

地球物理聯(lián)合反演問題使用Tikhonov反演方法,該方法將反演問題轉(zhuǎn)化為一個最優(yōu)化問題,其中描述巖石物性參數(shù)分布的高斯混合模型以負對數(shù)形式作為Tikhonov反演中的誤差估計[9],該誤差可以通過當前模型mi和參考模型mref之間的最小二乘誤差來近似,該項誤差和最小系數(shù)矩陣Ws需要在每次迭代時更新。動態(tài)參考模型和系數(shù)矩陣根據(jù)當前地球物理模型以及地質(zhì)結(jié)構(gòu)及巖石物理先驗信息來確定[9,17]。

3.1 多地球物理巖性約束反演目標函數(shù)

多地球物理屬性下巖性約束反演的目標函數(shù)用如下的最小二乘形式表示擬合誤差:

(9)

其中:m為地球物理模型;β、αs、αν,p為正則化參數(shù);φd(m)為包含多個地球物理屬性的擬合誤差;φs(m)為針對各模型的高斯混合模型擬合誤差;φν,p(m)是對每個方向和物理屬性的平滑項,其表達式分別為:

(10)

(11)

(12)

3.2 正則化參數(shù)的計算方式

反演目標函數(shù)中包含3個正則化參數(shù)β、αs、αν,p,計算中采取保持αν,p不變,僅更新β、αs來求解式(9)。在平滑項中,αν,p作為尺度參數(shù)控制著空間導數(shù)項,αν,p的數(shù)值采用對每個物理屬性進行加權(quán)獲得,用于均衡不同尺度(如密度、電導率和磁化率等)時,均衡平滑項對目標函數(shù)值的貢獻,具體加權(quán)系數(shù)根據(jù)經(jīng)驗給定。

正則化參數(shù)β、αs的計算采用Astic等[9]的計算策略,對于多物理場情況,本次研究通過以下方式對其進行改變:當所有的地球物理數(shù)據(jù)擬合誤差大于或者等于其目標值時,減小β的數(shù)值,當全部地球物理數(shù)據(jù)擬合誤差都等于或低于其目標值時,則增加參數(shù)αs的數(shù)值。

則參數(shù)αs按照式(13)增加其數(shù)值:

(13)

4 溫度預測方法

本文選用與溫度最為關(guān)聯(lián)的電導率參數(shù)來進行溫度場計算,溫度與電導率之間的關(guān)系根據(jù)Arrhenius定律[18]表示為:

σ=σ0e-E0/kT,

(14)

其中:σ為電導率(巖石物理試驗實測電導率),S/m;k為波爾茲曼常數(shù)0.8617×10-4ev/°K;E0為活化能eV;T為開氏溫度;σ0為前指因子,不同溫度、不同巖石其值不同;e為自然常數(shù)。這些參數(shù)可由實驗室實驗測得[19],一些巖石的活化能E0和前指因子σ0可參見表1。

表1 巖漿巖的活化能E0值與常數(shù)系數(shù)logσ0試驗室測量結(jié)果

5 模型試算

基于以上反演理論及算法,采用C語言聯(lián)合Python語言實現(xiàn)了基于巖性引導的聯(lián)合反演計算,為檢驗算法有效性,本文進行了模型的單獨反演和聯(lián)合反演模型試驗,設計模型如下:在電阻率為1 000 Ω·m、相對密度為0 g/cm3、異常體相對磁化率為0 SI的均勻半空間介質(zhì)中,設置長方體異常體,長寬高為1 km×1 km×0.4 km,頂面距離地面1 km,相對密度為-0.2 g/cm3,相對磁化率為0.2 SI,電阻率為10 Ω·m,異常體的空間分布見圖1所示;該理論模型由地下兩部分組成,即塊狀體和背景值,兩部分屬性值的分布符合高斯混合模型分布,其中塊狀體密度為20%,背景值的分布密度為80%,密度/磁化率和密度/電阻率二維分布及高斯混合模型分布見圖2。

圖1 正演模型參數(shù)

大地電磁正演觀測數(shù)據(jù)頻率為100～0.01 Hz之間的10個對數(shù)間隔頻率,在正演計算結(jié)果中加入2%隨機誤差的合成數(shù)據(jù);重磁反演的觀測數(shù)據(jù)來自于在重磁三維理論正演計算結(jié)果中加入2%隨機誤差的合成數(shù)據(jù),通過對重、磁、電合成數(shù)據(jù)進行基于高斯混合模型先驗信息約束下的聯(lián)合反演,反演的初始模型均使用均勻半空間模型,其半空間物性參數(shù)為:1e-4g/cm3、1e-4SI和1 000 Ω·m。磁法正反演過程所假設的磁性參數(shù)為:地球磁場強度50 000 nT與剖面夾角0°、剖面內(nèi)磁傾角90°、地磁場傾角90°。

由以上3種方法的單獨反演結(jié)果與聯(lián)合反演結(jié)果對比(圖3所示),可以在電阻率模型單獨反演結(jié)果中看出異常面積大于真實模型,異常上下邊界均偏差較大;在密度和磁化率單獨反演結(jié)果中,無法準確識別異常體上下界面,異常體面積偏大,范圍較理論模型有較大差距,異常體下部出現(xiàn)了一定范圍的假異常,表明該類方法垂向分辨率較弱。通過聯(lián)合反演得到的電阻率、密度、磁化率結(jié)果(圖3b、d、f)顯示與理論模型數(shù)值相比較,異常體邊界位置、幾何尺寸和物性參數(shù)值均更加接近理論模型,密度和磁化率聯(lián)合反演結(jié)果上下界面位置比單獨反演定位結(jié)果更加準確。通過以上分析,證明聯(lián)合反演效果好于單獨反演。

a—重力單獨反演結(jié)果;b—長方體模型聯(lián)合反演密度結(jié)果;c—磁法單獨反演結(jié)果;d—長方體模型聯(lián)合反演磁化率結(jié)果;e—大地電磁單獨反演結(jié)果圖;f—長方體模型聯(lián)合反演電阻率結(jié)果

6 實際案例

為了驗證該項技術(shù)的應用效果,在青海共和盆地進行了應用。青海共和盆地位于青藏高原東北緣,為高溫地熱異常盆地[20],現(xiàn)已建成我國干熱巖勘探與開發(fā)示范研究區(qū)[21-22]。研究區(qū)位于青海共和盆地二級構(gòu)造單元切吉凹陷東緣,切吉凹陷基底以三疊紀地層和印支期花崗巖為主,基底具東淺西深的斜坡帶性質(zhì)[20]。

研究區(qū)內(nèi)進行了大地電磁、重力、二維地震等多種地球物理勘查,以位于切吉凹陷中部1條NE向綜合測線為例介紹應用效果。該測線地震勘查結(jié)果(圖4)顯示覆蓋層西厚東薄,多條斷裂切穿基底,斷裂多為西傾逆斷層。共和組—咸水河組地層為砂巖和泥巖,最厚處達2 900 m,具有較好的保溫效果。

圖4 共和盆地典型地震勘查結(jié)果

通過地震劃分的地層作為模型結(jié)構(gòu)約束,聯(lián)合二維重力、大地電磁法兩種數(shù)據(jù)進行聯(lián)合反演,得到該條剖面的電阻率分布(圖5)和密度分布(圖6),由于該條測線磁法數(shù)據(jù)質(zhì)量較差,故磁法數(shù)據(jù)沒有加入該案例進行聯(lián)合反演。綜合分析認為共和組地層為次高阻、低密度層,具有一定的橫向連續(xù)性,局部電阻率橫向不均勻,推斷為因斷層發(fā)育導致的電阻率變化,橫向西部厚度大、東部厚度略小,與地震勘探顯示結(jié)構(gòu)一致,該層重力反演結(jié)果顯示為低密度響應;臨夏組和咸水河組地層為低阻、低密度響應,西部埋深大、東部埋深小,電阻率橫向連續(xù)性較好;下覆花崗巖為高阻、高密度特征。

圖5 研究區(qū)聯(lián)合反演電阻率結(jié)果

圖6 研究區(qū)聯(lián)合反演密度結(jié)果

從電阻率和密度交會(圖7)及密度分布圖中可見,密度由2個高斯分布組成,其均值分別為3.7 g/cm3和0.5 g/cm3,電阻率同樣由2個高斯分布組成,其均值分別為20 Ω·m和70 Ω·m。

圖7 密度和電阻率聯(lián)合反演結(jié)果交會(a)及密度分布(b)

將聯(lián)合反演后剖面的電阻率值利用式(13)計算得到剖面的溫度場分布結(jié)果,由于上覆地層巖性主要為泥巖和砂巖,為區(qū)域內(nèi)的保溫層,其電阻率和溫度對應關(guān)系不夠明確,因此本項工作僅計算了該剖面覆蓋層之下花崗巖地層的溫度場分布(圖8),其中活化能和前指因子σ0按照表1中實驗室利用花崗巖測定的數(shù)值來進行計算,其中E0為0.9 eV,logσ0為-2.4 Ω-1cm-1。由于該溫度場是通過電阻率剖面計算得來,其溫度場分布特征規(guī)律與電阻率剖面特征基本一致,花崗巖中3 000～4 000 m深度范圍內(nèi)溫度基本均勻,溫度范圍在200 ℃左右,等溫線呈西低東高的特征,隨著深度加深溫度逐漸升高,到5 000 m深剖面西段溫度達到280 ℃以上;但由于淺部地層電阻率較低,花崗巖頂部接近覆蓋層的區(qū)域所計算溫度場不符合地溫場特征,造成該種形態(tài)的原因是由于溫度場是通過電阻率直接計算得到,而該段處于電阻率從低向高變化的過渡帶區(qū)域,需要進行針對的校正研究才能符合真實地溫場特征。

圖8 溫度—深度剖面

7 結(jié)論

1)利用高斯混合模型能夠表征模型參數(shù)的多峰分布特點,將其用于聯(lián)合反演連續(xù)地球物理變量和離散巖性。相比傳統(tǒng)的反演算法,將巖石物性測量結(jié)果中同一巖性不同參數(shù)之間的統(tǒng)計學關(guān)系引入反演過程中,更加符合實際規(guī)律。

2)聯(lián)合反演相對于單一反演不論在數(shù)值恢復能力還是在異常邊界的刻畫能力上都有較大程度的提高,使得反演結(jié)果具有較好的空間連續(xù)性。

3)基于約束的聯(lián)合反演方案,雖然可以較好地將已知的結(jié)構(gòu)信息和巖石物性分布規(guī)律置于反演過程中,反演結(jié)果與先驗信息更加匹配,但所求得的反演最優(yōu)解與單獨反演相比,模型參數(shù)的小尺度變化被“平滑”,使得垂向分辨率變差。

4)直接將電阻率與溫度之間的實驗室測量參數(shù)應用在巖性簡單的大尺度溫度場,分布結(jié)果與鉆孔測溫結(jié)果基本相符,但在巖性過渡帶位置,溫度場與實際溫度測量結(jié)果相差較大,需要進行校正研究,進行針對性校正后可以用于溫度預測使用。