基于交叉梯度函數(shù)的三維重震聯(lián)合反演

張濟(jì)東，鄭 浩，王鳴謙，李 陽

(中國石化石油物探技術(shù)研究院，江蘇 南京 211103)

1 引 言

對一種地球物理數(shù)據(jù)反演的多解性不可避免，往往很難準(zhǔn)確地勘探出地下目標(biāo)體。國內(nèi)外學(xué)者經(jīng)過長期探索，提出了多種削弱多解性的方法，其中最有效的方法是在反演中加入更多信息，將各種類型的地質(zhì)信息整合到地球物理反演中，并利用多種物探方法對同一研究區(qū)域進(jìn)行觀測，生成同時符合地質(zhì)信息和地球物理數(shù)據(jù)的模型[1,2]，即聯(lián)合反演。聯(lián)合反演方法借助數(shù)據(jù)的互補(bǔ)性來減小反演的多解性，使得反演結(jié)果與實(shí)際地質(zhì)情況更加一致。

截至目前國內(nèi)外學(xué)者所研究的聯(lián)合反演可劃分為兩大類：基于巖性關(guān)系的聯(lián)合反演方法和基于結(jié)構(gòu)相似性的聯(lián)合反演方法。

基于巖性關(guān)系的聯(lián)合反演方法依賴于巖石不同物性之間的函數(shù)關(guān)系或統(tǒng)計關(guān)系。Colombo等利用密度、電阻率和速度三者間的經(jīng)驗(yàn)函數(shù)，進(jìn)行了二維電磁重震數(shù)據(jù)的聯(lián)合反演[3]，取得了較好效果；Jiajia Sun等[4]利用模糊C均值(FCM,Fuzzy C-means)聚類方法，實(shí)現(xiàn)了巖石物理約束的同步反演。該類方法仍然具有一定的局限性，得到多種巖石的物性信息需要一定的成本和地質(zhì)信息的支撐，不同地區(qū)不同巖性的物性關(guān)系函數(shù)差別較大，無法獲得統(tǒng)一的物性關(guān)系。因此，基于巖石物性關(guān)系的聯(lián)合反演發(fā)展受到了制約。

基于結(jié)構(gòu)相似性的聯(lián)合反演方法利用物性模型的結(jié)構(gòu)相似性進(jìn)行相互約束，因此不必獲得具體的巖石物性關(guān)系，在物性關(guān)系復(fù)雜的地質(zhì)條件下適用性更好。Haber最早提出了通過模型曲率來約束結(jié)構(gòu)相似性[5]，該方法未考慮模型梯度方向問題，只適合應(yīng)用于邊緣檢測。隨后， Gallardo和Meju利用交叉梯度函數(shù)定量衡量模型間的結(jié)構(gòu)相似性[6]，并進(jìn)行了二維地震和直流電法的聯(lián)合反演[7]，該方法通過對不同物性模型梯度之間叉乘來識別結(jié)構(gòu)邊界，此后基于交叉梯度函數(shù)的聯(lián)合反演以其靈活穩(wěn)定的優(yōu)點(diǎn)成為研究的熱點(diǎn)[7-9]。二維交叉梯度聯(lián)合反演已經(jīng)日趨成熟，并廣泛應(yīng)用于實(shí)際，高級和張海江等在Um等[10]和Bennington等[11]的基礎(chǔ)上對交叉梯度約束的三維聯(lián)合反演算法做了進(jìn)一步的改進(jìn)[12,13]，使得算法能夠更加方便地運(yùn)用于三維多種數(shù)據(jù)的聯(lián)合反演中。

重力勘探的優(yōu)勢在于受自然條件限制少、野外勘探效率高、成本低，在勘探初期應(yīng)用廣泛。但受到體積效應(yīng)的影響，重力勘探的分辨率較低，反演多解性強(qiáng)，難以精準(zhǔn)地恢復(fù)出物性分布情況。地震勘探作為最重要的物探方法被廣泛地應(yīng)用于石油勘探中，其分辨率較高，但同時成本高，應(yīng)用條件有限，在低降速帶較厚地區(qū)和復(fù)雜山區(qū)等應(yīng)用效果較差?？赏ㄟ^聯(lián)合地震數(shù)據(jù)來提高重力反演的分辨率，也可以對缺少地震波覆蓋的區(qū)域增加重力數(shù)據(jù)的約束，同時由于地層速度和密度的變化有一定的相關(guān)性，因此本文選用重力數(shù)據(jù)和地震初至走時數(shù)據(jù)進(jìn)行聯(lián)合反演研究。

