基于聚類和多元地質(zhì)統(tǒng)計(jì)學(xué)的電—震聯(lián)合建模約束反演技術(shù)及應(yīng)用

楊 博 張祥國 劉 展* 徐凱軍

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ②中國石油大慶油田公司勘探事業(yè)部，黑龍江大慶 163453)

0 引言

近年來，盡管單一的三維重、磁、電反演，尤其是三維大地電磁(MT)反演技術(shù)正漸趨成熟，但由于單一方法的局限性、觀測數(shù)據(jù)有限且難免存在觀測誤差等原因，使反演的多解性問題較嚴(yán)重，增加了地質(zhì)解釋的難度。不同地球物理方法對(duì)地下介質(zhì)不同物性參數(shù)的敏感度不同，綜合利用多種地球物理方法，可從不同角度對(duì)同一地質(zhì)體進(jìn)行研究。綜合多種地球物理方法進(jìn)行反演與解釋，減小多解性，已成為當(dāng)前地球物理反演研究的熱點(diǎn)問題[1-2]。地震資料通常對(duì)中、淺層具有較高分辨率，但在某些情況下，單一地震方法所得深層構(gòu)造成像并不清晰，重、磁、電勘探可作為其補(bǔ)充。其中大地電磁方法由于具有不受高阻層屏蔽、橫向分辨能力較強(qiáng)且探測深度大等優(yōu)點(diǎn)，在油氣勘探領(lǐng)域成為地震方法的重要輔助。同時(shí)，利用已知的地震、地質(zhì)、測井等資料進(jìn)行重、磁、電法數(shù)據(jù)約束反演，是提高重、磁、電勘探分辨率的重要途徑。

不同物性參數(shù)的耦合是聯(lián)合反演的關(guān)鍵，這些物性參數(shù)耦合約束主要分為巖石物理關(guān)系(經(jīng)驗(yàn)性的或統(tǒng)計(jì)性的)約束和結(jié)構(gòu)相似性約束兩大類。結(jié)構(gòu)相似耦合是從幾何關(guān)系上進(jìn)行約束，一般約束性相對(duì)較弱。Gallardo等[3]在進(jìn)行電—震聯(lián)合反演時(shí)利用交叉梯度算子識(shí)別結(jié)構(gòu)邊界的相似性，目前交叉梯度算子已成為應(yīng)用最廣泛的結(jié)構(gòu)相似度量方式之一，在多種聯(lián)合反演中得到應(yīng)用。

基于巖石物理關(guān)系的聯(lián)合反演的物性參數(shù)耦合性強(qiáng)，反演結(jié)果精度高。許多學(xué)者已針對(duì)構(gòu)建巖石物理關(guān)系這一問題展開研究。在傳統(tǒng)的基于經(jīng)驗(yàn)關(guān)系的巖石物理耦合方法中，先定義模型之間的關(guān)系函數(shù)，并將該函數(shù)應(yīng)用于整個(gè)區(qū)域，該類方法要求先驗(yàn)的巖石物理關(guān)系具有較高可靠性。Colombo等[4]建立了密度、電阻率與速度的經(jīng)驗(yàn)關(guān)系，并通過重、電、震數(shù)據(jù)聯(lián)合反演圈定了隱伏鹽丘； De Stefano等[5]利用Gardner公式[6]建立了地層密度與地震縱波速度的關(guān)系式，開展重力與地震旅行時(shí)數(shù)據(jù)聯(lián)合反演，并通過理論模型測試表明聯(lián)合反演能提高成像精度。

Zhdanov等[7]研究了一種廣義的耦合方式，即基于 Gramian 約束的多種地球物理數(shù)據(jù)聯(lián)合反演，該耦合方式不需先定義模型之間具體的關(guān)系函數(shù)。近年來，人們將模糊C均值聚類(FCM)思想引入聯(lián)合反演，通過聚類思想建立不同巖石物性之間的關(guān)系，該類方法可解決傳統(tǒng)的統(tǒng)計(jì)方法難以將多個(gè)巖石物理關(guān)系應(yīng)用于特定區(qū)域的問題。Sun等[8-9]提出基于引導(dǎo)FCM技術(shù)的聯(lián)合反演，并將其與巖性分類相結(jié)合，提高了反演過程的穩(wěn)定性及反演結(jié)果的可靠性； Carter-Mcauslan等[10]在重—震聯(lián)合反演目標(biāo)函數(shù)中引入FCM聚類，利用聚類約束實(shí)現(xiàn)速度與密度模型的巖性耦合。

