王曉， 周小鵬， 張新彥， 白志明， 滕吉文

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所，巖石圈演化國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 中國(guó)科學(xué)院大學(xué)， 北京 100049

?

上地殼縱橫波速度結(jié)構(gòu)相關(guān)反演成像方法

王曉1,2， 周小鵬1,2， 張新彥1,2， 白志明1*， 滕吉文1

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所，巖石圈演化國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 中國(guó)科學(xué)院大學(xué)， 北京 100049

基于縱橫波初至走時(shí)數(shù)據(jù)的層析成像方法越來(lái)越廣泛地被應(yīng)用于揭示不同構(gòu)造域殼幔速度結(jié)構(gòu)特征.我們從同一地質(zhì)體的縱橫波速度屬性相關(guān)這一基本思想出發(fā)，提出一種相關(guān)反演成像的方法：縱橫波速度反演交替進(jìn)行，在迭代反演過(guò)程中每通過(guò)一次反演獲得相應(yīng)的縱波速度(或橫波速度)結(jié)構(gòu)后，更新相應(yīng)的縱橫波速度比模型以及相應(yīng)的橫波(或縱波)速度反演的初始模型，然后繼續(xù)開展后續(xù)橫波(或縱波)速度反演工作.在反演過(guò)程中依據(jù)縱橫波速度的相關(guān)性信息和射線路徑長(zhǎng)度將走時(shí)殘差以不同權(quán)重分配到射線路徑經(jīng)過(guò)的單元，依據(jù)網(wǎng)格節(jié)點(diǎn)周圍平均的慢度擾動(dòng)更新速度模型.正反演過(guò)程分別基于有限差分走時(shí)計(jì)算方法和反投影成像方法.兩種典型模型試驗(yàn)表明，該技術(shù)應(yīng)用于上地殼速度結(jié)構(gòu)反演成像過(guò)程，可有效提高反演結(jié)果的可靠性，在很大程度上避免了常規(guī)單獨(dú)反演縱波和橫波速度過(guò)程容易帶來(lái)的畸變和失真.該方法應(yīng)用于重建青藏高原西部札達(dá)—泉水溝深地震測(cè)深(DSS)剖面下方的上地殼速度結(jié)構(gòu)，揭示出與青藏高原西緣板塊碰撞相關(guān)的上地殼速度結(jié)構(gòu)特征.

上部地殼； 縱橫波速度結(jié)構(gòu)； 相關(guān)反演； 有限差分； 反投影

1 引言

自20世紀(jì)80年代以來(lái)，基于深地震測(cè)深縱橫波初至走時(shí)數(shù)據(jù)的層析成像結(jié)果越來(lái)越廣泛地被應(yīng)用于揭示盆嶺、裂谷、斷陷等不同構(gòu)造的深部形態(tài)及上地殼內(nèi)高速巖體、巖漿囊和斷層等結(jié)構(gòu)的精細(xì)特征(Tarantola and Valette, 1982; Thurber, 1983; Vidale, 1988; Bregman et al., 1989; Saito, 1989; Hole and Zelt, 1995; 張先康，1996；Zhang and Toksozf, 1998; 丁志峰等，1999; Xu et al., 2006,2010,2014; Teng et al., 2013).Vidale(1988)最先提出了二維介質(zhì)的有限差分正演模擬算法，Hole (1992)將其改進(jìn)，推廣應(yīng)用于速度橫向變化劇烈的三維介質(zhì)，并提出了一種利用初至波走時(shí)數(shù)據(jù)Pg和/或Sg開展層析反演的方法.由于有限差分方法的通用性，該方法自提出以來(lái)不斷得到應(yīng)用和改進(jìn).經(jīng)過(guò)Ammon和Vidale (1993)改進(jìn)后，該技術(shù)已被廣泛地應(yīng)用于利用深地震測(cè)深(DSS)數(shù)據(jù)重建上地殼速度結(jié)構(gòu)領(lǐng)域.王椿鏞等(1997)應(yīng)用該方法得到大別造山帶人工地震測(cè)深剖面上部地殼橫向不均勻結(jié)構(gòu)圖像，為大別造山帶演化過(guò)程的研究提供了重要的深部地球物理證據(jù)；段永紅等(1999)詳細(xì)論述了二、三維層析成像的基本原理，并以唐山震區(qū)的實(shí)際反演結(jié)果驗(yàn)證了方法的有效性.之后，越來(lái)越多的學(xué)者將該方法用于研究地殼結(jié)構(gòu)和構(gòu)造特征，并取得了一系列的成果(Zelt, 1999; Bai et al., 2007; 王夫運(yùn)等，2008；Thybo and Nielson, 2009; Zhang et al., 2009b,2011;Karplus et al., 2011).但是，該方法仍然存在一定的局限性.

一般情況下由于受到地震資料信噪比限制，有效的S波到時(shí)數(shù)據(jù)相比P波數(shù)據(jù)數(shù)據(jù)量少、質(zhì)量差，單獨(dú)反演所得橫波速度結(jié)構(gòu)與單獨(dú)反演所得縱波速度結(jié)構(gòu)經(jīng)常存在顯著區(qū)別，由此計(jì)算得到的VP/VS圖像往往會(huì)產(chǎn)生畸變現(xiàn)象，這種不合理性就不可避免地給地質(zhì)解釋帶來(lái)巨大挑戰(zhàn)(Eberhart-Phillips, 1990).

針對(duì)以上缺陷，Thurber和Atre (1993)首先提出了一種基于P波到時(shí)數(shù)據(jù)及S-P時(shí)間差反演VP和VS的算法，這種算法假定P波和S波具有基本相同的射線路徑.顯然這一假定并不適合復(fù)雜三維速度結(jié)構(gòu)：此時(shí)由于P波射線路徑和S波射線路徑差異較大，基于反演得到的P波速度模型和S波初始速度模型并不能可靠地反演獲得S波速度結(jié)構(gòu)(Wagner et al., 2005)，就像由S波速度模型和P波初始速度模型并不能可靠地反演獲得P波速度結(jié)構(gòu)一樣(Eberhart-Phillips, 1990).

為此，Conder和Wiens (2006)提出采用奇異值分解法(SVD)聯(lián)合反演VP和VP/VS，該方法不受P波路徑和S波路徑相同這一基本假設(shè)的限制.其射線追蹤過(guò)程考慮了地球的曲率，采用梯形網(wǎng)格，用非線性方程表征S波走時(shí)和P波慢度、VP/VS比值之間的關(guān)系，并且將平滑約束加入偏導(dǎo)數(shù)矩陣中，幾乎利用了所有的SVD解，提高了反演的穩(wěn)定性，同時(shí)有效減小了速度模型的方差.雖然該方法給出了分辨率矩陣的定量描述，但是并不能給出相對(duì)分辨率.Zhang等(2009b)提出了另一種改進(jìn)方法：一是在每次迭代后檢查P波和S波射線路徑，在后續(xù)反演計(jì)算VP/VS時(shí)排除射線路徑差異較大的走時(shí)數(shù)據(jù)；二是使用了具有共同臺(tái)站事件對(duì)的差分到時(shí)數(shù)據(jù)以確定震源區(qū)VP/VS模型.該算法避免了一些假象的產(chǎn)生，例如直接依據(jù)VP和VS得到地殼VP/VS模型等，但在實(shí)踐中仍然存在一些不足，如需要人為確定射線路徑差異的門檻值，以及某些情況下反演得到的VP與VS圖像差異較大，地質(zhì)解釋存在困難等.為了有效解決P波路徑與S波路徑差異問(wèn)題，Wang (2014)采用同一炮檢對(duì)的P波和S波走時(shí)殘差，直接將VP/VS異常作為反演的模型參數(shù)，計(jì)算P波射線路徑上每一小段的相對(duì)P波慢度擾動(dòng)q，將q用于迭代反演沿S波射線路徑上每個(gè)格點(diǎn)的相對(duì)S波慢度擾動(dòng)r，直到每個(gè)網(wǎng)格點(diǎn)的r達(dá)到誤差允許范圍內(nèi).相比以往常規(guī)反演中假設(shè)P波、S波射線路徑相同而言，該方法反演得到的VP/VS比值誤差更小.但是，該方法的VS模型由S波數(shù)據(jù)單獨(dú)反演得到，其結(jié)果在很大程度上受到S波數(shù)據(jù)質(zhì)量、反演參數(shù)及迭代次數(shù)等限制.

基于此，本文提出了一種綜合利用深地震測(cè)深Pg和Sg走時(shí)數(shù)據(jù)反演上地殼縱橫波速度結(jié)構(gòu)的相關(guān)反演成像方法(圖4).與以上方法相比，該方法從假定Pg和Sg數(shù)據(jù)的走時(shí)殘差來(lái)自地下同一地質(zhì)異常體、且該地質(zhì)體的縱橫波速度異常圖像相關(guān)的這一基本思想出發(fā)，在交叉迭代反演縱橫波速度過(guò)程中使用縱橫波速度異常的相關(guān)性作為約束，將走時(shí)殘差更準(zhǔn)確可靠地反投影到目標(biāo)地質(zhì)異常體，從而有效保證了成像結(jié)果可靠性.