目前基于交叉梯度函數(shù)的聯(lián)合反演仍存在靈敏度矩陣數(shù)值差異較大，反演系統(tǒng)差異較大以及大型稀疏矩陣儲存和計算等問題。針對這些問題，本文利用一種基于模型更新量及交叉梯度結(jié)構(gòu)約束結(jié)合的聯(lián)合反演框架，分別交替進(jìn)行數(shù)據(jù)擬合和交叉梯度結(jié)構(gòu)約束，并提出正反演網(wǎng)格相分離的網(wǎng)格劃分策略，平衡了計算效率和精度。最后采用三維地震波初至和重力數(shù)據(jù)，設(shè)計了理論模型，驗(yàn)證了算法的有效性。

2 單一類型數(shù)據(jù)反演相關(guān)問題

對于類似地震初至走時和重力數(shù)據(jù)的反演問題，包含了數(shù)據(jù)擬合項以及模型約束和模型光滑約束正則化項的目標(biāo)函數(shù)表達(dá)式為：

(1)

式中,φd為數(shù)據(jù)擬合項，φm為模型約束項，φl為模型光滑約束項，α和β分別為模型約束項和模型光滑約束項的權(quán)重，g為正演算子，重力正演采用基于幾何格架的快速算法[14]，地震初至層析采用快速行進(jìn)三維有限差分算法(FMM,Fast Marching Method)[15,16]，Cd為觀測數(shù)據(jù)d的協(xié)方差矩陣，Cm為模型參數(shù)m的協(xié)方差矩陣，m0為參考模型，Clm為光滑模型的協(xié)方差矩陣，D為光滑度算子，本文采用一階導(dǎo)數(shù)進(jìn)行度量。

對于地震波初至層析問題，采用阻尼LSQR(LSQR,Least Square QR)方法迭代求解[17]，該方法在數(shù)據(jù)誤差較大的情況下仍能保證求解過程穩(wěn)定高效，并且其算法容易實(shí)現(xiàn)。對于重力反演問題，采用高斯牛頓法進(jìn)行迭代，該方法收斂速度更快、對線性反演問題適用性更好。同時，反演過程中引入懲罰函數(shù)、深度加權(quán)函數(shù)來進(jìn)一步提高反演的分辨率，降低反演的多解性。

3 聯(lián)合反演框架

Gallardo等(2004)提出利用交叉梯度函數(shù)度量模型間的結(jié)構(gòu)相似性[6]，兩模型m1、m2的交叉梯度表示為：

τ=?m1×?m2

(2)

該函數(shù)具有以下性質(zhì)：

1)當(dāng)兩種物性參數(shù)分布結(jié)構(gòu)相似，即兩者的等值面平行或一種物性均勻分布，此時交叉梯度值為零；

2)當(dāng)兩種物性參數(shù)分布結(jié)構(gòu)不相似，即兩者等值面不平行，此時交叉梯度值不為零。

Bennington等(2015)已經(jīng)給出了三維情況下的交叉梯度及其偏導(dǎo)數(shù)的計算公式[14]，不再贅述，本文采用中心差分格式進(jìn)行求解。

首先對非線性問題進(jìn)行線性化，即對于交叉梯度函數(shù)τ，在初始模型(ms0，mg0)處進(jìn)行一階泰勒展開，寫成矩陣形式為：

(3)

式中，ms、mg分別為兩參數(shù)模型，即重震聯(lián)合反演中分別為速度模型和密度模型。當(dāng)兩模型結(jié)構(gòu)相似時，交叉梯度值τ(ms,mg)為零，得到矩陣形式：

(4)

式中，Δms、Δmg分別為速度模型和密度模型的更新量。