陳曉等[11]提出基于寬范圍巖石物性約束的大地電磁與地震聯(lián)合反演，該類方法可降低由于先驗(yàn)巖石物理關(guān)系不夠精確帶來的影響。胡祖志等[12]利用已知的測井、地震剖面等先驗(yàn)信息進(jìn)行約束建模，并基于井—震約束的MT和重力數(shù)據(jù)實(shí)現(xiàn)了人工魚群聯(lián)合反演。相鵬等[13]提出一種變密度—速度關(guān)系的重力與地震同步聯(lián)合反演方法。

基于巖石物理關(guān)系的聯(lián)合反演的前提是可靠的巖石物理關(guān)系，而可靠的巖石物理關(guān)系一般通過實(shí)測巖心、測井曲線等資料獲得，實(shí)際應(yīng)用中并不一定能獲得足夠的上述資料。再者，目前建立的巖石物理關(guān)系是一種點(diǎn)對(duì)點(diǎn)的映射，這與實(shí)際的空間對(duì)空間的巖石物性關(guān)系是不相符的。

地質(zhì)統(tǒng)計(jì)學(xué)中的變差函數(shù)能定量表征地質(zhì)模型空間特征，反映物性值縱、橫向變化規(guī)律。若尚未知曉地下地質(zhì)體精確的物性統(tǒng)計(jì)特征，則可間接從地震信息、試井解釋、地質(zhì)推理與解釋等獲取變差函數(shù)。

本文在Yang等[14]研究的基礎(chǔ)上，提出一種基于聚類和多元地質(zhì)統(tǒng)計(jì)學(xué)的電—震聯(lián)合建模約束反演方法： 首先根據(jù)地震反演得到的速度模型、無約束MT反演得到的電阻率模型和已知鉆井信息建立初始電阻率模型； 再利用多元地質(zhì)統(tǒng)計(jì)學(xué)的交叉—變差函數(shù)建立速度與電阻率之間的巖石物理關(guān)系； 同時(shí)利用機(jī)器學(xué)習(xí)中的引導(dǎo)FCM算法進(jìn)行基于巖石物理關(guān)系的多重約束反演，實(shí)現(xiàn)電—震聯(lián)合建模； 將描述地質(zhì)模型特征表征為定量化數(shù)據(jù)，并與地球物理數(shù)據(jù)相結(jié)合，得到既能擬合地球物理數(shù)據(jù)，又符合地質(zhì)模型特征的最優(yōu)解。 將本文方法應(yīng)用于大楊樹盆地南部坳陷實(shí)際數(shù)據(jù)，取得了令人滿意的處理效果及最終綜合解釋成果。

1 電—震聯(lián)合建模約束反演方法

1.1 技術(shù)流程

各種反演方法都存在多解性，對(duì)反演結(jié)果直接進(jìn)行地質(zhì)解釋會(huì)產(chǎn)生較大誤差。解決該問題的常用方法是綜合利用各種已知信息進(jìn)行反演和解釋，以提高可靠性。本文提出采取基于機(jī)器學(xué)習(xí)和多元地質(zhì)統(tǒng)計(jì)學(xué)的聯(lián)合建模約束反演方法，分四個(gè)部分：

(1)對(duì)由疊前深度偏移(PSDM)數(shù)據(jù)反演得到速度模型降維，得到大地電磁尺度的速度模型，并對(duì)其進(jìn)行聚類分析；

(2)利用根據(jù)已知電阻率測井?dāng)?shù)據(jù)和無約束MT反演得到的電阻率模型(無測井?dāng)?shù)據(jù)時(shí))統(tǒng)計(jì)的平均電阻率值對(duì)聚類后的速度模型賦值，得到新的MT反演初始模型，更新模型的協(xié)方差模型，引入地震結(jié)構(gòu)約束[15]；