在資料處理和解釋過(guò)程中，如果測(cè)線彎曲程度不大，一般按照直線方式來(lái)處理；如果測(cè)線彎曲程度較大，最好將測(cè)線進(jìn)行分段直線劃分后開展后續(xù)工作.本文提出的相關(guān)反演方法是針對(duì)三維介質(zhì)進(jìn)行反演的，在反演開始的數(shù)據(jù)輸入過(guò)程中包含了接收點(diǎn)的x坐標(biāo)(包括測(cè)點(diǎn)相對(duì)原點(diǎn)的坐標(biāo)和測(cè)點(diǎn)相對(duì)炮點(diǎn)的坐標(biāo))、y坐標(biāo)、z坐標(biāo)(高程)等有效信息，有效解決了檢波器的高程差異問(wèn)題.幾種典型模型的對(duì)比試驗(yàn)驗(yàn)證了該方法的有效性和優(yōu)勢(shì)，作為示例本文將該方法應(yīng)用于青藏高原西部真實(shí)的DSS數(shù)據(jù)，揭示了上地殼內(nèi)印度-歐亞板塊陸陸碰撞的證據(jù).

2 上地殼縱橫波速度結(jié)構(gòu)相關(guān)反演成像方法

2.1 構(gòu)建一維速度模型

利用初至走時(shí)數(shù)據(jù)進(jìn)行有限差分反演成像的方法具有震相拾取簡(jiǎn)單、成像速度快、容許對(duì)模型進(jìn)行精細(xì)剖分、適應(yīng)性強(qiáng)等優(yōu)點(diǎn)(蘭海強(qiáng)等，2012a,2012b；趙烽帆等，2014)，因此近二十年以來(lái)該方法在重建上部地殼結(jié)構(gòu)，研究構(gòu)造精細(xì)特征方面得到了廣泛應(yīng)用(Vidale, 1990; 王椿鏞等，1997，2002；段永紅等，1999；Zhao et al., 2004; 趙愛華和丁志峰，2005；Zhang et al., 2009b；劉一峰和蘭海強(qiáng)，2012；Teng et al., 2013).Hole (1992)提出了利用初至走時(shí)數(shù)據(jù)反演地殼結(jié)構(gòu)的反投影重建方法.該方法基本思想是，以小的矩形網(wǎng)格單元?jiǎng)澐直碚髂Ｐ涂臻g,近似認(rèn)為地震波在網(wǎng)格單元上以平面波形式傳播.在此基礎(chǔ)上通過(guò)程函方程的有限差分算子計(jì)算出地震波從震源出發(fā)到所有網(wǎng)格點(diǎn)的初至走時(shí)場(chǎng)，根據(jù)計(jì)算得到的走時(shí)場(chǎng)沿最陡下降的方向由檢波器到炮點(diǎn)快速追蹤獲得射線路徑.反演過(guò)程采用反投影算法，即將走時(shí)殘差平均分配到射線路徑上，獲得各網(wǎng)格單元內(nèi)的慢度擾動(dòng)值，進(jìn)而獲得更新的速度模型.

反演過(guò)程中對(duì)走時(shí)擾動(dòng)與慢度擾動(dòng)的關(guān)系作線性化處理，即

(1)

通過(guò)迭代方法以簡(jiǎn)單反投影近似求解該方程，最終給出慢度擾動(dòng)模型：

(2)

其中，δtj和lj分別是第j條射線的走時(shí)殘差和總長(zhǎng)度，gj(r)是數(shù)據(jù)核，A是以射線為中心的橫截面積.假設(shè)射線路徑穿過(guò)K個(gè)節(jié)點(diǎn)，則每個(gè)節(jié)點(diǎn)上的慢度擾動(dòng)為：

(3)

其中，δtk和lk分別是第k個(gè)網(wǎng)格節(jié)點(diǎn)的走時(shí)殘差和射線長(zhǎng)度.

這種反演方法的計(jì)算效率允許對(duì)模型采取密集采樣，提供了空間分辨率較好的三維層析圖像.鑒于該方法在反演過(guò)程使用了擬線性最小二乘反演，因此選擇恰當(dāng)?shù)某跏寄Ｐ褪种匾?Kissling等(1994,1995)在解決天然源地震的精定位和一維速度模型問(wèn)題時(shí)提出確定最佳一維速度模型的方法, 最小一維速度模型方法計(jì)算的速度模型可使定位結(jié)果的走時(shí)殘差均方根最小, 定位精度更高.理論上最佳一維模型應(yīng)逼近真實(shí)的地球模型，然而在寬角反射/折射地震資料處理解釋領(lǐng)域，作為利用Pg或Sg走時(shí)數(shù)據(jù)反演上地殼速度結(jié)構(gòu)全過(guò)程的第一步，如何構(gòu)建上地殼速度結(jié)構(gòu)的最佳一維模型長(zhǎng)期以來(lái)國(guó)內(nèi)外文獻(xiàn)一直鮮有論述且未見得到應(yīng)用.一種變通的方法是嘗試比較若干不同一維初始模型，及不同參數(shù)進(jìn)行反演，然后比較不同反演結(jié)果對(duì)應(yīng)的射線覆蓋和走時(shí)擬合效果(徐濤等，2004；李飛等，2013)，從中選擇較為合理的速度模型作為真實(shí)地球的近似模型，這在某種程度上的確有助于逼近真實(shí)地球模型，然而初始模型選擇的主觀隨意性大大減少了結(jié)果的可信度.

為得到可靠的一維模型，基于Hole(1992)，Ammon和Vidale (1993)的反演代碼，我們?cè)诔跏挤囱葸^(guò)程中將反演參數(shù)的橫向網(wǎng)格單元和節(jié)點(diǎn)數(shù)設(shè)定為最大，覆蓋全剖面長(zhǎng)度，橫向光滑長(zhǎng)度也設(shè)定為剖面整個(gè)長(zhǎng)度，從而通過(guò)10次迭代反演巧妙地得到一維速度模型.為驗(yàn)證該方法的有效性，我們構(gòu)建了兩個(gè)理論速度模型，并計(jì)算相應(yīng)的初至走時(shí)數(shù)據(jù)(圖1，2)，目的是為了驗(yàn)證針對(duì)不同的異常體(簡(jiǎn)單規(guī)則異常和復(fù)雜不規(guī)則異常)，我們的反演方法是否均行之有效.剖面長(zhǎng)度均為380 km，地表有9炮和380個(gè)檢波器，間隔為1 km.需要注意的是，每張圖中的(a)和(b)分別表示擾動(dòng)的P波和S波速度模型，(c)和(d)是它們的理論折合走時(shí)曲線.這兩種速度模型具有相同的初始一維P波和S波背景(圖1a, 1b的插圖).與此同時(shí)，我們簡(jiǎn)單地假設(shè)兩個(gè)模型中VP/VS比值為恒定值1.80.

圖3a和3b分別表示利用10次迭代后得到的VP、VS模型最終反演得到的一維速度-深度曲線.模型1(圖1a，1b)的灰色曲線和模型2(圖2a，2b)的虛線基本上接近原始理論曲線(黑色曲線)，這表明即使在初始一維模型與理論一維模型差別較大的情況下，反演得到的一維速度模型仍然能很好地逼近理論模型，因此我們的方法是穩(wěn)定有效的，它可以廣泛地用于合理構(gòu)建與Pg和Sg數(shù)據(jù)對(duì)應(yīng)的一維縱橫波速度模型，從而穩(wěn)定有效地開展后續(xù)二維反演工作.

2.2 相關(guān)反演成像方法

盡管獲得了可靠的一維速度模型，但是通過(guò)反演構(gòu)建可靠的二維VP、VS和VP/VS模型仍然是一個(gè)很大的挑戰(zhàn).圖5和圖6的右圖分別顯示了利用Ammon和Vidale(1993)的反演代碼最終得到的模型1和模型2的VP、VS和VP/VS誤差.這兩個(gè)圖表明盡管重建的VP與VS模型顯示了與初始模型相似的特性，但是VS模型已經(jīng)在很大程度上發(fā)生了扭曲，最終得到的VP/VS模型已經(jīng)嚴(yán)重偏離了給定的初始值1.80.這些結(jié)果清楚地表明，分別單獨(dú)反演VP和VS不能可靠地推斷出最終的VP/VS模型，類似地，不能簡(jiǎn)單地依據(jù)VP和VP/VS的值推斷VS(Eberhart-Phillips,1990；Wagner et al., 2005；Conder and Wiens,2006; Wang, 2014).

為構(gòu)建可靠的二維模型，嘗試在反演過(guò)程中添加緊約束：縱橫波速度反演交叉進(jìn)行，且以每次迭代后生成的VP或VS模型以及VP/VS比值作為下一步反演的VS或VP初始速度模型.反演方法見圖4，“相關(guān)反演”的概念描述如下：

