起伏地形條件下瑞雷面波傳播特性研究

李 宏,楊心超,朱海波,張 偉,曲壽利

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.南方科技大學(xué),廣東深圳518055)

李宏,楊心超,朱海波,等.起伏地形條件下瑞雷面波傳播特性研究[J].石油物探,2017,56(6):-791

LI Hong, YANG Xinchao, ZHU Haibo,et al.Rayleigh wave propagation with undulating topography[J].Geophysical Prospecting for Petroleum,2017,56(6):-791

起伏地形條件下瑞雷面波傳播特性研究

李 宏1,楊心超1,朱海波1,張 偉2,曲壽利1

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.南方科技大學(xué),廣東深圳518055)

采用曲線網(wǎng)格有限差分方法研究了瑞雷面波通過起伏地形時(shí)的傳播特性。討論了地形的凹凸、高(深)度、跨度等因素對面波波形、能量、頻散等參數(shù)的影響。平緩地形對面波影響較小,而地形起伏較大時(shí),則會(huì)產(chǎn)生嚴(yán)重的散射現(xiàn)象。地形起伏對面波有濾波效應(yīng),在頻譜上會(huì)出現(xiàn)陷波點(diǎn),陷波點(diǎn)頻率隨地形的高度與跨度變化。另外,曲線網(wǎng)格有限差分方法適應(yīng)較復(fù)雜地形模型的地震波場數(shù)值模擬,可以應(yīng)用于實(shí)際復(fù)雜地形的地震波場響應(yīng)特征模擬與分析。

曲線網(wǎng)格;有限差分方法;起伏地形;面波模擬

復(fù)雜的近地表環(huán)境對地震波傳播特性影響顯著,地形特征對地震波能量的放大與衰減具有重要作用[1]。實(shí)際觀察與數(shù)值模擬分析發(fā)現(xiàn),2008年汶川地震中山脊對地震波能量的放大作用是造成一些地區(qū)受災(zāi)嚴(yán)重的主要原因[2-3]。地震面波作為近地表波場中重要組成成分,在大地震中對地表造成破壞的同時(shí),也提供了一種了解地球內(nèi)部結(jié)構(gòu)的重要手段。了解地震面波在復(fù)雜地形條件下的傳播特性,不僅可以為防災(zāi)減損提供建議,也可以指導(dǎo)面波勘探觀測系統(tǒng)設(shè)計(jì),盡量避開不利地形[4-5],或者進(jìn)行面波噪聲壓制[6]。

關(guān)于面波地形效應(yīng)的研究已有很多,KNOPOFF[7]和FUJII[8]通過實(shí)驗(yàn)觀測研究了瑞雷面波通過楔形地形時(shí)反射系數(shù)和透射系數(shù)的影響因素;HUDSON等[9]與MAL等[10]從理論上研究了楔形地形對瑞雷面波反射系數(shù)和透射系數(shù)的影響;陳偉等[4]采用勢函數(shù)方法模擬并討論了正、負(fù)楔形體角度對瑞雷面波傳播特征的影響;周紅等[11]采用局域離散波數(shù)法探討了瑞雷面波穿過凹陷地形時(shí)能量、頻率等參數(shù)的變化;汪利民[12]采用交錯(cuò)網(wǎng)格有限差分法并結(jié)合聲學(xué)/彈性界面近似與階梯近似起伏地形方法對瑞雷面波傳播特征進(jìn)行模擬分析,并且討論了不同陡峭程度的起伏地形對瑞雷面波傳播特性的影響[13];巴振寧等[14-15]采用間接邊界元方法對瑞雷面波入射含凸起與凹陷層狀地層模型的響應(yīng)進(jìn)行分析;劉中憲等[16]采用有限元-間接邊界積分耦合方法,對盆山耦合地形二維P波、P-SV波和瑞雷面波的散射波場進(jìn)行了模擬研究。