將不同類型數(shù)據(jù)單獨(dú)反演的模型更新量與交叉梯度約束相結(jié)合，得到加入結(jié)構(gòu)約束的線性反演方程組，聯(lián)合反演系統(tǒng)可寫為：

(5)

式中,I為單位矩陣，Δms、Δmg分別為加入交叉梯度約束后的兩個模型更新量，Δms0、Δmg0分別為兩個模型單獨(dú)反演的模型更新量，αs、αg分別為平衡不同模型更新量數(shù)量級差異的參數(shù)，可通過事先評估單獨(dú)反演模型更新量的數(shù)量級來確定，βt為平衡數(shù)據(jù)擬合和結(jié)構(gòu)約束的參數(shù)，其值越大，結(jié)構(gòu)約束的權(quán)重越大，通過調(diào)整αs、αg、βt三個參數(shù)來實(shí)現(xiàn)反演過程中各項約束的平衡。

圖1為該反演策略的流程圖，具體可通過以下幾個步驟來實(shí)現(xiàn)：

1)給出初始速度模型與密度模型；

2)根據(jù)阻尼LSQR方法和高斯牛頓方法分別求解得到單獨(dú)反演的速度模型更新量與密度模型更新量；

3)求解出更新前模型間的交叉梯度值和交叉梯度偏導(dǎo)數(shù)，帶入公式(5)進(jìn)行聯(lián)合反演，求解這樣的線性方程組，得到結(jié)構(gòu)約束下的模型更新量；

4)更新速度與密度模型，計算出模型之間的交叉梯度以及地震走時數(shù)據(jù)和重力數(shù)據(jù)的殘差；

5)判斷結(jié)果是否滿足標(biāo)準(zhǔn)，收斂標(biāo)準(zhǔn)一般為數(shù)據(jù)擬合情況、模型更新量大小、交叉梯度值大小等，若滿足則得到結(jié)果，不滿足則重新回到第2)步。

4 反演網(wǎng)格劃分

4.1 網(wǎng)格劃分策略

反演過程是對模型進(jìn)行數(shù)值求解的過程，需要對模型進(jìn)行數(shù)值化，數(shù)值化過程中存在反演精度和計算效率平衡的問題。若將網(wǎng)格劃分較密，通常情況下可以獲得分辨率較高的反演結(jié)果，但在反演中，對欠定問題的求解會存在困難，影響計算效率。若網(wǎng)格劃分較稀疏，反演的求解會變得更加容易，但恢復(fù)得到的模型會更加粗糙，分辨率較低，可能分辨不出需要獲得的構(gòu)造和異常體。因此在反演過程中，對模型網(wǎng)格的剖分需要謹(jǐn)慎處理，在保證分辨能力的同時，實(shí)現(xiàn)計算效率的最大化。

在聯(lián)合反演中，上述問題被放大，如果將網(wǎng)格劃分過密，則反演中交叉梯度偏導(dǎo)數(shù)的規(guī)模巨大，計算和儲存較為困難；如果將網(wǎng)格劃分過于稀疏，則地震波走時正演的精度較低，同樣會影響反演的效果。為了解決聯(lián)合反演中存在的正演精確度和偏導(dǎo)數(shù)矩陣儲存運(yùn)算的問題，提出了正演網(wǎng)格與反演網(wǎng)格相分離的網(wǎng)格劃分策略。在反演過程中，每個反演網(wǎng)格內(nèi)劃分多個正演網(wǎng)格，正演計算在更加密集的網(wǎng)格上進(jìn)行，既保證了正演的準(zhǔn)確性，又解決了偏導(dǎo)數(shù)矩陣計算和儲存的困難。

圖2 反演網(wǎng)格示意Fig.2 Inversion grid

以一個反演網(wǎng)格為例說明劃分原理，圖2中，紅色方框?yàn)榉囱菥W(wǎng)格，將每個反演網(wǎng)格再劃分為5×5×5的網(wǎng)格，形成了更加密集的正演網(wǎng)格(黑色網(wǎng)格線)，用于地震波走時正演計算，將計算得到的核矩陣經(jīng)過相應(yīng)網(wǎng)格上的累加得到反演網(wǎng)格上的核矩陣。模型數(shù)值化在粗糙網(wǎng)格上進(jìn)行，這樣的策略使得利用阻尼LSQR迭代過程的計算量減小，聯(lián)合反演中交叉梯度偏導(dǎo)數(shù)的計算和儲存也更加方便。

