基于多震相聯合反演江蘇地區一維速度模型研究

李細兵 趙啟光 朱峰 洪海春 王恒知

江蘇省地震局，南京市衛崗3號 210014

0 引言

傳統的利用天然地震走時獲取一維速度模型主要有Herglotz-Wiechert方法、線性擬合方法、τ（p）方法等（Shearer，2009），這些方法使用的前提是震源位置較為準確，這是因為其主要是利用臺站與震源之間的走時資料來反演一維速度結構。除了單獨反演速度結構外，現在有很多種方法能同時反演速度結構和震源參數，其中，使用較為廣泛的是Kissling等（1994）提出的最小一維速度模型方法，該方法已廣泛應用于國內外許多地區的速度模型反演中，并取得了豐富的研究結果（于湘偉等，2003；趙旭等，2007；李鉑等，2012；王小娜等，2015；柯乃琛等，2016）。此外，田玥等（2006）采用擬牛頓法與信賴域法結合的方法對北京地區地震震源和速度模型進行聯合反演，取得了較好的效果。但是，同時反演震源位置和速度結構需要在地震的定位誤差和速度的擾動足夠小的條件下，才能使得包括震源和速度在內的未知數聯合方程是線性的；否則，定位的誤差將會直接影響速度結構的準確性，進而造成速度結構的不準確，因此，這是進行震源和速度同時聯合反演的首要條件。實際上，由于地球介質存在橫向非均一性，大多數一維速度模型的構建過程都是各層介質速度平均化的過程。

準確的三維速度模型是諸如大尺度構造應力場、地球動力學等研究的基礎性資料，而較好的一維速度模型不僅可以提高地震定位的精度，還可以為各種三維地震成像結構提供可靠的初始模型。地震定位時很多地區的一維速度模型是建立在深反射地震剖面基礎上的，通過深反射地震剖面能得到較為準確的二維地殼速度結構和莫霍面深度信息，楊文采等（1999b）對穿過江蘇北部的郯城-漣水綜合地球物理剖面進行了深地震反射和大地電磁調查，但并未給出該反射剖面的具體深度和速度特征；在另一條從諸城到日照方向的剖面給出了地殼到地幔的速度結構（楊文采等，1999a），并將各深度的速度結構與巖性相對應起來，該速度結構可為江蘇北部地區的速度結構提供一種參考。王有學等（2004）利用布置在蘇北的中國大陸科學鉆探井的2條折射／反射地震剖面研究認為，該鉆探場址區的地殼速度結構可分為上地殼、中地殼、下地殼3層，各層厚度均為10km左右，整個地殼的平均速度為6.3 km／s。從符離集（安徽）到奉賢（上海）的HQ-13寬角反射地震剖面穿越整個江蘇中部地區，該速度結構在江蘇中南部具有一定的代表性，白志明等（2006）對該剖面進行了新的處理和解釋，揭示了該剖面上地殼深度為10～13km，速度為5.9～6.0km／s；中地殼深度為21～25km，速度為6.3～6.4km／s；下地殼分為2層，上層深度為25～29km，速度為6.6～7.0km／s，下層深度為30～36km，速度為6.9～7.0km／s。在層析成像方面，很多研究者使用的初始一維速度模型參考了周邊區域的深地震反射剖面結果（徐佩芬等，2000；徐紀人等，2003；黃耘等，2011）。對于江蘇地區，黃耘等（2011）綜合各種深地震反射資料結果，將該區域的地殼及上地幔分為7層速度模型，以此作為地震定位和區域三維速度反演的初始模型。盡管建立在深地震反射剖面基礎上的初始模型有較好的約束能力，具有一定的代表性，但往往成像覆蓋的地區范圍遠遠超過深地震反射所能覆蓋的區域，因此，存在一定的不足。如果能采用具有較好定位水平的地震數據本身構建出一維P波速度模型，由于地震射線覆蓋整個研究區，則能有效代表研究區的平均速度結構。在地震定位方面，目前江蘇臺網常規地震定位使用的一維速度模型主要是全球iasp91模型。這種2層模型具有分層簡單、能追蹤到多種震相等特點。但是，該模型是建立在整個全球尺度基礎上的，對于具體地區的地震定位水平還有較大的提升空間。

鑒于江蘇地震臺網記錄的震相資料中含有大量的續至Pg波震相，這些震相在傳統的走時反演速度結構中很少利用，本文同時利用初至波、續至波震相，采用Pg、Pn波震相走時聯合反演江蘇地區一維P波速度模型，以期為后續江蘇地區三維速度結構以及地震定位研究提供較可靠的一維速度模型。

