本發(fā)明涉及一種基于各向異性馬爾科夫隨機域的疊前地震參數(shù)反演方法，屬于地震資料反演的技術(shù)領(lǐng)域。

背景技術(shù)：

地震反演是定量解釋地震資料和預測地層彈性參數(shù)的重要手段，在油氣資源的勘探和開發(fā)過程中應用廣泛。根據(jù)地震反演所采用的正演模型，可以分為基于波動理論的波動方程反演和基于robinson褶積模型反演兩大類。在實際工作中主要使用基于褶積模型的反演方法，該類方法易于實現(xiàn)，并能夠滿足精度要求。疊前地震反演是一種基于褶積模型的反演，它基于振幅隨炮檢距變化理論，直接利用信息豐富的疊前地震道集數(shù)據(jù)，可以同步獲得縱波速度、橫波速度和密度等多種彈性參數(shù)，在彈性參數(shù)反演和儲層流體識別中發(fā)揮著不可替代的作用。

但疊前地震反演尚存諸多問題，如關(guān)于計算疊前道集反射系數(shù)，現(xiàn)有的文獻及商業(yè)化軟件多數(shù)使用精確左普利茲方程的線性近似式，但這些近似式需要一定的假設(shè)條件，會造成計算誤差。此外，疊前地震的多參數(shù)反演是一個高度非線性問題，需要解決反演病態(tài)問題、地震帶限問題。有關(guān)反演病態(tài)問題可以通過正則化方法得到解決，尤其是基于馬爾科夫隨機域的邊界保護正則化，而帶限問題可以使用測井數(shù)據(jù)的約束反演。在疊前地震反演方面，馬爾科夫隨機域技術(shù)沒有得到足夠的重視，現(xiàn)有研究中主要使用一些各向同性的馬爾科夫隨機域(或者各向同性的平滑因子)，這忽略了馬爾科夫隨機域的各向異性特點，尤其是對于層狀和復雜介質(zhì)。在層狀和復雜介質(zhì)中，水平、垂直和對角方向上模型梯度值的先驗作用差異很大。因此，對傳統(tǒng)的馬爾科夫隨機域技術(shù)需要進行改進和完善，以適應對于層狀或復雜地層的疊前地震反演，在解決疊前多參數(shù)反演病態(tài)問題的同時，獲取符合真實地質(zhì)構(gòu)造的彈性參數(shù)模型，以更好的指導復雜地質(zhì)構(gòu)造下油氣資源的勘探、開發(fā)和利用。

技術(shù)實現(xiàn)要素：

本發(fā)明所要解決的技術(shù)問題在于克服現(xiàn)有技術(shù)的不足，提供一種基于各向異性馬爾科夫隨機域的疊前地震參數(shù)反演方法，解決現(xiàn)有方法忽略了馬爾科夫隨機域的各向異性特點，尤其是對于層狀和復雜介質(zhì)，存在疊前多參數(shù)反演的非線性和不穩(wěn)定問題。

本發(fā)明具體采用以下技術(shù)方案解決上述技術(shù)問題：

基于各向異性馬爾科夫隨機域的疊前地震參數(shù)反演方法，包括以下步驟：

步驟一、建立疊前地震參數(shù)的目標函數(shù)，所述目標函數(shù)包含數(shù)據(jù)項和先驗約束項；

步驟二、利用佐普利茲方程計算地層的縱波反射系數(shù)，得到的縱波反射系數(shù)與零相位ricker子波褶積后生成合成地震記錄；根據(jù)測量所得地震記錄與合成地震記錄之間誤差的二范數(shù)獲得目標函數(shù)的數(shù)據(jù)項；

步驟三、利用各向異性擴散法中的邊界終止函數(shù)計算各向異性馬爾科夫隨機域權(quán)系數(shù)，以得到數(shù)據(jù)點所在不同方向的權(quán)系數(shù)值；

步驟四、基于各向異性馬爾科夫隨機域?qū)⑺脭?shù)據(jù)點所在不同方向的權(quán)系數(shù)值和縱波速度、橫波速度先和密度先驗約束項結(jié)合，獲得目標函數(shù)的先驗約束項；