4.2 網(wǎng)格劃分測試

筆者利用三維地震VSP(VSP，Vertical Seismic Profile)走時層析測試網(wǎng)格劃分策略的效果，建立三維速度模型如圖3(a)所示，模型大小為20 km×20 km×10 km，背景速度設(shè)為2 km/s，在模型中設(shè)立速度為3 km/s的高速異常，異常體的大小為4 km×4 km×3 km，頂面距地表2 km。分別將模型劃分為兩種網(wǎng)格，網(wǎng)格一：20×20×10個大小相等的網(wǎng)格，網(wǎng)格間距為1 km；網(wǎng)格二：100×100×50個大小相等的網(wǎng)格，網(wǎng)格間距為0.2 km。觀測系統(tǒng)如圖3(b)所示，在地面設(shè)36個炮點(diǎn)(圖中紅點(diǎn)所示)，分別在四口井中設(shè)24個接收點(diǎn)(圖中藍(lán)點(diǎn)所示)，四口井位置坐標(biāo)分別為(5，5)、(5，15)、(15，5)、(15，15)，井深均為5 km，在地面激發(fā)，井中接收，將100×100×50網(wǎng)格下正演得到的結(jié)果加入3 %高斯噪聲作為觀測數(shù)據(jù)。利用阻尼LSQR方法分別在兩種不同的網(wǎng)格下進(jìn)行反演，以及利用正反演網(wǎng)格分離的策略進(jìn)行反演，正演網(wǎng)格的間距為0.2 km，反演網(wǎng)格的間距為1 km，初始模型設(shè)為均勻背景速度，經(jīng)過同樣20次迭代計算后，分別對不同結(jié)果作出Y=10 km剖面和Z=4 km剖面，對比相應(yīng)的反演結(jié)果。

圖3 實(shí)際模型和觀測系統(tǒng)Fig.3 Actual model and observation system

圖4(a)和圖4(d)為網(wǎng)格一進(jìn)行反演的結(jié)果剖面圖，黑色虛線框表示異常的真實(shí)位置，結(jié)果中對異常的恢復(fù)不夠準(zhǔn)確，其原因可能為在此網(wǎng)格下正演的過程不夠準(zhǔn)確。且反演過程中為保證穩(wěn)定，阻尼因子選取了較大的值，導(dǎo)致了結(jié)果的擬合差較大。圖4(b)和圖4(e)為網(wǎng)格二下的反演結(jié)果，黑色虛線框表示異常的真實(shí)位置，相比于網(wǎng)格一，對于異常位置的恢復(fù)較好，大致上可以得出高速體的位置，在高速體的邊緣產(chǎn)生了部分拖尾異常，對異常形態(tài)的恢復(fù)有一定偏差。當(dāng)數(shù)據(jù)量不變，模型參數(shù)數(shù)量變大后，反演問題變?yōu)楦鼜?fù)雜的欠定問題，其運(yùn)算量大大增加。將反演結(jié)果與用正反演網(wǎng)格分離策略的反演結(jié)果(圖4(c)和圖4(f))對比，經(jīng)過同樣次數(shù)的迭代后，正反演網(wǎng)格分離網(wǎng)格下的反演結(jié)果最好，對高速體的恢復(fù)最為準(zhǔn)確，在這樣的網(wǎng)格劃分策略下，正演精確度高，反演計算量較小。

圖4 不同網(wǎng)格反演結(jié)果對比Fig.4 Comparison of inversion results of different grids

通過該模型試驗(yàn)，驗(yàn)證了這樣的網(wǎng)格劃分策略能夠有效地平衡正演精度和反演計算量，在后續(xù)的聯(lián)合反演模型測試中將應(yīng)用這種策略。

5 理論模型測試

