向 坤,陳 科,段心標,鄭棟宇

(1.中國石油化工股份有限公司石油物探技術研究院,江蘇南京 211103;2.成都理工大學,四川成都 610059)

波阻抗反演是LINDSETH[1]于1979年首次提出的方法,通過測井數(shù)據(jù)合成的地震記錄進行地震資料解釋,并且把波阻抗改寫成反射系數(shù)的積分形式。經過數(shù)十年的發(fā)展,波阻抗反演已成為巖石屬性、儲層表征和沉積相分析等領域中最重要的工具[2]。反演方法從確定性反演發(fā)展為隨機性反演,從線性反演發(fā)展為非線性反演,從疊后反演發(fā)展為疊前反演[3-4]。其中,疊前反演是根據(jù)縱橫波速度和密度計算彈性參數(shù),達到劃分巖性和檢測流體的目的。傳統(tǒng)的疊前反演是在確定性反演的基礎上,根據(jù)已有的測井信息建立初始模型,利用Zoeppritz方程的線性近似公式,結合疊前地震記錄估算縱橫波阻抗或彈性阻抗[5-6]。早期的確定性反演存在抗噪性較差和過度依賴初始模型等問題[7],而基于貝葉斯框架下的疊前三參數(shù)同步反演可以在地震數(shù)據(jù)和測井數(shù)據(jù)等先驗信息約束下對反演的結果進行不確定性評估[8-9]。在反演過程中,允許先驗信息加入貝葉斯函數(shù),然后得到最大似然分布以及后驗分布曲線。將得到的最大似然分布和后驗分布曲線應用于概率分析和儲層描述[10-11]。通常的貝葉斯反演在優(yōu)化過程中,采用共軛梯度或者擬牛頓法在空間中線性地尋找全局最優(yōu)位置。這種方法很容易陷入局部最小值,使得在優(yōu)化過程中不能跳出局部最小值尋找全局最小,因而得不到正確結果。基于隨機搜索的疊前同步三參數(shù)反演可以克服共軛梯度法的缺陷,充分利用先驗信息,輸出若干符合地震地質信息的結果[12]。反演目標函數(shù)中的正演算子,早期多采用線性近似方程,而目前多采用Zoeppritz方程,進而保證疊前反演的精度。在隨機反演中,關鍵問題之一就是優(yōu)化目標函數(shù)的效率。盡管人們提出了許多非線性隨機算法解決優(yōu)化問題,但這些方法效率較低,不能滿足實際生產需求。

自適應粒子群優(yōu)化算法(adaptive particle swarm optimization,APSO)是一種改進的非線性優(yōu)化方法,已經被證明比標準粒子群優(yōu)化更具成本效益[13]。在保證反演精度的同時,針對疊前反演中非線性優(yōu)化的效率問題提出了利用自適應粒子群優(yōu)化算法。該方法可以提供與地震反演最佳解一致的有限模型空間中全局最小值解,并且具有良好的抗噪性?？朔舜_定性反演的問題,即過度擬合和良好初始模型的需求。此外,該方法允許貝葉斯統(tǒng)計和馬爾可夫鏈法在數(shù)據(jù)域每一次迭代中更新粒子群。由于在寬角度區(qū)域,傳統(tǒng)的線性近似公式與Zoeppritz公式存在誤差,在反演中不能夠滿足精度條件。因此,在正演過程中選擇了Zoeppritz方程而非其近似公式。最后,利用NS工區(qū)的二維海洋地震資料和測井資料,估算縱橫波速度和密度,并且與傳統(tǒng)的隨機反演方法對比,驗證了該方法的有效性和可行性。

1 方法技術

1.1 疊前三參數(shù)同步反演