步驟五、提取待反演區(qū)域的測井數(shù)據(jù)，對測井數(shù)據(jù)進行反演參數(shù)的對數(shù)線性擬合，得到擬合系數(shù)及擬合誤差，用以控制目標函數(shù)極小尋優(yōu)過程的穩(wěn)定性；

步驟六、根據(jù)所得數(shù)據(jù)項和先驗約束項代入所建立的疊前地震參數(shù)目標函數(shù)，及確定目標函數(shù)，并基于快速模擬退火算法和擬合系數(shù)及擬合誤差進行目標函數(shù)的極小尋優(yōu)；

步驟七、完成對目標函數(shù)的迭代尋優(yōu)，輸出反演結(jié)果。

進一步地，作為本發(fā)明的一種優(yōu)選技術(shù)方案，所述步驟一建立疊前地震參數(shù)的目標函數(shù)為：

j(vp,vs,ρ)＝j(luò)1(vp,vs,ρ)+j2(vp,vs,ρ)

其中，j1為數(shù)據(jù)項，j2為先驗約束項，vp、vs和ρ分別為待反演地層的縱波速度、橫波速度和密度參數(shù)。

進一步地，作為本發(fā)明的一種優(yōu)選技術(shù)方案，所述步驟二獲得目標函數(shù)的數(shù)據(jù)項為：

其中，d為測量所得地震記錄；w為零相位ricker子波；r為縱波反射系數(shù)；θ為入射角；t為地震記錄的采樣時間。

進一步地，作為本發(fā)明的一種優(yōu)選技術(shù)方案，所述步驟三計算各向異性馬爾科夫隨機域權(quán)系數(shù)α采用公式:

其中，g(▽v)為邊界終止函數(shù)；▽v為梯度值；q為擴散系數(shù)，為預設(shè)定的常數(shù)。

進一步地，作為本發(fā)明的一種優(yōu)選技術(shù)方案，所述步驟四獲得目標函數(shù)的先驗約束項為：

其中，c為數(shù)據(jù)點鄰域內(nèi)點集；αc為各向異性馬爾科夫隨機域權(quán)系數(shù)；d為相鄰兩個數(shù)據(jù)點的一階差分；φ(t)為勢函數(shù)；vp、vs和ρ分別為待反演地層的縱波速度、橫波速度和密度參數(shù)，δp、δs和δd分別為與之對應的正則化參數(shù)。

進一步地，作為本發(fā)明的一種優(yōu)選技術(shù)方案，所述步驟五對測井數(shù)據(jù)進行反演參數(shù)的對數(shù)線性擬合采用公式：

lnvs＝k·lnvp+kc+δls

lnρ＝m·lnvp+mc+δld

其中，vp為待反演的縱波速度；vs為待反演的橫波速度；ρ為待反演的密度參數(shù)；δls和δld為擬合誤差；kc和mc為擬合系數(shù)。

本發(fā)明采用上述技術(shù)方案，能產(chǎn)生如下技術(shù)效果：

本發(fā)明是一種基于各向異性馬爾科夫隨機域的疊前地震多參數(shù)的反演方法，該方法通過各向異性擴散法計算各向異性馬爾科夫權(quán)系數(shù)，用以校正層狀或復雜地層對反演目標函數(shù)先驗約束項的影響，并基于精確佐普里茲方程進行正演計算，并采用測井數(shù)據(jù)的對數(shù)線性擬合關(guān)系來控制目標函數(shù)全局優(yōu)化反演的穩(wěn)定性。分別對合成地震數(shù)據(jù)及實際地震資料進行疊前反演，并取得準確的反演結(jié)果。該方法有效解決疊前多參數(shù)反演的非線性和不穩(wěn)定問題，可以有效的保護反演地層參數(shù)的橫向連續(xù)性和縱向分層性，特別適用于復雜地質(zhì)條件下的彈性參數(shù)反演和儲層參數(shù)預測，對油氣資源的勘探和開發(fā)具有重要指導和現(xiàn)實意義。

