本發(fā)明涉及一種混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法，具體涉及一種結(jié)合多維學(xué)習(xí)粒子群算法和快速模擬退火算法的混合全局優(yōu)化算法，屬于油氣地球物理勘探中的地震資料反演技術(shù)領(lǐng)域。

背景技術(shù)：

疊前地震反演是一種基于褶積模型的反演，它基于振幅隨炮檢距變化理論，可以同步獲得縱波速度、橫波速度和密度等多種彈性參數(shù)，可以有效識(shí)別儲(chǔ)層流體類型和地質(zhì)構(gòu)造特征。因此，疊前地震反演是目前應(yīng)用最廣泛、發(fā)展最成熟的地震反演技術(shù)之一，在油氣資源的勘探和開(kāi)發(fā)過(guò)程中發(fā)揮重要作用。

地震反演是一個(gè)非線性優(yōu)化問(wèn)題，即使得估算的地層參數(shù)所對(duì)應(yīng)的合成地震記錄與觀測(cè)地震記錄間的誤差達(dá)到最小。目前，地震反演的最優(yōu)化算法可以分為兩大類：第一類是局部梯度類算法，如最速下降法、共軛梯度法和牛頓類方法等；第二類是全局優(yōu)化算法，如模擬退火算法、遺傳算法和粒子群算法等。梯度類算法計(jì)算效率較高，對(duì)非線性反演問(wèn)題可以快速收斂獲得最優(yōu)解，但此類算法對(duì)初始模型依賴較高，并且對(duì)復(fù)雜的非線性反演問(wèn)題效果不理想。而全局最優(yōu)化方法可以解決非線性和多峰值反演問(wèn)題，并且不依賴于初始模型，特別適用于疊前地震多參數(shù)反演這類復(fù)雜的非線性反演問(wèn)題。

模擬退火算法和粒子群算法是目前發(fā)展較成熟的全局優(yōu)化算法，已有不少應(yīng)用于疊前地震反演的成功案例。但模擬退火算法計(jì)算量較大，收斂速度較慢，雖然粒子群算法運(yùn)算效率有顯著提高，但易不成熟收斂而陷入局部最小值。因此，結(jié)合粒子群算法和模擬退火算法的混合全局最優(yōu)化算法對(duì)兩類算法進(jìn)行優(yōu)勢(shì)互補(bǔ)，使得粒子群算法具有模擬退火算法的“跳躍”機(jī)制，擴(kuò)大了粒子搜索和擴(kuò)展的能力，并且不易陷落局部最優(yōu)解；此外，疊前地震反演是一種多參數(shù)的非線性反演問(wèn)題，多參數(shù)反演結(jié)果具有不穩(wěn)定性，而常規(guī)反演算法本身并未有針對(duì)該問(wèn)題的解決方法，考慮到多參數(shù)間存在基于地質(zhì)統(tǒng)計(jì)學(xué)的約束關(guān)系，在粒子群算法中添加多維學(xué)習(xí)項(xiàng)，可以有效解決疊前地震反演多參數(shù)的不穩(wěn)定問(wèn)題。

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

本發(fā)明所要解決的技術(shù)問(wèn)題是：提供一種混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法，克服疊前地震多參數(shù)反演的不穩(wěn)定問(wèn)題，并解決傳統(tǒng)粒子群算法收斂不成熟的問(wèn)題，可以同步而準(zhǔn)確的獲取地震三參數(shù)的反演結(jié)果。

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

一種混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法，包括如下步驟：

步驟1，確定粒子待反演參數(shù)的取值范圍，設(shè)定粒子的數(shù)量、最大迭代次數(shù)、模擬退火冷卻進(jìn)度表、模擬退火學(xué)習(xí)項(xiàng)加權(quán)系數(shù)和多維學(xué)習(xí)項(xiàng)加權(quán)系數(shù)；

步驟2，設(shè)置各個(gè)粒子的初始位置和初始速度，位置包括四個(gè)維度：縱波速度、橫波速度、密度、前述三個(gè)參數(shù)的聯(lián)合概率密度；