1 數據與方法

本文收集了2008～2016年間江蘇地震臺網記錄的地震P波震相數據，其包含有Pn、Pg震相走時數據，所有地震都給出了震源深度，同時具有相對較好的定位精度。為了進一步提高數據的精度，同時考慮到震源向下沿莫霍面折射滑行產生Pn波的特性，對數據進行了篩選，設置震源深度小于30km，每個地震的P波記錄不少于10個，這可使得震源具有相對較高的定位精度，從而提高數據走時的精度（圖1）。經過篩選后，共有76個臺站、661個地震，總的走時數據約14000條，其中，來自地殼的Pg震相走時數據約10000條，其余為Pn震相走時數據。這些地震的震源深度小于15km的占總數70%，說明江蘇地區絕大多數地震發生在中上地殼以內。由圖1可見，江蘇境內絕大部分地區都有地震記錄，在附近海域也有大量的地震分布；沿郯廬斷裂帶山東境內有較多的地震分布，而在江蘇境內的宿遷段地震數量明顯偏少，遠不如江蘇其他地區，在其右側的徐州和其左側的鹽城地區都有較多的地震發生，但震級都較小。從篩選后的時距曲線（圖2（a））可見，臺網記錄中含有大量來自地殼的續至Pg波，包含初至波和續至波在內的地殼Pg波與上地幔頂部Pn波震相區分較為明顯，當震中距大于300km時，續至波與Pn波逐漸分離開來，但是由于走時拾取等原因，時距曲線具有十分明顯的走時誤差，我們將采用新方法對其進行處理。

圖1 研究區2008～2016年地震空間分布與震源深度統計

本文采用Xu等（2010）提出的聯合Pg和Pn方法來反演江蘇地區地殼及上地幔一維P波速度結構。該方法對地殼Pg波震相和上地幔頂部Pn波震相走時分別采用不同方法計算，對前者采用直接利用球面坐標系下偽彎曲射線追蹤方法（Koketsu et al，1998），對后者利用Pn波震相的折射波性質分段計算走時。

對于地殼Pg波射線，從第j個震源到第k個臺站的走時殘差表示為

圖2 不同震相數據篩選前后對比

行的走時和從莫霍面到臺站的走時

式（2）、（3）中的Δuk為莫霍面深度變化引起的走時變化；Δtk為震源延遲項，主要包括事件深度和開始時刻等引起的震源誤差；Δhk為莫霍面深度的變化量；vc、vm分別為地殼底部的速度和上地幔頂部的速度。

初始模型選用的地殼速度為6.3km／s，是由Pg震相時距曲線進行擬合得到的，經過多次試算，最終選取上地幔頂部的速度為8.0km／s、深度為35km的反演結果較為理想。初始模型一共分為8層（含莫霍面），間隔為5km。在進行走時計算時，為了正確區分震相并去除拾取錯誤的震相，采用了一種自動識別震相的方法。考慮到江蘇地區的莫霍面深度在35km內，因此，設置初始距離為0.8°，當震中距小于該值時，只考慮式（1）中地殼Pg波走時的計算；而當震中距大于0.8°時，分別利用式（1）、（2）計算Pg、Pn走時。當單個走時殘差大于4s時，則該射線走時不參與反演，同時計算Pg、Pn波震相的走時殘差，當兩者的殘差差值小于0.5s時，也不考慮該走時，最后根據各自走時殘差的大小來確定該走時屬于哪種震相。

反演方法采用阻尼最小平方法（LSQR）（Paige et al，1982），對式（1）、（2）進行聯合求解，使下式達到最小，另外，為了克服方程求解的病態和不穩定性，加入了阻尼因子λ和平滑因子μ

其中，第1項為方程的最小平方殘差；第2項為阻尼項，用于控制整個解的穩定性；第3項為平滑項，用于控制相鄰網格之間解的平滑性。

圖3 反演前后殘差分布

2 結果與討論

如圖2（a）所示，目前由于各種原因很多震相數據的走時數據存在各種偏差，這很容易出現震相標識錯誤，將屬于地殼Pg波的震相標識為Pn波震相。由圖2（b）可見，對于Pn、Pg震相走時十分接近的臨界距離（震中距120～180km），難以進行有效的篩選；而對于震中距大于300km的震相，利用本文方法能對Pn、Pg震相進行十分有效的篩選識別，這種篩選提高了數據的使用精度，從而有助于提高后續層析成像結果的可靠性。