本發(fā)明的方法，可以有效解決疊前多參數(shù)反演的非線性和不穩(wěn)定問題，對于層狀或復雜地層，各向異性馬爾科夫權(quán)系數(shù)可以校正地層的各向異性對反演結(jié)果的影響，反演的模型參數(shù)可以準確反映地層的層狀特征，并對斷層和邊界起到很好的保護效果。

附圖說明

圖1為本發(fā)明方法的流程圖。

圖2為各向異性馬爾科夫隨機域權(quán)系數(shù)的示意圖。

圖3(a)為合成地震數(shù)據(jù)測試的縱波速度理論模型；圖3(b)為合成地震數(shù)據(jù)測試的橫波速度理論模型；圖3(c)為合成地震數(shù)據(jù)測試的密度速度理論模型。

圖4(a)為采用縱波速度模型計算的水平方向各向異性馬爾科夫權(quán)系數(shù)；圖4(b)為采用縱波速度模型計算的垂直方向各向異性馬爾科夫權(quán)系數(shù)。

圖5(a)為合成地震數(shù)據(jù)測試的縱波速度反演結(jié)果；圖5(b)為合成地震數(shù)據(jù)測試的橫波速度反演結(jié)果；圖5(c)為合成地震數(shù)據(jù)測試的密度速度反演結(jié)果。

圖6(a)實際地震資料的縱波速度反演結(jié)果；圖6(b)實際地震資料的橫波速度反演結(jié)果。圖6(c)實際地震資料的密度反演結(jié)果。

具體實施方式

下面結(jié)合說明書附圖對本發(fā)明的實施方式進行描述。

本發(fā)明提出的基于各向異性馬爾科夫隨機域的疊前地震多參數(shù)的反演方法，該方法如圖1所示，具體包括以下步驟：

步驟一、設(shè)置疊前地震參數(shù)的目標函數(shù)。本反演方法的目標函數(shù)包含數(shù)據(jù)項和先驗約束項：

j(vp,vs,ρ)＝j(luò)1(vp,vs,ρ)+j2(vp,vs,ρ)(1)

其中，j1為數(shù)據(jù)項，j2為先驗約束項，vp、vs和ρ分別為待反演地層的縱波速度、橫波速度和密度參數(shù)。

步驟二、使用佐普利茲方程計算地層的縱波反射系數(shù)，得到的縱波反射系數(shù)與零相位ricker子波褶積后生成合成地震記錄。目標函數(shù)的數(shù)據(jù)項為測量所得地震記錄與合成地震記錄之間誤差的二范數(shù)，具體表現(xiàn)形式如下：

其中，d為測量所得地震記錄；w為震源子波，這里采用零相位ricker子波近似代替；r為地層的縱波反射系數(shù)；θ為入射角；t為地震記錄的采樣時間。

步驟三、計算各向異性馬爾科夫隨機域權(quán)系數(shù)。本發(fā)明使用各向異性擴散法中的邊界終止函數(shù)g(▽v)計算層狀或復雜地層中的各向異性馬爾科夫隨機域權(quán)系數(shù)α，其具體表達形式如下：

其中，▽表示梯度算子，計算所需要的梯度值▽v采用地層的縱波速度參數(shù)。q為擴散系數(shù)，為預設(shè)定的常數(shù)。

本發(fā)明采用一階馬爾科夫鄰域計算各向異性馬爾科夫隨機域權(quán)系數(shù)，一階馬爾科夫鄰域?qū)挠?個不同方向的▽v，具體計算形式如下：

其中，下標(i,j)表示縱波速度參數(shù)的坐標；下標s和n表示2個垂直方向，e和w表示2個水平方向，wn、ws、en和es表示4個對角方向。

因此，可得到數(shù)據(jù)點周圍8個不同方向的權(quán)系數(shù)值，如圖2所示，即αn，αs，αe，αw，αwn，αws，αen，αns。

步驟四、計算先驗約束項，先驗約束項由縱波速度、橫波速度先和密度先驗約束項3部分組成：

j2(vp,vs,ρ)＝j(luò)2p(vp)+j2s(vs)+j2d(ρ)(5)

針對層狀地層的各向異性，基于各向異性馬爾科夫隨機域的先驗約束項為：