假定Pg和Sg數(shù)據(jù)的走時(shí)殘差完全來(lái)自地下同一地質(zhì)異常體，且該地質(zhì)體的縱橫波速度相對(duì)于圍巖同步地表現(xiàn)出速度高或速度低的特點(diǎn)，則該地質(zhì)體附近的VP和VS屬性的圖像特征存在明顯的正相關(guān)或負(fù)相關(guān).據(jù)此每次迭代后計(jì)算VP和VS模型的相關(guān)系數(shù)，分析每個(gè)網(wǎng)格點(diǎn)存在真實(shí)地質(zhì)體的概率.三維空間圍繞一個(gè)特定格點(diǎn)共有7個(gè)相鄰網(wǎng)格點(diǎn)(包括其本身)，因而每次迭代后網(wǎng)格點(diǎn)的系數(shù)被定義為：

圖1 模型1：具有矩形異常的上地殼速度模型((a) P波，(b)S波)以及對(duì)應(yīng)的理論折合走時(shí)曲線((c)，(d))Fig.1 Model 1: upper crust velocity models with rectangular anomaly ((a) for P-wave and (b) for S-wave),and their corresponding theoretical reduced traveltime curve ((c) and (d))

圖2 模型2：正弦速度模型((a) P波，(b)S波)以及對(duì)應(yīng)的理論折合走時(shí)曲線((c)，(d))加入的異常延伸至上地殼底部，該模型與圖1具有相同的一維速度背景.插圖：(a)和(b)右下角的說(shuō)明分別表示初始一維VP或VS模型.Fig.2 Model 2: sine-type velocity models ((a) for P-wave and (b) for S-wave) and their corresponding theoretical reduced traveltime curve ((c) and (d)) The added anomalies extend downward to the bottom of the upper crust, and this model has the same 1-D velocity background as that of Fig.1.Inset: illustrations in (a) and (b) at the bottom right corner respectively denotes the initial 1-D VP or VS model.

圖3 反演獲得的P波和S波一維速度曲線黑色實(shí)線：初始理論一維深度-速度曲線；點(diǎn)劃線：用于反演得到更加可靠的一維曲線的初始一維曲線，其適合于接下來(lái)的二維反演；灰色實(shí)線：模型1 10次迭代后的反演結(jié)果；虛線：模型2 10次迭代后的反演結(jié)果.Fig.3 Derived 1-D velocity curve of P-wave (a) and S-wave (b) using inversion methodBlack solid line: original theoretical 1-D depth-velocity curve; Dash dot line: initial 1-D curve for inversion to get a more reliable 1-D curve which is appropriate for following 2-D inversion; Gray solid line: inversion results after 10 times iteration for model 1; Dotted line: inversion results after 10 times iteration for model 2.

×(|ai|+|bi|)0.125,

(4)

其中，

ai=(VPnew-VPori)/VPori，

bi=(VSnew-VSori)/VSori.

事實(shí)上，一般固定C0為0.00001，VPnew和VSnew表示每次迭代后所考慮點(diǎn)的新速度值，VPori和VSori是初始一維速度模型的速度值.(4)式表示圍繞特定格點(diǎn)平均的相關(guān)系數(shù).用這個(gè)公式，經(jīng)過(guò)一次迭代后，如果鄰近點(diǎn)的VP和VS沒有發(fā)生變化，即ai=bi=0.0，系數(shù)將被保持在最小值C0；但是比如ai=bi=0.2時(shí)，我們將得到一個(gè)較大的系數(shù)，其值等于0.8917+C0.這意味著，如果每一步反演后圍繞一個(gè)網(wǎng)格點(diǎn)存在相似的速度變化，那么計(jì)算出的相關(guān)系數(shù)將是一個(gè)較大的值.另一方面可以使用上述相關(guān)系數(shù)來(lái)描述現(xiàn)有的地質(zhì)異常，計(jì)算的相關(guān)系數(shù)越大，對(duì)應(yīng)地質(zhì)體存在的概率越高.基于網(wǎng)格點(diǎn)的系數(shù)Cor1，可以很容易地通過(guò)線性插值計(jì)算出一個(gè)單元格的相關(guān)系數(shù)Cor2.然后利用這樣一種“相關(guān)反演”的想法對(duì)Hole (1992)的反演過(guò)程進(jìn)行了如下修改：

假設(shè)單元格的慢度擾動(dòng)與計(jì)算得到的相關(guān)系數(shù)Cor2成比例，即模型空間某一格點(diǎn)周圍縱橫波速度圖像特征的相關(guān)性越大，相應(yīng)地質(zhì)異常存在的可能性越大，該潛在異常體越可能導(dǎo)致較大的慢度擾動(dòng).上述可以被描述為：

(5)

(5)式假定射線路徑穿過(guò)n個(gè)單元格，單元格k的長(zhǎng)度為lk，Δtj是第j條射線的走時(shí)殘差，Δui和Cori分別是第i個(gè)單元格的慢度擾動(dòng)和相關(guān)系數(shù).(5)式表示針對(duì)第j條射線的走時(shí)殘差Δtj，根據(jù)第i個(gè)單元格的相關(guān)系數(shù)Cori，更新第i個(gè)單元格的慢度擾動(dòng)Δui.如此將相關(guān)系數(shù)成功引入到慢度擾動(dòng)計(jì)算中，進(jìn)而更新速度模型.如果相關(guān)系數(shù)Cori恒等于1.0，該種情況下，(5)式可簡(jiǎn)化為Hole (1992)的常規(guī)表達(dá)式，所計(jì)算的慢度擾動(dòng)變?yōu)椋?/p>

(6)

其中L是射線總長(zhǎng)度.

由以上分析可知，Hole(1992)有限差分反演的基本思想是，射線路徑上每個(gè)網(wǎng)格單元均可能引起走時(shí)殘差，且它們的幾率是相同的.VP與VS模型單獨(dú)反演，彼此之間沒有任何約束.與此相比，我們?cè)诜囱葸^(guò)程中通過(guò)引入P波和S波的相關(guān)系數(shù)Cor，將P波速度與S波速度緊密結(jié)合起來(lái)：利用每次迭代后的VP反演結(jié)果來(lái)約束VS反演過(guò)程，類似地，用VS的反演結(jié)果來(lái)約束VP反演，然后利用VP和VS模型之間的相關(guān)性來(lái)指導(dǎo)反演過(guò)程，從而得到更好的反演結(jié)果.

相應(yīng)地我們改進(jìn)了Hole(1992)，Ammon和Vidale (1993)的反演代碼，得到改進(jìn)的慢度擾動(dòng).射線追蹤和慢度擾動(dòng)計(jì)算后，取每個(gè)網(wǎng)格點(diǎn)慢度平均值，速度模型和VP/VS比值模型將被更新.圖5和圖6的右圖分別表示使用該方法最終反演得到模型1和2的VP、VS和VP/VS比值誤差的模型.單元格或網(wǎng)格點(diǎn)的反演參數(shù)重采樣尺寸和平滑尺寸保持最小.這兩張圖表明，與依據(jù)Hole (1992)等算法反演得到的結(jié)果相比，由我們改進(jìn)方法反演得到的VP和VS模型有明顯改善，顯然更接近初始理論模型(圖1和圖2)，而依據(jù)Hole (1992)算法分別反演縱橫波速度結(jié)構(gòu)，進(jìn)而計(jì)算得到的速度比模型卻嚴(yán)重失真；依據(jù)我們方法反演得到的VP/VS比值誤差最大值僅為0.02～0.03，相比之下依據(jù)傳統(tǒng)方法算得的VP/VS比值誤差最大值高達(dá)0.2.以上試驗(yàn)表明，采用不同反演方法得到的最終模型可能具有相同水準(zhǔn)的走時(shí)擬合效果，但這并不能保證反演結(jié)果的可靠性.即使在速度模型橫向變化較大情形下，我們的二維反演方法對(duì)于不同類型的合成走時(shí)數(shù)據(jù)集也取得滿意效果(圖7).

圖4 利用Pg和Sg數(shù)據(jù)獲取上地殼結(jié)構(gòu)的改進(jìn)二維反演流程圖Fig.4 Flow-chart of improved 2-D inversion procedure for upper crust structure using Pg and Sg data

圖5 利用合成Pg和Sg走時(shí)數(shù)據(jù)反演獲得的速度結(jié)構(gòu)(模型1)(a) 利用Ammon和Vidale(1993)反演代碼得到的上地殼VP結(jié)構(gòu)； (b) 利用我們的改進(jìn)方法得到的上地殼VP結(jié)構(gòu)； (c) 利用Ammon和Vidale(1993)反演代碼得到的上地殼VS結(jié)構(gòu)； (d) 利用我們的改進(jìn)方法得到的上地殼VS結(jié)構(gòu)； (e) 利用Ammon和Vidale(1993)反演代碼得到的VP/VS誤差； (f) 利用我們的改進(jìn)方法得到的VP/VS誤差.所有反演結(jié)果均經(jīng)過(guò)10次迭代獲得.Fig.5 Final 2-D inversion results for Model 1 using synthetic Pg and Sg traveltime data(a) Upper crust VP structure derived by the inversion code of Ammon and Vidale (1993); (b) Upper crust VP structure derived by our improved approach; (c) Upper crust VS structure derived by the inversion code of Ammon and Vidale (1993); (d) Upper crust VS structure derived by our improved approach; (e) VP/VS error derived by the inversion code of Ammon and Vidale (1993); (f) VP/VS error derived by our improved approach. All the inversion results were acquired after 10 times iteration.

