一種基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法
【專利摘要】本發(fā)明提供一種基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，該基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法包括：1）讀取野外觀測(cè)記錄及預(yù)設(shè)參數(shù)合成超道集；2）采用當(dāng)前反射系數(shù)模型，通過多震源激發(fā)正演模擬超道集，計(jì)算數(shù)據(jù)殘差；3）根據(jù)數(shù)據(jù)殘差計(jì)算更新梯度；4）通過隨機(jī)最優(yōu)化思想修改梯度并計(jì)算更新步長(zhǎng)；5）由梯度及更新步長(zhǎng)更新反射系數(shù)模型。本發(fā)明提供的方法提出了將隨機(jī)最優(yōu)化思想推廣到相位編碼的粘聲最小二乘逆時(shí)偏移中，通過加權(quán)平均之前的梯度，減小了梯度的隨機(jī)波動(dòng)，得到了較好的效果。
【專利說(shuō)明】
一種基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及油氣物探工程領(lǐng)域，特別是涉及到一種基于隨機(jī)最優(yōu)化的多震源粘聲 最小二乘逆時(shí)偏移方法。
【背景技術(shù)】
[0002] 實(shí)際地下介質(zhì)的粘滯性是普遍存在的，地震波在黏彈性介質(zhì)中的傳播主要表現(xiàn)為 速度頻散與振幅衰減。不考慮黏滯性的疊前成像算法不僅會(huì)使成像位置發(fā)生偏離，而且還 會(huì)引起成像振幅的欠估計(jì)，嚴(yán)重影響甚至誤導(dǎo)隨后的地震數(shù)據(jù)處理、解釋等工作。隨著油氣 勘探開發(fā)的深入，勘探精度要求逐漸提高，地震波成像也逐步從構(gòu)造成像向巖性成像發(fā)展。 然而，目前常規(guī)的偏移成像方法還無(wú)法滿足巖性油氣藏勘探開發(fā)的需求，究其原因是由于 常規(guī)偏移算子是正演算子的共輒轉(zhuǎn)置，而不是其逆算子。而基于反演思想的最小二乘偏移 相對(duì)于常規(guī)偏移來(lái)說(shuō)，具有更高的成像分辨率、振幅保真性和均衡性以及壓制偏移噪音等 優(yōu)勢(shì)，越來(lái)越受到學(xué)者的重視，然而計(jì)算量過于龐大限制了其進(jìn)一步推廣應(yīng)用。
[0003] 由于最小二乘逆時(shí)偏移(LSRTM)的計(jì)算量與炮數(shù)成線性關(guān)系，因而通過相位編碼 技術(shù)將多個(gè)炮集組合成一個(gè)超道集，可有效減小計(jì)算量。然而研究發(fā)現(xiàn)相位編碼算法的目 標(biāo)泛函是真實(shí)目標(biāo)泛函的隨機(jī)無(wú)偏估計(jì)，其梯度也是如此。由于相位編碼LSRTM的梯度是隨 機(jī)的，因而其步長(zhǎng)也應(yīng)該是隨機(jī)的。然而，目前相位編碼LSRTM求解時(shí)仍然采用與傳統(tǒng)LSRTM 算法相同的確定性最優(yōu)化解法，例如最速下降法、共輒梯度法等，忽略了梯度和步長(zhǎng)的隨機(jī) 性。為此我們發(fā)明了一種新的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，解決 了以上技術(shù)問題。

【發(fā)明內(nèi)容】

