王海燕

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

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

1 前言

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

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

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

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

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

（1）

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

（2）

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

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

（3）

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

圖1 射線路徑示意圖

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

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

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

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

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

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

3 層析成像的基本原理

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

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

（4）

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

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

（5）

從而有

（6）

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

（7）

寫成矩陣，有

D=CU （8）

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

（9）

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

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

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

根據(jù)眾多學(xué)者對(duì)三江地區(qū)的研究資料 [5-10] 可知，該區(qū)域地殼比較厚，莫霍面的深度達(dá)到了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進(jìn)行劃分。

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

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

圖3 P 波震相速度擾動(dòng)水平剖面

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

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

圖5 29°N緯度P 波震相速度擾動(dòng)垂直剖面

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

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

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

5 存在問題

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

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

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

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

參考文獻(xià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] 賀世杰，郭峰. 地幔柱的識(shí)別和演化綜合述評(píng) [J]. 地球科學(xué)進(jìn)展，2003，18（3）：433～439

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

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

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

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

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

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

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

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

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

[12] 曾融生，丁志峰， 吳慶舉. 青藏高原巖石圈構(gòu)造及動(dòng)力學(xué)過程研究[J]. 地球物理學(xué)報(bào)， 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