圖6 利用合成Pg和Sg走時(shí)數(shù)據(jù)反演獲得的速度結(jié)構(gòu)(模型2)(a) 利用Ammon和Vidale(1993)反演代碼得到的上地殼VP結(jié)構(gòu)； (b) 利用我們的改進(jìn)方法得到的上地殼VP結(jié)構(gòu)； (c) 利用Ammon和Vidale(1993)反演代碼得到的上地殼VS結(jié)構(gòu)； (d) 利用我們的改進(jìn)方法得到的上地殼VS結(jié)構(gòu)； (e) 利用Ammon和Vidale(1993)反演代碼得到的VP/VS誤差； (f) 利用我們的改進(jìn)方法得到的VP/VS誤差.所有反演結(jié)果均經(jīng)過(guò)10次迭代獲得.Fig.6 Final 2-D inversion results for Model 2 using synthetic Pg and Sg traveltime data(a) Upper crust VP structure derived by the inversion code of Ammon and Vidale (1993); (b) Upper crust VP structure derived by our improved approach; (c) Upper crust VS structure derived by the inversion code of Ammon and Vidale (1993); (d) Upper crust VS structure derived by our improved approach; (e) VP/VS error derived by the inversion code of Ammon and Vidale (1993); (f) VP/VS error derived by our improved approach. All the inversion results were acquired after 10 times iteration.

圖7 二維反演結(jié)果走時(shí)殘差均方根(a) 對(duì)應(yīng)于圖1和圖2中的VP模型; (b) 對(duì)應(yīng)于圖1和圖2中的VS模型.Fig.7 Root Mean Square (RMS) of time residual of the 2-D inversion results(a) For the VP models in Figs.1 and 2; (b) For the VS models in Figs.1 and 2.

2.3 模型的定量評(píng)估

(7)

對(duì)模型空間每一格點(diǎn)掃描計(jì)算速度擾動(dòng)引起的目標(biāo)函數(shù)(走時(shí)殘差均方根)變化，得到的殘差均方根變化值ΔF定量描述了模型空間某點(diǎn)速度變化的重要性，或觀測(cè)系統(tǒng)對(duì)模型空間某一點(diǎn)速度變化的分辨能力：圖8和圖9同時(shí)也分別展示了模型1和2的射線覆蓋(圖8a和8b，圖9a和9b)和以上擾動(dòng)函數(shù)ΔF值.對(duì)比二圖可見，速度的可靠性在很大程度上取決于射線覆蓋，然而在某些情況下，密集射線覆蓋并不總意味著反演得到的速度變化特征的可靠性，一般情形下采用擾動(dòng)函數(shù)ΔF的方法可以更簡(jiǎn)單直觀地描述模型的可靠性：某格點(diǎn)附近擾動(dòng)函數(shù)ΔF值越大，該點(diǎn)附近反演得到的速度圖像特征越可靠(圖8c和8d，圖9c和9d).

盡管如此，即使采用上述評(píng)估方法，我們?nèi)匀徊恢朗欠衲軌蛘_揭示潛在的地質(zhì)體.幸運(yùn)的是，由于我們采用了相關(guān)反演的辦法，反演所得縱橫波速度圖像特征的相關(guān)性在很大程度上指示了目標(biāo)地質(zhì)體的存在與否：圖8和9清楚顯示，采用相關(guān)反演所得縱橫波速度結(jié)構(gòu)的相關(guān)性圖像能夠準(zhǔn)確揭示出速度異常體在剖面中的位置和大致形態(tài)，這是其他評(píng)估辦法無(wú)法做到的，由此也充分體現(xiàn)了我們采用反演方法的優(yōu)越性.因此，計(jì)算VP和VS圖像特征相關(guān)性的方法可以幫助我們揭示VP與VS模型可能的同步速度異常，從而進(jìn)一步驗(yàn)證反演結(jié)果的可靠性.

3 札達(dá)—泉水溝深地震測(cè)深剖面上地殼速度結(jié)構(gòu)

為了驗(yàn)證我們的反演成像方法及解評(píng)價(jià)方法的適用性和優(yōu)越性，我們嘗試基于以上方法利用實(shí)際深地震測(cè)深(DSS)數(shù)據(jù)反演青藏高原西部札達(dá)—泉水溝剖面下方的上地殼速度結(jié)構(gòu).

印度板塊向歐亞大陸的俯沖碰撞對(duì)區(qū)域地形和環(huán)境帶來(lái)了顯著影響，如一千多公里的地殼縮短，西藏地殼的東向擠壓以及中亞和喜馬拉雅周邊許多著名的活動(dòng)高山.青藏高原西緣濃縮了青藏高原古生代和新生代的造山歷史，作為研究陸陸碰撞的理想窗口正日益成為國(guó)內(nèi)外地學(xué)領(lǐng)域關(guān)注的焦點(diǎn)(Murphy et al., 2000; Lacassin et al., 2004; Wittlinger et al., 2004; Valli et al., 2007; 李海兵等，2007,2008; Caldwell et al., 2009;許志琴等，2011).研究區(qū)位于西構(gòu)造結(jié)東側(cè)，處于喜馬拉雅、拉薩、羌塘、塔里木、拉達(dá)克、甜水海及西昆侖等不同地體疊合銜接部位，且發(fā)育巨大的喀喇昆侖走滑斷層.為探討研究區(qū)新生代以來(lái)的深部造山過(guò)程及其對(duì)環(huán)境的影響，中國(guó)科學(xué)院地質(zhì)與地球物理研究所從2011年9月至11月在該區(qū)域開展了380 km長(zhǎng)的深地震測(cè)深(DSS)項(xiàng)目(圖10).南北走向的深地震測(cè)深(DSS)剖面始于札達(dá)盆地，然后穿越阿伊拉日居山、喀喇昆侖斷裂(KF)、拉薩地塊、班公湖—怒江縫合帶(BNS)、羌塘地塊、郭扎—龍木錯(cuò)斷裂，并最終延伸到北端的甜水海地塊.本文基于前文所述方法利用獲得的DSS數(shù)據(jù)揭示青藏高原西部與俯沖碰撞在剖面上地殼留下的構(gòu)造遺跡.

圖8 模型1二維反演結(jié)果評(píng)估(a) 所得VP模型的射線覆蓋； (b) 所得VS模型的射線覆蓋； (c) VP模型每一格點(diǎn)速度擾動(dòng)引起的目標(biāo)函數(shù)(走時(shí)殘差均方根)變化值； (d) VS模型每一格點(diǎn)速度擾動(dòng)引起的目標(biāo)函數(shù)(走時(shí)殘差均方根)變化值； (e) 計(jì)算得到的VP與VS模型之間的相關(guān)系數(shù).Fig.8 Assessment of the 2-D inversion results as for Model 1 using our tomographic system(a) Ray-coverage of derived VP model; (b) Ray-coverage of derived VS model; (c) RMS variation on VP model that corresponds to the velocity perturbation at each point; (d) RMS variation on VS model that corresponds to the velocity perturbation at each point; (e) Computed correlation coefficient between VP and VS model.

圖9 模型2二維反演結(jié)果的評(píng)估(a) 所得VP模型的射線覆蓋； (b) 所得VS模型的射線覆蓋； (c) VP模型每一格點(diǎn)速度擾動(dòng)引起的目標(biāo)函數(shù)(走時(shí)殘差均方根)變化； (d) VS模型每一格點(diǎn)速度擾動(dòng)引起的目標(biāo)函數(shù)(走時(shí)殘差均方根)變化； (e) 計(jì)算得到的VP與VS模型之間的相關(guān)系數(shù).Fig.9 Assessment of the 2-D inversion results as for Model 2 using our tomographic system(a) Ray-coverage of derived VP model; (b) Ray-coverage of derived VS model; (c) RMS variation on VP model that corresponds to the velocity perturbation at each point; (d) RMS variation on VS model that corresponds to the velocity perturbation at each point; (e) Computed correlation coefficient between VP and VS model.

圖10 青藏高原西部札達(dá)—泉水溝深地震測(cè)深剖面及其構(gòu)造背景主要斷裂：MBT，主邊界沖斷裂；MCT，主中央沖斷裂；STDS，藏南拆離系；IYS，印度—雅魯藏布江縫合帶；KKF，喀喇昆侖斷裂；SF，獅泉河斷裂；BNS，班公—怒江縫合帶；DT，多瑪沖斷裂；LCF，龍木錯(cuò)斷裂. 紅色五角星表示炮點(diǎn)，藍(lán)色三角形表示接收點(diǎn).Fig.10 Tectonic setting and the deep seismic sounding profile in Western Tibetan Plateau which was carried out in September and October, 2011Main Faults: MBT, Main Boundary Thrust; MCT, Main Central/Panjal Thrust; STDS: Southern Tibetan Detachment System; IYS, Indus-Yarlung Zangbo River Suture Belt; KKF, Karakorum Fault; SF, Shiquanhe Fault; BNS, Bangong-Nujiang Suture Belt; DT: Domar Thrust; LCF, Longmu Co Fault. The red star denotes the shot point, and the blue triangle shows the receiver.

