天津—北京高分辨地震折射剖面初至波走時成像①

賈宇鵬，王夫運，田曉峰，段永紅，莘海亮，劉寶峰

（1．中國地震局蘭州地震研究所，甘肅 蘭州 730000；2．中國地震局地球物理勘探中心，河南 鄭州 450002）

0 引言

精細的地殼速度結構研究有助于了解基底和莫霍面形態(tài)，結合地質(zhì)情況可對活動斷層等作出合理的解釋，為研究地震的孕育環(huán)境和發(fā)生機理提供科學的依據(jù)。首都圈地區(qū)位于NE向的太行山隆起，EW向的燕山隆起和華北平原的交匯之處，地質(zhì)情況復雜。已有眾多學者對該區(qū)域做了大量研究工作，得到了很多成果。王夫運等［1］利用三維不分塊層析成像技術得到北京地區(qū)P波S波慢度圖和波速比圖，發(fā)現(xiàn)慢度和速度比的分布有呈NE向和NW向展布的特征。段永紅等［2］利用華北地區(qū)的30多條探測剖面得到了上部地殼的三維速度結構，結果表明華北盆地由兩個大的斷陷帶中間夾一個隆起組成，兩個大的斷陷帶又分為幾個小的斷陷盆地，區(qū)域內(nèi)的斷層多為鏟式正斷層，結晶基底深度變化較大。嘉世旭等［3－4］利用該高分辨剖面數(shù)據(jù)研究了華北東部裂陷盆地與燕山隆起區(qū)的地殼結構和構造耦合，給出了地殼速度模型和斷裂的位置，結果表明剖面西側燕山隆起區(qū)地殼結構高速穩(wěn)定但東側裂陷盆地低速松散。王帥軍等［5］得到了該剖面的二維地殼速度結構，可以清楚看到基底面、殼內(nèi)的反射界面和莫霍面的位置和形態(tài)。劉保全等［6］用深地震反射剖面探測了北京地區(qū)的地殼精細結構。

為了進一步探究該地區(qū)的精細地殼結構和活動斷層構造環(huán)境，2006年初中國地震局地球物理勘探中心在首都圈地區(qū)完成了一條高分辨折射地震探測剖面，取得了優(yōu)質(zhì)的爆破地震資料。本文使用Colin Zelt的正則化反演方法［7－8］，對反映該區(qū)域上地殼結構的Pg波走時數(shù)據(jù)進行成像，研究天津塘沽—北京順義高分辨折射探測剖面的上地殼二維精細速度結構，并探討該地區(qū)的潛伏斷層。

1 觀測系統(tǒng)和數(shù)據(jù)

高分辨地震折射剖面東起渤海之濱的天津塘沽，通過華北坳陷區(qū)北部，到達燕山臺褶帶南部邊緣的北京東北郊順義區(qū)，剖面全長185km。布設觀測點267個，最小觀測點距0．35km，最大觀測點距0．70km，平均觀測點距0．51km；炮點18個，平均炮間距12．15km。構成了較為完整的多重追逐和多重相遇觀測系統(tǒng)。位于測線東南部的華北坳陷區(qū)為新生代強烈坳陷區(qū)，其構造及塊體延伸方向大致為NNE－NE向；位于測線西北部的是北京平原西北部，測線最西北端到達燕山褶皺帶南端（圖1）。

采用全部18炮所拾取1 454個Pg波走時數(shù)據(jù)反演上地殼精細速度結構。如圖2（a）、2（b）、2（c）分別為剖面西北端樁號為297（第18炮），中間樁號225（第10炮）和東南端樁號118（第1炮）炮點的折合速度記錄圖。Pg波到時主要受基底速度和其上的沉積蓋層厚度及蓋層速度的影響。剖面西北段（圖2（a））的北京平原區(qū)由于靠近燕山隆起的南緣，沉積蓋層較薄，基底速度較高，Pg波折合到時約1～2s；而剖面中段（圖2（b））和東南段（圖2（c））的華北平原區(qū)由于受到冀中坳陷和黃驊坳陷沉積層的影響，Pg波的折合到時約1～4s。

圖1 天津—北京高分辨折射／反射測線位置圖Fig．1 Location Map of High Resolution Refraction／Wideangle－reflection Profile from Tianjin to Beijing．