步驟3，對(duì)第k次迭代，根據(jù)粒子的位置計(jì)算粒子適應(yīng)度，其中粒子適應(yīng)度由觀測(cè)地震記錄與合成地震記錄的誤差和步驟2所述三個(gè)參數(shù)的先驗(yàn)約束項(xiàng)構(gòu)成；

步驟4，根據(jù)步驟3的粒子適應(yīng)度，計(jì)算各粒子位置被選擇作為模擬退火學(xué)習(xí)項(xiàng)引導(dǎo)粒子的接收概率，并基于各粒子的接收概率，通過(guò)輪盤選擇法則確定模擬退火學(xué)習(xí)項(xiàng)引導(dǎo)粒子；

步驟5，對(duì)每個(gè)粒子，固定該粒子的橫波速度和密度，用所有粒子的縱波速度依次替換該粒子的縱波速度，并利用聯(lián)合概率密度公式計(jì)算所有粒子的縱波速度對(duì)應(yīng)的聯(lián)合概率密度，找出最大聯(lián)合概率密度；從每個(gè)粒子對(duì)應(yīng)的最大聯(lián)合概率密度中找出最大值及該最大值對(duì)應(yīng)的縱波速度vpmax；按上述同樣的方法，找到橫波速度vsmax和密度ρmax；根據(jù)vpmax、vsmax和ρmax，得到多維學(xué)習(xí)項(xiàng)引導(dǎo)粒子；

步驟6，根據(jù)模擬退火學(xué)習(xí)項(xiàng)引導(dǎo)粒子和多維學(xué)習(xí)項(xiàng)引導(dǎo)粒子更新各粒子的速度和位置；

步驟7，進(jìn)入第k+1次迭代，重復(fù)步驟3至步驟6，直至達(dá)到最大迭代次數(shù)，且到達(dá)模擬退火冷卻進(jìn)度表中的終止退火溫度，輸出粒子適應(yīng)度最優(yōu)的粒子。

作為本發(fā)明的一種優(yōu)選方案，步驟1所述待反演參數(shù)包括縱波速度、橫波速度和密度。

作為本發(fā)明的一種優(yōu)選方案，步驟1所述模擬退火冷卻進(jìn)度表包括初始退火溫度、終止退火溫度，溫度從初始退火溫度開(kāi)始逐漸降低直至終止退火溫度，所有退火溫度的個(gè)數(shù)與最大迭代次數(shù)相同，且一次迭代對(duì)應(yīng)一個(gè)退火溫度。

作為本發(fā)明的一種優(yōu)選方案，步驟2所述聯(lián)合概率密度表達(dá)式為：

其中，fi為第i個(gè)粒子縱波速度vpi、橫波速度vsi和密度ρi這三個(gè)參數(shù)的聯(lián)合概率密度，σ為縱波速度vpi、橫波速度vsi和密度ρi這三個(gè)參數(shù)的協(xié)方差矩陣，pi為第i個(gè)粒子位置，e(pi)為pi的期望，上標(biāo)t表示矩陣的轉(zhuǎn)置。

作為本發(fā)明的一種優(yōu)選方案，步驟3所述粒子適應(yīng)度表達(dá)式為：

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

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

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

作為本發(fā)明的一種優(yōu)選方案，所述步驟6更新公式為：

粒子速度更新公式：

粒子位置更新公式：

其中，d表示粒子的維度，vi為第i個(gè)粒子的速度，k、k-1分別為第k、k-1次迭代，c1、c2分別為模擬退火學(xué)習(xí)項(xiàng)加權(quán)系數(shù)、多維學(xué)習(xí)項(xiàng)加權(quán)系數(shù)，rand為[0,1]的隨機(jī)函數(shù)，pmin為模擬退火學(xué)習(xí)項(xiàng)引導(dǎo)粒子，pmax為多維學(xué)習(xí)項(xiàng)引導(dǎo)粒子，pi為第i個(gè)粒子位置。

本發(fā)明采用以上技術(shù)方案與現(xiàn)有技術(shù)相比，具有以下技術(shù)效果：

