簡介本資源是一個面向地球物理勘探與地震資料處理初學(xué)者的MATLAB入門級AVO正演建模工具包聚焦于振幅隨偏移距變化AVO理論的編程實現(xiàn)與可視化驗證。壓縮包為RAR格式僅含1個核心文件——avoMODING.m腳本大小925B結(jié)構(gòu)精簡但功能完整涵蓋AVO參數(shù)輸入、Shuey或Aki-Richards等經(jīng)典正解模型調(diào)用、角度域振幅計算及AVO響應(yīng)曲線繪制等關(guān)鍵環(huán)節(jié)。已有135人學(xué)習(xí)下載適用于高校地質(zhì)工程/地球物理學(xué)專業(yè)課程實踐、科研入門訓(xùn)練或地震解釋方法自學(xué)。讀者可直接運行該腳本通過修改速度模型、密度、泊松比等巖石物理參數(shù)實時觀察不同巖性組合下的AVO響應(yīng)特征快速建立理論公式與實際地震表現(xiàn)之間的映射關(guān)系并為后續(xù)流體識別與儲層預(yù)測打下編程與建?；A(chǔ)。 你有沒有過這種經(jīng)歷從某個網(wǎng)盤或者U盤里扒下一個名為 avoMODING.rar 的壓縮包解壓出來一坨 MATLAB 的 .m 文件文件名倒是挺規(guī)整——zoeppritz_avo.m、shuey_approx.m、fluid_replace.m——但當(dāng)你雙擊運行主腳本屏幕上要么報出一串矩陣維度的紅色錯誤要么畫出來的圖和你在地震教科書的PPT里看到的那張AVO道集完全對不上。這個壓縮包我在實驗室?guī)腿伺胚^好幾次了今天干脆把它徹底講明白AVO正演模擬到底在模擬什么、里面的 MATLAB 例程每行代碼在干嘛、以及你拿到這種來路不明的 rar 之后應(yīng)該怎么最快跑出第一張可用的道集。這篇文章適合三類人剛接手疊前道集解釋的勘探地球物理方向?qū)W生、需要在項目中快速搭一套AVO正演驗證流程的工程師、以及單純想搞懂“一個反射系數(shù)是怎么隨入射角變化的”的MATLAB使用者。我會從物理原理講到代碼實現(xiàn)再講到排查經(jīng)驗和擴展思路保證你讀完能自己動手改參數(shù)、畫出有意義的結(jié)果。1. avoMODING到底是干什么的一次AVO正演模擬的完整需求拆解1.1 壓縮包背后的物理問題AVO 全稱 Amplitude Variation with Offset中文一般叫“振幅隨偏移距變化”或“振幅隨入射角變化”。它要回答的問題非常直接當(dāng)一束地震波以不同角度打到地下某個巖性分界面上時反射回來的能量大小會不會變?nèi)绻麜冏兓姆绞胶蛶r層里的流體油、氣、水有什么關(guān)系這個問題的工程背景是傳統(tǒng)的地震剖面只能看到反射界面的“亮點”或“暗點”但亮點不一定是油氣可能是煤層、火成巖或者鈣質(zhì)夾層。而 AVO 引入了一個額外的維度——入射角。你可以把它想象成用不同角度的燈光照同一個物體如果這個物體是啞光的正面照和斜著照亮度差別不會太大如果它是鏡面的稍微換個角度反射光強度就劇烈變化。地下巖層里的流體種類恰恰會改變反射系數(shù)對入射角的“敏感度”。avoMODING 這個 rar 包里的 MATLAB 例程核心任務(wù)就是把這種“反射系數(shù)隨入射角變化”的曲線、道集、交會圖給算出來。它是疊前地震解釋的最前端工具后面接的 AVO 屬性分析、流體因子反演、彈性波阻抗反演全都建立在這個正演模擬的基礎(chǔ)上。1.2 一個例程至少應(yīng)該包含哪些模塊拿到一個 avoMODING.rar我建議你先別急著運行打開文件夾看看它有沒有這幾類文件。一個像樣的 AVO 正演例程至少應(yīng)該包含模塊對應(yīng)文件常見命名作用精確反射系數(shù)計算zoeppritz_avo.m / solve_zoeppritz.m用 Zoeppritz 方程求解四個反射/透射系數(shù)近似公式計算shuey_approx.m / aki_richards_approx.m用線性近似公式快速計算 R(θ)用于對比和屬性分析模型參數(shù)設(shè)置model_parameters.m / define_model.m定義上下層的 Vp、Vs、密度和入射角范圍道集生成與繪圖plot_avo_gather.m / wiggle_trace.m把反射系數(shù)顯示成道集或曲線流體替換fluid_replace.m / gassmann.m利用 Gassmann 方程在含水、含油、含氣之間切換看 AVO 響應(yīng)差異如果你的 rar 里只有前兩個文件那多半是個閹割版建議自己補一個參數(shù)設(shè)置腳本和繪圖函數(shù)不然沒法直觀看到結(jié)果。如果文件特別多而且互相亂調(diào)用也別慌先用matlab的依賴分析工具或者手動grep一下函數(shù)名理清調(diào)用關(guān)系。我在實際折騰這個例程包的時候發(fā)現(xiàn)一個規(guī)律大部分人拿它跑不出結(jié)果不是代碼本身的問題而是他根本不知道“正演的前提是先定義模型”。AVO 正演不是從地震數(shù)據(jù)里提取什么東西而是先假設(shè)“地下有一個含氣砂巖它的 Vp、Vs、密度是這樣”然后基于彈性波動理論算出這個模型應(yīng)該產(chǎn)生什么樣的反射振幅。所以參數(shù)的合理性直接決定結(jié)果的可用性。2. Zoeppritz方程和它的三個近似avomod的數(shù)學(xué)骨架2.1 精確解到底怎么求Zoeppritz 方程是 1919 年提出的它基于界面兩側(cè)位移連續(xù)和應(yīng)力連續(xù)的邊界條件聯(lián)立四個方程同時求解入射縱波在界面上產(chǎn)生的反射縱波PP、反射橫波PS、透射縱波TP、透射橫波TS四個振幅系數(shù)。在 MATLAB 里實現(xiàn)這個方程核心就是一個 4×4 矩陣的求解問題。我貼一段在例程包里最常見的實現(xiàn)方式單位統(tǒng)一用 m/s 和 kg/m3入射角用度內(nèi)部轉(zhuǎn)弧度function [Rpp, Rps, Tpp, Tps] zoeppritz_avo(vp1, vs1, rho1, vp2, vs2, rho2, theta1) th1 theta1 * pi / 180; p sin(th1) / vp1; % 射線參數(shù)Snell 定理 th2 asin(p * vp2); % 透射縱波角 ph1 asin(p * vs1); % 反射橫波角 ph2 asin(p * vs2); % 透射橫波角 M [ sin(th1) cos(ph1) -sin(th2) cos(ph2) cos(th1) -sin(ph1) cos(th2) sin(ph2) sin(2*th1) (vp1/vs1)*cos(2*ph1) (rho2*vp2*vs2)/(rho1*vp1*vs1)*sin(2*th2) -(rho2*vp2*vs2)/(rho1*vp1*vs1)*cos(2*ph2) cos(2*ph1) -(vs1/vp1)*sin(2*ph1) -(rho2*vp2)/(rho1*vp1)*cos(2*ph2) -(rho2*vs2)/(rho1*vp1)*sin(2*ph2) ]; B [ -sin(th1) cos(th1) sin(2*th1) -cos(2*ph1) ]; X M \ B; Rpp X(1); Rps X(2); Tpp X(3); Tps X(4); end這里最容易踩的坑是不同教材對 Zoeppritz 矩陣的符號約定不一樣有的把應(yīng)力的正方向定義成朝下有的把位移分量取正方向定義成朝上導(dǎo)致最終結(jié)果看起來差一個負號。所以寫完矩陣先別急著往下接先用垂直入射θ0驗證此時反射系數(shù)應(yīng)該約等于 (Z2-Z1)/(Z2Z1)ZρVp 是波阻抗。如果對不上優(yōu)先檢查第三行、第四行的符號而不是去改入射角。2.2 Shuey近似為什么是實際項目里最常用的Zoeppritz 的精確解雖然理論完整但公式復(fù)雜物理直覺差。1985 年 Shuey 在 Aki-Richards 線性近似的基礎(chǔ)上把反射系數(shù)改寫成關(guān)于入射角的顯式表達式R(θ) R0 G·sin2θ K·(tan2θ - sin2θ)其中R0 是法向入射反射系數(shù)也叫 AVO 截距Intercept反映垂直入射時的振幅強度G 是 AVO 梯度Gradient控制振幅隨入射角變化的速度是整個 AVO 分析里最核心的屬性K 與縱波速度相對變化率有關(guān)在入射角小于 30 度時第三項貢獻很小通常省略例程包里的 shuey_approx.m 實現(xiàn)通常長這樣function R shuey_approx(vp1, vs1, rho1, vp2, vs2, rho2, theta) th theta * pi / 180; dvp vp2 - vp1; drho rho2 - rho1; vp (vp1 vp2) / 2; rho (rho1 rho2) / 2; vs (vs1 vs2) / 2; R0 0.5 * (dvp/vp drho/rho); G R0 - (dvp/vp) * 4*(vs/vp)^2 - (drho/rho) * 2*(vs/vp)^2; K 0.5 * dvp/vp; R R0 G * sin(th).^2 K * (tan(th).^2 - sin(th).^2); end別看這公式簡單它把復(fù)雜的彈性波傳播問題壓縮成了三個參數(shù)和兩個三角函數(shù)項直接讓后續(xù)的截距-梯度分析成為可能。實際解釋的流程是把實際地震道集上每個反射界面的振幅隨角度的變化趨勢擬合出來得到截距 P 和梯度 G然后看 P×G 的異常。含氣砂巖的 P×G 通常會出現(xiàn)明顯負異常而含水砂巖雖然有負的 P但 G 不會顯著變負。2.3 誤差邊界什么時候不能再用兩項近似我見過不少人把 Shuey 近似當(dāng)萬能公式用入射角都采到 45 度了還在拿兩項近似做擬合結(jié)果梯度 G 被嚴重污染。實測下來當(dāng)入射角超過 30 度以后公式里的第三項 K·(tan2θ - sin2θ) 的貢獻會迅速增大如果你只取前兩項擬合出來的 R0 和 G 是有偏的。所以例程包里如果同時有精確解和近似解我建議你在主程序里同時計算兩條曲線并輸出相對誤差或直接疊圖。常規(guī)的界限是最大入射角推薦方法小于 20 度兩項 Shuey 近似完全夠用20~30 度三項 Shuey 近似注意密度項精度大于 30 度優(yōu)先用 Zoeppritz 精確解Shuey 只用于趨勢分析這個“先看角度范圍再選公式”的習(xí)慣能幫你避開很多解讀階段的假象。正演里算錯的反射系數(shù)到了反演階段就是地震資料上的假亮點。3. 跑通例程的完整流程從解壓到畫出第一張AVO道集3.1 解壓后第一件事檢查文件結(jié)構(gòu)與依賴關(guān)系把 avoMODING.rar 解壓到本地之后我強烈建議第一件事不是雙擊運行而是把文件夾放到一個純英文路徑下比如D:\codes\avoModing\。Windows 下 MATLAB 對中文路徑的支持時好時壞尤其是當(dāng)你后面要調(diào)用 MEX 文件、第三方工具箱或者寫入文件時中文路徑會帶來一堆莫名其妙的報錯。這個習(xí)慣花十秒鐘就能養(yǎng)成但它能幫你省下一整晚的排查時間。然后打開 MATLAB用cd切到該目錄運行depfun(main_avo.m) % 查看主腳本依賴的所有函數(shù)或者直接在編輯器里打開主腳本逐個點一下函數(shù)名看能否跳轉(zhuǎn)到對應(yīng)文件。如果發(fā)現(xiàn)有函數(shù)名標(biāo)紅找不到定義優(yōu)先檢查是不是子文件夾沒有加進路徑。用addpath(genpath(pwd))一次性把當(dāng)前目錄及所有子目錄加進搜索路徑是解決“函數(shù)未定義”最粗暴也最有效的辦法。3.2 主程序參數(shù)表哪些參數(shù)必須提前想清楚跑正演之前先把這個模型的“地質(zhì)身份”定好。一個 AVO 正演模型至少需要四組參數(shù)上覆泥巖的 Vp、Vs、密度下伏砂巖的 Vp、Vs、密度入射角范圍從 0 度到多少度輸出方式曲線、道集、交會圖以最常見的含氣砂巖模型為例參數(shù)大致是這樣的參數(shù)上覆泥巖含水砂巖含氣砂巖Vpm/s280030002600Vsm/s120015001500ρkg/m3235023502050Vp/Vs 比2.332.001.73注意含氣砂巖的 Vp 明顯比含水砂巖低但 Vs 幾乎不變這就是氣層導(dǎo)致的“縱波速度下降、橫波速度基本不變”的經(jīng)典流體響應(yīng)也是 AVO 能夠識別流體的底層邏輯。如果你在例程里把這些參數(shù)替換進去直接就能看到含水砂巖頂面的反射振幅隨角度變化緩慢而含氣砂巖頂面的反射振幅隨角度明顯變負。3.3 運行與驗證得到的道集合理嗎主腳本運行后你通常會看到兩類圖一類是反射系數(shù)曲線 R(θ)另一類是合成的 AVO 道集。道集怎么看橫軸是入射角或偏移距縱軸是時間或深度顏色代表振幅。在某個反射界面上如果振幅從左到右小角度到大角度越來越“亮”或越來越“暗”說明這個界面的 AVO 響應(yīng)強烈。跑完第一步先做三件驗證工作零角度處的反射系數(shù)用手算一下波阻抗差確認和曲線起點一致看大角度方向的曲線是否出現(xiàn)異常跳動如果有考慮臨界角效應(yīng)后面專門講把精確解和 Shuey 近似的曲線疊在一起看偏差是否在可接受范圍內(nèi)我自己的習(xí)慣是直接在命令行里打幾個關(guān)鍵值對比一下[R0_zoe] zoeppritz_avo(2800,1200,2350,2600,1500,2050,0); [R0_shu] shuey_approx(2800,1200,2350,2600,1500,2050,0); fprintf(Zoeppritz R0 %.4f, Shuey R0 %.4f\n, R0_zoe, R0_shu);如果這兩個值差超過 0.005說明某個函數(shù)的參數(shù)順序或者符號約定有問題先修這個再往下走。4. 結(jié)果解讀4類AVO異常和截距-梯度交會圖4.1 含氣砂巖在道集上長什么樣跑出第一張 AVO 道集之后最想知道的當(dāng)然是這個結(jié)果到底能不能說明地下含氣這里需要引入 Rutherford and Williams1989提出的含氣砂巖 AVO 分類框架。這個分類雖然老但現(xiàn)在工業(yè)界解釋疊前道集時仍然天天在用類型含氣砂巖阻抗法向入射反射系數(shù)振幅隨角度變化特征1類高阻抗比泥巖硬正值振幅先減后增可能出現(xiàn)極性反轉(zhuǎn)2類近零阻抗接近零反射很弱極性反轉(zhuǎn)常見3類低阻抗比泥巖軟負值振幅絕對值隨角度增大4類低阻抗更特殊負值振幅絕對值隨角度減小用上面那組含氣砂巖參數(shù)算出來的是典型的第 3 類法向反射系數(shù)為負并且隨入射角增大振幅的絕對值越來越大。對應(yīng)的圖形特征是道集上這個反射軸的“亮度”從左到右越來越強而且是負極性先負后正或先黑后白取決于顯示約定。為什么第 3 類最常見因為絕大多數(shù)淺層、中深層含氣砂巖都比圍巖泥巖更“軟”——縱波速度低、密度低導(dǎo)致阻抗差本來就很大再加上泊松比降低橫波速度差異相對小于是角度項進一步把負振幅拉大。你如果看到自己的正演結(jié)果居然在 20 度以后振幅往回縮那要看是不是參數(shù)里給出了異常的 Vs 值或者密度壓得太低。4.2 從正演到AVO屬性P-G交會圖怎么用例程包如果夠完整里面多半還有一個函數(shù)用來擬合法向入射截距 P 和梯度 G。做法很簡單對反射系數(shù)序列做最小二乘擬合theta_deg 0:0.5:30; Rpp zoeppritz_avo(vp1,vs1,rho1,vp2,vs2,rho2,theta_deg); A [ones(length(theta_deg),1), sin(theta_deg*pi/180).^2]; coef A \ Rpp(:); P coef(1); % 截距 G coef(2); % 梯度得到 P 和 G 之后把不同模型含水、含油、含氣的正演結(jié)果放到同一個 P-G 交會圖里你會看到它們分布在不同的象限或區(qū)域。典型含氣砂巖的 P×G 為正的負值區(qū)域第三象限或沿著負 P 負 G 方向含水砂巖則更靠近坐標(biāo)原點或正向區(qū)域。這個交會圖是 AVO 解釋里最有名的“甜點探測器”正演的意義就在于你知道一個真實氣藏對應(yīng)的 P、G 應(yīng)該在哪個位置再看實際數(shù)據(jù)的 P、G 點是否落進來。5. 實際跑代碼時最容易翻車的三個地方5.1 DLL初始化失敗可能是路徑和運行庫的問題很多人在 MATLAB 里調(diào)用外部代碼或 MEX 文件時會撞見類似這樣的報錯OSError: [WinError 1114] 動態(tài)鏈接庫(DLL)初始化例程失敗。Error loading D:...\xxx.dll這個錯誤我見到太多次了它在 Windows MATLAB 環(huán)境下高發(fā)原因通常不是代碼邏輯而是系統(tǒng)層面的 DLL 加載問題。最常見的誘因有三個路徑里有中文或空格導(dǎo)致 DLL 依賴的本地資源找不到目標(biāo) DLL 依賴的 Visual C 運行庫缺失需要裝 vc_redist.x64.exe殺毒軟件把 DLL 隔離或攔截了加載時初始化函數(shù)無法執(zhí)行排查建議按順序來先把整個工程目錄挪到D:\codes\這種純英文路徑再確認 MATLAB 的位數(shù)matlab -arch和你調(diào)用的 DLL 位數(shù)一致最后用Dependencies之類的工具打開 DLL看缺失的依賴項。不要一上來就懷疑 MATLAB 安裝壞了大多數(shù) WinError 1114 都是環(huán)境問題。5.2 矩陣維度報錯與復(fù)數(shù)結(jié)果臨界角沒有處理Zoeppritz 求解里asin(p * vp2)可能算出復(fù)數(shù)因為射線參數(shù) p sin(θ1)/vp1 是固定值當(dāng)入射角增大到一定程度時p * vp2 1導(dǎo)致反正弦函數(shù)的定義域越界。這個入射角就是臨界角。超過臨界角后透射波會變成非均勻波折射回介質(zhì)內(nèi)部反射系數(shù)在臨界角附近會出現(xiàn)劇烈的振幅變化。如果你不加處理直接把這個復(fù)數(shù)結(jié)果拿去畫道集圖里就會出現(xiàn)一撮“毛刺”或者 NaN 空洞。解決辦法是在循環(huán)里檢查abs(p * vp2)超過 1 就做截斷或直接丟棄該角度同時在道集繪制時限制最大顯示角度。經(jīng)驗法則是最大入射角取臨界角的 80% 左右既能保證信息量又不會讓臨界角附近的噪聲干擾注意力。5.3 符號約定不統(tǒng)一先和解析解比對再往下走我前面提到過 Zoeppritz 矩陣符號亂的問題這里再展開。不同代碼庫、不同論文里對“位移正方向”“反射系數(shù)極性”的定義經(jīng)常不一致。最典型的例子是同一個地質(zhì)模型你在 A 例程里算出的 3 類 AVO 道集是“負黑正白”在 B 例程里可能完全反過來。如果你拿自己的結(jié)果和別人的圖對比發(fā)現(xiàn)極性反了先別懷疑地質(zhì)參數(shù)先檢查是不是符號約定不同。怎么快速自查用最簡單的地質(zhì)界面上覆是高速高密度下伏是低速低密度計算垂直入射反射系數(shù)。正常約定下它應(yīng)該是負值即反射波與入射波相位相反。如果你的代碼算出正號那么整個道集的顏色約定就要整體取反或者你在繪圖時有意反轉(zhuǎn)了極性。把這個驗證腳本寫進例程的頭部注釋里能救很多人的命。6. 從單道正演到合成道集擴展例程的進階思路6.1 用Ricker子波做褶積生成合成地震道只畫反射系數(shù)曲線在正演層面雖然夠用但地震解釋人員看的是“地震道集”也就是反射系數(shù)經(jīng)過子波褶積后的結(jié)果。擴展例程很自然的下一步就是把 Ricker 子波和反射系數(shù)做褶積生成更接近真實地震記錄的道集。t 0:0.001:0.8; w ricker_wavelet(30, 0.001); % 30Hz Ricker 子波需要自備或自己寫 r_trace zeros(size(t)); r_trace(200) Rpp(1); % 假設(shè)一個界面在 0.2s for i 2:length(theta_deg) r_trace_i zeros(size(t)); r_trace_i(200) Rpp(i); synth(:,i) conv(w, r_trace_i, same); end這里最關(guān)鍵的是“道集上每個角度的子波波形要保持一致”否則你觀察到的振幅變化可能只是子波旁瓣的干涉結(jié)果而不是真實的 AVO 響應(yīng)。實際地震道集在近偏移距和遠偏移距上的子波會因為動校正拉伸而有差異這是另一個處理環(huán)節(jié)的問題正演階段可以暫時忽略但心里要有數(shù)。6.2 從正演走向反演AVO屬性提取與流體識別正演例程跑熟之后你會自然想到一個應(yīng)用如果我從合成道集或?qū)嶋H道集里擬合出 P 和 G能不能反推下伏巖層的彈性參數(shù)這一步就是從正演到反演的橋梁。常用的套路是利用 P、G 組合出流體因子比如流體因子 F P G某些物性條件下與含氣飽和度相關(guān)性好泊松比變化率 Δσ 的近似公式λρ、μρ 彈性參數(shù)反演我在實際項目里比較常用的是把 P-G 交會圖和流體替換結(jié)果結(jié)合先用 Gassmann 方程把同一個砂巖分別替換成含水、含油、含氣三種狀態(tài)正演出三組 P-G 點再把實測數(shù)據(jù)的 P-G 點投影到圖上看落在哪個流體附近。這樣一來正演就不再是單純畫幾條曲線而是直接參與儲層流體判別的決策鏈。Gassmann 流體替換的簡化實現(xiàn)并不復(fù)雜核心是把巖石骨架的體積模量從含水狀態(tài)換算到目標(biāo)流體狀態(tài)再重新算縱波速度K_sat1 rho1 * (vp1^2 - 4/3 * vs1^2); K_sat2 ...; % 帶入目標(biāo)流體參數(shù) vp2 sqrt((K_sat2 4/3 * vs2^2) / rho2);注意這里的單位要統(tǒng)一密度用 kg/m3速度用 m/s模量單位就是 Pa。我踩過一次坑密度用 g/cm3、速度用 km/s算出來的 K 小了 10 的 6 次方倍所有速度更新全錯。建議在腳本開頭強制做單位轉(zhuǎn)換把所有參數(shù)統(tǒng)一成國際單位制再計算。最后再分享一個我在實際使用中的體會拿到 avoMODING 這類例程包別急著貪多求全先把 Zoeppritz 精確解、Shuey 近似、單界面道集這三樣?xùn)|西跑明白比下載十個擴展包都管用。我剛接觸 AVO 那會兒曾在臨界角處理上栽過跟頭畫出過一條“振幅先增后減又暴增”的道集后來發(fā)現(xiàn)就是 p×vs2 越界導(dǎo)致的復(fù)數(shù)傳播?，F(xiàn)在我的習(xí)慣是每次修改參數(shù)后都固定輸出一組與解析解對比的驗證數(shù)值一旦結(jié)果偏離預(yù)期立刻回溯是參數(shù)問題還是代碼問題而不是埋頭在圖上找原因。希望這篇拆解能幫你省掉那些我已經(jīng)替你踩過的坑。本文還有配套的精品資源點擊獲取