疊前三參數(shù)同步反演,即利用疊前地震道集中包含的振幅和角度信息同時反演縱波速度、橫波速度和密度3個參數(shù),進而用于指示儲層的物性和含油氣性[14]?；谕暾问降腪oeppritz方程疊前三參數(shù)反演,可以精確地描述振幅隨入射角變化,并且在入射角接近臨界角時依然可以保證反演的精度[15-16]。由于波阻抗反演問題總是不適定的,并且地震數(shù)據(jù)的帶寬有限以及缺少低頻分量,因此幾乎不可能從地震道中直接估算絕對阻抗[17-18]。因此,在貝葉斯框架下的疊前三參數(shù)同步反演,可以根據(jù)生成的概率分布,減少反演結果的不確定性[19-20]。在貝葉斯框架的基礎上,結合隨機反演的特點,建立了疊前三參數(shù)同步反演流程,如圖1所示。首先,對疊前地震數(shù)據(jù)進行處理得到角道集,對測井數(shù)據(jù)進行處理得到低頻背景模型。此外,利用測井信息建立縱、橫波速度以及密度的擾動先驗分布。利用先驗分布結合Metropolis-Hasting方法實現(xiàn)隨機抽樣,建立初始粒子群。初始粒子群經過不斷優(yōu)化,最終得到3個參數(shù)的全局最優(yōu)值和后驗分布。

圖1 疊前三參數(shù)同步反演流程

在貝葉斯框架下的疊前三參數(shù)同步反演中,后驗分布由先驗分布表示,根據(jù)給定的疊前地震數(shù)據(jù)Dobs,那么三參數(shù)的后驗分布p(m|Dobs)可以表示為:

(1)

式中:p(m)是擾動模型m的先驗概率分布。在三參數(shù)疊前反演的問題中,m為縱橫波速度和密度。假定縱橫波速度和密度服從高斯分布,那么模型的先驗分布表達式為:

(2)

式中:σm是模型的協(xié)方差;μm是模型的期望值。在貝葉斯框架中,構造似然函數(shù)是檢驗模型選擇和評估模型不確定性考慮的另一種方法。在此函數(shù)中,如果假設噪聲均值為零,且具有噪聲協(xié)方差的高斯噪聲分布,那么噪聲的分布表達式為:

(3)

因此,根據(jù)疊前角道集Dobs和噪聲協(xié)方差σN,那么似然函數(shù)p(Dobs|m)的表達式為:

(4)

式中:R(m)是Zoeppritz縱波反射系數(shù)方程;Sw是

利用角道集提取的各個入射角子波。

1.2 自適應粒子群

在隨機反演的優(yōu)化過程中,粒子群優(yōu)化算法(particle swarm optimization,PSO)是被廣泛應用于目標函數(shù)的優(yōu)化和機器學習的一種方法[21]。該算法通過研究群體中個體之前的信息共享和協(xié)作在限定的空間中尋找最優(yōu)解。雖然粒子群算法已被廣泛應用于解決各種復雜的工程問題,但在實際問題中仍需研究克服粒子群算法計算量大和精度低等缺點[22]。尤其在地球物理反演中,針對高維度參數(shù)空間的反演,需要更為高效的算法。其中,自適應粒子群(APSO)算法加入了自適應參數(shù)和精英學習策略(elite learning strategy,ELS)。與傳統(tǒng)粒子群優(yōu)化算法相比,自適應粒子群算法在求解精度、收斂速度和算法可靠性等方面更具優(yōu)勢。由于收斂速度比傳統(tǒng)的粒子群算法更快,因而自適應粒子群算法在全局搜索上具有更高的計算效率。在疊前三參數(shù)同步反演中,自適應粒子群算法在優(yōu)化目標函數(shù)和尋找地震反演問題的全局最優(yōu)解方面比粒子群優(yōu)化算法有更好的計算能力[23]。

在粒子群優(yōu)化的過程中,將粒子的當前最優(yōu)值和全局最優(yōu)值分別表示為pBest[i]和pBest[g]。假設粒子群中粒子的個數(shù)為N,那么將第i次迭代中某個粒子的當前空間位置表示為xi。圖2顯示了空間Ω中第i次迭代粒子的當前位置和更新后的位置。由圖2可見,更新速度取決于粒子的當前最優(yōu)值和全局最優(yōu)值以及上次的更新速度。在傳統(tǒng)粒子群優(yōu)化算法的迭代中,粒子更新速度和更新后的新位置可以寫成:

