李緒宣 于更新 符力耘 溫書亮 管西竹

（1.中海油研究總院； 2.中國科學(xué)院地質(zhì)與地球物理研究所）

應(yīng)用邊界元模擬方法分析復(fù)雜海底地震散射特征＊

李緒宣1于更新2符力耘2溫書亮1管西竹2

（1.中海油研究總院； 2.中國科學(xué)院地質(zhì)與地球物理研究所）

邊界元法對隨機起伏的復(fù)雜海底界面具有良好的適應(yīng)性。比較了邊界元法與有限差分法對復(fù)雜斷層模型的模擬精度，并驗證了邊界元法的有效性。利用邊界元法對復(fù)雜海底模型進行波場模擬，反映起伏海底界面對地震波傳播的影響；利用統(tǒng)計參數(shù)描述復(fù)雜海底地貌特征，將崎嶇海底界面劃分為快、慢變化和強、弱起伏等4種特征類型。根據(jù)不同統(tǒng)計參數(shù)的選擇建立崎嶇海底理論模型，利用邊界元法對不同類型的崎嶇海底理論模型進行模擬研究，同時與實際海底資料相對比，分析了復(fù)雜海底地震散射特征。此項研究成果可為復(fù)雜海底地區(qū)目標(biāo)導(dǎo)向地震觀測系統(tǒng)設(shè)計和采集參數(shù)優(yōu)化提供理論依據(jù)。

邊界元法 復(fù)雜海底 地震波傳播 散射特征

深水海底地貌和結(jié)構(gòu)復(fù)雜，而地震勘探觀測一般是在水面進行，因此深刻理解地震波在復(fù)雜海底地貌和結(jié)構(gòu)上的散射特征對于地震采集觀測具有重要的應(yīng)用價值。邊界元法是繼有限元法發(fā)展起來的一種數(shù)值模擬方法，已在過去幾十年里得到了廣泛的研究和發(fā)展[1－11]。由于邊界元法只在研究區(qū)域的邊界上剖分單元，因而將求解問題的維數(shù)降低，大大減少了所要求解的節(jié)點數(shù)，顯著提高了計算效率。該方法的特點是：①自動滿足無窮遠處邊界條件，適于無限域問題；②在幾何上精確描述地下不規(guī)則界面；③顯式應(yīng)用地下界面的連續(xù)邊界條件；④可分別模擬來自地下不同部位的波。因此，筆者認為邊界元法對隨機起伏的復(fù)雜界面有著良好的適應(yīng)性，非常適用于復(fù)雜海底地貌和構(gòu)造的地震波正演模擬。

本文利用邊界元法對復(fù)雜海底模型進行波場模擬，反映起伏海底界面對地震波傳播的影響；并將邊界元法與有限差分法相對比，驗證了邊界元法的有效性；使用統(tǒng)計參數(shù)描述復(fù)雜海底地貌的統(tǒng)計特征，將崎嶇海底界面劃分為快、慢變化和強、弱起伏等4種特征類型；根據(jù)不同統(tǒng)計參數(shù)的選擇建立崎嶇海底理論模型，應(yīng)用邊界元法對不同類型的崎嶇海底理論模型進行模擬研究，同時與實際海底資料相對比，定量分析了復(fù)雜海底地震散射特征。此項研究可為復(fù)雜海底地區(qū)目標(biāo)導(dǎo)向地震觀測系統(tǒng)設(shè)計和采集參數(shù)優(yōu)化提供理論依據(jù)。

1 邊界元模擬方法

1.1 邊界積分方程及其離散

圖1所示為一起伏海底模型，由Ω0和Ω12個子域組成，子域Ω0的上下邊界分別以Γ0和Γ1表示，子域Ω1由邊界Γ1和無窮遠處邊界表示，且假定其兩端趨于無窮。界面Γ0為自由水面，震源位于子域Ω0內(nèi)（為簡化起見，這里只討論二維聲波問題），則子域Ω0內(nèi)的地震波響應(yīng)u（r）滿足邊界積分方程

式（1）中：Γ為邊界，Γ＝Γ0＋Γ1；u0（r）為入射場，u0（r）＝S（ω）G（r，r0）；?/?n 表示邊界 Γ 外法向?qū)?shù)：r為源節(jié)點的位置矢量；r′為易節(jié)點的位置矢量；G（r，r′）表示均勻背景介質(zhì)空間的格林函數(shù)，對于二維問題，可表示為

圖1 起伏海底模型示意圖

對式（1）應(yīng)用自由水面應(yīng)力為零的條件，則得到

對于子域Ω1，無窮遠處邊界上位移和應(yīng)力滿足遠場輻射條件，即

