簡介這是一套面向地球物理、工程波動模擬初學者的初步虛譜法偽譜法MATLAB程序用于在復雜介質(zhì)中模擬彈性波傳播兼顧譜方法的高精度與有限差分式的直接求解適合地震波、聲波和地下結(jié)構(gòu)探測等應用場景。壓縮包內(nèi)共2個m文件整體僅3KB均為可直接運行的MATLAB源代碼包含計算網(wǎng)格建立、材料參數(shù)設置、初始波場與邊界條件配置、波動方程求解及結(jié)果可視化等基礎功能模塊。程序基于快速傅里葉變換FFT實現(xiàn)用戶可按需調(diào)整網(wǎng)格密度、時間步長與物性參數(shù)從而適配不同研究目標。目前已有179人學習下載適合需要快速入手彈性波數(shù)值模擬的科研人員和工程師通過閱讀和修改源碼可進一步結(jié)合具體模型開展地震波傳播、地下探測等深入模擬研究。1. 初步虛譜法程序彈性波模擬選偽譜法而不是差分法的關(guān)鍵理由做彈性波正演模擬時大多數(shù)人第一步會想到有限差分成熟、資料多、隨手就能找到全套代碼。但模型稍微大一點差分法的代價立刻顯形——每個最小波長要放10到15個網(wǎng)格點三維模型一跑就是幾天起步。偽譜法也叫虛譜法改用FFT在波數(shù)域里對空間求導一個正弦分量理論上兩個網(wǎng)格點就能表示實際取4到5個點波場干凈程度就能超過八階差分這是它在彈性波模擬里最值錢的地方。這個“初步虛譜法程序”壓縮包就是一條偽譜法彈性波正演的完整落地路徑。下面按“原理→跑通→調(diào)參→避坑→驗證”的順序把這條路線講透適合想用粗網(wǎng)格換高精度、又不想反復調(diào)數(shù)值頻散的從業(yè)者。2. 偽譜法原理與彈性波方程離散為什么粗網(wǎng)格能換來高精度2.1 有限差分的分辨率瓶頸與偽譜法的替代思路偽譜法的本質(zhì)是把空間導數(shù)的計算從網(wǎng)格局部挪到波數(shù)域全局。有限差分算子無論階數(shù)多高本質(zhì)上是對Taylor展開的截斷。八階差分在波數(shù)較低時接近理想導數(shù)一旦波數(shù)逼近Nyquist它的振幅響應就會明顯偏離理想的ik——體現(xiàn)到波場里就是數(shù)值頻散高頻分量速度變慢或變快波前面出現(xiàn)拖著尾巴的振蕩。要壓住這種頻散只有加密網(wǎng)格這一條路而加密網(wǎng)格意味著內(nèi)存和計算量按模型維度的次方增長。偽譜法繞開了這個限制。它的做法是對波場做FFT正變換在波數(shù)域把每個譜分量乘上ik或者所需的任意階導數(shù)算子再反變換回空間域。FFT對正弦分量是全精度的最大可表示波數(shù)就是Nyquist波數(shù)π/dx所以理論上每個波長兩個網(wǎng)格點就能精確表示一個正弦波。實際模擬中取4到5個點/波長是為了照顧震源附近的奇異性和時間離散誤差但已經(jīng)比差分法少一半以上的網(wǎng)格。彈性波模擬尤其吃這個紅利。模型里P波和S波速度差異明顯Vp/Vs通常在根號二到根號三之間S波波長只有P波的一半左右。差分法為了保證S波不出頻散整個網(wǎng)格都要按S波最短波長加密而偽譜法在最稀疏的網(wǎng)格上也能同時分辨兩種波這是它在彈性波模擬里一直被保留的原因。對只需要做二維兩層模型驗證的場景來說這個優(yōu)勢更直接網(wǎng)格從300×300降到150×150內(nèi)存少了四倍單步耗時也大幅下降。空間離散方式每波長網(wǎng)格點最大精確波數(shù)頻散特征單步計算量二階差分20~30有限強頻散需極密網(wǎng)格小八階差分10~15較高輕微頻散中偽譜法4~5Nyquist無空間頻散每次求導兩次FFT順帶說一個檢索層面的坑偽譜法還有個別名叫虛譜法二者都是pseudo-spectral的不同譯法代碼結(jié)構(gòu)完全一致。看到“虛譜”別以為是另一個技術(shù)家族在文獻和程序包里兩個詞混用的情況非常普遍。2.2 彈性波方程用一階速度-應力形式寫比二階位移形式更順手偽譜法可以作用在二階位移方程上但工程上我更推薦一階速度-應力方程組。原因有三個二階方程里出現(xiàn)對x和z的混合二階偏導偽譜法雖然也能算但邊界條件和震源加載的物理意義不如一階直觀一階方程里每個空間導數(shù)都是對單軸的代碼結(jié)構(gòu)規(guī)整不容易寫錯時間上可以直接用二階中心差分做跳蛙遞推存儲量只有五個變量。方程寫出來是下面這樣五個未知量分別是水平速度vx、垂直速度vz以及三個應力分量σxx、σzz、σxzrho ?vx/?t ?σxx/?x ?σxz/?z rho ?vz/?t ?σxz/?x ?σzz/?z ?σxx/?t (λ2μ) ?vx/?x λ ?vz/?z ?σzz/?t λ ?vx/?x (λ2μ) ?vz/?z ?σxz/?t μ ?vx/?z μ ?vz/?xλ和μ是拉梅參數(shù)由Vp、Vs和密度換算λρ(Vp2?2Vs2)μρVs2。網(wǎng)格模型只要給每個點填上Vp、Vs、ρ三個量再逐點換算成λ和μ遞推里需要的所有系數(shù)就齊了。這里有個容易踩的換算細節(jié)有些初步程序直接以λ2μ和μ的形式存參數(shù)省去每步除法有的則是每步都算。前者快很多后者代碼易讀但耗時。模擬前先確認參數(shù)文件里的“vp”“vs”“rho”是模型數(shù)組還是標量以及有沒有做速度到拉梅參數(shù)的換算很多結(jié)果怪異的問題都出在這一步。時間遞推用跳蛙格式即速度在n1/2時刻、應力在n時刻交錯更新。它是二階精度的空間誤差由偽譜法控制在幾乎為零時間誤差就成了總誤差的主要來源。如果要做長時間模擬可以換四階Runge-Kutta但每步要算四次導數(shù)場成本高很多初步程序保持二階中心差分即可。2.3 波數(shù)域求導算子整個偽譜法程序的核心就這一段把空間導數(shù)封裝成一個函數(shù)后續(xù)所有遞推都復用它。Python實現(xiàn)如下import numpy as np def spectral_derivative(field, dx, axis0): 沿指定軸對場做波數(shù)域一階求導。 以二維波場形狀 (nz, nx) 為準 axis0 對應 z 方向間距為 dzaxis1 對應 x 方向間距為 dx。 nx field.shape[axis] # 角波數(shù)向量fftfreq 返回頻率索引乘 2*pi 后是角波數(shù)單位 rad/m k 2.0 * np.pi * np.fft.fftfreq(nx, ddx) # 把波數(shù)向量廣播到 field 的目標軸 shape [1] * field.ndim shape[axis] nx k k.reshape(shape) # 正變換、在波數(shù)域乘 i*k、反變換取實部 derivative np.fft.ifft( np.fft.fft(field, axisaxis) * (1j * k), axisaxis ).real return derivative這段的要點有三個。第一fftfreq(nx, ddx)返回的頻率索引從0到nx/2再到負半軸乘2π之后正好是角波數(shù)如果程序里FFT庫返回的是循環(huán)頻率而非角頻率乘的因子要相應調(diào)整。第二乘的是1jk這是頻域求導的傅里葉變換性質(zhì)如果要求二階導改成(1jk)**2即可偽譜法求高階導數(shù)就是一次FFT的事這也是它區(qū)別于差分法的重要特性。第三反變換后必須取實部——由于浮點誤差ifft會帶回極小的虛部直接參與遞推會被逐時間步放大最終污染整個波場。如果你拿到的是Fortran版本核心邏輯一模一樣先調(diào)用FFT庫做正變換把實數(shù)組轉(zhuǎn)成復數(shù)譜乘上虛數(shù)單位乘波數(shù)再逆變換取實部。區(qū)別只在于FFT庫的布局約定比如某些庫返回的是物理排列的實部虛部需要先做fftshift數(shù)值實現(xiàn)不復雜但移植時最容易在這些地方翻車。3. 把初步虛譜法程序跑起來文件確認、環(huán)境準備與最小兩層算例3.1 解壓之后先確認四類文件缺了別急著跑一個典型的初步偽譜法程序包解壓后通常包含四類東西主程序源碼可能是Fortran的.f90、Python的.py或Matlab的.m參數(shù)定義要么是獨立的文本/配置塊要么寫在主程序開頭的常量區(qū)輸出與繪圖腳本把模擬結(jié)果寫成二進制或文本的地震記錄以及一個模型/算例目錄。如果壓縮包里帶README先看README的“運行方式”一節(jié)那里會寫明預期的輸出文件名和物理單位。沒有README是常態(tài)。我拿到這類包一般先按文件大小排個序最大的多半是結(jié)果或模型數(shù)據(jù)文件最小且能直接讀的才是可執(zhí)行入口。用編輯器打開主程序先搜“main”或“program”找到時間遞推主循環(huán)的位置再搜“parameter”或“const”把網(wǎng)格尺寸、時間步長、震源位置這幾組常量抄出來。這一步花十分鐘后面能省下幾小時的翻車排查。環(huán)境方面最常出現(xiàn)的坑是終端直接報“gfortran不是內(nèi)部或外部命令”“conda不是內(nèi)部或外部命令”這類信息。它的本質(zhì)是編譯器或Python解釋器的路徑?jīng)]加入系統(tǒng)PATH而不是程序本身有問題。Windows下我建議統(tǒng)一裝Anaconda并創(chuàng)建一個專門環(huán)境裝好numpy和scipyFortran代碼則用gfortran編譯確保編譯器和運行時庫都是64位。32位和64位混用鏈接階段大概率會報“無法定位程序輸入點getcurrentpackagefullname”之類的動態(tài)庫錯誤這類報錯基本都和位數(shù)不匹配有關(guān)。3.2 最小兩層模型一套立刻能用的參數(shù)為了驗證程序能跑不用上來就上一個真模型我用一個兩層介質(zhì)模型上層2000m/s下層3000m/s橫波速度按根號三比例對應。網(wǎng)格200×200網(wǎng)格間距10米震源用20Hz的Ricker子波、垂直集中力放在深度500米處。記錄時長1.5秒時間步長0.5毫秒。參數(shù)值選取理由網(wǎng)格 nx×nz200×200兩層模型只驗證物理過程夠用即可dxdz10 mS波最短波長約57.8m約5.8點/波長上層 Vp/Vs/ρ2000 / 1155 / 2000 kg/m3Vp/Vs√3接近真實沉積巖比例下層 Vp/Vs/ρ3000 / 1732 / 2200 kg/m3界面反射系數(shù)適中便于觀察界面深度1000 m給反射波留出清晰的走時窗口震源Ricker20 Hz垂直集中力集中力同時激發(fā)P波和S波震源位置x1000 mz500 m離頂面和邊界都足夠遠dt0.5 ms約為二維穩(wěn)定極限的1/3偏保守記錄長度1.5 s反射波有足夠時間回到地表這里的關(guān)鍵是網(wǎng)格間距和震源主頻的匹配。20Hz主頻對應上層橫波波長約57.8m10m網(wǎng)格每波長約5.8個點滿足偽譜法4到5點的經(jīng)驗要求。如果把主頻提到40Hz最短波長降一半網(wǎng)格間距就要縮到5m左右計算量翻四倍這個權(quán)衡在第4章還會展開。3.3 主循環(huán)跳蛙遞推的順序不能寫反拿到程序后主循環(huán)通常是這樣的結(jié)構(gòu)我把它重寫成一個盡量貼近各類初步程序的Python版本# 偽譜法彈性波模擬主循環(huán)跳蛙格式二階時間差分 # 數(shù)組形狀統(tǒng)一為 (nz, nx)axis0 是深度 zaxis1 是水平 x for it in range(nt): # 第一步由應力更新速度分量 vx dt / rho * ( spectral_derivative(sxx, dx, axis1) # ?σxx/?x spectral_derivative(sxz, dz, axis0) # ?σxz/?z ) vz dt / rho * ( spectral_derivative(sxz, dx, axis1) # ?σxz/?x spectral_derivative(szz, dz, axis0) # ?σzz/?z ) # 在震源位置加載垂直集中力源只加在 vz 分量 vz[nsz, nsx] dt / rho[nsz, nsx] * wavelet[it] # 第二步由速度更新應力分量 sxx dt * ( (lam 2.0 * mu) * spectral_derivative(vx, dx, axis1) lam * spectral_derivative(vz, dz, axis0) ) szz dt * ( lam * spectral_derivative(vx, dx, axis1) (lam 2.0 * mu) * spectral_derivative(vz, dz, axis0) ) sxz dt * mu * ( spectral_derivative(vx, dz, axis0) # ?vx/?z spectral_derivative(vz, dx, axis1) # ?vz/?x ) # 第三步應用吸收邊界第4章展開 # 第四步在接收點處把 vx/vz 寫入記錄道注意這里的存儲細節(jié)。vx代表水平振動速度vz代表垂直振動速度nsz是深度索引nsx是水平索引。加載垂直集中力時改的是vz而不是vx否則輻射圖會繞著一個錯誤的軸轉(zhuǎn)。如果震源是爆炸源則應該同時往sxx、szz、sxz上加各向同性壓力而不是直接改速度分量——很多初步程序把爆炸源實現(xiàn)成“往所有點加同一個速度擾動”得到的結(jié)果看著有波但波型比例完全錯誤。時間遞推的順序是先更新速度再更新應力還是反過來其實可以互換只要震源加在正確的位置、并保持交錯時刻的一致性。但每個時間步內(nèi)部順序要統(tǒng)一先算完所有速度分量再算所有應力分量不能混著來否則時間同步被打破高頻成分會迅速失穩(wěn)。上面的寫法重在清晰效率不是最優(yōu)。spectral_derivative每調(diào)用一次就是一次FFT加一次逆FFT這個循環(huán)里一共調(diào)用了12次其中對vx的x方向?qū)?shù)和vz的z方向?qū)?shù)在速度更新和應力更新里重復算了。優(yōu)化時可以先把六個一階導數(shù)場一次性算好再組裝應力更新整體能省掉約1/3的FFT開銷。初步程序不追求性能但這個邏輯值得記著后續(xù)做三維擴展時會用到。3.4 跑通后的第一道驗收直達波與反射波的到達時間跑完之后先看接收器輸出的兩組記錄。vz記錄上第一個到達的是直達P波初走時約等于震源到接收點的距離除以上層縱波速度隨后會看到來自界面的反射P波和反射轉(zhuǎn)換波。如果vz上和vx上除了直達波外什么都沒有檢查震源類型和界面兩側(cè)波阻抗差——速度差太小也會讓反射系數(shù)低到看不見這時加大兩層速度比再試。一個快速的手工驗算是把震源到界面的垂直距離和接收點的水平距離代入初等幾何關(guān)系算出反射P波的走時再與程序輸出的記錄道對比。以第3.2節(jié)的參數(shù)為例震源深500m、界面在1000m、接收點水平距離100m時反射P波路徑長約1503m按上層Vp2000m/s算走時約0.75秒直達P波走時約0.255秒。誤差在1到2毫秒以內(nèi)說明程序核心邏輯基本正確超過這個量就要回去檢查網(wǎng)格方向或介質(zhì)參數(shù)是否裝反了。4. 三個必調(diào)參數(shù)時間步長、吸收邊界與震源子波改錯了就翻車4.1 時間步長偽譜法的穩(wěn)定極限不是差分法那個公式偽譜法的空間導數(shù)沒有頻散誤差但這不意味著可以無腦用大時間步長。如果時間差分仍然是二階中心差分穩(wěn)定性條件來自最大可表示的波數(shù)k_maxπ/dx與介質(zhì)最大波速vmax的乘積。一維情況下理論極限約為0.637·dx/vmax二維時波數(shù)向量可以沿對角方向疊加k_max變?yōu)棣小?/dx極限步長縮到約0.45·dx/vmax三維更嚴約0.37·dx/vmax。偽譜法能精確表示到Nyquist波數(shù)而差分法在高波數(shù)部分的振幅響應實際上是衰減的相當于天然濾掉了一部分不穩(wěn)定成分所以偽譜法對時間步長更敏感。我一般不會頂著極限值用而是取二維極限的一半左右dt 0.3·dx/vmax。這樣既留出安全余量又不會因為步長太小讓長時程模擬的步數(shù)猛增。以第3章那個兩層模型為例vmax取下層縱波速3000m/sdx10m二維穩(wěn)定極限約1.5毫秒取0.5毫秒是極限的1/3屬于穩(wěn)妥選擇。如果壓縮包代碼里時間步長是寫死的先按這個公式重新算一遍再跑。判斷步長是否過大不一定要等波場爆炸。最快的診斷方法是打印每一時間步的總能量在均勻無吸收模型里總能量應當基本守恒。如果看到某個分量能量隨步數(shù)單調(diào)上升比如從1e-2漲到1e0基本可以斷定步長越過穩(wěn)定極限。把dt縮小到原來的1/4再跑能量曲線趨于平穩(wěn)就說明問題出在此處而非程序邏輯。提示步長的大小對偽譜法的影響是“全有或全無”的越界一步就會在幾十步內(nèi)爆掉。養(yǎng)成每個新模型先跑50步看能量的習慣比跑完整個記錄才發(fā)現(xiàn)翻車要省時得多。4.2 吸收邊界阻尼帶的厚度和衰減系數(shù)要一起調(diào)初步程序很少帶PML最常見的是在計算域四周加一層阻尼帶也叫海綿邊界或吸收層。它的原理很簡單每時間步對邊界區(qū)的波場乘一個小于1的衰減因子讓波在到達人工邊界前衰減到可忽略。實現(xiàn)不難但參數(shù)配不對時阻尼帶本身就會變成反射源效果比不加還糟。阻尼系數(shù)一般取成空間位置的函數(shù)例如σ(x)σ_max·(x/L)2其中L是阻尼帶的網(wǎng)格數(shù)x是該點到計算域邊界的歸一化距離。σ_max的經(jīng)驗范圍是2到3倍的vmax/(L·dx)。L的取值至少要覆蓋一個中心波長中心波長用震源主頻對應的波長來算λ_cvmax/f0。在20Hz主頻、3000m/s最大速度的模型里中心波長150米L建議取15到20個網(wǎng)格dx10m時。L太薄時波在阻尼帶內(nèi)還沒衰減到位就撞到硬邊界反射能量依舊可觀。給一段阻尼帶實現(xiàn)可以直接替換第3.3節(jié)主循環(huán)里的“第三步”# 生成二維阻尼衰減系數(shù)場四個邊界各加 L 個網(wǎng)格 def build_damper(nz, nx, L, vmax, dt): sig_max 3.0 * vmax / (L * dx) # 單位 1/sL*dx 是帶的總長度米 damp np.ones((nz, nx), dtypenp.float64) for i in range(L): factor sig_max * ((i 1) / L) ** 2 * dt damp[i, :] * np.exp(-factor) # 上邊界 damp[-(i 1), :] * np.exp(-factor) # 下邊界 damp[:, i] * np.exp(-factor) # 左邊界 damp[:, -(i 1)] * np.exp(-factor) # 右邊界 return damp # 每個時間步在遞推之后執(zhí)行 vx * damp vz * damp sxx * damp szz * damp sxz * damp注意角點區(qū)域會被重復衰減這個實現(xiàn)在角點的衰減系數(shù)比邊上大一倍實際影響不大如果要嚴格處理需要按到最近邊界的距離分別計算x和z方向的衰減因子再相乘。更重要的是阻尼帶內(nèi)最好保持常數(shù)速度模型不要放界面或強速度梯度否則波在帶內(nèi)產(chǎn)生反射這部分反射同樣會污染內(nèi)部波場。4.3 震源子波Ricker子波的主頻和網(wǎng)格間距是配對關(guān)系震源子波最常用Ricker表達式是f(t)(1?2π2f?2(t?t?)2)exp(?π2f?2(t?t?)2)其中t?一般取1.2到1.5個主頻周期讓子波初始時刻接近零避免在t0時刻給波場一個階躍激勵。實現(xiàn)如下# Ricker 子波f0 為主頻dt 為時間步長 t np.arange(nt) * dt t0 1.2 / f0 wavelet (1.0 - 2.0 * (np.pi * f0 * (t - t0)) ** 2) * \ np.exp(-(np.pi * f0 * (t - t0)) ** 2)主頻f?越高波場分辨率越高能分辨更薄的層但代價是S波最短波長同步變短需要更細的網(wǎng)格。經(jīng)驗約束是每個最短波長至少要有4到5個網(wǎng)格點即dx ≤ v_s_min/(4·f?)。這里速度取整個模型里最小的S波速度因為S波波長最短最容易頻散。以第3章模型為例上層Vs1155m/sf?20Hz時最短波長約57.8mdx10m相當于每波長約5.8個點處于安全區(qū)間。如果把主頻從20Hz提到40Hz最短波長降一半dx就必須縮到5m左右計算量漲四倍這就是主頻和網(wǎng)格步長的直接權(quán)衡。如果壓縮包默認震源是爆炸源而你需要同時看P波和S波換成垂直集中力源即可。爆炸源只會輻射純縱波無論后來怎么調(diào)吸收邊界和網(wǎng)格橫波分量始終是零這一點在驗證環(huán)節(jié)最容易把人帶偏。震源加載位置建議離邊界至少10個網(wǎng)格否則即使有阻尼帶源與人工邊界之間的多次反射也會干擾早期波場。5. 偽譜法程序避坑指南5個最常見的翻車現(xiàn)場與排查方法5.1 波場圖上一片棋盤格噪聲高頻Nyquist分量在作怪現(xiàn)象模擬幾步后波場圖出現(xiàn)顆粒狀交替亮暗的棋盤格尤其在震源附近最明顯振幅隨步數(shù)增長。原因單點加載震源在空間上是一個極窄的尖峰它的頻譜在Nyquist波數(shù)附近仍然有可觀的能量。偽譜法對這個分量是全精度放大的不像差分法有天然的抑制于是波場里出現(xiàn)以單個網(wǎng)格為周期的交替擾動視覺上就是棋盤格。解決把震源先做空間平滑再乘子波。常見做法是給震源區(qū)一個高斯半徑比如σ_source1.5倍的dx讓源在空間上分布到8到10個網(wǎng)格點同時檢查FFT后是否取了實部虛部殘留也會產(chǎn)生類似的高頻噪聲。如果程序本身沒有平滑函數(shù)可以在加載震源前對相鄰網(wǎng)格按高斯權(quán)重分配能量。5.2 邊界反射比預期早出現(xiàn)阻尼帶沒蓋住最大波長現(xiàn)象波場圖上在計算域邊界附近出現(xiàn)強反射弧反射波到達內(nèi)部接收點的時間明顯早于模型里真實界面的理論走時。原因阻尼帶厚度L沒有按最大中心波長設計。L太薄時長波長成分在帶內(nèi)衰減不夠振幅在到達硬邊界時仍然可觀邊界反射自然回傳。解決把L加大到至少一個中心波長。用vmax/f0算出中心波長后再換算成網(wǎng)格數(shù)如果程序里阻尼帶厚度寫死改參數(shù)或預處理速度模型時把邊界區(qū)擴展。驗證方法是給一個無反射界面的均勻模型跑一次把接收點能量畫成時間曲線觀察末段是否有明顯長時間拖尾的反射能量。阻尼帶的σ_max也要同步調(diào)到2到3倍vmax/(L·dx)薄帶配大衰減、厚帶配小衰減兩種組合效果不同需要交叉驗證。5.3 振幅隨時間指數(shù)增長直到NaN時間步長越過穩(wěn)定極限現(xiàn)象前面的波形看著正常到幾百步之后某個應力分量量級從1e-2跳到1e20甚至直接變成NaN程序掛掉。原因按照4.1節(jié)算出的單方向穩(wěn)定條件只是一維理論在二維模型里波動能量沿多個方向傳播實際允許的步長通常更小。很多初步程序的dt是作者用他的模型試出來的換到你自己的網(wǎng)格尺寸和速度模型后穩(wěn)定余量可能已經(jīng)不夠。解決把dt縮小到當前值的一半甚至1/4重跑看是否仍然發(fā)散。同時建議在時間循環(huán)里加一個能量檢測每50步打印一次波場總能量看到指數(shù)上升就立即終止避免跑完整個記錄長度才發(fā)現(xiàn)翻車、白燒算力。穩(wěn)定步長與dx、vmax的具體取值參考4.1的公式但最終以你的模型能量曲線為準這是這類程序最不可省的一步基本功。5.4 橫波分量離奇失蹤震源類型和參數(shù)化把S波滅掉了現(xiàn)象接收記錄上只有縱波初至之后全是微弱的低頻尾巴理論上應當明顯的反射轉(zhuǎn)換波消失vx分量尤其干凈。原因兩類常見誤操作。一是用爆炸源加載它只激發(fā)P波S波天然為零二是參數(shù)換算時把μ設成了0或很小的值導致S波速度接近0波場根本傳播不出去。解決換成垂直集中力源加載在vz分量上同時檢查拉梅參數(shù)換算μρVs2如果模型文件里Vs列填了0或沒填μ就會變成0。一張快速自檢圖是把Vp、Vs畫成按深度的曲線看Vs站點是否與Vp同步變化若Vs全程為0程序里再聰明也算不出S波。5.5 程序在Windows下報動態(tài)庫或命令找不到環(huán)境沒有對齊現(xiàn)象終端執(zhí)行編譯命令時報“gfortran不是內(nèi)部或外部命令”運行Python時報“numpy模塊不存在”或者程序啟動直接報“無法定位程序輸入點getcurrentpackagefullname于動態(tài)鏈接庫…”運行就中斷。原因三類問題混在一起——編譯器或解釋器的PATH沒有配好、Python環(huán)境不對、以及32位/64位運行時庫混用。后者在下載了舊版編譯好的現(xiàn)成程序包時最容易出現(xiàn)因為動態(tài)鏈接庫的位數(shù)和主程序不匹配系統(tǒng)加載時就報找不到入口點。解決Fortran源碼重新用本地gfortran編譯別直接用網(wǎng)上別人編好的exePython部分統(tǒng)一到Anaconda的64位環(huán)境建環(huán)境后執(zhí)行conda install numpy scipy別用系統(tǒng)自帶的Python。檢查位數(shù)的方法是打開終端分別敲gfortran --version和python --version確認輸出里有沒有帶32位字樣。這一類報錯的排查邏輯和網(wǎng)上常見的“conda不是內(nèi)部或外部命令”完全一樣先確認環(huán)境變量再確認位數(shù)最后才是代碼問題。6. 驗證偽譜法程序正確性解析解對比與網(wǎng)格收斂性檢查寫完代碼、跑通模擬不等于程序是對的。我驗證任何正演程序都走固定的三步解析解走時對比、網(wǎng)格收斂性檢查和能量守恒檢查。這三步能過濾掉九成以上的隱性錯誤。第一步用兩層介質(zhì)模型或均勻半空間模型把接收點的波場與解析走時對比。均勻半空間里直達P波走時是r/Vp直達S波走時是r/Vs兩層模型里反射P波走時按鏡像源法計算公式簡單手算即可。把程序輸出的單道記錄拆成vx和vz兩列找到初至時間誤差在1到2毫秒內(nèi)算通過。嚴格檢查可以再加一個垂直自由表面邊界對比Rayleigh波存在與否但初步程序一般不需要。第二步是網(wǎng)格收斂性檢驗。把dx、dz同時減半dt等比縮小重跑同一個模型對比同一接收點的波形。偽譜法如果實現(xiàn)正確兩次結(jié)果的波形差異應該在1%以內(nèi)且差值主要集中在高頻尾部。如果減半網(wǎng)格后波形明顯變化說明原網(wǎng)格本身就不滿足分辨率要求需要按第4章的公式重新選擇網(wǎng)格間距而不是程序邏輯有問題。第三步是能量監(jiān)測這個前面提過。在沒有阻尼帶和震源持續(xù)加載的均勻模型中總能量應該守恒在帶阻尼帶的模型中能量應單調(diào)衰減而不是振蕩上升。把每步總能量畫出來曲線形狀正常程序才算真正通過驗收。我拿到的每一個偽譜法程序都會先跑這三步再做物理實驗。走時對不上先查震源類型能量發(fā)散了先查時間步長波形不收斂先查網(wǎng)格間距順序不要倒過來。這個習慣幫我擋掉了大量“看起來正常其實參數(shù)錯位”的翻車現(xiàn)場。希望幫到你。本文還有配套的精品資源點擊獲取