由于沿剖面地質(zhì)構造不同，加之采集數(shù)據(jù)時的一些因素，數(shù)據(jù)在不同地段的拾取誤差就有所不同?？梢钥闯?，圖2三張記錄圖在炮點附近約30km內(nèi)記錄的信噪比還是不錯的，樁號297的炮點要好于樁號118的炮點。通過分析全部18炮的數(shù)據(jù)記錄圖，發(fā)現(xiàn)Pg的信噪比在燕山南緣過渡帶（300～245）較好，追蹤較遠，在冀中坳陷（245～180）和黃驊拗陷（140～115）一般，追蹤較短。第1，2，3，4炮Pg在SE方向追蹤到側線端點，NW方向追蹤到冀中坳陷東南端（樁號180）附近；第5，6，7炮SE方向受黃驊拗陷的影響，追蹤較短，到達樁號140附近，NW方向能穿越大半個冀中坳陷，追蹤到樁號220附近；冀中坳陷內(nèi)部的第8，9，10炮Pg可完全追蹤覆蓋冀中坳陷；第11～18炮東南方向追蹤到樁號190附近，基本穿過了冀中坳陷，NW方向由于是燕山南緣過渡帶，可清晰追蹤到剖面的西北端點。

2 速度成像

2．1 成像原理

傳統(tǒng)上采用有限差分求解程函方程的正演算法，已有較多文獻進行討論［9－14］。對于反演成像方法，John Hole在1992年提出了將走時殘差分布在整個射線上的反投影算法，由于反投影算法無需進行矩陣計算，反演效率很高［9］。Toomey指出，反投影算法將殘差分布在整條射線上進行計算，對于射線密集的高速區(qū)域會帶來較大計算誤差［12］，Colin Zelt提出了修正的反投影算法以避免高速區(qū)射線集中帶來的誤差累積［7］。

本文采用Colin Zelt的正則化反演方法（FAST）［7］對剖面的上地殼速度結構進行精細成像。該方法在解病態(tài)方程組時在數(shù)據(jù)之外引入一些約束條件去處理解的欠定部分或防止數(shù)據(jù)過度擬合，這些約束條件通常是對解的復雜性的與約束性措施，如成像中的平滑度控制。因此，正則化過程在一定程度上可視為反演走時曲線而非擬合單點走時［9］。由此，可構造一個包含速度模型平滑度和數(shù)據(jù)擬合度的目標函數(shù)［7－8］：

其中m是模型矢量；δt是數(shù)據(jù)殘差；Cd是數(shù)據(jù)協(xié)方差矩陣；Ch和Cv分別是水平和垂向平滑度矩陣；λ是數(shù)據(jù)擬合與模型平滑度協(xié)調(diào)因子；sz是水平與垂向平滑度權重因子。于是，每一次的線性迭代問題可歸結為求解一個δm，使得目標函數(shù)最?。?］：

式中L是目標函數(shù)偏導數(shù)矩陣；m0是當前模型；δm是待求模型擾動，新模型矢量m＝m0＋δm。方程（2）可用LSQR算法求解。

采用橫向和縱向不同寬度的網(wǎng)格大小，同時合理選取正則化因子可以對剖面上部分無射線覆蓋區(qū)可以較好地給出合理的平滑解，得出整個研究區(qū)在橫向和縱向較合理的最小格點最平滑解［7］。

2．2 成像過程

2．2．1 初始模型

初始速度模型的選取與研究區(qū)域的地質(zhì)構造密切相關，該區(qū)域坳陷和隆起區(qū)相間，構造復雜。為了得到可靠的反演結果，選取了四個十分不同的初始模型進行試驗，如圖3。

圖3 一維初始速度模型Fig．3 1DStarting velocity models．

2．2．2 正演、反演參數(shù)

初至波成像方法的橫向分辨很大程度取決于觀測點距和爆破點的密集程度，而縱向分辨除此之外還取決于研究區(qū)域殼內(nèi)的速度梯度。由于剖面是高分辨的觀測系統(tǒng)，最小觀測點距達到0．35km，最大觀測點距0．70km，平均觀測點距0．51km，平均炮間距12．15km，所以正演模型采用較密集0．25km×0．25km的格點間隔?？紤]到反演單元格橫向?qū)挾纫M可能接近觀測點距和縱向上該區(qū)域的沉積盆地導致從地表到基底有較大速度梯度，因此反演采用0．50km×0．25km的單元格。反演迭代次數(shù)可以由χ2到達1附近而且隨迭代變化很小而確定。經(jīng)試驗迭代20次就能夠達到此要求，所以選取20次迭代。根據(jù)以往的經(jīng)驗，水平與垂向平滑權重因子Sz暫定為0．20。由此來討論初始模型對反演結果的影響。

圖4 不同初始模型得到的速度圖Fig．4 Testing of different starting velocity models．