則其對應(yīng)邊界積分方程為

在Ω0與Ω1之間的公共界面上Γ1位移和應(yīng)力連續(xù)，即

式（6）中：u＋（r）為邊界Γ1上界面的位移；u－（r）為邊界Γ1下界面的位移。

為了求解積分方程，首先要對邊界Γ進行離散。將邊界Γ劃分為L個邊界元，記為Γe（e＝1，2，…，L），總節(jié)點數(shù)為N。假設(shè)小單元Γe由節(jié)點I1、I2所確定，利用插值形函數(shù)φ，變量r、u和?u/?n即可由節(jié)點上相應(yīng)變量值的線性組合來近似，以u為例：

其中ξ為一個單元的局部坐標(biāo)。

將方程（1）寫為算子形式：

式（8）中：H、G為邊界積分算子；f為入射波場；u為邊界上的位移；t為u在邊界上的法向?qū)?shù)，t＝?u/?n。邊界積分算子H和G的離散形式可分別表示為

其中δij為Kronecker delta函數(shù)。這些積分可由高斯積分數(shù)值算法計算得出。在每個子域按式（8）建立方程，則可得到對應(yīng)的總體矩陣方程組，利用高斯消去法求解此線性方程組，即可求得模型中所有節(jié)點的地震響應(yīng)u（r）。

1.2 改進算法

為提高算法的效率，我們采用了改進的分塊高斯消去法來求解矩陣方程。在海上勘探中，震源和檢波器均位于海水表層，所以僅需計算自由水面的地震響應(yīng)。采用由深至淺的方式逐層進行邊界系數(shù)矩陣的數(shù)值計算、組裝和消元，消去2個相鄰子域公共邊界上的耦合數(shù)據(jù)，如此計算至水層，得到最后的矩陣方程，其維數(shù)僅由水層的節(jié)點數(shù)確定。將此矩陣與多測線多震源相結(jié)合，可計算出水面任意處的觀測波場。在整個計算過程中，主矩陣需要的最大內(nèi)存僅由最大子域的總節(jié)點數(shù)確定，這樣既節(jié)省了內(nèi)存又縮短了計算時間。

為了提高計算效率，在程序執(zhí)行中采用了一種變單元尺寸的方法。由于單位波長3個單元采樣足以確保結(jié)果的精度，而波長是頻率的函數(shù)，所以程序可根據(jù)介質(zhì)速度和每一個計算頻率自動計算出相應(yīng)的單元尺寸并自動剖分單元，從而提高了計算速度并使參數(shù)輸入簡化。對于存在人工邊界的子域，引入無限元吸收邊界[12]可使模型截斷邊界的散射最小化。與傳統(tǒng)吸收邊界條件方法相比，此方法不僅有效改善了吸收效果，而且節(jié)省了存儲空間和運算時間。

1.3 有效性驗證

為了驗證邊界元法對復(fù)雜海底地貌模型的正確性和有效性，首先設(shè)計了一個復(fù)雜斷層模型，模型尺寸為3800 m×1500 m，各層層速度如圖2a所示，計算頻率為0～40 Hz。分別使用邊界元法與有限差分法對該模型進行數(shù)值模擬，對比邊界元法模擬結(jié)果（圖2b）與有限差分法模擬結(jié)果（圖2c），可以看出兩者具有較好的一致性；但邊界元法可以清晰地刻畫出斷層各斷點的繞射，而有限差分法的模擬結(jié)果中各斷點的繞射則相對較弱。

為了進一步驗證邊界元法的模擬精度，分析該方法在崎嶇海底模擬中的主要優(yōu)勢，又設(shè)計了一個起伏海底的二維速度模型進行模擬，各層層速度如圖3a所示。該復(fù)雜海底模型大小為3000 m×2000 m，炮點位于（1500 m，0 m）處，計算頻率為0～40 Hz。圖3b—d分別為550、750和900 ms時刻計算的波場快照，可以看出模擬結(jié)果比較清晰地反映出了起伏海底界面對地震波傳播的影響。

圖2 設(shè)計的復(fù)雜斷層模型及其自激自收剖面

圖3 復(fù)雜海底地貌二維速度模型及其不同時刻的波場快照

根據(jù)對上述2個模型的模擬可以看出，邊界元法能高效準(zhǔn)確地模擬地震波在復(fù)雜地質(zhì)構(gòu)造中的傳播，精確描述各層界面，準(zhǔn)確反映各層間的波阻抗差異及地層厚度的變化，有效處理起伏界面引起的地震波散射，因此對崎嶇海底模型的數(shù)值模擬具有明顯的優(yōu)勢。

2 復(fù)雜海底地震散射特征分析