每次聯合反演都將上次反演結果作為初始模型進入下一次的聯合反演，經過4次迭代，反演結果趨于穩定，走時數據殘差的均方根值由1.21s降為0.62s，殘差減少了52%，與初始模型殘差分布相比，反演后的殘差以零值為基線基本呈正態分布（圖3（c）），數據殘差減小十分明顯。從反演前后的殘差與震中距間的關系來看（圖3（a）、3（b）），初始模型的殘差對于小震中距震相而言，以正殘差為主，說明淺層的速度偏大；而對于大于500km的震中距震相（主要為Pn震相），則以負殘差為主，說明初始模型的莫霍面速度偏小；經過反演大、小震中距下的走時殘差都得到了有效校正，基本沿殘差零值兩側分布。圖3（d）中顯示出Pg、Pn波折合走時隨震中距的增加震相區分越加明顯。圖4為反演前后一維速度模型分布。由圖4可見，反演后的P波速度為5.14～6.41km／s，莫霍面厚度為33.86km，上地幔頂部速度為8.06km／s。其中，淺層0～5km的速度最小，5～15km之間速度變化很小，為5.88km／s；與下地殼相比，中、上地殼20km以內的速度較小，整個地殼的平均速度為6.0km／s，與全球速度模型iasp91相比，顯著偏低。

圖4 反演前后一維速度模型分布

由于阻尼最小平方法（LSQR）不能對反演結果給出誤差估計，因此，本文采用常規Bootstrap方法對所有數據進行隨機采樣并進行相同參數下的反演，重復100次得到本次反演模型的誤差估計值。地殼不同深度下的速度誤差分別為0.036km／s（0km）、0.028km／s（5km）、0.022km／s（10km）、0.025km／s（15km）、0.018km／s（20km）、0.013km／s（25km）、0.052km／s（30km）；上地幔頂部速度和莫霍面深度的誤差分別為0.0035km／s、0.114km。

為了探討該模型是否對江蘇地區具有一定的代表性，選取了目前江蘇及其鄰區常用的幾種模型（圖5），計算相同地震數據下不同模型的殘差結果并對其進行簡單比較。此次一維P波速度模型只包括地殼、上地幔頂部的速度以及莫霍面深度參數。選取的模型中用于日常定位的有3種，分別為iasp91全球模型、華北模型和華南模型等；用于層析成像的模型主要來自于HQ13測線（綜合模型）（黃耘等，2011）。定位模型中，iasp91模型是目前江蘇臺網日常使用的，而華南模型目前在相鄰的安徽地區作為地震定位速度模型使用（謝石文等，2016）。

幾種模型中，iasp91模型屬于全球速度模型，華北模型和華南模型主要是針對特定地區地震定位的模型，因此，這幾種模型的應用范圍并不能代表江蘇地區的速度模型。由幾種模型計算得到的所有殘差均方根值可見，本文反演的一維P波速度模型殘差均方根最小；深地震反射獲取的速度模型（綜合模型）的殘差最大，這可能說明江蘇地區具有較為復雜的地質構造分布，單一的反射剖面不能很好地約束整個區域的地震走時。我們將此次反演模型與其他模型的結果進行比較，結果如圖6所示。由圖6可見，深地震反射速度模型能擬合的數據最少且殘差最大，其次為華南模型，然后是華北模型，最好的是iasp91模型。iasp91模型的走時殘差均方根值為0.88s，比此次反演模型的多近30%。由此可見，本文采用的聯合反演獲取的一維P波速度模型對初至波、續至波震相數據擬合最好。

圖5 研究區4種一維P波速度模型與反演速度模型對比

3 結論

本文利用地震走時數據采用聯合反演方法獲取了江蘇地區的一維P波速度模型，得到以下認識。

（1）與僅采用初至波走時的傳統天然地震走時獲取方法相比，該方法同時使用初至波、續至波震相走時，充分利用了大量存在的續至波參與反演，能有效提高中下地殼的反演能力。針對常規地震震相目錄中常存在的震相標識錯誤問題，采用的自動判別篩選震相方法能最大限度提高數據走時的精度，可以對不同震相進行有效區分。

圖6 研究區4種一維P波速度模型與反演模型的殘差分布

（2）通過與其他速度模型進行簡單比較，反演模型對Pg、Pn震相走時擬合效果最佳，具有最小的殘差。當所用走時數據定位精度較高時，該聯合反演方法能為研究區三維速度結構成像和地震定位提供較可靠的一維速度模型。此外，目前各臺網使用諸如華南或華北等模型作為地震定位的一維速度模型，盡管它們是簡單的二層速度模型，但卻有較高的定位效率。而本方法速度模型是多層的，后續工作將對該模型作進一步改進，并考慮運用爆破地震來有效檢驗該模型的定位能力。