[發明專利]一種提高分辨率的孔隙度反演方法

有效

• 申請號: 201911288767.8
• 申請日: 2019-12-12
• 公開（公告）號: CN111077571B
• 公開（公告）日: 2020-11-06
• 發明（設計）人: 周東勇;文曉濤;秦子雨;李壘
• 申請（專利權）人: 成都理工大學
• 主分類號: G01V1/30
• 分類號: G01V1/30;G01V1/28;G01V1/40
• 代理公司: 北京元本知識產權代理事務所(普通合伙) 11308
• 代理人: 王紅霞

1.一種提高分辨率的孔隙度反演方法，其特征在于，該方法包括如下步驟：假定存在待反演地震數據、測井反射系數及測井孔隙度，利用模擬退火算法實現孔隙度反演：

(1)推導垂直入射至含流體孔隙介質分界面上的快縱波反射系數近似式，基于巖石物理理論建立孔隙度和地震記錄之間的直接定量關系；

(2)根據孔隙度和地震記錄之間的關系，以測井孔隙度和測井反射系數作為先驗約束，構建非線性反演目標函數；

(3)利用測井數據信息計算孔隙度調節參數，給定約束參數、最大反演迭代次數及孔隙度初始解等先驗信息；

(4)將上述已知參數帶入非線性反演目標函數，結合快速模擬退火算法，實現孔隙度的高精度直接定量反演；

所述步驟(1)具體包括：孔隙度與地震數據之間的定量關系：

基于Biot彈性孔隙介質理論，快縱波垂直入射到含飽和流體孔隙介質分界面上的快縱波反射系數近似式為：

其中，r為快縱波反射系數，K_l，μ_l，l＝1，2分別表示巖石骨架的體積模量和剪切模量，α_l，M_l，l＝1，2為Biot參數；vp_lj，l＝1，2，j＝1，2為縱波速度，γ_lj，l＝1，2，j＝1，2為孔隙流體相對于巖石骨架的位移位振幅和巖石骨架的位移位振幅之比；下標l＝1，2分別表示上、下層孔隙介質，j＝1，2分別表示快、慢縱波；

巖石骨架模量、巖石基質模量和孔隙度之間存在聯系，利用(8)和(9)式描述他們之間的定量關系，

K_d＝(1.0-f)K_s/(1.0+c_k·f)， (8)

μ_d＝(1.0-f)μ_s/(1.0+c_μ·f)， (9)

其中，f為孔隙度，K_d、μ_d分別為巖石骨架的體積模量和剪切模量，K_s、μ_s分別為巖石基質的體積模量和剪切模量，c_k、c_μ為孔隙介質調節參數；

含飽和流體孔隙介質中存在快、慢兩類縱波，慢縱波在孔隙介質中傳播時會迅速衰減消失，因此地面檢波器接收到的地震記錄事實上是快縱波的反射，孔隙度f是時間t的函數f(t)，而由(1)-(9)式可知快縱波反射系數r是孔隙度f(t)的函數r[f(t)]，根據褶積公式，與孔隙度直接相關的快縱波合成地震記錄d[f(t)]可以表示為，

d[f(t)]＝w(t)*r[f(t)]+n(t)， (10)將上式離散化，并寫成向量的形式為，

d_N×1＝G_N×N·r_N×1+n_N×1. .(11)

其中，w(t)為時間域子波，n(t)為噪聲，下標N×1表示N維列向量，N×N表示N維矩陣，N為一道地震數據的采樣點數；符號‘*’表示褶積算子；

所述步驟(2)中具體包括：

根據(11)式中的正演過程，以測井中的先驗信息為約束，含流體孔隙介質孔隙度反演目標函數F為，

其中，f為預測孔隙度，R為預測反射系數，d為合成地震記錄；f₀和R₀分別為由測井數據得到的孔隙度和反射系數，d_data為實際地震數據；β為反射系數約束參數，γ為孔隙度約束參數，目標函數中的所有向量均為N×1維列向量，表示向量的L2范數；

上式中，第一項為觀測地震數據和合成地震記錄的擬合程度，第二項表示由預測孔隙度計算的反射系數與測井反射系數的相似性，第三項表示預測孔隙度與測井孔隙度的相似性，后兩項是對反演結果的測井約束；

所述步驟(4)具體包括：

(41)Metropolis收斂準則：

反演算法通過不斷迭代更新孔隙度來獲得最優解，在第i次迭代時，首先計算與孔隙度f_i對應的目標函數F_i，再更新孔隙度為f_i+1，并計算其對應的目標函數值F_i+1，兩函數值之差ΔF為，

ΔF＝F_i+1-F_i，i≥1， (13)

若ΔF＜0說明孔隙度的更新方向使目標函數值減小，則接受孔隙度的修改；若ΔF≥0，則進一步判斷ΔF是否滿足(14)式，若滿足條件，則仍接受孔隙度的修改，否則拒絕修改；

0＜p＜K， (14)

其中，p是概率密度，K是一個隨機數，其值域為0＜K＜1；k_b為玻爾茲曼常數，實際應用中常取1，exp{·｝為指數函數，T為當前退火溫度，常表示成指數形式，

T_i＝T₀·(0.95)^i-1，i≥1， (16)

其中,T₀為初始溫度，T_i為第i次退火溫度；

(42)解的迭代與更新：

將測井孔隙度f₀經平滑后作為初始孔隙度f₁，給定初始解后，采用依賴于溫度的似Cauchy分布產生更新向量Δf_i，即

f_i+1＝f_i+ξ·Δf_i，i≥1， (17)

Δf_i＝T·sign(q-0.5)·[(1+1/T)^|2·q-1|-1]. (18)

上式中，ξ為更新系數，f_i為當前預測孔隙度，f_i+1為更新后的孔隙度，其值域為f_i+1∈[0，0.3]；q是服從[0，1]分布的N維隨機向量，sign為符號函數；

(43)非線性孔隙度反演流程：

假設存在待反演地震數據d_data、測井反射系數R₀、測井孔隙度f₀、初始孔隙度f₁、最大迭代次數L及測井約束參數β和γ，基于模擬退火算法實現孔隙度反演的流程如下：

431)對于第i次迭代的孔隙度f_i，利用(17)式計算新的孔隙度值f_i+1，再基于(12)式計算孔隙度更新前后反演目標函數值F_i和F_i+1；

432)結合(13)式計算兩目標函數值之差ΔF，并判斷是否滿足條件ΔF＜0，若條件滿足則接受孔隙度值的更新，否則繼續判斷是否滿足(14)式，若滿足則仍可接受孔隙度值的更新，否則拒絕更新；

433)判斷是否達到最大迭代次數L，若是則退出循環，獲得最優孔隙度值；否則利用(16)式迭代更新當前退火溫度，并返回431)繼續迭代更新孔隙度值f。