本發(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)度表達式為：

其中，為第i個粒子適應(yīng)度，pi為第i個粒子位置，e(pi)為pi的期望，l為觀測數(shù)據(jù)采樣長度，θ為入射角，t為地震記錄采樣時間，w為震源子波，r為縱波反射系數(shù)，d為觀測地震記錄，σ為縱波速度vpi、橫波速度vsi和密度ρi這三個參數(shù)的協(xié)方差矩陣，λ1、λ2均為預(yù)設(shè)定系數(shù)，上標t表示矩陣的轉(zhuǎn)置。

作為本發(fā)明的一種優(yōu)選方案，步驟4所述接收概率表達式為：

其中，piaccept為第i個粒子的接收概率，為第i個粒子適應(yīng)度，n為粒子的數(shù)量，tk為第k次迭代的退火溫度。

作為本發(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＝[pi¹,pi²,pi³,pi⁴]＝[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)度值并計算各粒子位置被選擇作為模擬退火學(xué)習(xí)項的引導(dǎo)粒子pmin的接收概率，接收概率的表達式為：

其中，piaccpet為位置pi被接收的概率，tk為第k次迭代的退火溫度；為最優(yōu)粒子的適應(yīng)度值。

基于各粒子的接收概率，使用輪盤選擇算法決定模擬退火學(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)定性。

步驟六，更新粒子速度和粒子位置。

粒子速度的更新公式為：

其中，為模擬退火學(xué)習(xí)項，其中的引導(dǎo)粒子pmin由步驟四得到；為多維學(xué)習(xí)項，其中的引導(dǎo)粒子pmax由步驟五得到，上標d代表粒子的維度，每個粒子pi包含四個維度(即縱波速度、橫波速度、密度三參數(shù)，以及三參數(shù)的聯(liá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)。