其中，c為數(shù)據(jù)點鄰域內(nèi)點集，這里使用一階馬爾科夫領(lǐng)域，對于某個數(shù)據(jù)點，對應與之相鄰的8個方向的點集；d表示相鄰兩個數(shù)據(jù)點的一階差分；δp、δs和δd分別為與之對應的正則化參數(shù)。此外，先驗約束項是通過施加正則化勢函數(shù)φ實現(xiàn)的，這里采用具有邊界保護性質(zhì)的勢函數(shù)：

各向異性馬爾科夫隨機域權(quán)系數(shù)αc由步驟三得到，用以校正層狀地層的各向異性。對層狀地層，通常水平方向的權(quán)系數(shù)大于垂直方向的權(quán)系數(shù)，確保橫向梯度值較小，以保護地層的橫向連續(xù)性；而層間較大梯度值得以保留，以體現(xiàn)地層的縱向分層性。

步驟五：提取待反演區(qū)域的測井數(shù)據(jù)，對測井數(shù)據(jù)分別進行反演參數(shù)vp和vs及vp和ρ的對數(shù)線性擬合，得到對應的擬合系數(shù)kc和mc及擬合誤差δls和δld，用以控制目標函數(shù)極小尋優(yōu)過程的穩(wěn)定性，其具體形式如下：

反演參數(shù)的對數(shù)線性擬合關(guān)系式控制步驟六中快速模擬退火過程的穩(wěn)定性。

步驟六：根據(jù)所得數(shù)據(jù)項和先驗約束項代入所建立的疊前地震參數(shù)目標函數(shù)，確定目標函數(shù)及并基于快速模擬退火算法進行目標函數(shù)的極小尋優(yōu)，實現(xiàn)最優(yōu)化。通過步驟二、步驟三和步驟四建立的目標函數(shù)的具體形式為：

使用基于全局優(yōu)化算法的快速模擬退火算法進行目標函數(shù)的極小尋優(yōu)，快速模擬退火算法結(jié)合了三參數(shù)的同步反演。此外，對于橫波速度的反演，實際反演的參數(shù)是δls，最終vs通過由步驟五建立的對數(shù)線性擬合關(guān)系式得到。

三參數(shù)擾動量的具體表現(xiàn)形式為：

lnvs^(m+1)＝k·lnvp^(m+1)+kc+δls^(m+1)(11)

密度參數(shù)擾動量需要滿足附加條件：

m·lnvp^(m+1)+mc-δld≤lnρ^(m+1)≤m·lnvp^(m+1)+mc+δld(12)

式中vp^(m)、δls^(m)和ρ^(m)是當前模型參數(shù)值；vp^(m+1)、δls^(m+1)和ρ^(m+1)是擾動后的參數(shù)值；t^(m)是當前的溫度值；[vpmin,vpmax]、[δlsmin,δlsmax]和[ρmin,ρmax]是3個參數(shù)取值的范圍；ξ為分布在[0，1]的隨機數(shù)；sign(·)為符號函數(shù)；k和kc是測井數(shù)據(jù)中vp和vs的擬合系數(shù)；m和mc是測井數(shù)據(jù)中vp和ρ的擬合系數(shù)；δls和δld分別是對應的擬合誤差。

步驟七、完成對目標函數(shù)的迭代尋優(yōu)，輸出反演結(jié)果。

為了驗證本發(fā)明的方法能夠有效解決疊前多參數(shù)反演的非線性和不穩(wěn)定問題，對于層狀或復雜地層，各向異性馬爾科夫權(quán)系數(shù)可以校正地層的各向異性對反演結(jié)果的影響，本發(fā)明列舉驗證例進行驗證說明。

圖2為各向異性馬爾科夫隨機域權(quán)系數(shù)示意圖。以待反演的縱波速度模型的數(shù)據(jù)點vi,j為例，各權(quán)系數(shù)αn，αs，αe，αw，αwn，αws，αen，αns對應梯度值的方位如虛線所示，不同灰度的背景表明地層速度的橫向連續(xù)性和縱向分層性。垂直方向上的權(quán)系數(shù)值較小，而水平方向上的權(quán)系數(shù)較大，以校正地層的各項異性對先驗約束項的影響。以兩個應用實例對本發(fā)明進行說明。

