相 鵬 王金鐸 譚紹泉 陳學國

(中國石化勝利油田分公司勘探開發(fā)研究院，山東東營 257000)

0 引言

隨著油氣勘探領(lǐng)域轉(zhuǎn)向更深、更復雜的目標，以地震為代表的單一地球物理勘探技術(shù)遇到了極大的挑戰(zhàn)，低信噪比、頻帶限制、復雜地表、復雜構(gòu)造、埋深大等問題影響了地震成像的分辨率和可靠性，構(gòu)造解釋與建模多解性強，制約了油氣勘探進程。為了更準確地獲取復雜目標的信息，不同地球物理資料的聯(lián)合反演方法已成為近年來研究的熱點[1-4]。

聯(lián)合反演從聯(lián)合方式上可以分為兩大類：一類是相對較早的順序聯(lián)合反演[5-9]，它將一種地球物理資料(如重力)反演的模型(密度模型)通過物性關(guān)系轉(zhuǎn)換成另一種反演方法(如地震反演)的初始模型(速度模型)，然后進行反演，如此迭代直至獲得滿意的最終模型，雖然在很多領(lǐng)域已經(jīng)取得了較好的應用效果[10-13]，但其本質(zhì)仍是單一方法反演。另一類是同步聯(lián)合反演，因為其理論的完備性，受到越來越多學者的重視。多種數(shù)據(jù)在反演過程中相互補充、相互制約，可進一步縮小穩(wěn)定可靠解的范圍，達到同時擬合多種觀測數(shù)據(jù)，降低反演多解性。

同步聯(lián)合反演可以分為基于物性關(guān)系和基于幾何形狀相似性的同步聯(lián)合反演。基于物性關(guān)系的同步聯(lián)合反演[14-17]是利用Gardener公式[18]和Archie公式[19]或者二者的近似多項式分別建立速度與密度、速度與電阻率之間的關(guān)系。基于幾何形狀相似性的同步聯(lián)合反演則是假設(shè)地下目標的物性邊界相同，不同物性模型(速度、密度、電阻率、磁化率等)具有相似的幾何形狀，其中Gallardo等[20-21]提出的交叉梯度約束的電磁法與地震同步聯(lián)合反演，已經(jīng)成為目前同步聯(lián)合反演方法的主流。上述兩類同步反演方法各有優(yōu)缺點，基于幾何形狀相似性約束的同步聯(lián)合反演方法避免了建立物性關(guān)系的復雜過程，適合物性資料不足的研究區(qū)域，但是其幾何相似性的假設(shè)過于生硬，尤其是對兩種以上地球物理資料進行聯(lián)合反演時，如MT、重力和地震，反演得到的電阻率、密度和速度模型具有相似的幾何形狀，不一定符合研究區(qū)真實情況?；谖镄躁P(guān)系的同步聯(lián)合反演方法具有較為嚴謹?shù)膸r石物理基礎(chǔ)，若研究區(qū)的物性資料豐富，能夠充分利用已知物性資料獲得符合已有認識的反演結(jié)果。然而目前大部分反演方法中物性關(guān)系是固定不變的，若物性關(guān)系不準確則會導致反演結(jié)果出現(xiàn)較大偏差，因而限制了此類反演方法的適用范圍。

針對現(xiàn)有同步聯(lián)合反演方法的不足，本文提出一種變密度—速度關(guān)系的重力與地震同步聯(lián)合反演方法。初始物性關(guān)系在反演迭代過程中根據(jù)當前及上輪迭代的密度和速度模型得到更新，可降低不準確的初始物性關(guān)系對反演結(jié)果的不利影響，從而提高了重震同步聯(lián)合反演在低勘探程度區(qū)的應用潛力。

1 方法原理

假設(shè)速度與密度之間存在未知的線性關(guān)系，這個線性關(guān)系在模型的不同區(qū)域可以是不同的，并且在反演過程中可以根據(jù)每次迭代的結(jié)果自動修正。這種方法對存在由于物性、埋深、溫度等原因引起密度—速度關(guān)系變化的研究區(qū)具有實際意義。本節(jié)給出目標函數(shù)并推導出反問題的解估計公式，然后給出物性關(guān)系更新算法。