該DSS測(cè)線海拔高程在4.0～5.8 km，大部分檢波器被部署在海拔4.6～5.1 km的高度.在不同構(gòu)造單元(圖10)總共觸發(fā)了七炮，工作區(qū)每一炮點(diǎn)深度～40 m，TNT炸藥量3000 kg.共有200臺(tái)便攜式三分量數(shù)字地震儀沿剖面以2.0 km左右的間距接收地震數(shù)據(jù).這些地震儀記錄的每炮信號(hào)長(zhǎng)度為5 min，采樣間隔5 ms.對(duì)P波和S波進(jìn)行帶通濾波保留1～10 Hz的頻率成分.圖11和圖12分別展示了折合速度為6.0 km·s-1和3.5 km·s-1的P波和S波的所有折合地震記錄.可以利用走時(shí)互換原理從這些記錄中拾取出清晰的P波和S波初至(Pg和Sg)，一般來(lái)說(shuō)Pg的信噪比(SNR)要比Sg的高得多.

Pg震相一般可以在偏移距0～100 km觀測(cè)到，某些炮可觀測(cè)到最大偏移距達(dá)120 km的地方，如第一炮(SP1)的北支和第五炮(SP5)的南支(圖11a，11e).在SP2的南支(圖11b)偏移距-50～30 km，SP3的北支(圖11c)10～50 km，SP6的北支(圖11f)5～40 km，Pg走時(shí)發(fā)生明顯延遲.這些特征也可以由拾取的Sg數(shù)據(jù)(圖12)觀察到.S波記錄顯示在以下偏移距范圍內(nèi)Sg的到時(shí)數(shù)據(jù)也發(fā)生了明顯延遲：SP2的南支(圖12b)-30～5 km，SP3(圖12c)的-90～-60 km和10～50 km，SP5的南支(圖12e)-80～-60 km，SP6的北支(圖12f)5～40 km.這些延遲的Pg和Sg走時(shí)數(shù)據(jù)彼此相關(guān)性明顯，炮點(diǎn)SP3和SP6等記錄剖面上縱橫波初至走時(shí)同步滯后表明這兩炮點(diǎn)附近存在對(duì)應(yīng)的沉積盆地.

從記錄剖面上共拾取了398個(gè)Pg和382個(gè)Sg走時(shí)數(shù)據(jù)用于反演.正演過(guò)程將模型劃分為1 km×1 km的矩形單元.首先，利用前文改進(jìn)算法分別反演縱橫波一維速度結(jié)構(gòu)，10次迭代后得到用于二維反演的一維初始速度模型(圖13).然后，使用前文所述相關(guān)反演成像方法交叉反演，10次迭代后得到穩(wěn)定的縱橫波速度模型和VP/VS比值模型(圖15).Pg數(shù)據(jù)走時(shí)殘差均方根從初次迭代的0.336 s下降到10次迭代后的0.119 s，Sg數(shù)據(jù)走時(shí)殘差均方根從初次迭代的0.442 s下降到10次迭代后的0.169 s.與最終速度模型相應(yīng)的計(jì)算走時(shí)數(shù)據(jù)也被標(biāo)注在地震記錄剖面上(圖11，圖12)，可以看出我們的合成走時(shí)數(shù)據(jù)與觀測(cè)數(shù)據(jù)有良好的匹配效果.

圖11 札達(dá)—泉水溝寬角地震剖面P波地震記錄記錄剖面折合速度為6.0 km·s-1，十字表示拾取的Pg數(shù)據(jù)，曲線表示利用我們的層析成像系統(tǒng)計(jì)算的對(duì)應(yīng)最終VP模型的折合走時(shí)數(shù)據(jù).Fig.11 P-wave seismic record of Zhada-Quanshuigou wide-angle seismic profileThe record sections were plotted with a reduced velocity of 6.0 km·s-1, and the cross denotes the picked Pg data, while the curve represents the computed reduced traveltime data correspond to our final VP model using our tomographic system.

圖12 札達(dá)—泉水溝寬角地震剖面S波地震記錄記錄剖面折合速度為3.5 km·s-1，十字表示拾取的Sg數(shù)據(jù)，曲線表示利用我們的層析成像系統(tǒng)計(jì)算的對(duì)應(yīng)最終VS模型的折合走時(shí)數(shù)據(jù).Fig.12 S-wave seismic record of Zhada-Quanshuigou wide-angle seismic profileThe record sections were plotted with a reduced velocity of 3.5 km·s-1, and the cross denotes the picked Sg data, while the curve represents the computed reduced traveltime data correspond to our final VS model using our tomographic system.

圖13 用于札達(dá)—泉水溝寬角地震剖面Pg和Sg數(shù)據(jù)反演的一維速度模型Fig.13 The 1-D velocity models for inversion as for the Pg and Sg data of Zhada-Quanshuigou wide-angle seismic profile

圖14和圖15分別展示了使用Hole(1992)，Ammon和Vidale(1993)的反演代碼以及我們的相關(guān)反演成像方法所獲得的VP、VS和VP/VS比值模型.可以看出，基于Hole(1992)，Ammon和Vidale(1993)的反演方法分別反演縱橫波速度結(jié)構(gòu)的過(guò)程，由于缺少相互關(guān)聯(lián)約束，反演得到的VP和VS模型之間圖像特征存在較大區(qū)別(圖14a，14b)，互相關(guān)性較差，這無(wú)疑會(huì)使得VP/VS比值產(chǎn)生偏差或扭曲.與之對(duì)比用我們的相關(guān)反演成像方法得到的VP和VS模型(圖15a，15b)相互關(guān)聯(lián)的非常好，圖像特征較相似，證明我們采用新反演方法得到的上地殼結(jié)構(gòu)更為合理.圖16顯示的射線路徑和走時(shí)殘差均方根擾動(dòng)表明我們模型大多數(shù)地方的速度變化是可靠的.圖16e表示計(jì)算得到的VP與VS模型之間的相關(guān)性，其表明在剖面的某些位置存在明顯的正異常，如橫向坐標(biāo)50 km，深度10 km的位置，以及橫向坐標(biāo)120～140 km的位置，存在著明顯的地質(zhì)異常.

圖15a、15b顯示，該剖面最南端札達(dá)盆地下方，及獅泉河盆地下方，均存在明顯的低VP、低VS及低VP/VS異常，且異常向下延伸的深度達(dá)到8～10 km.這可能分別與印度洋板塊的北向俯沖及古班公—怒江洋盆的閉合有關(guān)，其北向俯沖向下延伸的幾何形態(tài)表明它可能是不同地質(zhì)歷史時(shí)期板塊階段俯沖碰撞所遺留的構(gòu)造痕跡或“化石”.推測(cè)其基本原理為，伴隨俯沖碰撞前陸盆地堆積大量海洋沉積覆蓋.在山脈隆起地形抬升過(guò)程中，水體被逐漸排泄出高原，從而使得巖土較為干燥.巖石在干燥條件下橫波速度伴隨壓力下降減小的幅度較慢，而縱波速度伴隨壓力下降減小的幅度較快，這一巖石物理基礎(chǔ)可以很好地解釋上述縱橫波速度比偏低的現(xiàn)象.

Kapp等(2003)的研究表明，獅泉河逆沖斷裂以南的獅泉河盆地覆蓋有大量中生代沉積，且獅泉河北向俯沖向下延深到10 km以下，這與我們的以上分析是一致的.

阿伊拉日居山的隆升剝蝕、巖漿活動(dòng)與喀喇昆侖斷裂的活動(dòng)密切相關(guān).對(duì)喀喇昆侖東南段阿伊拉日居山地區(qū)同構(gòu)造片麻巖-花崗巖的年代學(xué)研究認(rèn)為喀喇昆侖斷裂形成時(shí)代為23～25 Ma,持續(xù)變形到約12 Ma,之后阿伊拉日居山快速隆升(李海兵等，2006,2007,2008).王世鋒等(2008)對(duì)札達(dá)盆地內(nèi)新生代沉積給出了全面而系統(tǒng)的磁性地層學(xué)的時(shí)代約束,并探討了盆地的控盆構(gòu)造問(wèn)題, 揭示盆地內(nèi)新生代地層的古地磁年齡為9.5～2.6 Ma，且控盆構(gòu)造為喀喇昆侖斷裂.我們的成像結(jié)果顯示，與札達(dá)盆地及獅泉河盆地8～10 km的厚度相比較，噶爾盆地沉積厚度僅3～4 km，這可能與盆地形成機(jī)制不同有關(guān)：噶爾盆地為走滑斷陷盆地，而前兩者為板塊俯沖碰撞的前陸盆地.新生代以來(lái)伴隨喀喇昆侖斷裂在不同階段的走滑拉張運(yùn)動(dòng)，可能同步產(chǎn)生了阿伊拉日居山逐步抬升和噶爾盆地逐步沉降、斷層切割越來(lái)越深等地質(zhì)現(xiàn)象.