當(dāng)模型簡單時(shí),可采用解析或半解析方法求解面波記錄。該方法計(jì)算量較小且可以獲得較準(zhǔn)確的解,但當(dāng)復(fù)雜模型時(shí),如存在起伏地表或低速體等情況則方法的適應(yīng)性較差,所以需要采用對模擬區(qū)域格點(diǎn)離散的數(shù)值模擬技術(shù)[17-18]。采用格點(diǎn)離散方式進(jìn)行近地表波場模擬關(guān)鍵是如何實(shí)現(xiàn)地表處的自由邊界條件。有限元方法網(wǎng)格劃分靈活,自由邊界條件實(shí)現(xiàn)也比較容易、直接,但其計(jì)算量較大;而有限差分法求解波動(dòng)方程實(shí)現(xiàn)較容易且計(jì)算量適中。但若采用矩形網(wǎng)格則需要對起伏地表進(jìn)行階梯狀近似,要用較大的網(wǎng)格密度來壓制地表處網(wǎng)格散射[19]。曲線網(wǎng)格為解決起伏地表問題提供了一個(gè)有效手段[20-21]。

本文采用曲線網(wǎng)格有限差分方法研究了瑞雷面波通過起伏地形時(shí)的傳播特性。首先介紹了曲線網(wǎng)格有限差分法結(jié)合牽引力鏡像法處理起伏自由地表的模擬方法;然后討論了地形的凹凸、高(深)度、跨度等對瑞雷面波波形、能量、頻散等參數(shù)的影響。

1 方法原理

起伏地形網(wǎng)格剖分?jǐn)M合有兩種方式:一是ROBERTSSON[19]所采用的階梯狀網(wǎng)格近似方式,汪利民[12]結(jié)合聲學(xué)/彈性界面近似與Robertsson方法對面波傳播特征進(jìn)行模擬分析;二是通過網(wǎng)格變換方式,采用曲線網(wǎng)格來擬合地形的起伏,從而避免了階梯狀近似方法需大量格點(diǎn)壓制虛假散射問題[22]。

本文采用張偉[20]提出的曲線網(wǎng)格(圖1),該網(wǎng)格可以很好地?cái)M合起伏地形。其思想是在物理空間內(nèi)采用與物理邊界重合的曲線網(wǎng)格進(jìn)行離散,然后將網(wǎng)格坐標(biāo)變換到均勻的計(jì)算網(wǎng)格空間,最后在計(jì)算網(wǎng)格空間對波動(dòng)方程進(jìn)行求解。

圖1 曲線網(wǎng)格坐標(biāo)變換示意[20]

在二維笛卡爾坐標(biāo)系中,一階速度-應(yīng)力P-SV波動(dòng)方程表達(dá)式為:

(1)

式中:U為速度-應(yīng)力矢量;F為震源項(xiàng);A,B為系數(shù)矩陣。

(5)

其中,vx,vy為速度分量;σxx,σyy,σxy為應(yīng)力分量;T表示矩陣轉(zhuǎn)置;fx,fy為體力項(xiàng);Mxx,Myy,Mxy分別為地震矩張量;ρ為密度;λ,μ為拉梅常數(shù)。

物理空間與計(jì)算空間的坐標(biāo)映射方程為:

(6)

式中:(x,y)為物理空間坐標(biāo);(ξ,η)為計(jì)算空間坐標(biāo)。通過鏈?zhǔn)角髮?dǎo)法則,可得到在計(jì)算空間坐標(biāo)系中P-SV波動(dòng)方程表達(dá)式:

(7)

交錯(cuò)網(wǎng)格是目前地震波數(shù)值模擬中比較流行的網(wǎng)格結(jié)構(gòu),但其在自由界面格點(diǎn)處應(yīng)力或速度分量缺失,需要在該格點(diǎn)處進(jìn)行插值等處理,因而降低了計(jì)算精度。而本文所采用的曲線網(wǎng)格是基于同位網(wǎng)格,可很好地避免這一點(diǎn),本文空間差分采用低頻散低耗散的DRP/opt MacCormack(Dispersion Relation Preserving/optimized MacCormack)格式[23]:

