張濟東，鄭 浩，王鳴謙，李 陽

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

1 引 言

對一種地球物理數(shù)據(jù)反演的多解性不可避免，往往很難準確地勘探出地下目標體。國內(nèi)外學(xué)者經(jīng)過長期探索，提出了多種削弱多解性的方法，其中最有效的方法是在反演中加入更多信息，將各種類型的地質(zhì)信息整合到地球物理反演中，并利用多種物探方法對同一研究區(qū)域進行觀測，生成同時符合地質(zhì)信息和地球物理數(shù)據(jù)的模型[1,2]，即聯(lián)合反演。聯(lián)合反演方法借助數(shù)據(jù)的互補性來減小反演的多解性，使得反演結(jié)果與實際地質(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)驗函數(shù)，進行了二維電磁重震數(shù)據(jù)的聯(lián)合反演[3]，取得了較好效果；Jiajia Sun等[4]利用模糊C均值(FCM,Fuzzy C-means)聚類方法，實現(xiàn)了巖石物理約束的同步反演。該類方法仍然具有一定的局限性，得到多種巖石的物性信息需要一定的成本和地質(zhì)信息的支撐，不同地區(qū)不同巖性的物性關(guān)系函數(shù)差別較大，無法獲得統(tǒng)一的物性關(guān)系。因此，基于巖石物性關(guān)系的聯(lián)合反演發(fā)展受到了制約。

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

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

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

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

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

(1)

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

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

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

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

τ=?m1×?m2

(2)

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

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

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

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

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

(3)

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

(4)

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

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

(5)

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

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

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

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

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

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

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

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

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

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

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

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

以一個反演網(wǎng)格為例說明劃分原理，圖2中，紅色方框為反演網(wǎng)格，將每個反演網(wǎng)格再劃分為5×5×5的網(wǎng)格，形成了更加密集的正演網(wǎng)格(黑色網(wǎng)格線)，用于地震波走時正演計算，將計算得到的核矩陣經(jīng)過相應(yīng)網(wǎng)格上的累加得到反演網(wǎng)格上的核矩陣。模型數(shù)值化在粗糙網(wǎng)格上進行，這樣的策略使得利用阻尼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個炮點(圖中紅點所示)，分別在四口井中設(shè)24個接收點(圖中藍點所示)，四口井位置坐標分別為(5，5)、(5，15)、(15，5)、(15，15)，井深均為5 km，在地面激發(fā)，井中接收，將100×100×50網(wǎng)格下正演得到的結(jié)果加入3 %高斯噪聲作為觀測數(shù)據(jù)。利用阻尼LSQR方法分別在兩種不同的網(wǎng)格下進行反演，以及利用正反演網(wǎng)格分離的策略進行反演，正演網(wǎng)格的間距為0.2 km，反演網(wǎng)格的間距為1 km，初始模型設(shè)為均勻背景速度，經(jīng)過同樣20次迭代計算后，分別對不同結(jié)果作出Y=10 km剖面和Z=4 km剖面，對比相應(yīng)的反演結(jié)果。

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

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

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

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

5 理論模型測試

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

在均勻模型中設(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，在精細網(wǎng)格上進行正演，保證正演精度，粗糙網(wǎng)格上進行反演，大大減小了反演中大型稀疏矩陣儲存和計算的困難。

圖5 理論模型Fig.5 Theoretical model

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

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

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

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

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

圖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個炮點和66個接收點，地震波由炮點向下傳播，經(jīng)過界面反射，傳播到接收點，不考慮在異常體界面發(fā)生的反射。重力觀測系統(tǒng)同以上模型一致。將正演計算的結(jié)果加入3 %高斯噪聲作為反演的觀測數(shù)據(jù)，反演初始模型設(shè)為物性參數(shù)為0 g/cm3和4 km/s的均勻模型，進行單獨反演和聯(lián)合反演。

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

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

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

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

6 結(jié) 論

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

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

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