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