圖2 粒子群更新方向示意

vi=ωvi-1+c1*rand1*(pBest[i]-xi)+
c2*rand2*(pBest[g]-xi)

(5)

xi+1=vi+xi

(6)

式中:xi為當前粒子的位置;xi+1為迭代更新速度后粒子的位置;v為更新速度;c1和c2為學習因子;rand1和rand2為[0,1]之間的隨機數(shù);ω為慣性因子。

而自適應粒子群算法在粒子群優(yōu)化算法的基礎上,修正了慣性因子,并且增加了精英學習策略。圖3 描述了在反演過程中,自適應粒子群算法的具體工作流程,在迭代過程中每個粒子會不斷移動,直到滿足迭代的終止條件。在疊前三參數(shù)同步反演中,實現(xiàn)自適應粒子群算法包括以下幾個步驟。

圖3 自適應粒子群優(yōu)化算法技術流程

首先是粒子的初始化階段。對于每一個粒子,里面包含有縱波速度、橫波速度以及密度3個分量。根據(jù)縱波速度、橫波速度以及密度的先驗分布,構建粒子的初始狀態(tài)。針對粒子中每一道中的采樣點,粒子的初始狀態(tài)可以表示為:

(7)

(8)

ρ0=Lρ+Δρ

(9)

(10)

在粒子群演化的過程中,用參數(shù)fd表示粒子群的演化狀態(tài)。在迭代中,計算粒子k到其它粒子j的距離可以表示為:

(11)

式中:N代表粒子群中總的粒子數(shù);D是單個地震道參與反演的采樣點數(shù)。之后比較所有粒子的距離,對距離di從大到小進行排序,計算演化因子fd,定義如下:

(12)

式中:dg為全局最佳位置所對應的距離。fd在初始階段的默認值為0,可以通過演化因子判斷當前粒子群的狀態(tài)。圖4為不同的演化因子所對應的演化狀態(tài),共4個狀態(tài),即跳出階段,探索階段,開發(fā)階段以及收斂階段。不同的演化因子對應著不同的探索階段。

圖4 演化因子對應的狀態(tài)示意

然后,計算自適應慣性參數(shù),該參數(shù)可以平衡搜索局部最佳和全局最佳位置的能力。

(13)

式中:ω(fd)在每次迭代中從0.9減小到0.4。因此,由于演化因子的存在,ω(fd)使得粒子群優(yōu)化算法變化更為有效。由于fd初始值為0,因此ω(fd)的初始值為0.4。

在迭代過程中,需要定義慣性系數(shù)c1和c2,如:

ci(i+1)=ci(i)±δc,i=1,2

(14)

每一次迭代中,c1和c2的值都在變化。為了收斂,c1越大,c2的值越小;反之,c1的值越小,c2的值越大。c1和c2在[1.5,2.5]范圍內變化,它們的總和應在[3.0,4.0]范圍內,以便粒子在局部最小值和全局最小值之間進行調整,δc是在[0.05,0.10]范圍內均勻生成的隨機值。而在傳統(tǒng)的粒子群算法中,c1和c2等于2,不具備調整的能力。在迭代過程中,采取精英學習策略(ELS)尋找全局最佳位置。標準的粒子群優(yōu)化算法如果陷入局部最優(yōu),就不能幫助全局最優(yōu)從當前最優(yōu)中跳出來。因此,提出利用精英學習策略來克服這一不足。當從每一次迭代中找出全局最佳時,它將增加輕微的擾動量,通過高斯擾動利用學習策略在空間中移動,如:

(15)

σ=σmax-(σmax-σmin)*(i/T)

(16)

式中:σmax=1.0和σmin=0.1是從實證研究測試得出的。i為當前的迭代次數(shù),T為總的迭代次數(shù)。如果新的pk比當前位置更適合,那么它將被接受。當滿足預定最小值或達到最大迭代次數(shù)時的收斂準則時,輸出全局最優(yōu)解和若干次優(yōu)解。