實例一，合成地震數(shù)據(jù)測試：

圖3(a)、圖3(b)、圖3(c)為合成地震數(shù)據(jù)測試的理論模型，其分別為縱波速度、橫波速度和密度三組參數(shù)模型。

使用精確佐普里茲方程求解該理論模型的縱波反射系數(shù)，與主頻為45hz的零相位理論ricker子波進行褶積，獲得觀測地震記錄。觀測地震記錄有61個角度道集，每個道集有5個角度道(角度間隔為3°，角度覆蓋范圍為0-15°)。

按公式(2)建立目標函數(shù)的數(shù)據(jù)項。使用精確佐普里茲方程對初始模型進行正演計算，獲得合成地震記錄。

按公式(3)和(4)計算各向異性馬爾科夫隨機域的權(quán)系數(shù)，使用縱波速度模型作為梯度值，擴散系數(shù)q取850。除模型邊界外，數(shù)據(jù)點處均可得8個對應不同方向的權(quán)系數(shù)。

圖4(a)和圖4(b)分別為水平方向和垂直方向的各向異性馬爾科夫隨機域權(quán)系數(shù)值?？梢钥吹?，在地層分界面處，垂直方向上的權(quán)系數(shù)值較小，而水平方向上的權(quán)系數(shù)較大。

按公式(5)、(6)和(7)建立目標函數(shù)的先驗約束項。

將15和45道作為虛擬井，并對測井數(shù)據(jù)按公式(8)進行vp和vs及vp和ρ的對數(shù)線性擬合，得到對應的擬合系數(shù)和擬合誤差。

在本測試中，λ1和λ2的初值為0.3和0.6，vp、δld和密度參數(shù)對應的δ的初值分別為250.0、20.0和0.18。

快速模擬退火的初始溫度為0.05，終止溫度為0.00001，溫度衰減系數(shù)為0.9，模型參數(shù)擾動按公式(10)、(11)和(12)。

對所有數(shù)據(jù)點依次進行迭代尋優(yōu)，降溫后重復以上步驟，直到達到終止溫度，輸出最終的反演結(jié)果。圖5(a)、圖5(b)、圖5(c)為合成地震數(shù)據(jù)測試的反演結(jié)果，其分別為縱波速度反演結(jié)果，橫波速度反演結(jié)果，密度速度反演結(jié)果；可以看到反演結(jié)果與理論模型非常吻合，特別是地層的層狀特征得到很好保護。

實例二，實際地震資料反演：

實際地震資料中為海上地震資料的一條二維任意線。該任意線共有1981個角度道集，每個道集有15個角度道，角度范圍為3–45°，角度間隔為3°，道間距為12.5m，時間采樣率為2ms。使用工區(qū)內(nèi)幾個主要地層的參數(shù)值建立vp和密度的初始模型，δls的初始模型為0。

采用工區(qū)內(nèi)主要地層的層速度計算各向異性馬爾科夫權(quán)系數(shù)。其他參數(shù)設(shè)置與實例一相同。

圖6(a)、圖6(b)、圖6(c)為實際地震資料的反演結(jié)果，其分別為縱波速度反演結(jié)果，橫波速度反演結(jié)果，密度速度反演結(jié)果?？梢钥吹?，反演結(jié)果較好的體現(xiàn)地層的層狀特征，目的層的連續(xù)性很好。

綜上，本發(fā)明的方法有效解決疊前多參數(shù)反演的非線性和不穩(wěn)定問題，可以有效的保護反演地層參數(shù)的橫向連續(xù)性和縱向分層性，特別適用于復雜地質(zhì)條件下的彈性參數(shù)反演和儲層參數(shù)預測。

上面結(jié)合附圖對本發(fā)明的實施方式作了詳細說明，但是本發(fā)明并不限于上述實施方式，在本領(lǐng)域普通技術(shù)人員所具備的知識范圍內(nèi)，還可以在不脫離本發(fā)明宗旨的前提下做出各種變化。