1.1 同步聯(lián)合反演目標函數(shù)及解估計

重力與地震同步聯(lián)合反演的目標函數(shù)定義為

(1)

式中：m表示模型參數(shù),包括速度參數(shù)向量V和密度參數(shù)向量ρ；gobs是實測重力異常；gcal是正演重力異常；σg是重力數(shù)據(jù)方差；Sobs是實測地震數(shù)據(jù)；Scal是正演地震數(shù)據(jù)；σS是地震數(shù)據(jù)方差；B是密度—速度關(guān)系矩陣；w1、w2、w3是權(quán)重。分別對ρ、V求偏導數(shù)，得到

(2)

(3)

式中：Ag和AS分別是重力和地震的靈敏度矩陣。

令式(2)和式(3)等于0，可得

(4)

(5)

聯(lián)立式(4)和式(5)，可得

ρ=P-QR-1Q)-1×

(6)

(7)

式中

(8)

Q=w3BT

(9)

(10)

求解上述方程組可采用非線性共軛梯度法[22]，在保證反演精度的前提下，提高計算速度，降低內(nèi)存需求。

1.2 物性關(guān)系更新算法

Delyon等[23]提出的EM隨機近似算法，簡稱SAEM(Stochastic Approximation of EM)算法，主要包括三步：

(1)求不完全數(shù)據(jù)樣本x的第k次迭代樣本值x(k)；

(2)根據(jù)下式更新充分統(tǒng)計量T(k)的近似值

T(k)=T(k-1)+γ(k)[T(x(k),y)-T(k-1)]

(11)

式中：y是觀測數(shù)據(jù)； {γ(k)}k≥0是正值單調(diào)遞減步長序列，取γ(k)=1/k；

(3)最大化目標函數(shù)，求取統(tǒng)計參數(shù)集θ(k+1)。

已經(jīng)證明該算法在絕大多數(shù)情況下是收斂的，參見文獻[23]。

本文按照SAEM算法的步驟(2)計算每次迭代密度模型和速度模型的充分統(tǒng)計量，進而更新密度—速度關(guān)系。在重震聯(lián)合反演中，密度模型ρ和速度模型V對應SAEM算法中的樣本x，物性關(guān)系B為對角陣，若將所有地層分隔成J個具有不同速度—密度線性關(guān)系的區(qū)域，那么B的對角線上有J個不同的值，則SAEM算法的統(tǒng)計參數(shù)集θ=B=[b(1),b(2),…,b(J)]，其中b(j)(j=1,2,…,J)表示第j個區(qū)域內(nèi)的物性關(guān)系。物性關(guān)系更新算法步驟如下：

1)根據(jù)式(6)～式(10)求得當前迭代ρ(k)和V(k);

(12)

(13)

(3)計算θ=B=[b(1),b(2),…,b(J)]

(14)

2 理論模型與實際資料試算

采用SEG/EAGE常用的Overthrust模型對本文方法進行測試。速度網(wǎng)格橫向剖分數(shù)目為801，垂向剖分數(shù)目為187，縱、橫向間距均為25m。水平坐標范圍0～20000m。圖1a為Overthrust速度模型，速度范圍為2360～6000m/s。將速度模型通過關(guān)系式ρ=v/2000轉(zhuǎn)換成密度模型，進行重力正演。

設(shè)地震炮點埋深為25m，共放199炮，炮間距為100m，第一炮位置在100m處。利用主頻為10Hz的零相位雷克子波激發(fā)，地表全孔徑接收，道間距為100m，第一道位置在50m處，共200道。時間采樣率為4ms，采樣點數(shù)為8192。圖1b是炮點位于1000m處的單炮記錄，是在頻率域利用有限差分算法正演獲得的。重力觀測點位于地表，第一個觀測點位于50m處，共100個觀測點，間距為200m，正演重力曲線見圖1c?？梢钥闯觯亓Ξ惓δＰ椭胁康寞B瓦狀沖斷構(gòu)造以及中深部的背斜(圖1c中方框區(qū)域)均有明顯的響應。

2.1 重/震單獨反演