(3)以測井?dāng)?shù)據(jù)和MT反演結(jié)果建立反映空間統(tǒng)計(jì)巖石物理關(guān)系的交叉—變差函數(shù)(垂向交叉—變差函數(shù)由測井?dāng)?shù)據(jù)求得，由反演結(jié)果求得垂向與主方向和次方向交叉—變差函數(shù)的比例關(guān)系，進(jìn)而確定先驗(yàn)交叉—變差函數(shù)[16])；

(4)將交叉—變差函數(shù)整合到傳統(tǒng)的目標(biāo)函數(shù)中，并利用模糊約束反演(FCI)[17]進(jìn)行基于巖石物理關(guān)系的多重約束反演[14]。

上述具體技術(shù)流程如圖1，易見基于巖石物理關(guān)系的多重約束MT反演是其中關(guān)鍵。

1.2 基于巖石物理關(guān)系的多重約束MT反演

Kelbert等[18]提出了基于協(xié)方差矩陣的MT反演目標(biāo)函數(shù)

(1)

式中：m是模型參數(shù)；d是實(shí)測MT數(shù)據(jù)；Cd是數(shù)據(jù)誤差的協(xié)方差；f(m)是正演的MT數(shù)據(jù)；m0是初始模型；β是拉格朗日算子(正則化參數(shù))；Cm是模型協(xié)方差，也稱為平滑算子，表征反演期間模型元素之間的平滑度。從理論上講，Cm的每個(gè)元素都可針對(duì)m中的所有成對(duì)元素進(jìn)行獨(dú)立配置，從而可對(duì)相鄰單元之間的平滑度進(jìn)行相當(dāng)大的微調(diào)。為了降低反演的多解性，擬采用綜合不同地球物理數(shù)據(jù)的聯(lián)合反演方法。對(duì)于基于巖石物理關(guān)系的聯(lián)合反演而言，構(gòu)建合理的巖石物理關(guān)系是關(guān)鍵。本文將多元地質(zhì)統(tǒng)計(jì)學(xué)中的交叉—變差函數(shù)整合到MT反演的目標(biāo)函數(shù)中，以表征速度與電阻率之間關(guān)系。

圖1 本文電—震聯(lián)合建模約束反演技術(shù)流程

1.2.1 交叉—變差函數(shù)求取技術(shù)

上已述及，地質(zhì)統(tǒng)計(jì)學(xué)中變差函數(shù)可定量表征地質(zhì)模型空間特征，反映物性值縱、橫向變化規(guī)律。 同時(shí)，變差函數(shù)也是地質(zhì)統(tǒng)計(jì)學(xué)中最基本和最重要的模擬工具，用于描述數(shù)據(jù)值的空間互相關(guān)性，即數(shù)據(jù)點(diǎn)在空間上相距越遠(yuǎn)，相關(guān)性越小。其數(shù)值可表示為區(qū)域化變量Z(x)在x和x+h兩點(diǎn)之差的方差的一半，此處h表示滯后距，或稱距離、間隔、步長。變差函數(shù)與數(shù)據(jù)點(diǎn)位置無關(guān)，只依賴于h。

此處電阻率ρ被h分割為相隔h的點(diǎn)xi與xi+h上的N(h)對(duì)電阻率{ρ(xi)，ρ(xi+h)}(i=1，2，…，N(h))，可根據(jù)下式計(jì)算其實(shí)驗(yàn)變差函數(shù)

(2)

式中：i為樣本序號(hào)；N(h)為數(shù)據(jù)對(duì){ρ(xi)，ρ(xi+h)}的個(gè)數(shù)，記為點(diǎn)對(duì)數(shù)。

以h為橫坐標(biāo)、γ(h)為縱坐標(biāo)作出圖2所示的變差圖。

γ(h)可用四個(gè)參數(shù)來描述[19]，即變差函數(shù)類型(如高斯模型、球狀模型、指數(shù)模型)、變程a(表示超過該距離后，數(shù)據(jù)點(diǎn)之間就不再有明顯相關(guān)性，也稱作影響距離)、塊金值C0(表示相距很近兩點(diǎn)的樣品變化)和基臺(tái)值C0+C(反映變量的變化幅度)。