(10)

(11)

在時(shí)間迭代過程中,采用了四階精度的Runge-Kutta格式[24]。為避免偏心差分算子帶來的頻散和不穩(wěn)定性,在Runge-Kutta格式中交替使用前向與后向差分格式:

(12)

其中,系數(shù)α2=0.5,α3=0.5,α4=1.0;β1=1/6,β2=1/3,β3=1/3,β4=1/6。為方便討論,將(12)式簡寫為:

(13)

對前向和后向差分的各種選擇,(13)式有多種組合方式,本文在時(shí)間迭代過程中采用了4個(gè)計(jì)算時(shí)間步長為一個(gè)循環(huán)的處理方式:

(14)

LEVANDER[25]首先在交錯(cuò)網(wǎng)格中采用應(yīng)力鏡像法對水平自由地表進(jìn)行處理,但該方法并不適應(yīng)于任意曲線網(wǎng)格。ZHANG等[26]將該思想引入到貼體網(wǎng)格中并提出牽引力鏡像法,解決了曲線網(wǎng)格中起伏地表自由邊界條件處理問題。當(dāng)網(wǎng)格在自由地表處正交時(shí),牽引力鏡像法可以退化為應(yīng)力鏡像法[20-21]。本文采用ZHANG等[26]提出的牽引力鏡像法來實(shí)現(xiàn)自由邊界條件。

1) 速度導(dǎo)數(shù)求取。在地表自由邊界處,要求地表外法向牽引力T為零,即:

(15)

式中:n為地表某點(diǎn)的法向矢量;σ為點(diǎn)的應(yīng)力張量。將(7)式代入(15)式整理后可得:

(16)

矩陣X,Y的具體表達(dá)式為:

(17)

(18)

其中,χ=λ+2μ。

由(16)式可以看出,自由地表網(wǎng)格層速度的η方向?qū)?shù)可以通過速度ξ方向?qū)?shù)獲得。對于地表下兩層網(wǎng)格,為避免使用地表以上虛擬格點(diǎn)上的速度值,本文采用4/4-MacCormack緊致差分格式[27]。

2) 應(yīng)力分量導(dǎo)數(shù)的求取。由于在自由地表處格點(diǎn)上牽引力為零:

(19)

利用牽引力鏡像法,自由地表以上3層格點(diǎn)上的應(yīng)力值可以通過(20)式得到:

(20)

式中:jf為自由地表處對應(yīng)的η方向格點(diǎn)坐標(biāo),j=1,2,3。在文獻(xiàn)[26]中詳細(xì)討論對比了采用牽引力鏡像法處理起伏自由地表的結(jié)果與邊界元等方法得到的結(jié)果,驗(yàn)證了該方法的準(zhǔn)確性,本文采用該方法研究地形效應(yīng)對面波傳播過程的影響。

另外,吸收邊界處理方式是影響長時(shí)程面波模擬穩(wěn)定性的一個(gè)重要因素。這里我們采用的是ZHANG等[28]提出的ADE CFS-PML(Auxiliary Differential Equation Complex Frequency Shifted Perfectly Matched Layer)吸收邊界。該方法分析了PML吸收層內(nèi)的自由表面情況,修改PML方程以滿足吸收層的自由表面條件,消除了PML層和自由地表?xiàng)l件不相容導(dǎo)致的不穩(wěn)定問題。

2 地形效應(yīng)模擬與分析

2.1 模型設(shè)計(jì)

