3D主動源地震折射／反射聯合快速波前成像研究

劉 妍 徐朝繁

1 中國地震局蘭州地震研究所，蘭州市東崗西路450號，730000

2 中國地震局地球物理勘探中心，鄭州市文化路75號，450002

在地殼深淺部細結構及物性特征的探測與研究中，寬角折射／反射地震剖面得到廣泛應用［1－2］。深淺部3D 精細結構是深入理解與認識地殼構造運動演化過程及震害評估最重要的依據，然而由于受人工地震測深專用儀器設備實用性、自動化程度及數量的限制，以前所實施的深地震測深剖面只能給出2D 地殼結構圖像，其結果的不足是顯而易見的。隨著主動源地震觀測設備的快速發展，目前已具備采集3D 折射／反射深地震測深數據體的設備條件，相應的技術儲備已顯得十分必要。一些學者基于在一些局部區域已實施的較為密集的地震測線網或專門設計的人工地震觀測系統資料，利用單一的Pg或PmP 震相構建3D 地殼結構模型［3－4］，取得了一些進展，但缺乏對所反演的3D 地殼結構模型的有效約束，直接影響到反演結果的分辨能力與精度。雖然國內已有學者對多種震相聯合反演進行研究［5－6］，但為數不多。此外，有限差分反演在強速度對比結構模型中，波場的一階導數可能不存在［7］，即程函方程不成立，將導致正演計算的不穩定。本文試圖基于快速波前追蹤算法FMM（fast marching method）［8］，利用合成3D 折射／反射數據體，對強速度對比結構及反射界面橫向起伏變化模型進行聯合波前成像研究，為3D 折射／反射地震聯合試驗探測奠定了基礎。

1 主動源地震折射／反射震相記錄特征

人工地震記錄截面圖上的震相類別及其運動學、動力學特征與地殼的結構特征密切相關，對所識別的震相資料進行分析處理，可建立探測區較精細的地殼物性結構及構造模型。在寬角反射折射記錄截面圖上，可追蹤到十分清晰的殼內波組Pg、PmP，這兩個震相具有全球普遍性，若探測剖面足夠長且激發能量足夠大，還可觀測到低頻的上地幔頂部弱折射震相Pn。一般來說，在一張觀測系統較為完備的共炮點寬角反射／折射記錄圖上，都可追蹤到十分清晰的Pg、PmP，有時還能觀測到能量較弱的Pn震相，這3 個震相的運動學和動力學特征與地殼結晶基底、Moho及上地幔頂部結構密切相關，其他的殼內反射震相Pi，其能量一般很弱且不穩定，只在一些局部的區段被追蹤到，反映了地殼的小區域局部小尺度的結構特征。圖1為華南某地的一張典型寬角折射／反射地震探測記錄截面圖，可以發現，來自于中上地殼的折射震相Pg 及Moho 面的寬角反射震相PmP十分清晰。我們將基于目前可用的人工地震觀測設備，設計3D 觀測系統及強速度對比結構模型，合成折射／反射數據體，利用Pg、PmP聯合成像方法，對3D 主動源地震折射／反射聯合成像的有效性及分辨能力進行數值模擬計算分析，為即將開展的3D 深部細結構探查提供參考。

圖1 人工地震寬角折射／反射探測典型記錄截面Fig.1 Typical seismic record section of wide angle refraction／reflection investigation

2 數值模擬試驗

2.1 理論模型及觀測系統

大量的深地震測深結果表明，大陸造山帶及火山區地殼表現為十分復雜的結構特征，如川西的龍門山及周邊地區，不但地殼厚度變化大，而且在殼內存在高速的雜巖體及低速帶，在其他一些熱點地區還存在地幔柱及火成巖墻，這些地區的地殼往往表現為較強的速度對比結構特征。為此，本文計算3D折射／反射數據體理論模型的設計重點將放在強速度對比結構及傾斜Moho的起伏變化上，由非均勻層、高低速圓柱、直立十字墻及較為復雜的Moho結構界面組成（圖2）。圖2（a）、（b）、（c）分別為該模型在0、25km 深度的切片、Moho深度等值線及疊置在其上的觀測系統示意圖。

在該模型中，非均勻層的頂界為一平界面，P波速度為5.8km／s，底界由從東40km 向西過渡到55km 的傾斜界面疊加峰值±2km 的隨機擾動組成，用來模擬Moho界面的起伏變化，其底界P波速度為6.3km／s，非均勻層下P 波速度為8.1km／s；直立墻（圖2中的亮綠色）相交于模型中心，墻芯P波速度為6.5km／s；在模型的東北、西北、西南、東南4個方位各有一高、低速直立圓柱，高速圓柱（圖2中的深藍色）柱芯的波速最高，為7.5km／s；低速圓柱（圖2中的紅色）柱芯的波速最低，為3.5km／s。觀測系統由986個位于模型頂界（相當于地表）的接收點（圖2中的藍色三角形）和42個激發點（圖2中的紅色五角星）組成。