利用理論數(shù)據(jù)對上述的重震聯(lián)合反演算法進(jìn)行驗(yàn)證，通過比較聯(lián)合反演和單獨(dú)反演的結(jié)果來驗(yàn)證算法的有效性。

在均勻模型中設(shè)立兩異常體，分別為高速高密度體和低速低密度體，如圖5所示，高速高密度體的速度和剩余密度值分別為6 km/s和1 g/cm3，低速低密度體的速度和剩余密度值分別為2 km/s和-1 g/cm3，背景速度值為4 km/s。將模型劃分為100×100×50的正演網(wǎng)格，網(wǎng)格間距為0.2 km，反演網(wǎng)格數(shù)量為20×20×10，在精細(xì)網(wǎng)格上進(jìn)行正演，保證正演精度，粗糙網(wǎng)格上進(jìn)行反演，大大減小了反演中大型稀疏矩陣儲存和計算的困難。

圖5 理論模型Fig.5 Theoretical model

在Z=0觀測面上均勻設(shè)立400個重力數(shù)據(jù)觀測點(diǎn)，測點(diǎn)間隔為1 km，地震初至觀測系統(tǒng)與圖3(b)一致，對模型進(jìn)行正演計算，得到的結(jié)果加入3 %高斯噪聲作為觀測數(shù)據(jù)。利用上述算法進(jìn)行反演，設(shè)初始模型為速度4 km/s、剩余密度0的均勻模型。慢度更新量和密度更新量在數(shù)量級上相同，因此取αs為1，αg為1，參數(shù)βt從10-5至105以間隔10倍關(guān)系分別計算地震走時、重力觀測數(shù)據(jù)的擬合差與交叉梯度值，經(jīng)試驗(yàn)獲得合適的參數(shù)βt值為10。對單獨(dú)反演結(jié)果和聯(lián)合反演結(jié)果分別做Y=10 km的垂直剖面和Z=4 km的水平剖面，如圖6、圖7所示，圖中黑色虛線框表示異常的真實(shí)位置。

圖6 單獨(dú)反演結(jié)果剖面Fig.6 Single inversion result

圖7 聯(lián)合反演結(jié)果剖面Fig.7 Joint inversion result

由反演結(jié)果看出，地震單獨(dú)反演結(jié)果對異常的恢復(fù)不夠準(zhǔn)確，高速異常位置的恢復(fù)有所偏離，低速異常的幅值和位置與真實(shí)情況有較大差別，且在模型左側(cè)產(chǎn)生了大面積的虛假低速異常。這是由于地震波傳播過程中，穿過低速體的射線較少，導(dǎo)致了僅利用地震波走時數(shù)據(jù)難以通過反演準(zhǔn)確識別低速體。重力單獨(dú)反演結(jié)果可以大致反應(yīng)地下異常的形態(tài)，但異常的幅值小于真實(shí)值，通過反演結(jié)果的切片圖可以看出，橫向分辨率較好，可以較好地刻畫出異常在水平方向上的邊界，同時重力反演能夠有效地識別出負(fù)異常，通過重力數(shù)據(jù)的補(bǔ)充，地震聯(lián)合反演結(jié)果得到了明顯的改善，圖中可識別出正負(fù)兩異常體，且位置、幅值和理論模型更加接近，位于模型左側(cè)的虛假異常也得到了一定程度的削弱，同時聯(lián)合反演得到的速度和密度模型在結(jié)構(gòu)上更加相似。通過地震數(shù)據(jù)和重力數(shù)據(jù)相互約束，反演的分辨率得到了提高，對異常的恢復(fù)更加準(zhǔn)確。

為了進(jìn)一步分析聯(lián)合反演和單獨(dú)反演對物性參數(shù)的恢復(fù)情況，分別做出反演結(jié)果網(wǎng)格點(diǎn)的物性值分布圖，如圖8所示。從圖8的交會圖中可以看出，單獨(dú)反演結(jié)果的物性參數(shù)集中在(4，0)點(diǎn)附近，呈不規(guī)則分布，且大部分網(wǎng)格單元的物性點(diǎn)都集中在零點(diǎn)附近，與真實(shí)的物性值有較大差距，這樣的分布是光滑約束在物性值上的體現(xiàn)。而聯(lián)合反演結(jié)果交會圖顯示出兩條明顯的特征線指向真實(shí)物性點(diǎn)，雖然仍未到達(dá)真實(shí)物性點(diǎn)，其指向表明了分辨物性異常的能力。