每對(duì)區(qū)域化變量Z1(x)和Z2(x)間的相關(guān)范圍與方向差異問題定義為交叉—變差函數(shù)。針對(duì)實(shí)際情形，同一地質(zhì)體可通過不同物性參數(shù)表征，如速度v和電阻率ρ，其交叉—變差函數(shù)可表示

圖2 變差函數(shù)示意圖

[v(xi)-v(xi+h)]

(3)

1.2.2 基于巖石物理關(guān)系的約束技術(shù)

將多元地質(zhì)統(tǒng)計(jì)學(xué)中的交叉—變差函數(shù)整合到式(1)，以表征速度與電阻率之間的關(guān)系。交叉—變差函數(shù)項(xiàng)為

(4)

(5)

基于巖石物理關(guān)系的約束技術(shù)使得MT反演結(jié)果不僅可擬合觀測MT數(shù)據(jù)，同時(shí)滿足先驗(yàn)的電阻率—速度關(guān)系。

1.2.3 模糊約束技術(shù)

許多學(xué)者已將FCM聚類技術(shù)集成到確定性反演框架中，以進(jìn)一步約束地球物理反演[8-10，17，20-21]。

為使反演得到的物性值聚類為已知的巖石物性值，Sun等[8-9]提出引導(dǎo)FCM技術(shù)

(6)

式中：M是待聚類數(shù)據(jù)的數(shù)量；c是聚類中心的數(shù)量；mk是第k個(gè)數(shù)據(jù)；μki是隸屬度值，用于度量第k個(gè)數(shù)據(jù)屬于第i個(gè)聚類的程度；q是模糊化參數(shù)；ui是通過FCM算法自動(dòng)更新的第i個(gè)聚類中心；vi是根據(jù)巖石樣本測量確定的先驗(yàn)物性值。

該策略的一個(gè)優(yōu)勢是不會(huì)損害FCM算法的良好收斂性，同時(shí)還能根據(jù)先驗(yàn)的巖石物理信息將搜索聚類中心引導(dǎo)到所需位置。

對(duì)式(6)分別對(duì)聚類中心ui和隸屬度μki求微分，并使其等于零，得到更新后的對(duì)應(yīng)值

(7)

(8)

本文所使用FCI技術(shù)的基本思想是在反演的每一次迭代過程中利用引導(dǎo)FCM技術(shù)引入已知的電阻率信息來約束模型參數(shù)。

FCI技術(shù)的額外輸入是聚類中心的數(shù)量和聚類中心的值(比如每種巖性的電阻率平均值)。

(9)

通過這種方式實(shí)施反演，一起迭代、優(yōu)化了三個(gè)參數(shù)，即模型m、聚類中心u和隸屬度μ。

基于引導(dǎo)FCM分析技術(shù)的FCI技術(shù)使反演過程中實(shí)現(xiàn)反演與地質(zhì)分類的互相改進(jìn)，得到的電阻率模型不僅可擬合觀測到的MT數(shù)據(jù)，同時(shí)更接近真實(shí)電阻率值，且更利于地質(zhì)解釋。

2 應(yīng)用效果分析

為了驗(yàn)證本文方法的實(shí)際應(yīng)用效果，選取大楊樹盆地南部坳陷實(shí)測數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。

大楊樹盆地主要目的層由火山巖、火山碎屑巖和沉積巖組成。由于該區(qū)火山巖覆蓋程度較高，使之對(duì)于地震波具有低通效應(yīng)且對(duì)反射能量屏蔽、吸收衰減嚴(yán)重，導(dǎo)致九峰山組地震資料品質(zhì)差，無法開展有效的沉積巖預(yù)測和構(gòu)造解釋[22-24]。通常從MT數(shù)據(jù)可獲取深層結(jié)構(gòu)，但其分辨率不及地震數(shù)據(jù)。 因此，為了獲得更可靠的地下結(jié)構(gòu)，實(shí)施了電—震聯(lián)合建模約束反演。

研究區(qū)是一北東向傾斜的面積為13km×10km的長方形區(qū)域(圖3)，測網(wǎng)密度為0.50km×0.25km，布設(shè)26條測線，共計(jì)858個(gè)測點(diǎn)，每個(gè)測點(diǎn)MT數(shù)據(jù)頻點(diǎn)數(shù)不少于38。MT40線與地震S3線一致并穿過鉆孔，故此次選取40線為例進(jìn)行分析。