由表1可以看出，使用十分不同的初始速度模型進行反演計算，初始均方根殘差區(qū)別很大，但是經(jīng)過20次迭代后，最終均方根殘差都達到65ms左右。圖4中，兩條黑線為3．5km／s和4．5km／s的速度等值線，較密集的三條白色線為5．8km／s，5．9 km／s，6．0km／s的速度等值線?？梢钥闯觯╝）速度圖和其余三個（b），（c），（d）速度圖有一定的差別，但是差別主要體現(xiàn)在射線覆蓋不夠的底部區(qū)域，四張速度圖的上部區(qū)域大致形態(tài)還是相似的。這說明反演結果對初始模型的依賴性不大，走時數(shù)據(jù)起決定性作用。說明該方法（FAST）比較可靠，可以用來處理高分辨觀測系統(tǒng)的資料。

圖5 不同平滑權重因子Sz得到的速度圖Fig．5 Testing of different vertical versus horizontal model smoothness factor Sz．

表1 不同初始模型對結果的影響

表2 不同平滑權重因子Sz對結果的影響

參考已有的研究成果，依據(jù)該地區(qū)地殼模型的先驗信息，從中選取初始速度模型B來進一步討論水平與垂向平滑權重因子Sz對反演結果的影響，反演依然迭代20次。

表2中不同的Sz值對模型的初始、最終均方根殘差影響比較抽象，而圖5就比較直觀了?？梢悦黠@的看出Sz對兩條黑線為3．5km／s和4．5km／s的速度等值線影響較小，但對可能反映基底形態(tài)的三條白色線5．8km／s，5．9km／s，6．0km／s的速度等值線的影響較大。考慮到射線的覆蓋密集程度及未采樣區(qū)域，每一組閉合白線的上部更能夠反映基底形態(tài)。參考已有的研究成果，最終選取Sz為0．20。

2．2．3 成像結果

圖6（a）為初始速度模型B和Sz為0．20的最終P波速度模型圖，上下兩個圖分別為深度方向擴大比例為4和1。其初始均方根殘差為664．60ms，最終均方根殘差達到64．90ms。圖中五角星為爆炮點，深色虛線為潛伏斷裂，淺色虛線為地層尖滅線。

2．3 檢測板分析

采用檢測板方法對最終模型（圖6（a））的分辨進行評估。在最終速度模型加入5%的速度擾動后，使用該速度模型作為檢測板輸入模型，由此產(chǎn)生正演數(shù)據(jù)并加入5%的噪音進行反演恢復，經(jīng)過20次迭代得到檢測板恢復模型。圖7是網(wǎng)格大小為2．5km×0．5km，5km×1km，10km×2km和20 km×5km的恢復模型。

可以看出，20km×5km和10km×2km網(wǎng)格的測試在整個模型上恢復的很好，分辨完全可以達到網(wǎng)格的大??；5km×1km的網(wǎng)格能較好的恢復0～4km深度；2．5km×0．5km的網(wǎng)格在模型的淺部都不能恢復，恢復情況不好。圖7（b）顯示模型分辨在燕山南緣過渡帶要優(yōu)于幾個拗陷的內(nèi)部。

圖6 最終速度模型及射線和殘差分布Fig．6 The final velocity and distribution of ray paths，traveling time residues．

3 結論和討論

圖6（a）中，上地殼速度分層明顯，基本反映出了該地區(qū)不同構造單元的沉積厚度和基底形態(tài)。白堊紀沉積，第三系沉積和基底的埋深極為相似，這說明該地區(qū)受到第三系以來的地殼運動的影響較大。

圖7 不同網(wǎng)格大小的檢測板試驗恢復模型Fig．7 Traveltime checkboard resolution tests using difference grid sizes．

兩條黑線為3．5km／s和4．5km／s的速度等值線，可能分別對應于第三系以來沉積層厚度和早第三系與白堊紀沉積巖分界，較密集的三條白色線為5．8km／s，5．9km／s，6．0km／s的速度等值線，結合圖6（b）的射線在基底為滑行波，或許反映了沿剖面基底界面的大致形態(tài)?；鬃顪\處約2km，最深處約達到8km。剖面樁號290～275的北京凹陷基底埋深明顯較深；剖面樁號275～245通縣凸起區(qū)基底埋深約為2km；剖面樁號245～225的大廠凹陷基底埋深約5km；剖面樁號225～180的武清凹陷，基底埋深約4～8km；剖面樁號180～140的滄縣隆起基地埋深約4km；剖面樁號140～115的黃驊拗陷，基底埋深約6km。