1.3 基于粒子群優(yōu)化-馬爾科夫鏈蒙特卡洛的隨機模擬反演實現(xiàn)過程

在疊前三參數(shù)同步反演中,已經假設模型的擾動分布服從高斯分布。根據(jù)三參數(shù)擾動先驗分布抽樣獲取初始粒子群[24]。首先,利用測井曲線中的三參數(shù)信息,去掉低頻成分,生成三參數(shù)擾動曲線。通過對擾動進行統(tǒng)計,建立三參數(shù)樣本點x和樣本集X。

典型的蒙特卡洛方法以相對較小的步長圍繞著分布移動,而這些步長沒有具備朝著相同方向進行的趨勢。盡管這些方法易于實現(xiàn)和分析,但可能需要很長時間才能探索空間中的最佳位置[25]。因此為了加速馬爾科夫鏈蒙特卡羅算法(MCMC),提出了與自適應粒子群優(yōu)化算法相結合的馬爾科夫鏈蒙特卡羅算法,即自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法(APSO-MCMC)。這樣在保證隨機性的同時,進一步提高全局搜索的計算效率。在隨機模擬的過程中,馬爾科夫鏈蒙特卡羅算法需要一個實用的解決方案來解決生成躍遷矩陣的問題。躍遷矩陣與自適應粒子群優(yōu)化算法相關聯(lián),在粒子群中根據(jù)每個粒子的目標函數(shù)值優(yōu)選后,將優(yōu)選得到的粒子按照概率p1繼續(xù)在空間中搜索,而剩下的粒子則以一定概率p2加入隨機分量,如果在加入隨機分量后目標函數(shù)值小于之前,則保留當前結果;否則停留在先前結果。按照這個方法構建躍遷矩陣,在加速搜索過程的同時,保持了其自身的收斂性。在本方法中,躍遷矩陣的表達式為:

(17)

式中:p1和p2為自身轉移的概率,其大小可以根據(jù)實際情況定義。當以p1的概率搜索時,表達式為:

xi=xi+p1*rand(0,1)*Vi

(18)

式中的Vi由自適應粒子群優(yōu)化算法在每次迭代過程中獲得。圖5描述了自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法的簡要技術流程?；谀繕撕瘮?shù)和過渡狀態(tài),可以計算新概率,然后確定過程是否繼續(xù)。主要步驟為:

圖5 自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡洛算法在第M次迭代過程中的流程示意

1) 從三參數(shù)的擾動分布中隨機采樣,生成初始狀態(tài)X0,并對初始粒子群進行優(yōu)選;

2) 在迭代時,對每個粒子在轉化過程中應用躍遷矩陣P,然后計算目標函數(shù)的值。按目標函數(shù)值降序對粒子排序,然后輸出最佳粒子;

3) 粒子群應用自適應粒子群優(yōu)化算法方法,獲得局部最優(yōu)pBest[i]和全局最優(yōu)pBest[g]。計算自適應慣性系數(shù)以確定演化狀態(tài)ω(fd),然后更新粒子迭代后的位置XM+1和速度VM+1;

4) 根據(jù)目標函數(shù)的值確定每個粒子移動到下一個位置或停止移動;

5) 確定收斂準則。如果全局最佳pBest[g]小于預設定值ε或達到最大迭代次數(shù)后,則終止迭代并轉到步驟7)。如果未達到收斂準則,進行步驟6);

6) 計算更新粒子的分布,增加迭代次數(shù),然后返回步驟2);

7) 輸出最佳粒子,然后獲得最大似然結果以及后驗分布。

2 應用實例

2.1 NS工區(qū)概述

