王海燕

摘 要：本文利用二維介質(zhì)中的射線追蹤完成在地球內(nèi)部地震波的三維的射線追蹤，根據(jù)走時殘差，建立地震波旅行時層析成像系統(tǒng)方程組。通過川滇地區(qū)寬頻帶地震儀所記錄的遠震P波走時殘差，獲得該區(qū)域的三維P波速度結(jié)構(gòu)，同時結(jié)合當?shù)氐刭|(zhì)構(gòu)造和已有的層析成像結(jié)果，對本論文的成像結(jié)果進行分析，并對青藏高原東緣地區(qū)深部速度結(jié)構(gòu)和殼慢動力學特點進行討論與研究。

關(guān)鍵詞：三維地震層析成像；射線追蹤；地殼結(jié)構(gòu)

1 前言

地球內(nèi)部蘊含著豐富的資源的同時其內(nèi)部的構(gòu)造活動又給人類帶來巨大的災(zāi)難，而利用天然地震所產(chǎn)生的地震波信息，已經(jīng)獲得了許多對地球科學產(chǎn)生深遠影響的發(fā)現(xiàn)，比如地球內(nèi)部成層結(jié)構(gòu)的建立、深源地震的發(fā)現(xiàn)等[1-5]。由于地震波的速度取決于介質(zhì)的密度、物質(zhì)組成、彈性模量、溫度等因素，因而地震波場所攜帶信息的差異反應(yīng)了地球內(nèi)部介質(zhì)特性的變化，地震層析成像技術(shù)就是利用地震波的到時、波形或頻散特征等，建立地球內(nèi)部的三維速度結(jié)構(gòu)影像。但是由于研究對象的特殊性，如震源定位不是很準確，地震波射線的透射角度有限，在非均勻介質(zhì)中地震波傳播的復(fù)雜性等種種原因，致使地震層析成像技術(shù)很難達到醫(yī)學CT的效果。

2 IASP91地球模型中的射線追蹤

如果我們把橢圓的地球近似看做球形，且內(nèi)部成層狀構(gòu)造，可以認為地球是由具有對稱性的無限多個均勻的同心圈層構(gòu)成，地震波在地球介質(zhì)內(nèi)部傳播，因而，可以利用二維介質(zhì)中的射線追蹤完成地震波在地球內(nèi)部的射線追蹤，包括地震波在圈層介質(zhì)分界面上的反射波等震相的追蹤。

根據(jù)IASP91地球模型的速度參數(shù)表可知，其內(nèi)部速度是地球半徑的函數(shù)，并且每個地層的速度梯度不等，為保證地層速度變化的連續(xù)性，網(wǎng)格結(jié)點之間的速度值利用三次樣條插值獲得，進而對介質(zhì)的縱向非均勻特征進行描述。

針對上述地球模型，我們可以將地震波在地球內(nèi)部傳播的縱波及橫波速度分別表示為

（1）

在二維介質(zhì)的情況下，地震波傳播的射線軌跡滿足程函方程

（2）

可以推導(dǎo)出二維情形下射線軌跡滿足的一階微分方程組。

在數(shù)值計算中，利用數(shù)值微分可以寫為

（3）

其中（x0，z0）為射線起始點E的坐標（圖1），為射線的離源角，及為速度場在x、z兩個坐標軸方向上的速度梯度，為給定的時間步長。射線經(jīng)過時間后，射線將到達A點，其空間坐標為（x，z），離源角為。當A點位于某一界面時，需要根據(jù)斯奈爾定律來計算經(jīng)過折射后的離源角。

圖1 射線路徑示意圖

（1）當射線由震源E以離源角α0出發(fā)，按規(guī)定的步長dx逐步前進，到達A點，計算A點的坐標，以及走時TA，并判斷點A是否到達或者超過第一層下界面。

（2）若A點沒有達到下一層界面，則根據(jù)點所在地層的層序數(shù)，點到地心的距離R1，獲取射線當前所在位置的速度值以及震相。

（3）根據(jù)層序數(shù)判斷是射線是否到達地層的最底層，并對地層或塊體內(nèi)部進行追蹤，確定射線傳播的路徑及旅行時間，并返回射線終點的坐標及離源角α等參數(shù)。