模型設(shè)計(jì)的目的是為了討論地形不同的深(高)度和跨度對面波傳播特性的影響。本文所采用的介質(zhì)為泊松體,其縱波速度vP=4000m/s,橫波速度vS=2310m/s,密度ρ=2300kg/m3,瑞雷波速度vR≈2123m/s;采用炸藥震源,地震子波為Ricker子波,主頻為20Hz。為了充分激發(fā)面波,震源位于地表附近(-500m,-80m)處,距地形起伏處有一定距離。這樣可以使瑞雷波與P波能夠充分分離,方便時(shí)間窗口截取瑞雷波進(jìn)行后續(xù)討論;地表檢波器排列間距為20m,如圖2所示。網(wǎng)格離散空間步長Δx≈5m,離散最小波長所用格點(diǎn)數(shù)約為8.5。但由于曲線網(wǎng)格存在一定扭曲,最大空間步長約為5.75m,最小空間步長約為4.25m;時(shí)間步長Δt=0.5ms。

模型的水平距離為7km,深度為1km。地形采用高斯形山(谷)模型,其計(jì)算公式為:

(21)

式中:D0為控制起伏程度(高度或深度);w為控制水平寬度參數(shù)(跨度);x0=3km為起伏地形中心位置;“±”中的“+”表示山形地形,“-”表示谷形地形。以D0=150m,w=150m為例,其形態(tài)如圖2所示。

采用高斯形山(谷)模型主要是為了使自由表面從水平地形到山(谷)體平滑過渡,避免采用圓形或多邊形(三角形或梯形)模型在過渡區(qū)域的拐點(diǎn)。當(dāng)面波通過拐點(diǎn)時(shí),會(huì)產(chǎn)生很強(qiáng)的散射現(xiàn)象,模擬結(jié)果可能會(huì)影響到面波傳播規(guī)律的討論。

圖2 高斯形谷模型(a)與高斯形山模型(b)

2.2 凹陷模型結(jié)果

對不同深度(25～150m,間隔25m)和跨度(50～150m,間隔25m)共計(jì)30個(gè)凹陷模型進(jìn)行模擬。圖3顯示了其中兩種情形(D0=50m,w=150m;D0=150m,w=50m)的波場快照,可以看出,不同的深度和跨度凹陷地形對面波波場傳播的影響。當(dāng)?shù)匦屋^為平緩時(shí)(圖3a),面波穿過地形畸變較小,但會(huì)產(chǎn)生轉(zhuǎn)換波(P波和S波)向下傳播;當(dāng)?shù)匦屋^為陡峭時(shí)(圖3b),面波在地形凹陷處產(chǎn)生很強(qiáng)的散射現(xiàn)象,并且出現(xiàn)反射面波,透射面波能量出現(xiàn)顯著衰減,轉(zhuǎn)換波能量較強(qiáng)。在實(shí)際地震勘探中,如果地表起伏較為劇烈,向下方傳播的面波轉(zhuǎn)換波也可能是造成波場復(fù)雜的一個(gè)重要原因。圖4為不同的深度和跨度凹陷地形附近(1～5km)截取的面波地震記錄(1.5～3.0s)。當(dāng)?shù)匦屋^為陡峭時(shí),可以看到清晰的反射面波和透射面波以及轉(zhuǎn)換波(散波P波和S波)(圖4b)。

在凹陷地形兩側(cè)面波與散射波傳播較為穩(wěn)定處取兩點(diǎn)(1km與5km處),分別計(jì)算反射和透射面波能量,來討論地形對面波能量分配的影響。面波能量計(jì)算公式為:

(22)

式中:vx(t)和vy(t)為速度分量,面波波形位于時(shí)間[t1,t2]區(qū)間內(nèi),區(qū)間長度略大于面波波形長度;A0為歸一化系數(shù);E為所求面波能量。圖5中所示紅框?yàn)橛?jì)算直達(dá)瑞雷面波和反射瑞雷面波的時(shí)間區(qū)間,其中二次瑞雷面波為來自與P波經(jīng)過地形時(shí)產(chǎn)生的散射瑞雷面波。

圖3 不同深度和跨度凹陷地形模型波場快照(vx分量) a D0=50m,w=150m; b D0=150m,w=50m