日松背斜及炮點(diǎn)SP6附近存在明顯高速異常體.炮點(diǎn)SP4南側(cè)的異常體具有較高的VP、VS和VP/VS比值，上述性質(zhì)表明近地表巖石主要由基性或超基性材料構(gòu)成，這也由部分暴露于地表的超基性鐵鎂質(zhì)巖石得到驗(yàn)證.

我們的VP與VS模型還顯示，從SP6到SP7的剖面下方存在另一個(gè)淺盆地，其南緣可能被一個(gè)相對(duì)較深的斷裂(深度大約6 km)所切割.除此之外，在著名的班公湖—怒江縫合帶北側(cè)還有可能存在著最大深度達(dá)8 km的淺斷裂.與此同時(shí)，盡管龍木錯(cuò)斷裂被認(rèn)為是羌塘地體和甜水海地體之間的主要邊界斷裂，我們的VP與VS模型并沒有顯示出與之相應(yīng)的顯著速度變化，這可能意味著它新生代的活動(dòng)并不活躍.

圖14 利用Ammon和Vidale (1993)反演代碼得到的札達(dá)—泉水溝剖面下的上地殼VP、VS和VP/VS模型對(duì)比(a) VP速度模型； (b) VS速度模型； (c) VP/VS模型.Fig.14 The upper crustal VP, VS and VP/VS model beneath the Zhada-Quanshuigou for comparison which were derived using the inversion codes of Ammon and Vidale (1993)(a) The derived VP model; (b) The derived VS model; (c) The resulted VP/VS model.

圖15 利用我們反演方法重建的札達(dá)—泉水溝寬角地震剖面下方的上地殼結(jié)構(gòu)(a) VP速度模型； (b) VS速度模型； (c) VP/VS模型.圖頂說(shuō)明表示沿剖面的地勢(shì)與構(gòu)造背景.Fig.15 The upper crustal structure beneath the Zhada-Quanshuigou wide-angle seismic profile, which were inverted using our tomographic system(a) The derived VP model; (b) The derived VS model; (c) The resulted VP/VS model.The illustration on the top in this figure shows the relief and the tectonic setting along the profile.

圖16 札達(dá)—泉水溝剖面下方上地殼速度模型二維反演結(jié)果評(píng)估(a) VP模型對(duì)應(yīng)的射線覆蓋； (b) VS模型對(duì)應(yīng)的射線覆蓋； (c) VP模型格點(diǎn)速度擾動(dòng)引起的走時(shí)殘差均方根變化值； (d) VS模型格點(diǎn)速度擾動(dòng)引起的走時(shí)殘差均方根變化值； (e) 計(jì)算得到VP與VS模型之間的相關(guān)系數(shù).Fig.16 Assessment of the 2-D inversion results as for our upper crustal velocity model beneath Zhada-Quanshuigou profile(a) Ray-coverage of derived VP model; (b) Ray-coverage of derived VS model; (c) RMS variation on VP model that corresponds to the velocity perturbation at each point; (d) RMS variation on VS model that corresponds to the velocity perturbation at each point；(e) The computed correlation coefficient between VP and VS model.

4 結(jié)論

當(dāng)前基于深地震測(cè)深縱橫波初至走時(shí)數(shù)據(jù)的層析成像結(jié)果越來(lái)越廣泛地被應(yīng)用于揭示造山帶、盆地等不同類型構(gòu)造域下方的上地殼結(jié)構(gòu)精細(xì)特征.本文基于縱橫波初至走時(shí)數(shù)據(jù)異常來(lái)源于同一地質(zhì)體，因而該地質(zhì)體附近縱橫波速度圖像特征密切相關(guān)這一思想，在Hole(1992)及Ammon和Vidale(1993)工作基礎(chǔ)上提出了反演上地殼縱橫波速度結(jié)構(gòu)的相關(guān)反演成像方法，這一方法可描述為：

(1)縱橫波速度反演交替進(jìn)行，在迭代反演過(guò)程中每通過(guò)一次反演獲得相應(yīng)的縱波速度(或橫波速度)結(jié)構(gòu)后，更新相應(yīng)的縱橫波速度比模型以及相應(yīng)的橫波(或縱波)速度反演的初始模型，然后開展下一步橫波(或縱波)速度反演工作.

(2)在反演過(guò)程中依據(jù)縱橫波速度結(jié)構(gòu)的相關(guān)性信息和射線路徑長(zhǎng)度將走時(shí)殘差以不同權(quán)重分配到射線路徑經(jīng)過(guò)的單元，依據(jù)網(wǎng)格節(jié)點(diǎn)周圍平均的慢度擾動(dòng)更新速度模型.

(3)正反演過(guò)程分別基于Vidale (1988)的有限差分走時(shí)計(jì)算方法和Hole (1992)的反投影成像方法.兩種典型模型試驗(yàn)表明，該方法應(yīng)用于上地殼速度結(jié)構(gòu)反演成像過(guò)程，可有效提高反演結(jié)果的可靠性，在很大程度上避免了常規(guī)單獨(dú)反演縱波和橫波速度過(guò)程容易帶來(lái)的畸變和失真.

此外我們發(fā)展了獲取初始一維速度模型及直接定量描述每個(gè)網(wǎng)格點(diǎn)速度值可靠性的解評(píng)價(jià)方法：通過(guò)計(jì)算每個(gè)網(wǎng)格點(diǎn)參數(shù)單獨(dú)變化引起的目標(biāo)函數(shù)——走時(shí)殘差均方根的改變量來(lái)評(píng)估觀測(cè)系統(tǒng)對(duì)特定格點(diǎn)速度變化的分辨能力.這種方法相比射線覆蓋數(shù)評(píng)價(jià)方法更為簡(jiǎn)潔直觀.相關(guān)反演最終結(jié)果所顯示縱橫波速度圖像特征之間的相關(guān)性，則可以直接用于估計(jì)地質(zhì)異常體存在與否及可能的位置.

以上相關(guān)反演成像方法被成功地應(yīng)用于重建青藏高原西部札達(dá)—泉水溝深地震測(cè)深(DSS)剖面下方的上地殼速度結(jié)構(gòu)，揭示出與青藏高原西緣板塊碰撞相關(guān)的上地殼速度結(jié)構(gòu)特征：札達(dá)盆地及獅泉河盆地下方存在明顯的VP和VS低速異常及低VP/VS異常，可能是印度板塊北向俯沖及班公湖—怒江特提斯洋閉合在上地殼遺留的地球物理“化石”和證據(jù).噶爾盆地沉積厚度較小，其形成與新生代以來(lái)喀喇昆侖斷裂的走滑拉張運(yùn)動(dòng)有關(guān).

致謝 感謝J.A.Hole教授、C.J.Ammon教授和J.E.Vidale教授提供用于該項(xiàng)研究的反演軟件，并衷心感謝來(lái)自田小波研究員、徐濤副研究員和梁曉峰副研究員的寶貴意見和建議.謹(jǐn)以此文紀(jì)念中國(guó)科學(xué)院地質(zhì)與地球物理研究所張忠杰研究員(1964—2013).

Ammon C J, Vidale J E. 1993. Tomography without rays.Bull.Seism.Soc.Am., 83(2): 509-528.

Bai Z M, Zhang Z J, Wang Y H. 2007. Crustal structure across the Dabie-Sulu orogenic belt revealed by seismic velocity profiles.J.Geophys.Eng., 4(4): 436-442.

Bregman N D, Bailey R C, Chapman C H. 1989. Crosshole seismic tomography.Geophysics, 54(2): 200-215, doi: 10.1190/1.1442644.

Caldwell W B, Klemperer S L, Rai S S, et al. 2009. Partial melt in the upper-middle crust of the northwest Himalaya revealed by Rayleigh wave dispersion.Tectonophysics, 477(1-2): 58-65.

Conder J A, Wiens D A. 2006. Seismic structure beneath the Tonga arc and Lau back-arc basin determined from jointVP,VP/VStomography.Geochem.Geophys.Geosyst., 7(3): Q03018, doi: 10.1029/2005GC001113.

Ding Z F, He Z Q, Sun W G, et al. 1999. 3-D crust and upper mantle velocity structure in eastern Tibetan plateau and its surrounding areas.ChineseJ.Geophys. (in Chinese), 42(2): 197-205.Duan Y H, Lai X L, Zhang X K, et al. 1999. Finite-difference 2-D and 3-D seismic traveltime velocity tomography.NorthChinaEarthquakeSciences(in Chinese), 17(4): 53-60.

Eberhart-Phillips D. 1990. Three-dimensional P and S velocity structure in the Coalinga Region, California.J.Geophys.Res.:SolidEarth, 95(B10): 15343-15363.