（4）從A點出發(fā)，繼續(xù)以步長dx前進，到達B點，計算B點的坐標以及走時TB，同樣進行判斷B點所處位置，若依然沒有到達或者超過界面，則重復(fù)以上計算，直到射線到達和超過界面為止；若射線超過界面，則采用“二分法”求射線與界面的交點，以及交點處界面傾角，并判斷是否進行廻折波射線追蹤。

（5）如果到達的不是目的界面，則要計算出入射線與法線的夾角、折射線與橫軸的夾角等相關(guān)參數(shù)；射線到達目的界面后，則進行反射波上行射線追蹤，直到追蹤到地面為止。

（6）改變震源處的離源角，重復(fù)上述步驟就得到了一條條經(jīng)過地下介質(zhì)到達地面的射線。

3 層析成像的基本原理

地震波旅行時層析成像是典型的地球物理反演問題，依據(jù)地震波穿過介質(zhì)內(nèi)部所需的時間來重建介質(zhì)的速度結(jié)構(gòu)。

通常情況下，射線的走時t 可以寫成：

（4）

其中，i是射線條數(shù)，l是射線元，v、s為分別介質(zhì)的速度、慢度函數(shù)。

如果假設(shè)介質(zhì)的慢度分布函數(shù)為，其中為介質(zhì)的慢度初始模型，為介質(zhì)的慢度擾動量，那么，根據(jù)（4）式有

（5）

從而有

（6）

所以，由上式可知，如果知道地震波的走時擾動量（走時殘差），就可以研究介質(zhì)中慢度（速度）的擾動分布特征。將積分離散化則層析成像系統(tǒng)的方程組具有如下形式

（7）

寫成矩陣，有

D=CU （8）

其中為走時殘差，由實際觀測值得到；為介質(zhì)分塊中的慢度分布，為待求向量；C為不同單元格內(nèi)所有線段的長度所形成的一個大型稀疏系數(shù)矩陣，由射線的傳播路徑確定。

（9）

4 地震層析成像的應(yīng)用

青藏高原的形成原因一直備受國內(nèi)外地學家的關(guān)注，尤其是青藏高原巖石圈的深部結(jié)構(gòu)。2003-2004年期間，美國麻省理工大學與國土資源部成都地質(zhì)礦產(chǎn)研究所合作，在三江地區(qū)進行了寬頻帶地震觀測，本文利用遠震P波走時殘差，獲得該區(qū)域的三維P波速度結(jié)構(gòu)。

我們研究的區(qū)域范圍為24.0°N～31.0°N，99.5°E～104.0°E，本文使用了25臺寬頻帶地震儀所記錄的地震數(shù)據(jù)，從中挑選出了信噪比較高的714個地震事件，震級大于4.5級，震中距為25.0°～180.0°，震源深度為0～641km，并且每一個遠震事件至少有10個臺站具有清晰的地震記錄，每個記錄的走時殘差絕對值小于5s。本文按照上文闡述的射線追蹤以及層析成像方法對研究區(qū)域?qū)掝l帶地震數(shù)據(jù)記錄進行地震波旅行時層析成像。

根據(jù)眾多學者對三江地區(qū)的研究資料 [5-10] 可知，該區(qū)域地殼比較厚，莫霍面的深度達到了70km左右，并且該地殼有上、中、下三層結(jié)構(gòu)，其厚度分別為10km、25km、35km，所以本文在IASP91地球模型基礎(chǔ)上稍微做了修改，修改后的IASP91地球模型參數(shù)我們稱為初始模型。網(wǎng)格在經(jīng)度方向和緯度方向上均以1.0°間距作水平方向的劃分，在深度方向上以10、20、30、40、50、65、80、100、120、150、200、250km進行劃分。

我們對三江地區(qū)觀測數(shù)據(jù)進行層析成像，則由0～125 km深度處的速度擾動水平剖面圖中可以看出，川西高原地殼和上地幔的速度存在明顯的橫向不均勻性，川中塊體的速度相對較高，而川青塊體以及川滇菱形塊體的速度相對較低，并且川滇菱形塊體相對較低的速度區(qū)域呈南北向分布；大型斷裂帶兩側(cè)存在明顯的速度差異，如龍門山斷裂帶、安寧河斷裂帶、金沙江斷裂帶等為研究區(qū)域中幾大塊體的速度分界線。

