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

劉 妍 徐朝繁

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

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

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

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

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

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

2 數(shù)值模擬試驗

2.1 理論模型及觀測系統(tǒng)

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

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

2.2 方法

數(shù)值模擬計算采用澳大利亞國立大學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），實現(xiàn)對多次折射和反射波場的快速追蹤，加之在模型優(yōu)化中采用了收斂較快的子空間反演技術(shù)［10］，使得快速波前追蹤成像方法不但速度快，而且穩(wěn)定和高效，目前已經(jīng)得到許多成功的應(yīng)用［11－13］。

圖2 3D 結(jié)構(gòu)模型及觀測系統(tǒng)Fig.2 3D model slices，Moho depth contours and observation system

2.3 模擬結(jié)果和討論

首先對于§2.1給出的3D 結(jié)構(gòu)模型及觀測系統(tǒng)（圖2），利用FMM 算法獲得折射Pg及反射PmP波走時數(shù)據(jù)；然后將這些數(shù)據(jù)作為觀測走時，從均勻?qū)觾A斜界面初始模型出發(fā)，對§2.1給出的強速度對比結(jié)構(gòu)及Moho界面起伏變化模型進行重建，并對其結(jié)果進行分析討論。

2.3.1 速度結(jié)構(gòu)成像

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

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

2.3.2 Moho界面成像

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

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

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

2.3.3 分析和討論

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

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

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

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

3 結(jié) 語

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

［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］張先康.阿尼瑪卿縫合帶東段上地殼結(jié)構(gòu)［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］楊卓欣.華北莫霍面構(gòu)造形態(tài)—深地震測深數(shù)據(jù)的三維反演［J］.地震地質(zhì)，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］段永紅.華北地區(qū)上部地殼結(jié)構(gòu)的三維有限差分層析成像［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］白超英，黃國嬌，李忠生.三維復(fù)雜層狀介質(zhì)中多震相走時聯(lián)合反演成像［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］黃國嬌，白超英.多震相走時聯(lián)合三參數(shù)同時反演成像［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