利用本方法對NS工區(qū)二維海洋地震數(shù)據(jù)進行反演。由于海洋地震數(shù)據(jù)包含大量多次波,因而在反演前需要對數(shù)據(jù)進行預處理。在預處理后,通過疊前時間偏移得到共反射點道集和含有角度信息的疊前數(shù)據(jù)體。測井解釋的結果表明,在2000m左右含有豐富的含氣顯示。反演之前,根據(jù)實際地下情況建立含氣飽和地層模型,然后利用Zoeppritz方程分析工區(qū)含氣儲層的AVO效應。根據(jù)測井信息將模型簡化成一個雙層模型,上層為泥巖,下層為含氣砂巖,并且假定雙層模型的各層是均勻介質。根據(jù)模型得到的正演模擬結果如圖6所示,其曲線特征符合第三類AVO響應,即地震波的負振幅隨著反射角的增加而顯著增加。將此特征關系與實際的疊前道集中的含氣儲層區(qū)段的振幅比較后,振幅特征與AVO曲線特征一致,說明工區(qū)內發(fā)育有豐富的含氣儲層。而后利用井旁道的低頻信息,結合疊前數(shù)據(jù),對三參數(shù)最大后驗估計進行計算,并和實際的測井曲線對比,估算反演的精度和誤差。完成井旁道反演后,再對整個剖面進行反演。此外,利用同樣的數(shù)據(jù),將本方法與傳統(tǒng)的粒子群算法(PSO)和遺傳算法(GA)對比,觀察同等迭代次數(shù)的情況下不同方法的計算效率。

圖6 AVO正演結果

2.2 反演結果分析

在反演過程中,首先對井旁道進行反演,以驗證反演結果的有效性。通過提取井旁道的道集,結合測井曲線中的低頻信息,建立擾動分布,生成初始粒子群后進行疊前三參數(shù)反演。圖7為縱波速度的擾動統(tǒng)計結果和先驗模型,利用統(tǒng)計結果建立針對縱波速度擾動的高斯分布。用同樣的方法建立橫波速度和密度的擾動分布。利用該高斯分布建立先驗模型,然后對該分布進行隨機采樣,生成若干初始擾動模型,形成初始粒子群。粒子群由40個粒子構成。通過對整個粒子群的反演,不斷估算演化因子,判定收斂狀態(tài)。圖8顯示了演化因子隨迭代次數(shù)變化的情況。由圖8可見,隨著迭代次數(shù)的增加,演化因子在初始階段逐步增加,而后逐漸變小。說明在反演初期階段,粒子群中的粒子在空間中不斷搜索全局最優(yōu)位置?？s小搜索的范圍后,狀態(tài)逐漸穩(wěn)定,演化因子變小,進入收斂狀態(tài)。

圖8 演化因子隨迭代次數(shù)變化的情況

從全部40個粒子中選擇5個最優(yōu)解。圖9為5個粒子的縱波反演結果。盡管粒子之間的反演結果不一致,但這些最優(yōu)解的粒子在搜索空間中彼此接近,并且反演結果都符合巖石物理模型。從輸出模型的后驗分布中選擇最大后驗估計(MAP),并將其與測井曲線進行比較。圖10為反演得到的最大后驗估計與實際真值和低頻模型的對比?？v波速度和密度顯示出與實際測井信息對比有良好的一致性,總體誤差小于5%,并且反演結果不過分依賴于低頻模型。雖然反演得到的橫波速度的誤差相比縱波速度和密度大,但是整體趨勢仍然與實際真值保持一致。將測井解釋結果與泊松比曲線對比發(fā)現(xiàn),三參數(shù)反演計算得到的泊松比,在氣層處顯示有低泊松比的特征。因此利用這一特征,可以幫助識別含氣儲層。此外,提取了粒子群中的某個粒子的目標函數(shù)值和全局最優(yōu)值,其對比情況如圖11所示。粒子在不斷收斂,最后粒子群中的所有粒子都會在搜索空間中趨于一個方向。在反演過程中追蹤粒子群中的第30個粒子可以發(fā)現(xiàn),在迭代初期階段目標函數(shù)值為34.12,隨著迭代次數(shù)的增加,目標函數(shù)的值受全局最小值的影響不斷減少,并且逐步收斂。而目標函數(shù)的全局最小值在反演初期從最初的23.91經過500次迭代后逐漸降低至7.36,進一步證實了該算法的有效性。圖12為反演后的合成地震記錄(藍色)與實際地震記錄(紅色)。對比發(fā)現(xiàn),相關系數(shù)達到了0.97,反演結果符合測井信息和地震信息。在驗證了井旁道反演的有效性后,而后針對工區(qū)的整條剖面進行反演。圖13、圖14和圖15分別為反演得到的二維剖面上的縱、橫波速度和密度。從二維剖面上可以看出,反演結果在保留低頻信息的同時,增加了很多細節(jié),并且與地震波形的變化特征保持一致。圖16為根據(jù)三參數(shù)反演結果計算的泊松比剖面,結合井上含氣儲層的泊松比特征,可以進一步利用泊松比指示目的層含氣儲層的存在。