縱觀圖2的4個深度剖面，我們發(fā)現(xiàn)，在金沙江斷裂帶和安寧河斷裂帶之間的研究區(qū)域，隨著深度的增加，特別是深度150～250km區(qū)間，低速度區(qū)域的面積逐漸增大。在川滇菱形塊體東側(cè)的低速區(qū)域，已穿過其速度分界線龍門山斷裂、安寧河斷裂帶延伸至四川盆地內(nèi)部。根據(jù)曾融生、Blackman、Sliver 等[11-14]提出青藏高原物質(zhì)有橫向流動的可能性，推斷出現(xiàn)這一現(xiàn)象的原因可能為：印度板塊向青藏高原的下部俯沖，使青藏高原的地殼抬升，在巨大的擠壓下，增多的地殼物質(zhì)向青藏高原東部地區(qū)的上地幔流動。其中東流的物質(zhì)在龍門山斷裂帶附近遇到了四川盆地的阻擋，部分改向東南方向流動，但是仍有部分向下侵入到四川盆地的西南區(qū)域，致使區(qū)域 P 波速度較低。并且在深部斷裂帶呈現(xiàn)負的速度異常，更有助于地殼塊體沿斷裂的側(cè)向擠出。

圖3 P 波震相速度擾動水平剖面

從102°E經(jīng)度P波震相的速度擾動垂直剖面（圖4）可以看出，該剖面自南往北依次穿過滇西塊體、滇中塊體、川西塊體、川東南塊體、川青塊體。在0～150 km深度上，23°N～31°N范圍內(nèi)的地殼和上地幔速度存在明顯的不均勻性，低速和高速相間。產(chǎn)生這一現(xiàn)象的原因可能是印度板塊在與歐亞板塊碰撞、俯沖、擠壓過程中，由于地殼應(yīng)力場的非均勻性和鮮水河斷裂帶的存在，在高速的四川盆地邊緣發(fā)生了斜上方的逆沖運動造成的。在深度150km以下，該剖面的速度值逐漸變小，且低速度區(qū)域的面積逐漸增大。這一結(jié)果也在一定程度上驗證了川西高原在下地殼至上地幔的范圍內(nèi)存在較軟物質(zhì)的可能性。

圖4 102°E經(jīng)度P 波震相速度擾動垂直剖面

圖5 29°N緯度P 波震相速度擾動垂直剖面

從29°N緯度P波震相的速度擾動垂直剖面（圖5）可以看出，在0～100 km深度范圍內(nèi)，地殼和上地幔的速度存在橫向的不均勻性。在川滇菱形塊體內(nèi)101°E～102°E范圍的80km深度以上，該處存在一個高速區(qū)域。100km以下，川滇菱形塊體低速度區(qū)域的面積逐漸增大。并且低速體在西部以金沙江斷裂為界，東部以小江斷裂為界，據(jù)此我們推測，在青藏高原受到印度板塊NNE方向的擠壓后，川滇地區(qū)受到連帶影響，由于斷裂帶的巖石張力的破碎，使四川盆地沿著小江斷裂逆向俯沖（圖5箭頭所示）。

綜上所述，我們推測：印度板塊以NNE方向運動與歐亞板塊碰撞，并向青藏高原下部俯沖，致使青藏高原地殼增厚，并向東強烈擠壓，由于地殼應(yīng)力場的非均勻性和斷裂帶的存在，在這種擠壓下，川滇地區(qū)的地殼發(fā)生變形，區(qū)內(nèi)的次級塊體及邊界產(chǎn)生逆沖或者扭轉(zhuǎn)，對川滇地區(qū)的地質(zhì)構(gòu)造造成影響。由研究區(qū)域大面積的低速體的存在，推測該現(xiàn)象可能是由地殼和上地幔存在高導(dǎo)層、高熱流值造成的。

縱觀整個研究區(qū)域，并結(jié)合前人的研究成果，我們得出川滇地區(qū)地殼結(jié)構(gòu)的總體特征是：在上地殼速度異常分布中，存在明顯的橫向不均勻性，但是總體上川滇地區(qū)地殼、上地幔的速度均速度較低；龍門山斷裂帶、安寧河斷裂帶，金沙江斷裂帶在地殼和上地幔一定深度內(nèi)的速度異常中顯示出構(gòu)造分界特征，并且深部斷裂帶呈現(xiàn)負速度異常；地殼厚度變化劇烈，地殼和上地幔存在高導(dǎo)層、高熱流值。

