本發(fā)明涉及一種混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法,具體涉及一種結(jié)合多維學(xué)習(xí)粒子群算法和快速模擬退火算法的混合全局優(yōu)化算法,屬于油氣地球物理勘探中的地震資料反演技術(shù)領(lǐng)域。
背景技術(shù):
疊前地震反演是一種基于褶積模型的反演,它基于振幅隨炮檢距變化理論,可以同步獲得縱波速度、橫波速度和密度等多種彈性參數(shù),可以有效識別儲層流體類型和地質(zhì)構(gòu)造特征。因此,疊前地震反演是目前應(yīng)用最廣泛、發(fā)展最成熟的地震反演技術(shù)之一,在油氣資源的勘探和開發(fā)過程中發(fā)揮重要作用。
地震反演是一個非線性優(yōu)化問題,即使得估算的地層參數(shù)所對應(yīng)的合成地震記錄與觀測地震記錄間的誤差達到最小。目前,地震反演的最優(yōu)化算法可以分為兩大類:第一類是局部梯度類算法,如最速下降法、共軛梯度法和牛頓類方法等;第二類是全局優(yōu)化算法,如模擬退火算法、遺傳算法和粒子群算法等。梯度類算法計算效率較高,對非線性反演問題可以快速收斂獲得最優(yōu)解,但此類算法對初始模型依賴較高,并且對復(fù)雜的非線性反演問題效果不理想。而全局最優(yōu)化方法可以解決非線性和多峰值反演問題,并且不依賴于初始模型,特別適用于疊前地震多參數(shù)反演這類復(fù)雜的非線性反演問題。
模擬退火算法和粒子群算法是目前發(fā)展較成熟的全局優(yōu)化算法,已有不少應(yīng)用于疊前地震反演的成功案例。但模擬退火算法計算量較大,收斂速度較慢,雖然粒子群算法運算效率有顯著提高,但易不成熟收斂而陷入局部最小值。因此,結(jié)合粒子群算法和模擬退火算法的混合全局最優(yōu)化算法對兩類算法進行優(yōu)勢互補,使得粒子群算法具有模擬退火算法的“跳躍”機制,擴大了粒子搜索和擴展的能力,并且不易陷落局部最優(yōu)解;此外,疊前地震反演是一種多參數(shù)的非線性反演問題,多參數(shù)反演結(jié)果具有不穩(wěn)定性,而常規(guī)反演算法本身并未有針對該問題的解決方法,考慮到多參數(shù)間存在基于地質(zhì)統(tǒng)計學(xué)的約束關(guān)系,在粒子群算法中添加多維學(xué)習(xí)項,可以有效解決疊前地震反演多參數(shù)的不穩(wěn)定問題。
技術(shù)實現(xiàn)要素:
本發(fā)明所要解決的技術(shù)問題是:提供一種混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法,克服疊前地震多參數(shù)反演的不穩(wěn)定問題,并解決傳統(tǒng)粒子群算法收斂不成熟的問題,可以同步而準確的獲取地震三參數(shù)的反演結(jié)果。
本發(fā)明為解決上述技術(shù)問題采用以下技術(shù)方案:
一種混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法,包括如下步驟:
步驟1,確定粒子待反演參數(shù)的取值范圍,設(shè)定粒子的數(shù)量、最大迭代次數(shù)、模擬退火冷卻進度表、模擬退火學(xué)習(xí)項加權(quán)系數(shù)和多維學(xué)習(xí)項加權(quán)系數(shù);
步驟2,設(shè)置各個粒子的初始位置和初始速度,位置包括四個維度:縱波速度、橫波速度、密度、前述三個參數(shù)的聯(lián)合概率密度;
步驟3,對第k次迭代,根據(jù)粒子的位置計算粒子適應(yīng)度,其中粒子適應(yīng)度由觀測地震記錄與合成地震記錄的誤差和步驟2所述三個參數(shù)的先驗約束項構(gòu)成;
步驟4,根據(jù)步驟3的粒子適應(yīng)度,計算各粒子位置被選擇作為模擬退火學(xué)習(xí)項引導(dǎo)粒子的接收概率,并基于各粒子的接收概率,通過輪盤選擇法則確定模擬退火學(xué)習(xí)項引導(dǎo)粒子;
步驟5,對每個粒子,固定該粒子的橫波速度和密度,用所有粒子的縱波速度依次替換該粒子的縱波速度,并利用聯(lián)合概率密度公式計算所有粒子的縱波速度對應(yīng)的聯(lián)合概率密度,找出最大聯(lián)合概率密度;從每個粒子對應(yīng)的最大聯(lián)合概率密度中找出最大值及該最大值對應(yīng)的縱波速度vpmax;按上述同樣的方法,找到橫波速度vsmax和密度ρmax;根據(jù)vpmax、vsmax和ρmax,得到多維學(xué)習(xí)項引導(dǎo)粒子;
步驟6,根據(jù)模擬退火學(xué)習(xí)項引導(dǎo)粒子和多維學(xué)習(xí)項引導(dǎo)粒子更新各粒子的速度和位置;
步驟7,進入第k+1次迭代,重復(fù)步驟3至步驟6,直至達到最大迭代次數(shù),且到達模擬退火冷卻進度表中的終止退火溫度,輸出粒子適應(yīng)度最優(yōu)的粒子。
作為本發(fā)明的一種優(yōu)選方案,步驟1所述待反演參數(shù)包括縱波速度、橫波速度和密度。
作為本發(fā)明的一種優(yōu)選方案,步驟1所述模擬退火冷卻進度表包括初始退火溫度、終止退火溫度,溫度從初始退火溫度開始逐漸降低直至終止退火溫度,所有退火溫度的個數(shù)與最大迭代次數(shù)相同,且一次迭代對應(yīng)一個退火溫度。
作為本發(fā)明的一種優(yōu)選方案,步驟2所述聯(lián)合概率密度表達式為:
其中,fi為第i個粒子縱波速度vpi、橫波速度vsi和密度ρi這三個參數(shù)的聯(lián)合概率密度,σ為縱波速度vpi、橫波速度vsi和密度ρi這三個參數(shù)的協(xié)方差矩陣,pi為第i個粒子位置,e(pi)為pi的期望,上標t表示矩陣的轉(zhuǎn)置。
作為本發(fā)明的一種優(yōu)選方案,步驟3所述粒子適應(yīng)度表達式為:
其中,
作為本發(fā)明的一種優(yōu)選方案,步驟4所述接收概率表達式為:
其中,piaccept為第i個粒子的接收概率,
作為本發(fā)明的一種優(yōu)選方案,所述步驟6更新公式為:
粒子速度更新公式:
粒子位置更新公式:
其中,d表示粒子的維度,vi為第i個粒子的速度,k、k-1分別為第k、k-1次迭代,c1、c2分別為模擬退火學(xué)習(xí)項加權(quán)系數(shù)、多維學(xué)習(xí)項加權(quán)系數(shù),rand為[0,1]的隨機函數(shù),pmin為模擬退火學(xué)習(xí)項引導(dǎo)粒子,pmax為多維學(xué)習(xí)項引導(dǎo)粒子,pi為第i個粒子位置。
本發(fā)明采用以上技術(shù)方案與現(xiàn)有技術(shù)相比,具有以下技術(shù)效果:
1、本發(fā)明方法將粒子群算法和快速模擬退火算法有效結(jié)合,解決傳統(tǒng)粒子群算法易不成熟收斂的問題,并在粒子群算法中添加基于三參數(shù)聯(lián)合概率密度擇優(yōu)組合的多維學(xué)習(xí)項,克服疊前地震多參數(shù)同步反演的不穩(wěn)定,可以同步而準確的獲取縱波速度、橫波速度和密度三參數(shù)反演結(jié)果。
2、本發(fā)明方法與傳統(tǒng)粒子群算法的反演結(jié)果相比反演效果改善明顯。
附圖說明
圖1是本發(fā)明混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法的流程圖。
圖2是本發(fā)明各反演參數(shù)的理論模型,其中,(a)為縱波速度,(b)為橫波速度,(c)為密度。
圖3是理論模型的合成地震數(shù)據(jù)(觀測地震數(shù)據(jù))。
圖4是使用本發(fā)明方法的各反演參數(shù)反演結(jié)果與理論模型的對比,其中,(a)為縱波速度,(b)為橫波速度,(c)為密度。
圖5是使用傳統(tǒng)粒子群算法的各反演參數(shù)反演結(jié)果與理論模型的對比,其中,(a)為縱波速度,(b)為橫波速度,(c)為密度。
圖6是使用本發(fā)明方法和傳統(tǒng)粒子群算法進行反演的收斂速度的對比。
具體實施方式
下面詳細描述本發(fā)明的實施方式,所述實施方式的示例在附圖中示出。下面通過參考附圖描述的實施方式是示例性的,僅用于解釋本發(fā)明,而不能解釋為對本發(fā)明的限制。
如圖1所示,為本發(fā)明反演方法的流程圖,具體步驟如下:
步驟一,初始化設(shè)置。待反演模型離散化,設(shè)定各位置待反演的縱波速度、橫波速度和密度的取值范圍;設(shè)定粒子的數(shù)量n、最大迭代次數(shù)k和模擬退火冷卻進度表;設(shè)置學(xué)習(xí)項加權(quán)系數(shù)c1和c2。
加權(quán)系數(shù)c1和c2分別對應(yīng)粒子速度更新的模擬退火學(xué)習(xí)項和多維學(xué)習(xí)項。對這兩種學(xué)習(xí)項的選擇通過加權(quán)系數(shù)c1和c2進行平衡。模擬退火學(xué)習(xí)項結(jié)合快速模擬退火算法,模擬退火算法的“跳躍”機制使得這種學(xué)習(xí)有一定的概率接受非最優(yōu)位置的粒子作為引導(dǎo)粒子,從而擴大粒子的更新范圍和搜索能力,不易陷入局部最優(yōu)解;多維學(xué)習(xí)項是粒子向其他維度粒子位置的學(xué)習(xí),多維學(xué)習(xí)項的引導(dǎo)粒子是基于三參數(shù)間的統(tǒng)計學(xué)關(guān)系而選擇的具有最大概率密度值的組合粒子,可以有效控制三參數(shù)同步反演的穩(wěn)定性。
模擬退火冷卻進度表包括初始溫度t1和終止溫度tend,設(shè)定溫度個數(shù)與最大迭代次數(shù)相同,即每個迭代k與一個溫度tk對應(yīng),隨著迭代進行,溫度逐漸降低直至終止溫度。溫度控制模擬退火學(xué)習(xí)項中各粒子位置被選擇作為引導(dǎo)粒子的接收概率。
步驟二,設(shè)置粒子初始位置和初始速度。每個粒子位置pi包含四個維度(即縱波速度、橫波速度和密度三參數(shù),以及三參數(shù)的聯(lián)合概率密度),即pi=[pi1,pi2,pi3,pi4]=[vpi,vsi,ρi,fi],i=1,2,…,n。
其中,三參數(shù)的聯(lián)合概率密度的具體表達式為:
其中,e為期望;σ為縱波速度、橫波速度和密度三參數(shù)的協(xié)方差矩陣,由測井數(shù)據(jù)獲取。
在取值范圍內(nèi),隨機設(shè)置粒子的初始位置和初始速度。粒子的初始位置組成[4×n]的矩陣,粒子的初始速度設(shè)為零。
步驟三,設(shè)置目標函數(shù),即粒子適應(yīng)度表達式。本反演方法的目標函數(shù)由觀測地震記錄與合成地震記錄的誤差和三參數(shù)的先驗約束項構(gòu)成,具體表現(xiàn)形式如下:
其中,f為粒子適應(yīng)度;r為使用精確佐普里茲方程計算的縱波反射系數(shù);w為震源子波;d為觀測地震記錄;l為觀測數(shù)據(jù)采樣長度;λ1和λ2為預(yù)設(shè)定系數(shù);θ為入射角;t為地震記錄采樣時間。
步驟四,設(shè)置模擬退火學(xué)習(xí)項。按公式(2)計算各粒子pi的適應(yīng)度值
其中,piaccpet為位置pi被接收的概率,tk為第k次迭代的退火溫度;
基于各粒子的接收概率,使用輪盤選擇算法決定模擬退火學(xué)習(xí)項的引導(dǎo)粒子pmin,顯然,適應(yīng)度最小的粒子具有最大的接收概率,但非最優(yōu)粒子也有一定概率“跳躍”成為引導(dǎo)粒子,且這種概率在高溫時最大,并隨著溫度降低而減小。這種“跳躍”機制增加了粒子的搜索擴展能力,避免因不成熟收斂而陷入局部最優(yōu)解。
步驟五,設(shè)置多維學(xué)習(xí)項。
固定橫波速度和密度,由公式(1)計算各粒子的縱波速度對應(yīng)的聯(lián)合概率密度值,取最大概率密度值對應(yīng)的縱波速度vpmax;按同樣方法,分別獲得最大概率密度對應(yīng)的橫波速度vsmax和密度ρmax,即:
多維學(xué)習(xí)項的引導(dǎo)粒子pmax是具有最大概率密度值的擇優(yōu)組合粒子,即
pmax=[vpmax,vsmax,ρmax,fmax]
多維學(xué)習(xí)項的引導(dǎo)粒子是基于三參數(shù)間的統(tǒng)計學(xué)關(guān)系而選擇的具有最大概率密度值的組合粒子,該引導(dǎo)粒子不局限于單一粒子所對應(yīng)的位置,而是多維度、多粒子間擇優(yōu)組合的結(jié)果,可以有效控制疊前地震三參數(shù)同步反演的穩(wěn)定性。
步驟六,更新粒子速度和粒子位置。
粒子速度的更新公式為:
其中,
粒子位置的更新公式為:
步驟七,進入下一次迭代k+1,并使溫度降低至tk+1,重復(fù)步驟三至步驟六。直至達到最大迭代數(shù)k,溫度降至終止溫度tend,輸出適應(yīng)度最優(yōu)的粒子。
下面以一個合成地震數(shù)據(jù)測試進行具體說明:
合成地震數(shù)據(jù)測試的理論模型如圖2所示,其中,(a)為縱波速度,(b)為橫波速度,(c)為密度。使用精確佐普里茲方程計算理論模型中地層的縱波反射系數(shù),計算的角度間隔為5°,角度覆蓋范圍為0-50°。分別將每層的11組反射系數(shù)與主頻為50hz的零相位理論ricker子波進行褶積,得到如圖3所示的合成地震記錄(角度道集),即觀測地震數(shù)據(jù),其中有11個角度道,時間采樣間隔2ms。
使用本發(fā)明方法對觀測地震數(shù)據(jù)進行反演,以獲得縱波速度、橫波速度和密度三參數(shù)的反演結(jié)果,使反演結(jié)果的合成地震數(shù)據(jù)與觀測地震數(shù)據(jù)之間誤差的二次范數(shù)極小,具體實現(xiàn)方式如下:
設(shè)置粒子位置的取值范圍,每個粒子位置包含四個維度,即縱波速度、橫波速度、密度和聯(lián)合概率密度,對應(yīng)的最小值為[1.2,0.6,1.2,0]和最大值為[5.6,3.5,3.2,1],粒子數(shù)量n為20,最大迭代次數(shù)k為300,初始溫度t1為0.9,終止溫度tend為0.003,學(xué)習(xí)項加權(quán)系數(shù)c1=c2=1。
在取值范圍內(nèi),隨機生成粒子的初始位置,并對粒子位置進行對數(shù)歸一化處理;設(shè)定各粒子的初始速度為零。
按步驟三計算各粒子的適應(yīng)度值,即各粒子對應(yīng)參數(shù)模型的合成地震記錄與觀測地震數(shù)據(jù)的誤差的二次范數(shù)。
按步驟四計算模擬退火學(xué)習(xí)項,基于各粒子的接收概率,按輪盤選擇算法確定該項的引導(dǎo)粒子pmin;按步驟五計算多維學(xué)習(xí)項,擇優(yōu)組合粒子位置,使其具有最大聯(lián)合概率密度值,構(gòu)成該項的引導(dǎo)粒子pmax,這里使用三參數(shù)模型的平滑結(jié)果作為測井數(shù)據(jù),計算三參數(shù)的協(xié)方差矩陣。
對各粒子按步驟六進行速度更新和位置更新。
降低溫度并進入下一次迭代,利用本發(fā)明方法對粒子不斷迭代更新,直至終止溫度和最大迭代次數(shù),輸出最優(yōu)適應(yīng)度的粒子,作為觀測地震數(shù)據(jù)的反演結(jié)果。
圖4為使用本發(fā)明方法的三參數(shù)反演結(jié)果與理論模型的對比,其中,(a)為縱波速度,(b)為橫波速度,(c)為密度;可以看到反演結(jié)果(虛線)與理論結(jié)果(實線)吻合很好,特別是一些高頻成分得到很好的恢復(fù);圖5為使用傳統(tǒng)粒子群算法的三參數(shù)反演結(jié)果與理論模型的對比,其中,(a)為縱波速度,(b)為橫波速度,(c)為密度;可以看到反演結(jié)果(虛線)與理論結(jié)果(實線)吻合較差。圖6為分別使用傳統(tǒng)粒子群算法與本發(fā)明方法進行反演的收斂速度的對比,可以看到使用本算法的收斂速度(實線)比使用傳統(tǒng)粒子群算法的收斂速度(虛線)明顯較快。
以上實施例僅為說明本發(fā)明的技術(shù)思想,不能以此限定本發(fā)明的保護范圍,凡是按照本發(fā)明提出的技術(shù)思想,在技術(shù)方案基礎(chǔ)上所做的任何改動,均落入本發(fā)明保護范圍之內(nèi)。