圖9 反演得到的縱波速度

圖10 反演得到的最大后驗估計和對應的真值以及低頻模型a 縱波速度; b 橫波速度; c 密度; d 泊松比

圖11 迭代中單個粒子和全局最優(yōu)的更新

圖13 反演得到的縱波速度剖面

圖14 反演得到的橫波速度剖面

圖15 反演得到的密度剖面

圖16 根據(jù)反演結果計算得到的泊松比剖面

此外,通過與傳統(tǒng)的粒子群算法和遺傳算法對比后發(fā)現(xiàn),該算法具有更快的收斂速率,表1對比了3種方法在同一目標函數(shù)下不同迭代次數(shù)的耗時情況,Gbest代表每次迭代過程中目標函數(shù)的值;t代表迭代所需要的時間。在迭代初期,3種方法差距不大,并且粒子群優(yōu)化算法的單次迭代時間要少于遺傳算法的單次迭代時間。隨著迭代次數(shù)的增加,可以發(fā)現(xiàn),自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法的優(yōu)化速率明顯快于其它兩種算法。雖然自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法的計算時間比粒子群優(yōu)化算法所需要的時間略長,但該方法能夠在固定的迭代次數(shù)內更加快速收斂,并且在固定的迭代次數(shù)下,自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法花費的時間比遺傳算法花費的時間少。與常規(guī)粒子群優(yōu)化算法相比,雖然每一步中增加了演化因子的計算,潛在地增加了時間,但是,與粒子群優(yōu)化算法相比,自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法的收斂更快,能夠在空間中更快地搜索全局最小值。針對高維度的反演問題,需要更加快速的收斂,因此自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法相比粒子群優(yōu)化算法更高效。

表1 同一目標函數(shù)下不同優(yōu)化方法的比較

3 結論

本文討論了傳統(tǒng)的疊前波阻抗地震反演方法存在的不足,在此基礎上分析了同步隨機反演的優(yōu)勢。在利用Zoeppritz方程保證隨機反演精度的同時,實現(xiàn)了基于自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅反演算法,該方法用于對NS地區(qū)實際地震資料的測試結果表明:

1) 在構建正演算子過程中利用Zoeppritz公式而非近似公式,保證了反演的精度,反演過程中不斷隨機尋找全局最優(yōu)值而非確定性方法,確保了反演的可靠性;

2) 基于自適應粒子群優(yōu)化-馬爾科夫鏈蒙特卡羅算法的反演算法能夠在有限的空間中尋找全局最優(yōu)位置,方法中包含的演化因子可以根據(jù)其大、小判定搜索狀態(tài),方法中的精英策略能夠幫助粒子跳出局部最優(yōu)位置;

3) NS地區(qū)二維疊前數(shù)據(jù)的應用結果表明該地區(qū)的含氣儲層符合第三類AVO特征響應,經過本方法的測試驗證了應用的有效性和可行性。將本方法與傳統(tǒng)的粒子群算法和遺傳算法進行對比測試后發(fā)現(xiàn),新的方法具有更高的計算效率。

未來將巖性、深度和地層反射系數(shù)等約束信息加入到方法中,可以使反演的結果與沉積相、地震相和巖相獲得的認知更為一致,進一步提高反演的可靠性和實用性。