圖4 凹陷地形地表處面波的地震記錄(vx分量) a D0=50m,w=150m; b D0=150m,w=50m

圖5 面波波形截取 D0=150m,w=50m模型vx分量地震記錄,檢波點(diǎn)位于(1km,0)處

圖6顯示了經(jīng)過不同的深度和跨度凹陷地形面波能量分配情況,可以看出,對于相同的高度,隨跨度的減小反射面波能量逐漸增強(qiáng)(圖6a),透射面波能量逐漸減弱(圖6b),轉(zhuǎn)換波能量呈整體逐漸增強(qiáng)趨勢(轉(zhuǎn)換波能量為入射面波減去反射面波和透射面波能量)(圖6c);而對于相同的跨度(圖6a,圖6b,圖6c),隨高度的增加,同樣出現(xiàn)反射面波能量增強(qiáng)、透射面波能量逐漸減弱和轉(zhuǎn)換波能量增強(qiáng)的現(xiàn)象。值得注意的是,當(dāng)深度大于某個(gè)數(shù)值時(shí)(100m,與波長相關(guān)),轉(zhuǎn)換波能量隨跨度減小而略微減小,主要原因是反射面波急劇增強(qiáng)分配了部分能量。圖7給出了1km與5km處反射與透射面波的頻譜,可以看出,地形對面波具有濾波效應(yīng),體現(xiàn)為反射和透射面波頻譜上出現(xiàn)陷波點(diǎn)。其陷波點(diǎn)頻率受深度和跨度的影響,當(dāng)深度與跨度增加時(shí),陷波點(diǎn)向低頻移動(dòng)。

圖6 經(jīng)過凹陷地形后面波能量變化 a 反射面波能量; b 透射面波能量; c 轉(zhuǎn)換波能量

圖7 反射面波與透射面波頻譜(凹陷地形) a 反射面波頻譜(跨度50m); b 透射面波頻譜(跨度50m); c反射面波頻譜(深度100m); d透射面波頻譜(深度100m)

2.3 凸起模型結(jié)果

為討論不同高度和跨度的凸起地形對面波傳播的影響,我們共設(shè)計(jì)30個(gè)凸起模型進(jìn)行模擬(高度為25～150m,跨度為50～150m,間隔25m),介質(zhì)參數(shù)和震源參數(shù)與凹陷地形模型相同,其中兩種情形的波場快照如圖8所示。圖9為凸起地形地表處面波的地震記錄。從圖9可以看出,當(dāng)?shù)匦屋^為平緩時(shí),對面波傳播影響較小(圖9a),但也會(huì)產(chǎn)生向下傳播的轉(zhuǎn)換波;當(dāng)?shù)匦巫兓^為劇烈時(shí),其波場變得較為復(fù)雜,并產(chǎn)生反射面波和透射面波以及轉(zhuǎn)換波(圖9b)。

同樣,在面波與其散射波傳播較為穩(wěn)定的1km與5km兩點(diǎn)處分析反射和透射面波能量的變化情況。圖10顯示了經(jīng)過凸起地形的面波能量變化情況,可以看出,影響能量分配的受控因素較為復(fù)雜。整體看,當(dāng)跨度減小時(shí),其反射面波與轉(zhuǎn)換能量增強(qiáng),透射面波能量減弱,但當(dāng)高度大于某個(gè)閾值(100m),跨度小于某個(gè)閾值(75～100m)時(shí),情況出現(xiàn)反轉(zhuǎn)。其原因是當(dāng)跨度小于面波橫向波長時(shí),面波可以不沿凸起地形表面?zhèn)鞑?直接穿過凸起地形,這與凹陷情形不同。