從正演頻率域地震數(shù)據(jù)中選取3.5～20.6Hz的7個頻率用于全波形反演，保留真實模型淺層100m深度范圍內(nèi)的速度值，深層(≥100m)經(jīng)平滑處理后作為初始模型(圖2a)。設(shè)置最大迭代次數(shù)為80， 3.5、 9.6、 20.6Hz三個頻點的反演結(jié)果分別見圖2b～圖2d。對比圖2a,可見反演結(jié)果邊界清晰，取得了很好的效果。

圖1 Overthrust模型(a)速度模型; (b)1000m處正演單炮記錄; (c)重力異常曲線

全波形反演(FWI)方法具有局部收斂性，嚴重依賴于初始速度。若初始模型較完整地保留了低頻信息時，F(xiàn)WI經(jīng)初步迭代后應當就能夠收斂在真實解附近。圖3為初始模型為常速度3500m/s時的反演結(jié)果?？梢姡诔跏寄Ｐ筒焕硐氲那闆r下，地震數(shù)據(jù)單獨反演的結(jié)果誤差很大，出現(xiàn)了更多的高頻擾動，沒有反映出從淺部到深部速度逐漸增大的趨勢。

采用重加權(quán)共軛梯度法[24]對重力異常數(shù)據(jù)單獨進行反演，得到密度模型。初始模型的密度值設(shè)為常數(shù)2.0g/cm3。在不加任何約束的情況下，反演結(jié)果如圖4所示?？梢钥闯?，分辨率較低，不能刻畫出疊瓦狀沖斷構(gòu)造的基本形態(tài)，但其反映的趨勢背景是正確的，基本上能反映出了密度從上到下逐漸增大的趨勢，模型中部的背斜輪廓亦可大致看出，左右兩邊的水平沉積地層反映得較為清晰。

2.2 重震同步聯(lián)合反演

采用本文方法對該模型進行重力和地震數(shù)據(jù)聯(lián)合反演。

速度初始模型設(shè)為常值3000m/s， 密度模型設(shè) 為常值2.0g/cm3， 則初始速度—密度關(guān)系為v=1500ρ。

首先開展固定物性關(guān)系的同步聯(lián)合反演試驗，以分析在物性關(guān)系變化的情況下，采用固定速度—密度關(guān)系反演的結(jié)果。展示3個頻點的反演結(jié)果，見圖5。可見反演效果較圖3地震單獨反演結(jié)果有所提高，反映出從淺部到深部速度逐漸增大的趨勢，但是高頻部分失真嚴重，難以準確解釋逆斷層。說明不準確的密度—速度關(guān)系會降低聯(lián)合反演的精度和可靠性。

圖2 理想速度初始模型的地震全波形反演結(jié)果(速度模型)(a)初始速度模型; (b)3.5Hz反演結(jié)果; (c)9.6Hz反演結(jié)果; (d)20.6Hz反演結(jié)果

圖3 常速度(3500m/s)初始模型的不同頻率地震全波形反演結(jié)果(速度模型)(a)3.5Hz； (b)9.6Hz； (c)20.6Hz

圖4 常密度(2.0g/cm3)初始模型重力反演剖面

用本文算法開展聯(lián)合反演試驗，檢驗在反演過程中修正物性關(guān)系對結(jié)果產(chǎn)生的影響。圖6為采用變物性關(guān)系進行同步聯(lián)合反演結(jié)果?？梢姕\部地層除兩側(cè)邊界受到邊界效應的影響，其他部分都得到了較好的恢復；深部逆掩構(gòu)造不很清晰，但是基本能夠分辨出構(gòu)造形態(tài)，整體效果優(yōu)于固定物性關(guān)系的聯(lián)合反演結(jié)果(圖5)。根據(jù)反演剖面，速度—密度的關(guān)系為Binv=1985.7，較準確地描述了模型準確的物性關(guān)系B=2000。

取得理想效果的主要原因是，反演過程中通過重力信息建立了低頻背景模型，從而將地震信息的反射分量與透射分量分離。反射分量可集中恢復速度模型中的高頻成分。而且在反演過程中不斷根據(jù)每次迭代模型的更新量調(diào)整不準確的初始密度—速度關(guān)系，使重、震之間的相互約束關(guān)系更合理。