2.2 方法

數值模擬計算采用澳大利亞國立大學Rawlinson 博士提供的快速波前追蹤成像程序包（http：／／rses.anu.edu.au／～nick／fmtomo.htm）［9］進行，該程序包在求解程函方程中引進迎風有限差分算子（upwind finite difference operators）來近似程函方程中的梯度項，利用窄帶擴展（narrow band expanding）模擬地震波前的傳播，窄帶中激活格點的搜索采用返回指針堆排序算法（Heap sort algorithm with back pointers），在遇到界面時，首先激活界面上的格點，并將其視為新的窄帶，采用分步快速波前追蹤法（Multi－stage FMM），實現對多次折射和反射波場的快速追蹤，加之在模型優化中采用了收斂較快的子空間反演技術［10］，使得快速波前追蹤成像方法不但速度快，而且穩定和高效，目前已經得到許多成功的應用［11－13］。

圖2 3D 結構模型及觀測系統Fig.2 3D model slices，Moho depth contours and observation system

2.3 模擬結果和討論

首先對于§2.1給出的3D 結構模型及觀測系統（圖2），利用FMM 算法獲得折射Pg及反射PmP波走時數據；然后將這些數據作為觀測走時，從均勻層傾斜界面初始模型出發，對§2.1給出的強速度對比結構及Moho界面起伏變化模型進行重建，并對其結果進行分析討論。

2.3.1 速度結構成像

成像的初始模型為一均勻結構層，層速度為Vp＝6.0km／s，反演模型水平網格大小為0.2°×0.2°，垂直間距為2km。分別利用Pg、PmP 和Pg＋PmP走時數據進行反演。圖3（a）、（b）、（c）分別為用Pg、PmP、Pg＋PmP重建速度結構的走時殘差均方根隨迭代次數的變化圖，一般在前4次迭代內，均方根殘差迅速減小，7 次迭代后，均方根殘差趨于穩定，其反演結果作為最終恢復模型。圖4 分別給出了單一震相Pg（圖4（a1）、（a2））、PmP（圖4（b1）、（b2））及折反聯合震相Pg＋PmP（圖4（c1）、（c2））對強速度對比結構成像結果在0、25km 深度的切片圖。

圖3 走時殘差均方根隨迭代次數變化圖Fig.3 Travel time RMS residuals with iteration times

2.3.2 Moho界面成像

由于Pg對Moho結構不能起到直接的約束作用，我們只利用PmP 對Moho界面進行重建。首先用0.2°×0.2°的網格把界面網格化，反演的初始模型為一深度由東部35km 過渡到西部50 km 的傾斜界面，迭代4 次后均方根殘差已不再減小，該模型作為反演的最終結果。圖5（a）、（b）、（c）分別為Moho成像結果的深度、深度誤差等值線圖及走時殘差均方根隨迭代次數變化圖。

圖4 折射（Pg）、寬角反射（PmP）及折反聯合（Pg＋PmP）成像結果Fig.4 Seismic imaging slices with Pg，PmP and Pg ＋PmP

圖5 Moho界面的PmP波成像結果、深度誤差及反演均方根殘差變化Fig.5 Moho imaging，depth errors and RMS residuals with iteration times

2.3.3 分析和討論

比較圖2 中理論模型在0、25km 深度的切片（圖2（a）、（b））與圖4中相應深度由單一震相Pg 對強速度對比結構模型的恢復結果（圖4（a1）、（a2））可以看出，在0km 深度的切片上（圖2（a）、圖4（a1）），Pg反演結果較好，不論是十字直立墻還是高、低速圓柱體以及背景速度結構都得到了很好的恢復，與理論模型極為接近。隨著深度的增加，各切片的分辨能力明顯降低，在25km深度切片上（圖2（b）、圖4（a2）），只能得到低速柱的大體輪廓，直立十字墻及模型外角區域的高速直立柱未能恢復。總體上，淺層的速度結構得到了較好的重建，而隨著深度的增加，其分辨能力越來越差，這是因為Pg射線主要集中在淺層且交叉較為密集，能對其細結構進行較好的約束，可得到很好的反演結果，而對深部結構的分辨能力很差，這與實際Pg波震相記錄所反映的地殼物性結構特征相一致。

圖4（b1）、（b2）給出了0、25km 深度的PmP反演結果，與相應深度的理論模型（圖2（a）、（b））對比可以發現，PmP對整個地殼的結構均有不同程度的分辨能力，雖對中地殼的分辨能力明顯高于地殼淺部和下地殼，但都能給出各深度高、低速柱及直立十字墻的結構輪廓。

圖4（c1）、（c2）為折射Pg和反射PmP（即Pg＋PmP）聯合成像在0、25km 深度的切片，將這些結果與真實模型（圖2）及單一折射震相Pg成像切片（圖4（a1）、（a2））和PmP 成像結果（圖4（b1）、（b2））對比可以看出，Pg＋PmP的成像結果明顯優于單個Pg或PmP，尤其是背景梯度層的恢復最為明顯，總體上給出了更為清晰的強速度對比及梯度層結構特征。