圖11顯示了1km與5km處反射面波和透射面波的頻譜,可以看出,其形態(tài)變化也較為復(fù)雜。隨著高度不同,透射面波頻譜出現(xiàn)了陷波現(xiàn)象,其陷波點(diǎn)頻率隨高度增加有減小趨勢(圖11b)。對于相同高度凸起地形(150m),當(dāng)跨度小于某個(gè)閾值時(shí)(50m),其頻譜形態(tài)產(chǎn)生本質(zhì)變化,其反射能量急劇減弱并出現(xiàn)陷波點(diǎn)(圖11c),透射能量增強(qiáng)(圖11d)。

圖8 不同深度和跨度凸起地形模型波場快照(vx分量) a D0=50m,w=150m; b D0=150m,w=50m

圖9 凸起地形地表處面波的地震記錄(vx分量) a D0=50m,w=150m; b D0=150m,w=50m

圖10 經(jīng)過凸起地形后面波能量變化 a 反射面波能量; b 透射面波能量; c 轉(zhuǎn)換波能量

圖11 反射面波與透射面波頻譜(凸起地形) a 反射面波頻譜(跨度150m); b 透射面波頻譜(跨度150m); c 反射面波頻譜(高度150m); d 透射面波頻譜(高度150m)

3 結(jié)論

本文采用貼體網(wǎng)格有限差分方法研究了瑞雷面波經(jīng)過起伏地形時(shí)的傳播特性,討論了地形的凹凸、高(深)度和跨度等因素對面波波形、能量、頻散等參數(shù)的影響。

對于凹陷地形,地表起伏平緩時(shí),面波可以平滑地繞過地形;地表起伏劇烈時(shí),面波產(chǎn)生嚴(yán)重的散射現(xiàn)象。地形不僅對面波的高頻部分有影響,對低頻部分還存在濾波效應(yīng),表現(xiàn)為在頻譜上出現(xiàn)陷波點(diǎn),陷波點(diǎn)的頻率與地形的高度、跨度以及面波橫向波長有一定關(guān)系。

對于凸起地形,地表起伏平緩時(shí)對面波傳播影響較小,起伏劇烈時(shí)波場變得復(fù)雜;當(dāng)跨度較小、起伏較大時(shí),面波通過凸起地形的機(jī)理與凹陷地形不同,而且表現(xiàn)形式更加復(fù)雜。凹陷地形表現(xiàn)為阻擋,產(chǎn)生較強(qiáng)反射面波,而凸起地形表現(xiàn)為面波直接穿過地形,此時(shí)透射面波能量增大,其能量分配規(guī)律與頻譜表現(xiàn)也更加復(fù)雜。

從模擬結(jié)果看,可能存在純地形與瑞雷波能量陷波的經(jīng)驗(yàn)關(guān)系,我們將在下一步研究工作中討論該問題。另外,貼體網(wǎng)格有限差分方法結(jié)合牽引力鏡像法可以很好地處理起伏地表問題,為研究更復(fù)雜的近地表?xiàng)l件下地震波傳播規(guī)律提供有效手段。

[1] WEN Y Y,MA K F,OGLESBY D D.Variations in rupture speed,slip amplitude and slip direction during the 2008 Mw 7.9 Wenchuan earthquake[J].Geophysical Journal International,2012,190(1):379-390

[2] SHEN Z K,SUN J,ZHANG P,et al.Slip maxima at fault junctions and rupturing of barriers during the 2008 Wenchuan earthquake[J].Nature Geoscience,2009,2(10):718-724

[3] WEN J,CHEN X.Variations in f_max along the Ruptured Fault during the Mw 7.9 Wenchuan earthquake of 12 May 2008[J].Bulletin of the Seismological Society of America,2012,102(3):991-998

[4] 陳偉,徐義賢,張煜.瑞雷面波在彈性楔狀體中的傳播特性研究[J].工程地球物理學(xué)報(bào),2009,6(2):143-149

CHEN W,XU Y X,ZHANG Y.Study of transmission characteristics of Rayleigh waves on an elastic wedge[J].Chinese Journal of Engineering Geophysics,2009,6(2):143-149