圖3 研究區(qū)MT測點(diǎn)、地震測線及鉆井位置 黑點(diǎn)/黑字分別表示MT測點(diǎn)位置/線號(hào)，粉字表示測點(diǎn)號(hào)； 藍(lán) 線/藍(lán)字分別表示地震測線/線號(hào)； 紅點(diǎn)/紅字分別表示井位/井名

2.1 無地震約束2D MT反演

首先對(duì)40線做1D MT blocky反演，然后針對(duì)其中MT13測點(diǎn)的反演結(jié)果(圖4)與實(shí)測視電阻率及阻抗相位數(shù)據(jù)進(jìn)行正演擬合(圖5)，效果很好。

以1D反演結(jié)果作為初始模型，進(jìn)行了無地震數(shù)據(jù)約束的2D MT反演，獲得的40線電阻率剖面如圖6a所示； 根據(jù)MT反演的電阻率數(shù)值的大小與高低變化進(jìn)行地層解釋，并依托S3地震數(shù)據(jù)進(jìn)行斷裂識(shí)別和描述(圖6b)。可見該方法反演的視電阻率(圖7)及阻抗相位(圖8)數(shù)據(jù)的正演響應(yīng)很好地?cái)M合了對(duì)應(yīng)的實(shí)測數(shù)據(jù)。

由圖6可知，甘河組(K1g)電阻率值相對(duì)較高，這主要是由于該層大面積發(fā)育玄武巖等火山巖；九峰山組(K1j)是低阻，因?yàn)樵搶又饕l(fā)育砂泥巖；龍江組(K1l)為電阻率較高，主要發(fā)育碎屑巖和酸性熔巖；基底為最高阻。反演的電阻率模型與已知的信息有較好一致性。通過鉆井標(biāo)定可看出，該電性異常剖面對(duì)目標(biāo)層九峰山組頂、底界面有所反映，但對(duì)九峰山組內(nèi)細(xì)節(jié)刻畫欠清晰。

圖4 MT13測點(diǎn)1D MT blocky反演結(jié)果

圖5 MT13測點(diǎn)1D MT blocky反演結(jié)果的正演擬合圖 (a)視電阻率； (b)阻抗相位

圖6 MT40線無約束MT反演電阻率 剖面(a)及其解釋剖面(b)

圖7 MT40線實(shí)測與無約束MT反演擬合視電阻率對(duì)比 測點(diǎn)位置如圖3

圖8 MT40線實(shí)測與無約束MT反演擬合阻抗相位對(duì)比 測點(diǎn)位置如圖3

2.2 電—震聯(lián)合建模約束反演

為了更精細(xì)地描述目的層九峰山組地層細(xì)節(jié)，在無地震約束MT反演結(jié)果的基礎(chǔ)上進(jìn)行電—震聯(lián)合建模約束反演。在應(yīng)用聚類方法之前，應(yīng)預(yù)知聚類中心的數(shù)量和聚類中心值。 經(jīng)綜合整理分析工區(qū)測井?dāng)?shù)據(jù)、巖心露頭及1D MT blocky反演結(jié)果，得到如表1所示的每個(gè)地層的電阻率及速度。據(jù)此確定聚類中心數(shù)為7，速度、電阻率的聚類中心值分別為[4215，4325，4405，4645，4830，5000，6000](m/s)、[25，50，62，93，182，345，1260](Ω·m)，以此作為聚類先驗(yàn)信息，優(yōu)化反演結(jié)果。

表1 各地層的速度及電阻率統(tǒng)計(jì)平均值

2.2.1 構(gòu)建先驗(yàn)電阻率模型

先驗(yàn)電阻率模型的構(gòu)建是MT反演的重要部分。本文從地震數(shù)據(jù)體提取與地電結(jié)構(gòu)相關(guān)的關(guān)鍵地震層構(gòu)建先驗(yàn)?zāi)Ｐ偷牡仉娍蚣埽鳛镸T反演的輸入。 以40線為例，構(gòu)建先驗(yàn)電阻率模型的步驟如下：