5 存在問題

地震波旅行時層析成像包括數(shù)據(jù)處理，正演模擬和數(shù)據(jù)反演等環(huán)節(jié)，而每個環(huán)節(jié)都有可能導(dǎo)致最終結(jié)果的不確定。該方法在應(yīng)用中主要存在問題具體如下：

（1）反演的數(shù)據(jù)量不足。由于臺站的數(shù)量有限，參與層析成像反演的數(shù)據(jù)量不足，從而導(dǎo)致成像結(jié)果的分辨率不高。

（2）震源定位不準確。由于臺網(wǎng)的稀疏，使得地震的定位存在很大的誤差，尤其是震源的深度，那么對地震層析成像反演結(jié)果有很大影響的的旅行時就會變的不準確

（3）初始參考模型引起的誤差。在進行射線追蹤時，由于不知道真實的射線路徑，我們首先引入了參考模型，然后通過逐步迭代對參數(shù)進行修正。如果成像反演時利用了折射或者反射震相，則初始參考模型中速度界面信息的準確性尤為重要，因為速度間斷面對這些震相的射線傳播路徑影響很大，所以可靠的速度間斷面信息對地震層析成像反演是至關(guān)重要的。

參考文獻

[1] Wang Y X， Jiang M， Nábelek J，et al. The new insight of the velocity structure beneath the Indian-Eurasian continental collision zone. 2006， AGU West. Pac. Geophys. Meet， Suppl.， Abstract S45A～04.

[2] Zhao Dapeng， Akira Hasegawa. P wave tomographic imaging of the crust and upper mantle

beneath the Japan islands [J]. J. Geophys. Res.， 1993， 98（B3）： 4333～4353

[3] 賀世杰，郭峰. 地幔柱的識別和演化綜合述評 [J]. 地球科學進展，2003，18（3）：433～439

[4] 劉建華，劉福田，吳華，等. 中國南北帶地殼和上地幔的三維速度圖像 [J]. 地球物理學報，1989，32（2）：143～151

[5] 王有學，姜枚，熊盛青，等. 西昆侖巖石圈的拆沉作用及其深部構(gòu)造含義—地震層

析成像及航磁異常證據(jù) [J]. 中國地質(zhì)，2006，33（2）：299～308

[6] 姚振興，李白基，梁尚鴻，等. 青藏高原地區(qū)瑞利波群速度和地殼構(gòu)造[J]. 地球物理學報，1981，24（3）：287～29

[7] 丁志峰，何正勤，吳建平，等. 青藏高原地震波三維速度結(jié)構(gòu)的研究[J]. 中國地震， 2001，17（2）：202～209

[8] 錢輝，姜枚， Chen Wangping，等. 青藏高原吉隆—魯谷（Hi-Climb）層析成像與引藏碰撞的消減作用[J]. 地球物理學報，2007，50（5）：1427～1436

[9] 吳慶舉，曹融生，趙文津. 喜馬拉雅-青藏高原的上地幔傾斜構(gòu)造與陸-陸碰撞過程[J]. 中國科學D輯，2004，34（10）：919～925

[10] 賀日政，趙大鵬，高銳，等. 西昆侖造山帶下巖石圈地幔速度結(jié)構(gòu)[J]. 地球物理學報，2006，49（3）：778～787

[11] 曾融生，朱介壽，周兵，等. 青藏高原及其東部鄰區(qū)的三維地震波速度結(jié)構(gòu)與大陸碰撞模型 [J]. 地震學報，1992，11（S）：520～533

[12] 曾融生，丁志峰， 吳慶舉. 青藏高原巖石圈構(gòu)造及動力學過程研究[J]. 地球物理學報， 1994，37（S）：99～116

[13] Blackman D K. Seicmic anisotrophy in the mantel beneath an oceanic spreading center[J]. Nature，1993， 366： 675～677

[14] Sliver P G， Chan W W. Implication for continental structure and evolution from seismic anisotropy [J]. Nature， 1998， 335（1）： 34～39