1、本發(fā)明方法將粒子群算法和快速模擬退火算法有效結(jié)合，解決傳統(tǒng)粒子群算法易不成熟收斂的問(wèn)題，并在粒子群算法中添加基于三參數(shù)聯(lián)合概率密度擇優(yōu)組合的多維學(xué)習(xí)項(xiàng)，克服疊前地震多參數(shù)同步反演的不穩(wěn)定，可以同步而準(zhǔn)確的獲取縱波速度、橫波速度和密度三參數(shù)反演結(jié)果。

2、本發(fā)明方法與傳統(tǒng)粒子群算法的反演結(jié)果相比反演效果改善明顯。

附圖說(shuō)明

圖1是本發(fā)明混合全局優(yōu)化算法的疊前地震多參數(shù)反演方法的流程圖。

圖2是本發(fā)明各反演參數(shù)的理論模型，其中，(a)為縱波速度，(b)為橫波速度，(c)為密度。

圖3是理論模型的合成地震數(shù)據(jù)(觀測(cè)地震數(shù)據(jù))。

圖4是使用本發(fā)明方法的各反演參數(shù)反演結(jié)果與理論模型的對(duì)比，其中，(a)為縱波速度，(b)為橫波速度，(c)為密度。

圖5是使用傳統(tǒng)粒子群算法的各反演參數(shù)反演結(jié)果與理論模型的對(duì)比，其中，(a)為縱波速度，(b)為橫波速度，(c)為密度。

圖6是使用本發(fā)明方法和傳統(tǒng)粒子群算法進(jìn)行反演的收斂速度的對(duì)比。

具體實(shí)施方式

下面詳細(xì)描述本發(fā)明的實(shí)施方式，所述實(shí)施方式的示例在附圖中示出。下面通過(guò)參考附圖描述的實(shí)施方式是示例性的，僅用于解釋本發(fā)明，而不能解釋為對(duì)本發(fā)明的限制。

如圖1所示，為本發(fā)明反演方法的流程圖，具體步驟如下：

步驟一，初始化設(shè)置。待反演模型離散化，設(shè)定各位置待反演的縱波速度、橫波速度和密度的取值范圍；設(shè)定粒子的數(shù)量n、最大迭代次數(shù)k和模擬退火冷卻進(jìn)度表；設(shè)置學(xué)習(xí)項(xiàng)加權(quán)系數(shù)c1和c2。

加權(quán)系數(shù)c1和c2分別對(duì)應(yīng)粒子速度更新的模擬退火學(xué)習(xí)項(xiàng)和多維學(xué)習(xí)項(xiàng)。對(duì)這兩種學(xué)習(xí)項(xiàng)的選擇通過(guò)加權(quán)系數(shù)c1和c2進(jìn)行平衡。模擬退火學(xué)習(xí)項(xiàng)結(jié)合快速模擬退火算法，模擬退火算法的“跳躍”機(jī)制使得這種學(xué)習(xí)有一定的概率接受非最優(yōu)位置的粒子作為引導(dǎo)粒子，從而擴(kuò)大粒子的更新范圍和搜索能力，不易陷入局部最優(yōu)解；多維學(xué)習(xí)項(xiàng)是粒子向其他維度粒子位置的學(xué)習(xí)，多維學(xué)習(xí)項(xiàng)的引導(dǎo)粒子是基于三參數(shù)間的統(tǒng)計(jì)學(xué)關(guān)系而選擇的具有最大概率密度值的組合粒子，可以有效控制三參數(shù)同步反演的穩(wěn)定性。

模擬退火冷卻進(jìn)度表包括初始溫度t1和終止溫度tend，設(shè)定溫度個(gè)數(shù)與最大迭代次數(shù)相同，即每個(gè)迭代k與一個(gè)溫度tk對(duì)應(yīng)，隨著迭代進(jìn)行，溫度逐漸降低直至終止溫度。溫度控制模擬退火學(xué)習(xí)項(xiàng)中各粒子位置被選擇作為引導(dǎo)粒子的接收概率。

步驟二，設(shè)置粒子初始位置和初始速度。每個(gè)粒子位置pi包含四個(gè)維度(即縱波速度、橫波速度和密度三參數(shù)，以及三參數(shù)的聯(lián)合概率密度)，即pi＝[pi¹,pi²,pi³,pi⁴]＝[vpi,vsi,ρi,fi],i＝1,2,…,n。