(1)利用測井信息進(jìn)行標(biāo)定，綜合PSDM數(shù)據(jù)與無約束MT反演得到的電阻率模型(圖9a)進(jìn)行綜合解釋，得到關(guān)鍵層位解釋結(jié)果(圖9b)；

(2)基于層位解釋結(jié)果進(jìn)行地質(zhì)統(tǒng)計(jì)學(xué)反演，得到速度模型(圖9c)；

(3)對(duì)速度模型進(jìn)行聚類分析(圖9d)；

(4)利用無約束MT反演得到的電阻率模型的平均電阻率值對(duì)聚類后的速度模型賦值，得到新一輪MT反演初始模型(圖9e)。

所得的先驗(yàn)電阻率模型不僅反映了相應(yīng)的地震構(gòu)造形態(tài)，同時(shí)保持合理的電阻率分布。

2.2.2 基于巖石物理關(guān)系的多重約束MT反演

利用測井?dāng)?shù)據(jù)與40線速度和電阻率剖面建立交叉—變差函數(shù)，進(jìn)行基于巖石物理關(guān)系多重約束反演，得到電阻率數(shù)據(jù)(圖10a)，然后對(duì)其進(jìn)行定量解釋(圖10b)。圖11和圖12分別表明基于巖石物理關(guān)系的多重約束MT反演的視電阻率及阻抗相位數(shù)據(jù)的正演響應(yīng)很好地?cái)M合了相應(yīng)實(shí)測數(shù)據(jù)。

圖9 基于聚類分析的先驗(yàn)電阻率模型構(gòu)建流程示意 (a)無約束MT反演的電阻率模型； (b)PSDM圖像及 綜合解釋得到關(guān)鍵層位； (c)MT尺度的速度模型； (d)聚類后速度模型； (e)生成的先驗(yàn)電阻率模型

本文反演策略與無約束反演結(jié)果(圖6a)的對(duì)比表明： ①本文反演方法所得電阻率剖面成層性更好，且與主要地震構(gòu)造形態(tài)一致，特別是對(duì)目的層九峰山組的劃分更精細(xì)、清晰； ②先驗(yàn)電阻率數(shù)據(jù)表明甘河組(普遍呈電阻率較高值)廣泛分布火山巖，多重約束反演結(jié)果與該物性特征一致； ③龍江組下部電阻率值應(yīng)大于其上部電阻率值，與實(shí)際情況相吻合。

圖13顯示電—震聯(lián)合建模約束反演所得交叉—變差函數(shù)與先驗(yàn)交叉—變差函數(shù)相一致(圖13a)。

圖10 本文反演策略所得40線電阻率 剖面(a)及其解釋剖面(b)

圖11 MT40線實(shí)測與約束MT反演擬合視電阻率對(duì)比 所選測點(diǎn)位置如圖3所示

圖12 MT40線實(shí)測與約束MT反演擬合阻抗相位對(duì)比圖 測點(diǎn)位置如圖3

圖13b顯示了電—震聯(lián)合建模約束反演結(jié)果與電阻率測井?dāng)?shù)據(jù)趨勢一致。證實(shí)了反演結(jié)果的可靠性。

圖13 約束反演結(jié)果與先驗(yàn)信息對(duì)比 (a)交叉—變差函數(shù)； (b)W1測井曲線

3 結(jié)論

單一地球物理勘探方法不可避免地具有片面性和局限性。為了提高對(duì)目標(biāo)地質(zhì)體的識(shí)別能力，利用聚類和多元地質(zhì)統(tǒng)計(jì)學(xué)方法進(jìn)行MT與地震數(shù)據(jù)聯(lián)合建模反演，實(shí)現(xiàn)從無約束MT反演方法到基于巖石物理關(guān)系多重約束反演的流程。由于通過構(gòu)建先驗(yàn)電阻率模型而引入地震結(jié)構(gòu)約束、以交叉—變差函數(shù)引入先驗(yàn)電阻率與速度關(guān)系約束、用FCI技術(shù)引入先驗(yàn)電阻率信息并結(jié)合地質(zhì)分類，所以多重約束反演得到的電阻率模型更符合先驗(yàn)信息且更利于地質(zhì)解釋。實(shí)際數(shù)據(jù)的應(yīng)用結(jié)果證實(shí)了將多種地球物理方法整合到地球物理勘探中的必要性。