Hole J A. 1992. Nonlinear high-resolution three-dimensional seismic travel time tomography.J.Geophys.Res., 97(B5): 6553-6562.

Hole J A, Zelt B C.1995. 3-D finite-difference reflection traveltimes.Geophys.J.Int., 121(2): 427-434.

Kapp P, Murphy M A, Yin A, et al. 2003. Mesozoic and Cenozoic tectonic evolution of the Shiquanhe area of western Tibet.Tectonics, 22(4), doi: 10.1029/2001TC001332.

Karplus M S, Zhao W, Klemperer S L, et al. 2011. Injection of Tibetan crust beneath the south Qaidam Basin: Evidence from INDEPTH IV wide-angle seismic data.J.Geophys.Res., 116(B7): B07301, doi: 10.1029/2010JB007911.

Kissling E, Ellsworth W L, Eberhart-Phillips D E, et al. 1994. Initial reference models in local earthquake tomography.J.Geophys.Res.:SolidEarth, 99(B10): 19635-19646.

Kissling E, Solarino S, Cattaneo M. 1995. Improved seismic velocity reference model from local earthquake data in Northwestern Italy.TerraNova, 7(5): 528-534.

Lacassin R, Valli F, Arnaud N, et al. 2004. Large-scale geometry, offset and kinematic evolution of the Karakorum fault, Tibet.EarthPlanet.Sci.Lett., 219(3-4): 255-269.

Lan H Q, Zhang Z, Xu T, et al. 2012a. A comparative study on the fast marching and fast sweeping methods in the calculation of first-arrival traveltime field.ProgressinGeophys. (in Chinese), 27(5): 1863-1870, doi: 10.6038/j.issn.1004-2903.2012.02.005.Lan H Q, Zhang Z, Xu T, et al. 2012b. Effects due to the anisotropic stretching of the surface-fitting grid on the traveltime computation for irregular surface by the coordinate transforming method.ChineseJ.Geophys. (in Chinese), 55(10): 3355-3369, doi: 10.6038/j.issn.0001-5733.2012.10.018.

Li F, Xu T, Wu Z B, et al. 2013. Segmentally iterative ray tracing in 3-D heterogeneous geological models.ChineseJ.Geophys. (in Chinese), 56(10): 3514-3522, doi: 10.6038/cjg20131026.Li H B, Valli F, Xu Z Q, et al. 2006. Deformation and tectonic evolution of the Karakorum fault, western Tibet.GeologyinChina(in Chinese), 33(2): 239-255.

Li H B, Valli F, Liu D Y, et al. 2007. The formation age of the Karakorum Fault in western Tibet: Constraints from SHRIMP U-Pb dating of zircons.ChineseScienceBulletin, 52(4):438-447.Li H B, Valli F, Arnaud N, et al. 2008. Rapid uplifting in the process of strike-slip along the Karakorum fault zone in western Tibet: Evidence from40Ar/39Ar thermochronology.ActaPetrologicaSinica(in Chinese), 24(7): 1552-1556.

Liu Y F, Lan H Q. 2012. Study on the numerical solutions of the eikonal equation in curvilinear coordinate system.ChineseJ.Geophys. (in Chinese), 55(6): 2014-2026, doi: 10.6038/j.issn.0001-5733.2012.06.023.

Murphy M A, Yin A, Kapp P, et al. 2000. Southward propagation of the Karakoram fault system, southwest Tibet: Timing and magnitude of slip.Geology, 28(5): 451-454.

Saito H. 1989. Traveltimes and raypaths of first-arrival seismic waves: Computation method based on Huygens′ Principle. ∥ 59thAnn. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstract, 244-247.Tarantola A, Valette B. 1982. Generalized nonlinear inverse problems solved using the least squares criterion.Rev.Geophys.SpacePhys., 20(2): 219-232.Teng J W, Zhang Z J, Zhang X K, et al. 2013. Investigation of the Moho discontinuity beneath the Chinese mainland using deep seismic sounding profiles.Tectonophysics, 609: 202-216, doi: 10.1016/j.tecto.2012.11.024.

Thurber C H. 1983. Earthquake locations and three-dimensional crustal structure in the Coyote Lake area, central California.J.Geophys.Res., 88(B10): 8226-8236.

Thurber C H, Atre S R. 1993. Three-dimensionalVP/VSvariations along the Loma Prieta rupture zone.Bull.Seism.Soc.Am., 83(3): 717-736.

Thybo H, Nielsen C A. 2009. Magma-compensated crustal thinning in continental rift zones.Nature, 457(7231): 873-876, doi: 10.1038/nature07688.

Valli F, Arnaud N, Leloup P H, et al. 2007. Twenty million years of continuous deformation along the Karakorum fault, western Tibet: a thermochronological analysis.Tectonics, 26(4), doi: 10.1029/2005TC001913.

Vidale J E. 1988. Finite-difference traveltime calculation.Bull.Seism.Soc.Am., 78: 2062-2076.

Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions.Geophysics, 55(5): 521-526.

Wagner L S, Beck S, Zandt G. 2005. Upper mantle structure in the south central Chilean subduction zone (30°to 36°S).J.Geophys.Res.:SolidEarth, 110(B1): B01308, doi: 10.1029/2004JB003238.

Wang C Y, Zhang X K, Ding Z F, et al. 1997. Finite-difference tomography of upper crustal structure in Dabieshan orogenic belt.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 40(4): 495-502.

Wang C Y, Mooney W D, Wang X L, et al. 2002. Study on 3-D velocity structure of crust and upper mantle in Sichuan-Yunnan region, China.ActaSeismologicaSinica(in Chinese), 24(1): 1-16, doi: 10.3321/j.issn:0253-3782.2002.01.001.Wang F Y, Duan Y H, Yang Z X, et al. 2008. Velocity structure and active fault of Yanyuan-Mabian seismic zone—The result of high-resolution seismic refraction experiment.ScienceinChinaSeriesD:EarthSciences, 51(9): 1284-1296.

Wang S F, Zhang W L, Fang X M, et al. 2008. Magnetic & stratigraphic characteristics and tectonic significance of Zhada Basin in Southwest Tibet.ScienceinChina(in Chinese), 53(6): 676-683.

Wang Z. 2014. Joint inversion of P-wave velocity andVP-VSratio: imaging the deep structure in NE Japan.AppliedGeophysics, 11(2): 119-127.Wittlinger G, Vergne J, Tapponnier P, et al. 2004. Teleseismic imaging of subducting lithosphere and Moho offsets beneath western Tibet.EarthPlanet.Sci.Lett., 221(1-4): 117-130.

Xu T, Xu G M, Gao E G, et al. 2004. Block modeling and shooting ray tracing in complex 3D media.ChineseJ.Geophys. (in Chinese), 47(6): 1118-1126, doi: 10.3321/j.issn.0001-5733. 2004.06.027.

Xu T, Xu G M, Gao E G, et al. 2006. Block modeling and segmentally iterative ray tracing in complex 3D media.Geophysics, 71(3): T41-T51.

Xu T, Zhang Z, Gao E, et al. 2010. Segmentally iterative ray tracing in complex 2D and 3D heterogeneous block models.Bull.Seismol.Soc.Am., 100(2): 841-850.

Xu T, Li F, Wu Z B, et al. 2014. A successive three-point perturbation method for fast ray tracing in complex 2D and 3D geological models.Tectonophysics, 627: 72-81, doi: 10.1016/j.tecto.2014.02.012.Xu Z Q, Li H B, Tang Z M, et al. 2011. The transformation of the terrain structures of the Tibet Plateau through large-scale strike-slip faults.ActaPetrologicaSinica(in Chinese), 27(11): 3157-3170.

Zelt C A. 1999. Modelling strategies and model assessment for wide-angle seismic traveltime data.Geophys.J.Int., 139(1): 183-204.Zhang H J, Thurber C, Bedrosian P. 2009a. Joint inversion forVP,VS, andVP/VSat SAFOD, Parkfield, California.Geochem.Geophys.Geosyst., 10(11), doi: 10.1029/2009GC002709.

Zhang J, Toksozf M N. 1998. Nonlinear refraction traveltime tomography.Geophysics, 63(5): 1726-1737.

Zhang X K, Wang C Y, Liu G D, et al. 1996. Fine crustal structure in Yanqing-Huailai region by deep seismic reflection profiling.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 39(3): 356-364.

Zhang Z J, Bai Z M, Mooney W, et al. 2009b. Crustal structure across the Three Gorges area of the Yangtze platform, central China, from seismic refraction/wide-angle reflection data.Tectonophysics, 475(3-4): 423-437, doi: 10.1016/j.tecto.2009.05.022.

Zhang Z J, Klemperer S, Bai Z M, et al. 2011. Crustal structure of the Paleozoic Kunlun orogeny from an active-source seismic profile between Moba and Guide in east Tibet, China.GondwanaResearch, 19(4): 994-1007.