其中，三參數(shù)的聯(lián)合概率密度的具體表達(dá)式為：

其中，e為期望；σ為縱波速度、橫波速度和密度三參數(shù)的協(xié)方差矩陣，由測(cè)井?dāng)?shù)據(jù)獲取。

在取值范圍內(nèi)，隨機(jī)設(shè)置粒子的初始位置和初始速度。粒子的初始位置組成[4×n]的矩陣，粒子的初始速度設(shè)為零。

步驟三，設(shè)置目標(biāo)函數(shù)，即粒子適應(yīng)度表達(dá)式。本反演方法的目標(biāo)函數(shù)由觀測(cè)地震記錄與合成地震記錄的誤差和三參數(shù)的先驗(yàn)約束項(xiàng)構(gòu)成，具體表現(xiàn)形式如下：

其中，f為粒子適應(yīng)度；r為使用精確佐普里茲方程計(jì)算的縱波反射系數(shù)；w為震源子波；d為觀測(cè)地震記錄；l為觀測(cè)數(shù)據(jù)采樣長(zhǎng)度；λ1和λ2為預(yù)設(shè)定系數(shù)；θ為入射角；t為地震記錄采樣時(shí)間。

步驟四，設(shè)置模擬退火學(xué)習(xí)項(xiàng)。按公式(2)計(jì)算各粒子pi的適應(yīng)度值并計(jì)算各粒子位置被選擇作為模擬退火學(xué)習(xí)項(xiàng)的引導(dǎo)粒子pmin的接收概率，接收概率的表達(dá)式為：

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

基于各粒子的接收概率，使用輪盤選擇算法決定模擬退火學(xué)習(xí)項(xiàng)的引導(dǎo)粒子pmin，顯然，適應(yīng)度最小的粒子具有最大的接收概率，但非最優(yōu)粒子也有一定概率“跳躍”成為引導(dǎo)粒子，且這種概率在高溫時(shí)最大，并隨著溫度降低而減小。這種“跳躍”機(jī)制增加了粒子的搜索擴(kuò)展能力，避免因不成熟收斂而陷入局部最優(yōu)解。

步驟五，設(shè)置多維學(xué)習(xí)項(xiàng)。

固定橫波速度和密度，由公式(1)計(jì)算各粒子的縱波速度對(duì)應(yīng)的聯(lián)合概率密度值，取最大概率密度值對(duì)應(yīng)的縱波速度vpmax；按同樣方法，分別獲得最大概率密度對(duì)應(yīng)的橫波速度vsmax和密度ρmax，即：

多維學(xué)習(xí)項(xiàng)的引導(dǎo)粒子pmax是具有最大概率密度值的擇優(yōu)組合粒子，即

pmax＝[vpmax,vsmax,ρmax,fmax]

多維學(xué)習(xí)項(xiàng)的引導(dǎo)粒子是基于三參數(shù)間的統(tǒng)計(jì)學(xué)關(guān)系而選擇的具有最大概率密度值的組合粒子，該引導(dǎo)粒子不局限于單一粒子所對(duì)應(yīng)的位置，而是多維度、多粒子間擇優(yōu)組合的結(jié)果，可以有效控制疊前地震三參數(shù)同步反演的穩(wěn)定性。

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

粒子速度的更新公式為：

其中，為模擬退火學(xué)習(xí)項(xiàng)，其中的引導(dǎo)粒子pmin由步驟四得到；為多維學(xué)習(xí)項(xiàng)，其中的引導(dǎo)粒子pmax由步驟五得到，上標(biāo)d代表粒子的維度，每個(gè)粒子pi包含四個(gè)維度(即縱波速度、橫波速度、密度三參數(shù)，以及三參數(shù)的聯(lián)合概率密度)。

粒子位置的更新公式為：

步驟七，進(jìn)入下一次迭代k+1，并使溫度降低至tk+1，重復(fù)步驟三至步驟六。直至達(dá)到最大迭代數(shù)k，溫度降至終止溫度tend，輸出適應(yīng)度最優(yōu)的粒子。