為了研究復(fù)雜海底地貌對地震波傳播的影響，建立了崎嶇海底理論模型和實際模型，并分別對2類模型進行邊界元地震數(shù)值模擬，對地震波的傳播散射特征進行了定量分析。

一個隨機起伏的海底界面可以由不同的統(tǒng)計參數(shù)來描述。均方根高用來描述一條曲線或者一個表面的平均起伏高度，可以用來表征一個界面的崎嶇程度。設(shè)界面高度分布函數(shù)為h（r），其中h為界面距離參考界面的高度，r為參考界面上點的位矢，則此界面的均方根高為相關(guān)長度用來描述一條曲線或者一個表面的起伏頻率，可以用來表征一個界面的崎嶇程度變化的頻率，其計算公式為C（r）＝〈h（r′）h（r′＋r）〉/σ2，式中〈h（r′）h（r′＋r）〉表示自相關(guān)，σ2為歸一化因子。利用均方根高和相關(guān)長度這2個統(tǒng)計參數(shù)可以得到不同的隨機起伏海底模型。圖4a為一組均方根高相同、相關(guān)長度不同的隨機起伏界面，其均方根高d為20 m，相關(guān)長度a分別為22 m（細線）和220 m（粗線）。圖4b為一組相關(guān)長度相同、均方根高不同產(chǎn)生的模型，其相關(guān)長度a為40 m，均方根高d分別為50 m（細線）和10 m（粗線）。從圖4中可以看出：均方根高越大，其界面隨機起伏程度越強；反之，均方根高越小，其界面隨機起伏程度越弱。而相關(guān)長度越小，其界面起伏頻率變化越快；反之，相關(guān)長度越大，其界面起伏頻率變化越慢。

圖4 隨機起伏海底理論模型

圖5 不同隨機起伏海底模型的反射能量曲線（a、b）與反散射能量曲線（c、d）

分別使用以上4條曲線作為海底隨機起伏界面進行數(shù)值模擬。設(shè)水層速度為1500 m/s，下伏地層速度為3000 m/s，炮點坐標(biāo)為（0 m，0 m），計算頻率為0～45 Hz，主頻為20 Hz，利用對復(fù)雜海底理論模型的地震模擬結(jié)果可計算得到反射能量曲線（圖5a、b）與反散射能量曲線（圖5c、d）。反射能量曲線是通過對單炮記錄所有頻率分量總和進行計算得到，而反散射能量曲線則是利用歸一化后的起伏海底模型的反射能量曲線與水平海底模型的反射能量曲線相減獲得。由圖5可知，能量曲線在水平距離大約為600 m處出現(xiàn)能量峰值，之后隨偏移距增大而逐漸衰減。根據(jù)反射能量曲線變化可以看出，快變海底（相關(guān)長度小）反散射損失與慢變海底（相關(guān)長度大）反散射損失基本持平；強起伏海底（均方根高大）反散射損失大，而弱起伏海底（均方根高?。┓瓷⑸鋼p失較小。類似地，我們也計算了透散射能量曲線，其變化規(guī)律與反散射情況比較相似。對于強起伏快變化海底模型（均方根高大、相關(guān)長度小），地震波變化劇烈且與其復(fù)雜程度表現(xiàn)出一定的相關(guān)性，其透散射能量在近偏移距變化很大且隨偏移距增大而快速衰減。為了分析起伏海底對不同頻率的響應(yīng)特征，我們將計算頻率劃分為低（0～25Hz）、中（26～50 Hz）、高（51～75 Hz）等3段分別計算其反射能量并進行比較（圖6），結(jié)果表明中低頻情況下快變海底在中遠偏移距下的反射能量均小于慢變海底情況，而在高頻情況下兩者相近。

此外，我們根據(jù)某地區(qū)實際海底起伏設(shè)計模型（圖7a）并進行了數(shù)值模擬。該實際海底模型的均方根高為84 m，相關(guān)長度為5350 m。設(shè)水層速度為1500 m/s，下伏地層速度為3000 m/s，炮點位于水面中心位置，計算頻率為0～40 Hz，主頻為15 Hz。對此海底實際模型進行地震模擬，可計算得到其對應(yīng)的反射能量曲線（圖7b）及其透射能量曲線（圖7c），可以看出，此實際模型的反射能量曲線變化與起伏海底界面復(fù)雜程度呈現(xiàn)出一定相關(guān)性，而其透射能量曲線與水平海底透射能量曲線相接近。

圖6 不同頻率下快變海底模型（a＝22 m，細線）與慢變海底模型（a＝220 m，粗線）的反射能量曲線

圖7 實際起伏海底模型及其散射能量分析

3 結(jié)論與建議