[5] 沈鴻雁,李慶春,嚴(yán)月英,等.多道瞬態(tài)面波相速度分析[J].石油物探,2016,55(5):692-702

SHEN H Y,LI Q C,YAN Y Y,et al.Phase velocity analysis of multi-channel transient surface wave[J].Geophysical Prospecting for Petroleum,2016,55(5):692-702

[6] 岳龍,劉懷山,尹燕欣,等.基于連續(xù)小波變換的面波衰減方法研究[J].石油物探,2016,55(2):214-222

YUE L,LIU H S,YIN Y X,et al.Attenuation of ground roll based on continuous wavelet transform[J].Geophysical Prospecting for Petroleum,2016,55(2):214-222

[7] KNOPOFF L.Transmission and reflection of Rayleigh waves by wedges[J].Geophysics,1960,25(6):1203-1214

[8] FUJII K.Rayleigh-wave scattering at various wedge corners:investigation in the wider range of wedge angles[J].Bulletin of the Seismological Society of America,1994,84(6):1916-1924

[9] HUDSON J A,KNOPOFF L.Transmission and reflection of surface waves at a corner,2,Rayleigh waves (theoretical)[J].Journal of Geophysical Research Atmospheres,1964,69(2):281-289

[10] MAL A K,KNOPOFF L.Transmission of Rayleigh waves at a corner[J].Bulletin of the Seismological Society of America,1966,56(2):455-466

[11] 周紅,陳曉非.凹陷地形對Rayleigh面波傳播影響的研究[J].地球物理學(xué)報(bào),2007,50(4):1182-1189

ZHOU H,CHEN X F.A study on the effect of depressed topography on Rayleigh surface wave[J].Chinese Journal of Geophysics,2007,50(4):1182-1189

[12] 汪利民.起伏地表三維高頻瑞雷面波傳播特性研究[D].武漢:中國地質(zhì)大學(xué),2013

WANG L M.Propagation of high-frequency Rayleigh waves in 3d medium with topography[D].Wuhan:China University of Geoscience,2013

[13] WANG L M,XU Y X,XIA J H,et al.Effect of near-surface topography on high-frequency Rayleigh-wave propagation[J].Journal of Applied Geophysics,2015,116:93-103

[14] 巴振寧,梁建文.Rayleigh波斜入射下層狀場地中凸起地形的三維響應(yīng)分析[J].中國科學(xué):技術(shù)科學(xué),2015,45(8):874-888

BA Z N,LIANG J W.Three dimensional responses of a hill in a layered half-space for obliquely incident Rayleigh waves[J].Scientia Sinica Techolgica,2015,45(8):874-888

[15] 巴振寧,梁建文.層狀場地中凹陷地形Rayleigh波斜入射下三維地震響應(yīng)分析[J].振動(dòng)工程學(xué)報(bào),2015,28(5):809-821

BA Z N,LIANG J W.3-D seismic responses for oblique incident Rayleigh waves of a canyon cut in a layered half-space[J].Journal of Vibration Engineering,2015,28(5):809-821

[16] 劉中憲,武鳳嬌,王冬.沉積盆地-山體耦合場地對平面P波、SV波和Rayleigh波的二維散射分析[J].工程力學(xué),2016,33(2):160-171

LIU Z X,WU F J,WANG D.Two dimensional scattering analysis of plane P,SV,and Rayleigh waves by coupled alluvial basin-mountain terrain[J].Engineering Mechanics,2016,33(2):160-171

[17] 朱多林,白超英.基于波動(dòng)方程理論的地震波場數(shù)值模擬方法綜述[J].地球物理學(xué)進(jìn)展,2011,26(5):1588-1599

ZHU D L,BAI C Y.Review on the seismic wavefield forward modelling[J].Progress in Geophysics,2011,26(5):1588-1599

[18] 張明財(cái),熊章強(qiáng),張大洲.基于MPI的三維瑞雷面波有限差分并行模擬[J].石油物探,2013,52(4):354-362