下面以一個(gè)合成地震數(shù)據(jù)測(cè)試進(jìn)行具體說(shuō)明：

合成地震數(shù)據(jù)測(cè)試的理論模型如圖2所示，其中，(a)為縱波速度，(b)為橫波速度，(c)為密度。使用精確佐普里茲方程計(jì)算理論模型中地層的縱波反射系數(shù)，計(jì)算的角度間隔為5°，角度覆蓋范圍為0-50°。分別將每層的11組反射系數(shù)與主頻為50hz的零相位理論ricker子波進(jìn)行褶積，得到如圖3所示的合成地震記錄(角度道集)，即觀測(cè)地震數(shù)據(jù)，其中有11個(gè)角度道，時(shí)間采樣間隔2ms。

使用本發(fā)明方法對(duì)觀測(cè)地震數(shù)據(jù)進(jìn)行反演，以獲得縱波速度、橫波速度和密度三參數(shù)的反演結(jié)果，使反演結(jié)果的合成地震數(shù)據(jù)與觀測(cè)地震數(shù)據(jù)之間誤差的二次范數(shù)極小，具體實(shí)現(xiàn)方式如下：

設(shè)置粒子位置的取值范圍，每個(gè)粒子位置包含四個(gè)維度，即縱波速度、橫波速度、密度和聯(lián)合概率密度，對(duì)應(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í)項(xiàng)加權(quán)系數(shù)c1＝c2＝1。

在取值范圍內(nèi)，隨機(jī)生成粒子的初始位置，并對(duì)粒子位置進(jìn)行對(duì)數(shù)歸一化處理；設(shè)定各粒子的初始速度為零。

按步驟三計(jì)算各粒子的適應(yīng)度值，即各粒子對(duì)應(yīng)參數(shù)模型的合成地震記錄與觀測(cè)地震數(shù)據(jù)的誤差的二次范數(shù)。

按步驟四計(jì)算模擬退火學(xué)習(xí)項(xiàng)，基于各粒子的接收概率，按輪盤選擇算法確定該項(xiàng)的引導(dǎo)粒子pmin；按步驟五計(jì)算多維學(xué)習(xí)項(xiàng)，擇優(yōu)組合粒子位置，使其具有最大聯(lián)合概率密度值，構(gòu)成該項(xiàng)的引導(dǎo)粒子pmax，這里使用三參數(shù)模型的平滑結(jié)果作為測(cè)井?dāng)?shù)據(jù)，計(jì)算三參數(shù)的協(xié)方差矩陣。

對(duì)各粒子按步驟六進(jìn)行速度更新和位置更新。

降低溫度并進(jìn)入下一次迭代，利用本發(fā)明方法對(duì)粒子不斷迭代更新，直至終止溫度和最大迭代次數(shù)，輸出最優(yōu)適應(yīng)度的粒子，作為觀測(cè)地震數(shù)據(jù)的反演結(jié)果。

圖4為使用本發(fā)明方法的三參數(shù)反演結(jié)果與理論模型的對(duì)比，其中，(a)為縱波速度，(b)為橫波速度，(c)為密度；可以看到反演結(jié)果(虛線)與理論結(jié)果(實(shí)線)吻合很好，特別是一些高頻成分得到很好的恢復(fù)；圖5為使用傳統(tǒng)粒子群算法的三參數(shù)反演結(jié)果與理論模型的對(duì)比，其中，(a)為縱波速度，(b)為橫波速度，(c)為密度；可以看到反演結(jié)果(虛線)與理論結(jié)果(實(shí)線)吻合較差。圖6為分別使用傳統(tǒng)粒子群算法與本發(fā)明方法進(jìn)行反演的收斂速度的對(duì)比，可以看到使用本算法的收斂速度(實(shí)線)比使用傳統(tǒng)粒子群算法的收斂速度(虛線)明顯較快。

以上實(shí)施例僅為說(shuō)明本發(fā)明的技術(shù)思想，不能以此限定本發(fā)明的保護(hù)范圍，凡是按照本發(fā)明提出的技術(shù)思想，在技術(shù)方案基礎(chǔ)上所做的任何改動(dòng)，均落入本發(fā)明保護(hù)范圍之內(nèi)。