相 鵬 王金鐸 譚紹泉 陳學國

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

0 引言

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

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

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

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

1 方法原理

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

1.1 同步聯合反演目標函數及解估計

重力與地震同步聯合反演的目標函數定義為

(1)

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

(2)

(3)

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

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

(4)

(5)

聯立式(4)和式(5)，可得

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

(6)

(7)

式中

(8)

Q=w3BT

(9)

(10)

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

1.2 物性關系更新算法

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

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

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

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

(11)

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

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

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

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

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

(12)

(13)

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

(14)

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

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

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

2.1 重/震單獨反演

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

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

全波形反演(FWI)方法具有局部收斂性，嚴重依賴于初始速度。若初始模型較完整地保留了低頻信息時，FWI經初步迭代后應當就能夠收斂在真實解附近。圖3為初始模型為常速度3500m/s時的反演結果?？梢?，在初始模型不理想的情況下，地震數據單獨反演的結果誤差很大，出現了更多的高頻擾動，沒有反映出從淺部到深部速度逐漸增大的趨勢。

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

2.2 重震同步聯合反演

采用本文方法對該模型進行重力和地震數據聯合反演。

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

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

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

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

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

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

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

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

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

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

2.3 實際資料試算

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

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

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

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

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

圖8 測線2重力異常剖面

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

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

3 結論

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

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

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