對比圖2（c）中理論模型的Moho深度等值線圖與圖5（a）中Moho界面的PmP 波成像結果可以看出，除邊界部分外，模型的界面形態得到了很好的恢復，在有效探測區域，數據擬合較好（圖5（b）、（c）），這是因為PmP 對于Moho界面的起伏變化比較敏感，對界面有較強的分辨能力。

3 結 語

本文基于最新發展起來的快速波前跟蹤成像方法、深地震測深震相記錄特點及中國地震局地球物理勘探中心已有的深地震測深專用儀器，對單個折射震相Pg、Moho寬角反射震相PmP成像及二者聯合成像對強速度對比結構模型的重建能力進行數值模擬計算，并對其分辨能力進行分析。折射震相Pg 對上地殼細結構有較高的分辨能力，寬角反射震相PmP可對中下地殼及Moho界面的細結構特征進行較好的約束，但對淺部結構的分辨能力較差，而利用Pg、PmP 聯合成像方法，可以有效恢復深、淺部3D 空間的強速度對比結構及Moho界面形態特征，分辨能力較好，可構建較為精細的3D 地殼結構模型。這些研究為開展3D 折射／反射地震聯合試驗探測奠定了必要的基礎，對于復雜構造區深部細結構的3D 主動源地震探查具有重要意義。本數值模型中使用的計算方法、調試的程序模塊及有關參數還可為3D深地震測深觀測系統的設計提供參考。

［1］Teng J W.Investigation of the Moho Discontinuity Beneath the Chinese Mainland Using Deep Seismic Sounding Profiles［J］.Tectonophysics，2013，609：202－216

［2］張先康.阿尼瑪卿縫合帶東段上地殼結構［J］.地震學報，2007，29（6）：592－604（Zhang Xiankang.Upper Crust Structure of Eastern Animaqen Suture Zone：Results of Barkam－Luqu－Gulang Deep Seismic Sounding Profile［J］.Acta Seismologica Sinica，2007，29（6）：592－604）

［3］楊卓欣.華北莫霍面構造形態—深地震測深數據的三維反演［J］.地震地質，2000，22（1）：74－80（Yang Zhuoxin.Structure of Moho in North China－3DInversion from Deep Seismic Sounding Data［J］.Geology and Seismology，2000，22（1）：74－80）

［4］段永紅.華北地區上部地殼結構的三維有限差分層析成像［J］.地球物理學報，2002，45（3）：362－369（Duan Yonghong.Three－Dimensional Finite－Difference Tomography of Velocity Structure of the Upper Crustal in North China［J］.Chinese J Geophys，2002，45（3）：362－369）

［5］白超英，黃國嬌，李忠生.三維復雜層狀介質中多震相走時聯合反演成像［J］.地球物理學報，2011，54（1）：182－192（Bai Chaoying，Huang Guojiao，Li Zhongsheng.Simultaneous Inversion Combining Multiple－Phase Travel Times within 3DComplex Layered Media［J］.Chinese J Geophys，2011，54（1）：182－192）

［6］黃國嬌，白超英.多震相走時聯合三參數同時反演成像［J］.地球物理學報，2013，56（12）：4 215－4 225（Huang Guojiao，Bai Chaoying.Simultaneous Inversion of Three Model Parameters with Multiple Classes of Arrival Times［J］.Chinese J Geophys，2013，56（12）：4 215－4 225）

［7］Sethian J A，Popovici A M.3－D Traveltime Computation Using the Fast Marching Method［J］.Geophysics，1999，64（2）：516－523

［8］Rawlinson N，Sambridge M.Wavefront Evolution in Strongly Heterogeneous Layered Media Using the Fast Marching Method［J］.Geophysical Journal International，2004，156：631－647

［9］Rawlinson N.Structure of the Tasmanian Lithosphere from 3－D Seismic Tomography［J］.Australian Journal of Earth Sciences，2010，57：381－394

［10］Kennett B L N.Subspace Methods for Large Scale Inverse Problems Involving Multiple Parameter Classes［J］.Geophys J，1988，94：237－247

［11］Rawlinson N，Kennett B LN.Teleseismic Tomography of the Upper Mantle Beneath the Southern Lachan Orogen，Australia［J］.Physics of the Earth and Planetary Interiors，2008，167，84－97

［12］Rawlinson N，Urvoy M.Simultaneous Inversion of Active and Passive Source Datasets for 3－D Seismic Structure with Application to Tasmania［J］.Geophysical Research Letters，2006，33：1－5

［13］Rawlinson N.Seismic Wavefront Tracking in 3DHeterogeneous Media：Applications with Multiple Data Classes［J］.Exploration Geophysics，2006，37：322－330