圖5 基于固定密度—速度關(guān)系的不同頻率下重震同步聯(lián)合反演結(jié)果(速度模型)(a)3.5Hz； (b)9.6Hz； (c)20.6Hz

圖6 基于變密度—速度關(guān)系的不同頻率重震同步聯(lián)合反演結(jié)果(速度模型)(a)3.5Hz；(b)9.6Hz；(c)20.6Hz

這個實驗證明，重震同步聯(lián)合反演在初始模型和初始密度—速度關(guān)系都不理想的情況下，仍能取得較滿意的反演效果。

2.3 實際資料試算

為驗證本文方法的實際應用效果，選取中國西部M區(qū)數(shù)據(jù)進行實驗。圖7所示為該區(qū)1∶50000區(qū)域重力異常圖，其中有四條重力—地震聯(lián)合攻關(guān)線(圖中測線1～測線4)。由圖可見，研究區(qū)重力異常整體表現(xiàn)出北東向重力高值異常帶，但沿測線2異常發(fā)生扭曲、錯動，測線兩側(cè)異常形態(tài)、規(guī)模存在較大差異，以測線2為界，東西兩側(cè)地質(zhì)構(gòu)造特征差異較大。

圖7 M區(qū)區(qū)域重力異常圖

選擇測線2開展聯(lián)合反演實驗。圖8為該線重力—地震聯(lián)合攻關(guān)所得重力異常，測線兩端對應沉積凹陷，在剖面中間構(gòu)造帶部位是較單一的高異常，形態(tài)、幅值均低于測線1的對應位置。

首先按照常規(guī)速度分析流程生成初始層速度模型(圖9a)?？梢钥闯鲈撈拭鎻臏\到深、從南向北速度逐漸增大，剖面中部為高速區(qū)，整體上速度場畸變比較嚴重，深層構(gòu)造刻畫不清。從偏移剖面(圖9b)上看，沉積地層反映清晰，但其中深層成像模糊，構(gòu)造解釋存在一定困難。

應用本文變密度—速度關(guān)系的重磁聯(lián)合反演得到的速度剖面見圖10a?？梢钥闯觯撍俣绕拭鎸Ω咚俚貙拥男螒B(tài)刻畫清晰，較好地展現(xiàn)了深層逆掩構(gòu)造。圖10b是基于圖10a的疊前深度偏移剖面?？梢娖拭娌ńM特征清楚，中深層信噪比和連續(xù)性也得到增強，斷點、斷面清楚，能夠較準確、清楚地反映構(gòu)造面貌。

圖8 測線2重力異常剖面

圖9 初始速度模型(a)和深度偏移剖面(b)

圖10 基于變密度—速度關(guān)系的聯(lián)合反演速度剖面(a)及深度偏移剖面(b)

3 結(jié)論

本文提出了一種變密度—速度關(guān)系的重力與地震同步聯(lián)合反演方法，給出包含密度—速度關(guān)系項的重震同步聯(lián)合反演目標函數(shù)，詳細推導了求解密度模型與速度模型的解估計公式。該方法利用重力資料補充地震資料中缺少的低頻信息，建立低頻背景模型，使地震高頻信息能夠更準確地恢復模型中的高頻成分。

方法的創(chuàng)新之處在于在反演迭代過程中不斷更新密度—速度關(guān)系，針對不確定的物性關(guān)系對反演產(chǎn)生不利影響的問題，提出了一種基于密度模型和速度模型的高階統(tǒng)計量的物性關(guān)系更新算法。模型試算證明，該方法可以有效提高初始模型和初始密度—速度關(guān)系都不確定情況下聯(lián)合反演的效果，較固定密度—速度關(guān)系的同步聯(lián)合反演方法，反演效果有顯著提高。實際資料試算檢驗了方法的實際應用效果，證明方法具有較好的實用性。

本文僅對重震同步聯(lián)合反演的方法原理做了詳細的闡述，但文中所述的反演思路與公式可經(jīng)簡單修改后推廣應用于電震、重電震同步聯(lián)合反演，有待后續(xù)深入研究。