斷裂和底層尖滅線一般出現(xiàn)于速度梯度變化較大的界面，參考其他研究所給出的斷裂和底層尖滅線位置［1－5，15］，結合本研究所得到的天津—北京的高分辨地震折射剖面的上地殼精細速度結構（圖6（a）），給出了該地區(qū)的潛伏斷層與地層尖滅線的大致位置。F1懷柔—涿縣斷裂和F2順義良鄉(xiāng)斷裂位于北京凹陷、通縣突起和大廠凹陷的交界轉(zhuǎn)換區(qū)域，這里速度梯度變化很大；F3夏埑斷裂和F4寶坻銅川斷裂依據(jù)為在3．5km／s和4．5km／s的速度等值線左右兩邊和沿著代表基底形態(tài)的白色的速度等值線迅速下降；F5底層尖滅線在地表兩端分別的速度為1．6km／s和2．3km／s，速度差異較大，在3．5 km／s和4．5km／s的速度等值線左右兩邊和基底的速度等值線處有明顯變化，也可認為F5是滄東斷裂；F6滄東斷裂依據(jù)為從地表到基底速度明顯變化的區(qū)域。

圖6（b）最終模型射線追蹤圖可以看出，射線分布在地表極為密集，基底以上較為密集，對所研究區(qū)域上地殼有較好覆蓋。在基底處，射線表現(xiàn)為滑行波，在此圖中有明顯表現(xiàn)。圖中一共1 454根射線，射線完全密集覆蓋所研究區(qū)域，成像結果可靠。

圖6（c）和（d）為初始和最終走時殘差。經(jīng)過20次迭代后，走時殘差分布迅速減少到約200ms以下，這說明最終模型數(shù)據(jù)與觀測數(shù)據(jù)獲得了較好的擬合。

總之，成像結果的精細程度與可靠性取決于觀測系統(tǒng)是否為高分辨，數(shù)據(jù)質(zhì)量較好和使用合適的成像方法。通過上文的分析，這三方面都達到要求。在使用Colin Zelt正則化反演方法（FAST）的過程中，分析討論了不同初始模型和參數(shù)變化對反演結果的影響，并采用檢測板方法對最終模型的分辨進行評估。因此可以說，本文中給出的天津－北京剖面的地殼上地殼精細速度模型（圖6（a））是高分辨的，并且是比較可靠的。

致謝：感謝天津城市活斷層試驗探測項目的野外采集和室內(nèi)處理人員的辛勤勞動。感謝Colin Zelt教授提供的FAST軟件包。

［1］ 王夫運，張先康，陳棋福，等．北京地區(qū)上地殼三維細結構層析成像［J］．地球物理學報，2005，48（2）：359－366．

［2］ 段永紅，張先康，方盛明．華北地區(qū)上部地殼結構的三維有限差分層析成像［J］．地球物理學報，2002，45（3）：362－369．

［3］ 嘉世旭，張成科，趙金仁，等．華北東北部裂陷盆地與燕山隆起地殼結構［J］．地球物理學報，2009，52（1）：99－110．

［4］ 嘉世旭，齊 誠，王夫運，等．首都圈地殼網(wǎng)格化三維結構［J］．地球物理學報，2005，48（6）：1316－1324．

［5］ 王帥軍，張先康，張成科，等．武清—北京—赤城二維地殼結構和構造［J］．地球物理學報，2007，50（6）：1769－1777．

［6］ 劉保金，胡平，孟勇奇，等．北京地區(qū)地殼精細結構的深地震反射剖面探測研究［J］．地球物理學報，2009，52（9）：2264－2272．

［7］ Colin A Zelt，P J Barton．3Dseismic refraction tomography：A comparison of two methods applied to data from the Faeroe Basin［J］．J．Geophys．Res．，1998，103，7187－7210．

［8］ Priyank Jaiswal，Colin A Zelt，et al．2－D traveltime and waveform inversion for improved seismic imaging：Naga Thrust and Fold Belt，India［J］．Geophys．J．Int．，2008，173，642－658．

［9］ Hole J A．Nonlinear High－Resolution Three－Dimensional Seismic Travel Time Tomography［J］．J．Geophys．Res．，1992，97（B5）：6553－6562．

［10］ Vidale John．Finite－difference calculation of travel times［J］．Bulletin of the Seismological Society of America，1988，78（6）：2062－2076．

［11］ Vidale John．Finite－difference calculation of traveltimes in three dimensions［J］．Geophysics，1990：55，521－526．

［12］ Toomey D R，Solomon S C，Purdy G M．Tomographic imaging of the shallow crustal structure of the East Pacific Rise at 9°30 3/4N［J］．J．geophys．Res．，1994，99，24135－24157．

［13］ 張元生，徐果明．聯(lián)合利用走時與波形反演技術研究地殼三維速度結構（Ⅰ）：理論與方法［J］．西北地震學報，1998，20（2）：8－15．

［14］ 張元生，徐果明．聯(lián)合利用走時與波形反演技術研究地殼三維速度結構（Ⅱ）：應用［J］．西北地震學報，1998，20（3）：44－51．

［15］ 王椿庸，張先康，丁志峰，等．大別造山帶上部地殼結構的有限差分層析成像［J］．地球物理學報．1997，40（4）：459－501．