（1）應(yīng)用邊界元法對復(fù)雜海底理論模型和實際模型的數(shù)值模擬與散射能量分析表明：快變化強起伏海底散射損失大，在進行海上勘探時建議變觀以避開局部強變化海底的影響，或改變激發(fā)角度減弱散射損失。

（2）快變化海底散射引起地震波能量隨偏移距起伏變化，變化規(guī)律與海底起伏具有一定相關(guān)性，據(jù)此可進行地表一致性校正。

（3）快變化強起伏海底，其散射能量在近偏移距變化很大且隨偏移距衰減很快；起伏海底透散射作用弱于反散射作用。

（4）不同起伏海底對頻率具有一定選擇性，據(jù)此可在地震采集中針對實際海底調(diào)整激發(fā)頻率以取得理想效果。

[1] HE Nanqun，MU Yongguang，LI Guibiao，et al.VSP migration using the boundary element method[C]∥Beijing（89）International Symposium on Exploration Geophysics.Beijing，1989.

[2] SáNCHEZ－SESMA F J，CAMPILLO M.Diffraction of P，SV and Rayleigh waves by topographic features：a boundary integral formulation[J].Bull.Seism.Soc.Am.，1991，81：2234－2253.

[3] 徐世浙，王慶乙，王軍.用邊界單元法模擬二維地形對大地電磁場的影響[J].地球物理學(xué)報，1992，35（3）：380－388.

[4] 徐世浙，阮百堯，周輝，等.大地電磁場三維地形影響的數(shù)值模擬[J].中國科學(xué)：D輯，1997，27（1）：15－20.

[5] 符力耘，牟永光.彈性波邊界元法正演模擬[J].地球物理學(xué)報，1994，37（4）：521－529.

[6] FU L Y，WU R S.A hybrid BE－GS method for modeling regional wave propagation[J].Pure and Applied Geophysics，2001，158：1251－1277.

[7] SCHUSTER G T.The application of boundary integral equations to interactive modeling[C]∥The 3rd MIDAS Meeting.New York，1982.

[8] KAWASE H.Time－domain response of a semi－circular canyon for incident SV，P，and Rayleigh waves calculated by the discrete wave number boundary element method[J].Bull.Seism.Soc.Am.，1988，78：1415－1437.

[9] BOUCHON M，CAMPILLO M，GAFFET S.A boundary integral equation－discrete wavenumber representation method to study wave propagation in multilayered media having irregular interfaces[J].Geophysics，1989，54：1134－1140.

[10] FU L Y.Seismogram synthesis for piecewise heterogeneous media[J].Geophys.J.int.，2002，150：800－808.

[11] FU L Y.Numerical study of generalized lippmann－schwinger integral equation including surface topography[J].Geophysics，2003，68（2）：665－671.

[12] FU L Y，WU R S.Infinite boundary element absorbing boundary for wave propagation simulations[J].Geophysics，2000，65（2）：596－602.

Analysing seismic scattering characteristics of complex seabed by using the boundary－element simulation method

Li Xuxuan1Yu Gengxin2Fu Liyun2Wen Shuliang1Guan Xizhu2
（1.CNOOC Research Institute，Beijing，100027；2.Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing，100029）

The boundary－element method （BEM）has a good adaptability for simulating irregularly rough and complex seabed.The simulation accuracy of BEM for a complex fault model was compared with that of the finite－difference method，and the effectiveness of BEM was confirmed.BEM can be used to conduct wave simulation of rough seabed models，reflecting the impacts of rough seabed on seismic wave propagation.The statistical parameters were used to describe complex seabed topography，and then four types of rough seabed interface can be identified，i.e.fast lateral change，slow lateral change，strong vertical relief and weak vertical relief.The theoretical models of rough seabed can be build by selecting various statistical parameters，and BEM was used to make simulation of different theoretical models of rough seabed.Simultaneously，some actual seabed data was compared and the seismic scattering characteristics of complex seabed were analyzed.These results will provide some theoretical foundations for the seismic acquisition design of complex seabed and the optimization of seismic acquisition parameters.

the boundary－element method；complex seabed；seismic wave propagation；seismic scattering characteristics

＊國家重點基礎(chǔ)研究計劃（973計劃）項目“南海深水區(qū)復(fù)雜地質(zhì)結(jié)構(gòu)地震采集基礎(chǔ)理論研究（2009CB219403）”資助部分研究成果。

李緒宣，男，高級工程師，現(xiàn)任中海油研究總院地球物理總師。地址：北京市東城區(qū)東直門外小街6號海油大廈（郵編：100027）。E－mail：lixx＠cnooc.com.cn。

2010－10－26改回日期：2011－06－12

（編輯：周雯雯）