Zhao A H, Zhang Z J, Teng J W. 2004. Minimum travel time tree algorithm for seismic ray tracing: Improvement in effciency.J.Geophys.Eng., 1(4): 245-251, doi:10.1088/1742-2132/1/4/001.

Zhao A H, Ding Z F. 2005. A double-grid algrithm for calculating traveltimes of wide-angle reflection waves.ChineseJ.Geophys. (in Chinese), 48(5) : 1141-1147.

Zhao F F, Ma T, Xu T. 2014. A review of the travel-time calculation methods of seismic first break.ProgressinGeophys. (in Chinese), 29(3): 1102-1113, doi: 10.6038/pg20140313.

附中文參考文獻(xiàn)

丁志峰, 何正勤, 孫為國(guó)等. 1999. 青藏高原東部及其邊緣地區(qū)的地殼上地幔三維速度結(jié)構(gòu). 地球物理學(xué)報(bào), 42(2): 197-205.

段永紅, 賴曉玲, 張先康等. 1999. 有限差分二、三維速度層析成象. 華北地震科學(xué), 17(4): 53-60.

蘭海強(qiáng), 張智, 徐濤等. 2012a. 地震波走時(shí)場(chǎng)模擬的快速推進(jìn)法和快速掃描法比較研究.地球物理學(xué)進(jìn)展, 27(5): 1863-1870, doi: 10.6038/j.issn.1004-2903.2012.02.005.

蘭海強(qiáng), 張智, 徐濤等. 2012b. 貼體網(wǎng)格各向異性對(duì)坐標(biāo)變換法求解起伏地表下地震初至波走時(shí)的影響.地球物理學(xué)報(bào), 55(10):3355-3369, doi: 10.6038/j.issn.0001-5733.2012.10.018.

李飛, 徐濤, 武振波等. 2013. 三維非均勻地質(zhì)模型中的逐段迭代射線追蹤. 地球物理學(xué)報(bào), 56(10): 3514-3522, doi: 10.6038/cjg20131026.

李海兵, Valli F, 許志琴等. 2006. 喀喇昆侖斷裂的變形特征及構(gòu)造演化. 中國(guó)地質(zhì), 33(2): 239-255.

李海兵, Valli F, 劉敦一等. 2007. 喀喇昆侖斷裂的形成時(shí)代: 鋯石SHRIMP U-Pb年齡的制約. 科學(xué)通報(bào), 52(4): 438-447.

李海兵, Valli F, Arnaud N等. 2008. 喀喇昆侖斷裂帶走滑過(guò)程中伴隨的快速隆升作用:熱年代學(xué)證據(jù). 巖石學(xué)報(bào), 24(7): 1552-1556.

劉一峰, 蘭海強(qiáng). 2012. 曲線坐標(biāo)系程函方程的求解方法研究．地球物理學(xué)報(bào), 55(6): 2014-2026, doi: 10.6038/j.issn.0001-5733.2012.06.023.

王椿鏞, 張先康, 丁志峰等. 1997. 大別造山帶上部地殼結(jié)構(gòu)的有限差分層析成像. 地球物理學(xué)報(bào), 40(4): 495-502.

王椿鏞, Mooney W D, 王溪莉等. 2002. 川滇地區(qū)地殼上地幔三維速度結(jié)構(gòu)研究. 地震學(xué)報(bào), 24(1): 1-16, doi: 10.3321/j.issn:0253-3782.2002.01.001.

王夫運(yùn), 段永紅, 楊卓欣等. 2008. 川西鹽源—馬邊地震帶上地殼速度結(jié)構(gòu)和活動(dòng)斷裂研究——高分辨率地震折射實(shí)驗(yàn)結(jié)果. 中國(guó)科學(xué), 38(5): 611-621.

王世鋒, 張偉林, 方小敏等. 2008. 藏西南札達(dá)盆地磁性地層學(xué)特征及其構(gòu)造意義. 中國(guó)科學(xué), 53(6): 676-683.

徐濤, 徐果明, 高爾根等. 2004. 三維復(fù)雜介質(zhì)的塊狀建模和試射射線追蹤. 地球物理學(xué)報(bào), 47(6): 1118-1126, doi: 10.3321/j.issn:0001-5733.2004.06.027.

許志琴, 李海兵, 唐哲民等. 2011. 大型走滑斷裂對(duì)青藏高原地體構(gòu)架的改造. 巖石學(xué)報(bào), 27(11): 3157-3170.

張先康, 王椿鏞, 劉國(guó)棟等. 1996. 延慶—懷來(lái)地區(qū)地殼細(xì)結(jié)構(gòu)——利用深地震反射剖面. 地球物理學(xué)報(bào), 39(3): 356-364.

趙愛華, 丁志峰. 2005. 寬角反射地震波走時(shí)模擬的雙重網(wǎng)格法. 地球物理學(xué)報(bào), 48(5): 1141-1147.

趙烽帆, 馬婷, 徐濤. 2014. 地震波初至走時(shí)的計(jì)算方法綜述. 地球物理學(xué)進(jìn)展, 29(3): 1102-1113, doi: 10.6038/pg20140313.

(本文編輯 何燕)

Tomographic imaging of velocity structure in upper crust based on correlated inversion ofVPandVS

WANG Xiao1,2, ZHOU Xiao-Peng1,2, ZHANG Xin-Yan1,2, BAI Zhi-Ming1*, TENG Ji-Wen1

1StateKeyLaboratoryofLithosphereEvolution,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China

Presently seismic tomographic methods based on first arrivals of deep seismic sounding traveltime data have been widely used to image the upper crustal velocity structure. However, generally there are some difficulties in the inversion process using conventional methods: (1) How to construct reliable upper crust models using Pg and Sg data; (2) How to quantitatively estimate the reliability of reconstructed models; and (3) How to decrease the uncertainty and distortion validly.Based on the idea that theVPandVSanomalies from a same geological body are closely related each other, we develop a correlated inversion method:VPandVSinversion are carried out alternately, in which the correspondingVP/VSmodel and initialVS(orVP) model are updated using the obtainedVP(orVS) model after each inversion. The slowness perturbation resulted from the traveltime residuals is weighted and projected to the grid-points on the ray trace based on the updated correlation information, thus the velocity models are updated. Forward and inversion processes are respectively based on finite-difference calculations of travel times and the back-projection inversion method.Correlated inversion ofVPandVSmodels and sequential inversion technique are incorporated into the tomographic process to ensure the reliability of inversion results and avoid the distortion of the resultingVP/VSmodel. Comparison with the conventional inversion results using the code of Ammon and Vidale (1993), both of the derivedVPandVSmodels by our approach have been obviously improved, which are clearly closer to the original theoretical ones. While those constructed using the code of Ammon and Vidale (1993) have been seriously distorted. Moreover, the maximumVP/VSerror based on the conventional approach is up to 0.2, while the maximumVP/VSerror with our approach is only 0.02～0.03.Tests using two typical models show that the technology used in inversion tomography of the velocity structure within upper crust, which can effectively improve the reliability of the inversion results, can avoid the uncertainty and distortion resulted from the conventional singleVPandVSinversion process on a large scale. This method has been successfully applied to the western Tibetan plateau to reconstruct the velocity structure of upper crust beneath the Zhada-Quanshuigou deep seismic sounding profile. The derived models reveal the remained upper crustal velocity structure features resulted from plate collision in the western Tibetan plateau. Two similar obvious low-value anomalies ofVP,VSandVP/VSwith the maximum depth 8～10 km are present beneath the Zhada and Shiquanhe basin, which are probably the effect of the subduction of Indian plate beneath the Eurasian continent since Cenozoic or the closure of the Mesozoic Bangong-Nujiang oceanic basin. The disclosed Gar basin has a sedimentary thickness of 3～4 km, where the sedimentation and deep processes may be due to the gradual uplift of the Ayilariju Mountains and the strike-slip movement of the Karakoram fault.

Upper crust;VPandVSstructure; Correlated inversion; Finite-difference; Back-projection

10.6038/cjg20151011.

Wang X, Zhou X P, Zhang X Y, et al. Tomographic imaging of velocity structure in upper crust based on correlated inversion ofVPandVS.ChineseJ.Geophys. (in Chinese),58(10):3553-3570,doi:10.6038/cjg20151011.

中國(guó)地震局公益性行業(yè)科研專項(xiàng)(201408023)和國(guó)家自然科學(xué)基金(41374062,41174075，41274070，41474068)聯(lián)合資助.

王曉，女，1989年生，博士研究生，主要從事地殼結(jié)構(gòu)及地震資料偏移成像研究.E-mail: wangxiao@mail.iggcas.ac.cn

*通訊作者 白志明，男，1971年生，副研究員，主要研究方向?yàn)榈貧そY(jié)構(gòu)成像及動(dòng)力學(xué)研究.E-mail: bbzzmm@mail.iggcas.ac.cn

10.6038/cjg20151011

P315

2014-12-29，2015-09-15收修定稿

王曉， 周小鵬， 張新彥等. 2015. 上地殼縱橫波速度結(jié)構(gòu)相關(guān)反演成像方法.地球物理學(xué)報(bào)，58(10):3553-3570,