圖8 反演結(jié)果物性交會圖Fig.8 Inversion result crossplot

設(shè)計了更加復(fù)雜的模型來測試聯(lián)合反演的效果。建立三個棱柱異常體模型，如圖9(a)所示，幾何參數(shù)分別為，棱柱體(小)長2 km，寬4 km，高2 km，頂面埋深3 km；棱柱體(中)長3 km，寬4 km，高3 km，頂面埋深2 km；棱柱體(大)長4 km，寬4 km，高4 km，頂面埋深1 km。三個棱柱異常體的物性參數(shù)相同，剩余密度為1 g/cm3，地震波速度為6 km/s，圍巖速度為4 km/s，設(shè)模型中異常體的下方存在地震反射界面(Z=7 km平面)。地震觀測系統(tǒng)為地面觀測，如圖9(b)，在地面(Z=0平面)均勻布設(shè)25個炮點(diǎn)和66個接收點(diǎn)，地震波由炮點(diǎn)向下傳播，經(jīng)過界面反射，傳播到接收點(diǎn)，不考慮在異常體界面發(fā)生的反射。重力觀測系統(tǒng)同以上模型一致。將正演計算的結(jié)果加入3 %高斯噪聲作為反演的觀測數(shù)據(jù)，反演初始模型設(shè)為物性參數(shù)為0 g/cm3和4 km/s的均勻模型，進(jìn)行單獨(dú)反演和聯(lián)合反演。

圖9 理論模型和觀測系統(tǒng)Fig.9 Actual model and observation system

得到單獨(dú)反演和聯(lián)合反演結(jié)果，對結(jié)果分別做出Y=10 km垂直剖面和Z=4 km水平剖面，如圖10所示，圖11所示。單獨(dú)反演結(jié)果中，密度模型剖面無法分辨出三個棱柱，只對體積較大的棱柱異常體恢復(fù)較好。速度模型能夠較準(zhǔn)確地恢復(fù)出三個棱柱體異常，但其形態(tài)、幅值和真實(shí)情況有一定偏差。相較于單獨(dú)反演，聯(lián)合反演結(jié)果對于異常的形態(tài)和物性值的恢復(fù)都有了很大程度上的提升，密度模型中能夠有效分辨出三個棱柱異常體，速度模型對異常形態(tài)和幅值的恢復(fù)都更加接近真實(shí)情況。通過地震數(shù)據(jù)的補(bǔ)充，重力反演的分辨率得到了顯著提高，同時，在重力數(shù)據(jù)的約束下，地震反演結(jié)果也更加符合實(shí)際情況。

圖10 單獨(dú)反演結(jié)果剖面Fig.10 Single inversion result

圖11 聯(lián)合反演結(jié)果剖面Fig.6 Joint inversion result

6 結(jié) 論

本文基于結(jié)構(gòu)相似性原理，對重力和地震初至走時數(shù)據(jù)的交叉梯度聯(lián)合反演方法進(jìn)行了研究，提出了相應(yīng)的反演策略，通過對比分析模型測試結(jié)果，得到如下結(jié)論：

1)聯(lián)合反演相較于單獨(dú)反演，獲得的物性分布更加準(zhǔn)確。通過地震波走時和重力數(shù)據(jù)間的相互約束，解決了地震方法射線在低速區(qū)分布少反演結(jié)果差和重力反演縱向分辨率低的問題，有效減小了反演的多解性，提高了反演的分辨率。

2)本文提出的新的網(wǎng)格劃分策略有效平衡了正演的精度和反演的計算效率。模型測試證明，新的網(wǎng)格劃分下反演能達(dá)到較好的效果，提高了計算效率的同時，也不會影響到反演的分辨率，有效解決了聯(lián)合反演中交叉梯度偏導(dǎo)數(shù)儲存和求解的難題。