[0004]本發(fā)明的目的是提供一種考慮實(shí)際地下粘滯性的尚效率和尚精度的基于隨機(jī)最 優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法。
[0005] 本發(fā)明的目的可通過如下技術(shù)措施來(lái)實(shí)現(xiàn):基于隨機(jī)最優(yōu)化的多震源粘聲最小二 乘逆時(shí)偏移方法，該基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法包括：1)讀取野 外觀測(cè)記錄及預(yù)設(shè)參數(shù)合成超道集;2)采用當(dāng)前反射系數(shù)模型，通過多震源激發(fā)正演模擬 超道集，計(jì)算數(shù)據(jù)殘差;3)根據(jù)數(shù)據(jù)殘差計(jì)算更新梯度;4)通過隨機(jī)最優(yōu)化思想修改梯度并 計(jì)算更新步長(zhǎng);5)由梯度及更新步長(zhǎng)更新反射系數(shù)模型。
[0006] 本發(fā)明的目的還可通過如下技術(shù)措施來(lái)實(shí)現(xiàn)：
[0007] 在步驟1，輸入初始反射系數(shù)模型、偏移速度場(chǎng)、觀測(cè)數(shù)據(jù)、品質(zhì)因子、迭代終止的 閾值及偏移參數(shù)，初始反射系數(shù)模型的值為〇,即第1次迭代與常規(guī)粘聲逆時(shí)偏移等價(jià)。 [0008]在步驟2,計(jì)算數(shù)據(jù)殘差時(shí)，基于標(biāo)準(zhǔn)線性固體模型的擾動(dòng)波場(chǎng)的控制方程為：
[0009]
[0010] 其中，ps為擾動(dòng)波場(chǎng)，VQ為背景速度，p為密度，I為標(biāo)準(zhǔn)線性固體的個(gè)數(shù)，τε：?，Li為 松弛時(shí)間;H(t)為單位階躍函數(shù)，▽為梯度算子，▽?為散度算子，*為時(shí)間上的卷積算子，m (X)為模型參數(shù)，即反射系數(shù)模型，P〇為背景波場(chǎng)，即在背景介質(zhì)中傳播的波，其控制方程 為：
[0011]
[0012] 其中，f為經(jīng)過編碼的震源項(xiàng)，
[0013] 公式（1)可用矩陣算子形式表示為：
[0014] ps = Lm (3)
[0015] 其中，L為粘聲介質(zhì)擾動(dòng)波場(chǎng)的線性正演算子。
[0016] 在步驟2,在計(jì)算擾動(dòng)波場(chǎng)和背景波場(chǎng)時(shí)，松弛時(shí)間τ<3、τε的計(jì)算公式如公式(4)所 示：
[0017]
[0018]其中，w為圓頻率，Q為品質(zhì)因子。
[0019] 在步驟3,通過數(shù)據(jù)殘差計(jì)算更新模型的梯度方向，計(jì)算公式如公式(5)所示：
[0020] g = L*(Lm-p〇bs) (5)
[0021] 其中，g為梯度，Pobs為觀測(cè)記錄，為正演算子的共輒轉(zhuǎn)置，即逆時(shí)偏移算子，L為 擾動(dòng)波場(chǎng)的線性正演算子，m為反射系數(shù)模型。
[0022] 該基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法還包括，在步驟3之后，判 斷梯度是否滿足迭代終止條件，即梯度的模小于預(yù)設(shè)的閾值，若滿足，輸出當(dāng)前反射系數(shù)模 型，流程結(jié)束;否則流程進(jìn)入到步驟4。
[0023] 在步驟4,將隨機(jī)最優(yōu)化思想推廣到相位編碼的最小二乘逆時(shí)偏移算法中修改梯 度，隨機(jī)最優(yōu)化方法需要加權(quán)平均之前的梯度，因此在前幾次迭代時(shí)不需要修改梯度，修改 后的梯度如公式(6)所示，
[0024]
[0025] 其中，/為修改后的梯度，g為梯度，k為當(dāng)前迭代次數(shù);j為加權(quán)的前期迭代次數(shù)， 綜合考慮效果和效率，令其等于l〇;e為自然常數(shù)，a為衰減因子，取為0.4。
[0026] 在步驟4,由修改后的梯度計(jì)算更新步長(zhǎng)，如公式(7)所示，
[0027]
[0028]其中，ak為第k次迭代的更新步長(zhǎng)，gk為第k次迭代的梯度，L為擾動(dòng)波場(chǎng)的線性正演 算子。
[0029] 在步驟5,由公式(8)更新反射系數(shù)模型，
[0030]
[0031] 其中，Pk為預(yù)處理算子，mk為第k次迭代的反射系數(shù)模型，ak為第k次迭代的更新步 長(zhǎng)，g k為第k次迭代的梯度，采用背景波場(chǎng)的能量來(lái)近似Hessian矩陣的對(duì)角元素，在減少計(jì) 算量的同時(shí)加速了收斂速度，更新反射系數(shù)模型后流程返回到步驟2。
[0032] 本發(fā)明中的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，是一種有效提 高計(jì)算效率、快速壓制串?dāng)_噪聲的，且補(bǔ)償粘滯性的真振幅最小二乘逆時(shí)偏移方法。該方法 將實(shí)際地下的粘滯性引入到反演成像的最小二乘逆時(shí)偏移中，不僅能很好的避免常規(guī)粘聲 介質(zhì)成像的不穩(wěn)定性，且有效補(bǔ)償了由粘滯性造成的能量的吸收衰減?？紤]到基于常規(guī)相 位編碼的多震源最小二乘逆時(shí)偏移方法壓制串?dāng)_較慢，將隨機(jī)最優(yōu)化思想推廣到多震源最 小二乘逆時(shí)偏移中，有效提高了計(jì)算效率和收斂速度。
【附圖說(shuō)明】
[0033] 圖1為本發(fā)明的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法的一具體實(shí) 施例的流程圖；
[0034]圖2為本發(fā)明的一具體實(shí)施例中Marmousi模型速度場(chǎng)、擾動(dòng)模型及品質(zhì)因子模型 的不意圖；
[0035] 圖3為本發(fā)明的一具體實(shí)施例中基于隨機(jī)最優(yōu)化的多震源粘聲LSRTM反演成像結(jié) 果的不意圖；
[0036] 圖4為本發(fā)明的一具體實(shí)施例中數(shù)據(jù)殘差曲線的示意圖。
【具體實(shí)施方式】
[0037] 為使本發(fā)明的上述和其他目的、特征和優(yōu)點(diǎn)能更明顯易懂，下文特舉出較佳實(shí)施 例，并配合附圖所示，作詳細(xì)說(shuō)明如下。
[0038] 如圖1所示，圖1為本發(fā)明的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法 的流程圖。
[0039] 在步驟101，讀取野外觀測(cè)記錄及預(yù)設(shè)參數(shù)。輸入初始反射系數(shù)模型、偏移速度場(chǎng)、 觀測(cè)數(shù)據(jù)、品質(zhì)因子、迭代終止的閾值及偏移參數(shù)。初始模型的值為〇,即第1次迭代與常規(guī) 粘聲逆時(shí)偏移等價(jià)。流程進(jìn)入到步驟102。
[0040] 在步驟102,通過相位編碼合成超道集。流程進(jìn)入到步驟103。
[0041 ]在步驟103，采用當(dāng)前反射系數(shù)模型，通過多震源激發(fā)正演模擬擾動(dòng)波場(chǎng)，并計(jì)算 數(shù)據(jù)殘差，基于標(biāo)準(zhǔn)線性固體(GSLS)模型的擾動(dòng)波場(chǎng)的控制方程為：
[0042]
[0043] 其中，ps為擾動(dòng)波場(chǎng)，νο為背景速度，P為密度，I為標(biāo)準(zhǔn)線性固體的個(gè)數(shù)，τε?，hi為 松弛時(shí)間;H(t)為單位階躍函數(shù)，▽為梯度算子，▽?為散度算子，*為時(shí)間上的卷積算子，m (X)為模型參數(shù)，即反射系數(shù)模型。P〇為背景波場(chǎng)，即在背景介質(zhì)中傳播的波，其控制方程 為：
[0044]
[0045] 其中，f為經(jīng)過編碼的震源項(xiàng)。
[0046] 公式(1)可用矩陣算子形式表示為：
[0047] ps = Lm (3)
[0048] 在計(jì)算擾動(dòng)波場(chǎng)和背景波場(chǎng)時(shí)，松弛時(shí)間的計(jì)算公式如公式(4)所示：
[0049]
[0050]其中，w為圓頻率，Q為品質(zhì)因子。
[0051 ] 流程進(jìn)入到步驟104。
[0052]在步驟104,通過數(shù)據(jù)殘差計(jì)算更新模型的梯度方向，計(jì)算公式如公式(5)所示。
[0053] g = L*(Lm-p〇bs) (5)
[0054] 其中，g為梯度，P-為觀測(cè)記錄，為正演算子的共輒轉(zhuǎn)置，即逆時(shí)偏移算子。流程 進(jìn)入到步驟105。
[0055] 在步驟105，判斷梯度是否滿足迭代終止條件，即梯度的模小于預(yù)設(shè)的閾值，若滿 足，輸出當(dāng)前反射系數(shù)模型，流程進(jìn)入到步驟109;否則流程進(jìn)入到步驟106。
[0056] 步驟106,將隨機(jī)最優(yōu)化思想推廣到相位編碼的最小二乘逆時(shí)偏移算法中修改梯 度，隨機(jī)最優(yōu)化方法需要加權(quán)平均之前的梯度，因此在前幾次迭代時(shí)不需要修改梯度。修改 后的梯度如公式(6)所示，
[0057]
[0058] 其中，為修改后的梯度，k為當(dāng)前迭代次數(shù);j為加權(quán)的前期迭代次數(shù)，綜合考慮 效果和效率，令其等于l〇;e為自然常數(shù)， a為衰減因子，取為0.4。流程進(jìn)入到步驟107。
[0059] 步驟107,由修改后的梯度計(jì)算更新步長(zhǎng)，如公式(7)所示，
[0060]
[0061 ]其中，ak為第k次迭代的更新步長(zhǎng)，gk為第k次迭代的梯度。流程進(jìn)入到步驟108。 [0062]步驟108,由公式(8)更新反射系數(shù)模型，
[0063] mk+1=mk-akPkgk (8)
[0064] 其中，Pk為預(yù)處理算子，最優(yōu)的預(yù)處理算子應(yīng)為Hessian矩陣的逆，但計(jì)算量巨大， 直接求取并不現(xiàn)實(shí)，本發(fā)明中用背景波場(chǎng)的能量來(lái)近似Hessian矩陣的對(duì)角元素，在減少計(jì) 算量的同時(shí)加速了收斂速度。流程返回到步驟103。
[0065]步驟109,輸出最終偏移成像結(jié)果。
[0066] 實(shí)施例一
[0067]以Marmousi模型為例測(cè)試本發(fā)明所述方法，速度場(chǎng)如圖2(a)所示，對(duì)應(yīng)的擾動(dòng)模 型如圖2(b)所示，品質(zhì)因子模型如圖2(c)所示。從圖2可以看出，Marmousi模型中斷塊發(fā)育， 速度變化劇烈，可用于檢驗(yàn)成像算法的優(yōu)劣。模型參數(shù)如下:橫向網(wǎng)格點(diǎn)數(shù)為737，網(wǎng)格間距 12.5m，縱向網(wǎng)格點(diǎn)數(shù)為375，網(wǎng)格間距8m。計(jì)算參數(shù)為:在地表以50m間隔均勻分布185炮，每 炮都是737個(gè)檢波點(diǎn)全接收，震源為主頻20Hz的雷克子波，時(shí)間采樣間隔0.6ms。
[0068]圖3為成像測(cè)試結(jié)果，其中，圖3(a)為采用傳統(tǒng)極性編碼得到的多震源粘聲LSRTM 結(jié)果，圖3(b)為采用隨機(jī)最優(yōu)化極性編碼的多震源粘聲LSRTM成像結(jié)果，圖3(c)采用隨機(jī)最 優(yōu)化極性編碼的多震源聲波LSRTM成像結(jié)果。對(duì)比圖3(a)和圖3(b)可知，隨機(jī)最優(yōu)化算法相 比傳統(tǒng)算法成像質(zhì)量有了明顯改善，低頻噪音及串?dāng)_噪音更弱，分辨率略高。對(duì)比圖3(b)和 圖3(c)可以看出，由于考慮了地下介質(zhì)的衰減特性，本發(fā)明算法相比隨機(jī)最優(yōu)化多震源聲 波LSRTM成像振幅更加均衡，尤其是在深部地區(qū)的照明補(bǔ)償效果及分辨率更高。
[0069]相應(yīng)的三種算法的數(shù)據(jù)殘差曲線如圖4所示（圖中縱坐標(biāo)為對(duì)數(shù)坐標(biāo)），其中"Γ線 為采用傳統(tǒng)極性編碼的多震源粘聲LSRTM收斂曲線，"II"線為本發(fā)明算法，"III"線為采用 隨機(jī)最優(yōu)化極性編碼的多震源聲波LSRTM收斂曲線。從圖4可以清楚發(fā)現(xiàn)，本發(fā)明算法的誤 差收斂更快，在第30次迭代時(shí)的數(shù)據(jù)殘差已經(jīng)優(yōu)于其他兩種算法迭代第60次的數(shù)據(jù)殘差， 即采用本算法所需迭代次數(shù)更少，從而可顯著節(jié)約計(jì)算量、提高計(jì)算效率。
【主權(quán)項(xiàng)】
1. 基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特征在于，該基于隨機(jī)最 優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法包括： 1) 讀取野外觀測(cè)記錄及預(yù)設(shè)參數(shù)合成超道集； 2) 采用當(dāng)前反射系數(shù)模型，通過多震源激發(fā)正演模擬超道集，計(jì)算數(shù)據(jù)殘差； 3) 根據(jù)數(shù)據(jù)殘差計(jì)算更新梯度； 4) 通過隨機(jī)最優(yōu)化思想修改梯度并計(jì)算更新步長(zhǎng)； 5) 由梯度及更新步長(zhǎng)更新反射系數(shù)模型。2. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，在步驟1，輸入初始反射系數(shù)模型、偏移速度場(chǎng)、觀測(cè)數(shù)據(jù)、品質(zhì)因子、迭代終止的闊 值及偏移參數(shù)，初始反射系數(shù)模型的值為0,即第1次迭代與常規(guī)粘聲逆時(shí)偏移等價(jià)。3. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，在步驟2,計(jì)算數(shù)據(jù)殘差時(shí)，基于標(biāo)準(zhǔn)線性固體模型的擾動(dòng)波場(chǎng)的控制方程為：其中，Ps為擾動(dòng)波場(chǎng)，VO為背景速度，P為密度，I為標(biāo)準(zhǔn)線性固體的個(gè)數(shù)，Tei，Toi為松弛 時(shí)間;H(t)為單位階躍函數(shù)，▽為梯度算子，▽?為散度算子，*為時(shí)間上的卷積算子，m(x) 為模型參數(shù)，即反射系數(shù)模型，PO為背景波場(chǎng)，即在背景介質(zhì)中傳播的波，其控制方程為：其中，f為經(jīng)過編碼的震源項(xiàng)， 公式(1 )可用矩陣算子形式表示為： Ps = Lm (3) 其中，L為粘聲介質(zhì)擾動(dòng)波場(chǎng)的線性正演算子。4. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，在步驟2,在計(jì)算擾動(dòng)波場(chǎng)和背景波場(chǎng)時(shí)，松弛時(shí)間Tn、TE的計(jì)算公式如公式(4)所 示：其中，W為圓頻率，Q為品質(zhì)因子。5. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，在步驟3,通過數(shù)據(jù)殘差計(jì)算更新模型的梯度方向，計(jì)算公式如公式(5)所示： g = L*(Lm-p〇bs) (5) 其中，g為梯度，Pebs為觀測(cè)記錄，正演算子的共輛轉(zhuǎn)置，即逆時(shí)偏移算子，L為擾動(dòng)波 場(chǎng)的線性正演算子，m為反射系數(shù)模型。6. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，該基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法還包括，在步驟3之后，判 斷梯度是否滿足迭代終止條件，即梯度的模小于預(yù)設(shè)的闊值，若滿足，輸出當(dāng)前反射系數(shù)模 型，流程結(jié)束;否則流程進(jìn)入到步驟4。7. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，在步驟4,將隨機(jī)最優(yōu)化思想推廣到相位編碼的最小二乘逆時(shí)偏移算法中修改梯 度，隨機(jī)最優(yōu)化方法需要加權(quán)平均之前的梯度，因此在前幾次迭代時(shí)不需要修改梯度，修改 后的梯度如公式(6)所示，其中，/為修改后的梯度，g為梯度，k為當(dāng)前迭代次數(shù)；j為加權(quán)的前期迭代次數(shù)，綜合 考慮效果和效率，令其等于l〇;e為自然常數(shù)，a為衰減因子，取為0.4。8. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，在步驟4,由修改后的梯度計(jì)算更新步長(zhǎng)，如公式(7)所示，其中，Qk為第k次迭代的更新步長(zhǎng)，gk為第k次迭代的梯度，L為擾動(dòng)波場(chǎng)的線性正演算 子。9. 根據(jù)權(quán)利要求1所述的基于隨機(jī)最優(yōu)化的多震源粘聲最小二乘逆時(shí)偏移方法，其特 征在于，在步驟5，由公式(8)更新反射系數(shù)模型，其中，pk為預(yù)處理算子，mk為第k次迭代的反射系數(shù)模型，Qk為第k次迭代的更新步長(zhǎng)，gk 為第k次迭代的梯度，采用背景波場(chǎng)的能量來(lái)近似化SSian矩陣的對(duì)角元素，在減少計(jì)算量 的同時(shí)加速了收斂速度，更新反射系數(shù)模型后流程返回到步驟2。
【文檔編號(hào)】G01V1/28GK106033124SQ201610494147
【公開日】2016年10月19日
【申請(qǐng)日】2016年6月29日
【發(fā)明人】張猛, 匡斌, 單聯(lián)瑜, 趙慶國(guó), 王慧, 隆文韜, 高麗, 王修銀, 廉西猛, 王賢真, 李振春, 張傳強(qiáng), 陳震林
【申請(qǐng)人】中國(guó)石油化工股份有限公司, 中國(guó)石油化工股份有限公司勝利油田分公司物探研究院