ZHANG M C,XIONG Z Q,ZHANG D Z.3D finite difference parallel simulation of Rayleigh wave based on message passing interface[J].Geophysical Prospecting for Petroleum,2013,52(4):354-362

[19] ROBERTSSON J O A.A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography[J].Geophysics,1996,61(6):1921-1934

[20] 張偉.含起伏地形的三維非均勻介質(zhì)中地震波傳播的有限差分算法及其在強(qiáng)地面震動(dòng)模擬中的應(yīng)用[D].北京:北京大學(xué),2006

ZHANG W.Finite difference seismic wave modelling in 3D heterogeneous media with surface topography and its implementation in strong ground motion study[D].Beijing:Peiking University,2006

[21] 丘磊.正交曲線坐標(biāo)系下的地震波數(shù)值模擬技術(shù)研究[D].杭州:浙江大學(xué),2011

QIU L.Study on seismic wave simulations in orthogonal curvilinear coordinate system[D].Hangzhou:Zhejiang University,2011

[22] 裴正林.任意起伏地表彈性波方程交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬[J].石油地球物理勘探,2004,39(6):629-634

PEI Z L.Numerical modeling using staggered-grid high-order finite-difference of elastic wave equation on arbitrary relief surface[J].Oil Geophysical Prospecting,2004,39(6):629-634

[23] HIXON R.On increasing the accuracy of MacCormack schemes for aeroacoustic applications[J/OB].AIAA,1997:1-33[2017-05-12].https://www.researchgate.net/publication/24298166

[24] HU F Q,HUSSAINI M Y,MANTHEY J L.Low-dissipation and low-dispersion runge-kutta schemes for computational acoustics[J].Journal of Computational Physics,1996,124(1):177-191

[25] LEVANDER A R.Fourth-order finite-difference P-SV seismograms[J].Geophysics,1988,53(11):1425-1436

[26] ZHANG W,CHEN X.Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J].Geophysical Journal International,2006,167(1):337-353

[27] LELE S K.Compact finite difference schemes with spectral-like resolution[J].Journal of Computational Physics,1992,103(1):16-42

[28] ZHANG W,SHEN Y.Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling[J].Geophysics,2010,75(4):T141-T154

(編輯：顧石慶)

Rayleighwavepropagationwithundulatingtopography

LI Hong1, YANG Xinchao1, ZHU Haibo1, ZHANG Wei2, QU Shouli1

(1.SinopecGeophysicalResearchInstitute,Nanjing211103,China;2.SouthernUniversityofScienceandTechnology,Shenzhen518055,China)

This study adopts the body-fitted grid finite-difference method to investigate Rayleigh wave propagation with undulating topography.Factors including concavity-convexity,height,and span of topography influence wave shape,energy,and frequency spectrum.It was found that when undulating topography changes smoothly,a Rayleigh wave can traverse the topography with no obvious distortion;conversely,strenuously undulating topography can result in clearly discernible frequency dispersion.It was found that undulating topography has a filtering effect on the frequency spectrum,which can cause a notch on the spectrum of the Rayleigh wave’s reflection/transmission,such that the notch frequency depends on the height/width of the topographic undulation.The curvilinear grid finite-difference method could be effectively applied to seismic modeling with undulating topography.

curvilinear grid,finite difference,undulating topography,surface wave modeling

2017-05-24;改回日期2017-06-13。

李宏(1986—),男,博士,主要從事水力壓裂微地震資料處理解釋及地震波數(shù)值模擬研究。

國家科技重大專項(xiàng)“大型油氣田及煤層氣開發(fā)及勘探評(píng)價(jià)-涪陵頁巖氣開發(fā)示范工程”(2016ZX05060-005)資助。

This research is financially supported by the National Science and Technology Major Project (Grant No. 2016ZX05060-005).

P631

A

1000-1441(2017)06-0782-10

10.3969/j.issn.1000-1441.2017.06.002