基于局部一維模型與走時標定的區(qū)域三維速度模型構建技術研究

張慧民 靳 平 劉文學 李 欣 張銳迎

1）中國西安710024西北核技術研究所

2）中國西安710049西安交通大學

3）中國西安710043北方光電股份有限公司

引言

準確定位對于《全面禁止核試驗條約》（comprehensive nuclear test ban treaty，簡寫為CTBT）的地震監(jiān)測工作具有重要意義.實踐中面臨的主要挑戰(zhàn)之一是稀疏臺網條件下對小震級地震事件進行精確定位.目前地震定位中采用的主要是基于IASP91（Kennett，Engdahl，1991）或AK135（Kennett et al，1995）一維地球模型的走時表.雖然這些一維模型很好地反映了全球的平均性質，但實際地球是三維的，穿過三維介質的實際震相走時與由一維模型導出的理論走時之間往往存在偏差，尤其是對區(qū)域震相，在很多地區(qū)這種偏差可能很嚴重.在我國的新疆地區(qū)，地震事件在部分區(qū)域臺站上的實際走時與理論走時之間就存在著較為明顯的偏差.這種區(qū)域震相走時的偏差可能會影響到地震臺網的定位精度.為解決這一問題，目前通常采用兩種方法：① 利用爆炸或精確定位的大地震進行區(qū)域走時標定；②構建監(jiān)測區(qū)域的地殼與上地幔三維模型，并據(jù)此得到高精度區(qū)域震相走時.相比之下，第二種方法更能反映問題的物理實質，并具有更強的實用性.

為得到能夠提高區(qū)域地震定位精度的三維速度模型，以往主要通過兩種方式進行模型構建：一種方式是直接反演法，即將模型區(qū)域劃分為三維網格并反演其速度結構，或反演模型區(qū)域內部分觀測點下方的一維模型后進行插值得到三維模型；另一種方式是利用已有的地震學研究成果，結合該地區(qū)的地質構造等資料，以震中與發(fā)震時間已知的地震事件（ground truth，簡寫為GT）為參考，對已有的介質模型進行選擇、組合及必要的修正，從而得到準確反映區(qū)域走時特點的三維介質模型（Russian Federation／United States Calibration Working Group，2002；Steck et al，2004；Cormier，2001；Ryaboy et al，2001）.對禁核試地震監(jiān)測工作而言，構建區(qū)域性地殼與上地幔三維速度模型的主要目的在于提高區(qū)域地震定位能力.出于此目的，Ryaboy等（2001）采用第二種方式構建了北美、北歐地區(qū)的三維介質模型，有效減小了區(qū)域地震定位誤差.其作法的核心為：根據(jù)標定事件繪制區(qū)域走時偏差圖（即從由標定事件獲取的走時圖上減去由IASP91模型得出的走時圖），將其與現(xiàn)有模型的理論走時偏差圖進行對比，找出兩者的主要差異，然后結合附近的一維模型、有關地質構造信息等，確定如何修正現(xiàn)有模型以減小兩走時偏差圖的差異，并據(jù)此得到新的模型.如此不斷修正現(xiàn)有模型以減小理論走時偏差圖與標定事件走時偏差圖之間的差異，最終即可獲得滿意的區(qū)域三維模型.與直接反演法不同，這一構建方法工作量很大且對人員的專業(yè)經驗要求較高，實現(xiàn)難度較大.為此我們對其進行了改進，以類似直接反演的方法為構建過程提供數(shù)值參考，從而降低了此類構建過程的難度并減小了最終模型對人員主觀因素的依賴.

1 方法

在已往的此類模型構建方法中，如何對已有模型進行修正，主要是由構建人員根據(jù)走時偏差圖和自身的經驗確定的.在實踐中，這種做法需要構建人員具有非常豐富的經驗，而且其主觀因素會對模型的形成產生比較大的影響，這使得模型構建工作實施起來比較困難，且結果可能不夠客觀.為此，本文提出了借助于體波走時反演方法，為構建過程中如何修正模型提供較為客觀的、半定量的參考意見，據(jù)此協(xié)助構建人員實現(xiàn)區(qū)域模型的構建，以降低構建難度并減少對構建人員個人經驗的依賴程度.與一般直接反演方法不同的是，我們只是利用GT事件和走時反演方法，為構建人員修正現(xiàn)有模型提供參考數(shù)據(jù)，而并非直接進行介質速度結構反演.

本文的三維模型構建工作主要按以下步驟進行：① 根據(jù)已有的模型數(shù)據(jù)和地質資料，給出區(qū)域初始模型；② 計算現(xiàn)有模型走時偏差圖，并以走時反演法對該模型進行計算，記錄此方法對現(xiàn)有模型的修正量；③觀察走時偏差圖，找出主要的走時偏差異常，對比走時反演法對現(xiàn)有模型的修正數(shù)據(jù)，找出其中與主要走時偏差異常相關的部分，根據(jù)現(xiàn)有資料判斷是否需要對此修正數(shù)據(jù)再進行必要的校正（例如，反演法所得修正量過大，會導致模型與實際介質不符等），之后以校正后的修正數(shù)據(jù)對現(xiàn)有模型進行調整，得出一個新模型；④反復進行上述②、③兩步驟，直至走時偏差圖上不再有明顯的異常；⑤以走時反演法對已無明顯走時偏差異常的新模型進行微調，得到最終的區(qū)域三維介質模型；⑥ 檢驗所得三維模型在區(qū)域地震定位中的作用.

為實現(xiàn)上述模型構建步驟，需解決如下幾個問題：

1）三維介質模型中的射線追蹤.我們實現(xiàn)了二維格點化介質中計算初至波走時的初至波逐列遞推算法（Schneider et al，1992），通過做剖面將其推廣到三維介質的初至波走時計算中，并針對介質速度結構反演、地震定位等不同應用對其進行了相應的改進.實驗表明，此算法計算速度快，且相對計算誤差小于1‰，完全可以滿足后續(xù)工作的需要.

2）理論走時偏差圖計算.實踐中直接逐點計算走時偏差圖計算量較大.為此，我們結合數(shù)值插值設計了合理的格點計算順序，使區(qū)域走時偏差圖的計算量降低了約一個量級.

3）走時法介質三維速度模型反演.以射線追蹤為基礎，采用聯(lián)合迭代重建方法進行.為提高射線追蹤算法效率，本文設計了基于逐列遞推算法的反方向射線追蹤算法，即首先由逐列遞推算法計算出地震射線在接收點的出射角；然后由接收點出發(fā)，以出射角的反方向作為入射角，計算出射線傳播路徑；之后，以此射線為參考，采用線性修正的方式修改反方向入射角，直至所得射線路徑滿足預設要求.

4）三維介質中的區(qū)域地震定位.本文采用Powell共軛方向算法進行.

2 數(shù)據(jù)

本文所用數(shù)據(jù)主要包括GT事件、部分臺站下方的一維速度模型、相關地震地質資料等.

GT事件數(shù)據(jù)引自新疆維吾爾自治區(qū)地震局2006—2008年的地震公報.在其標稱定位精度為1級（誤差小于5km）的地震事件中，選取滿足數(shù)據(jù)質量與臺網評定質量均為1、定位臺站數(shù)不少于8、臺站最大張角小于90°、定位臺站最小震中距小于80km等條件的286個地震事件作為GT事件.后面的工作中我們認為這些事件的震中參數(shù)是準確的.我們將挑選出的GT事件分為兩組：其中一組包含146個事件（圖1），含P波走時數(shù)據(jù)2 004條，S波走時數(shù)據(jù)1 883條，用于三維介質模型構建；另外一組的140個GT事件則作為震中參數(shù)已知的參考事件用于區(qū)域地震定位實驗，以檢驗所構建三維模型用于定位時的效果.模型構建過程中使用的地震臺站均為新疆地震臺網臺站.

圖1 模型構建所用數(shù)據(jù)及臺站分布圖藍色圓點表示GT事件，紅色三角表示地震臺站，黃色方塊則表示該位置已有一維模型Fig.1 Data and stations used in the construction of 3D velocity model Blue dots denote GT events，red triangles represent stations，yellow squares stand for locations of 1Dpartial velocity model

已有的地震地質資料包括Crustal 2.0全球地殼模型數(shù)據(jù)，新疆地區(qū)莫霍面深度等高線圖，以及新疆及周邊地區(qū)介質結構的相關研究成果（胥頤等，2006；姜枚等，1999；李順成等，2005；米寧等，2005；李海鷗等，2006；錢輝等，2006；高銳等，2001，2002；賀日政等，2001；趙俊猛等，2003，2008；馮梅，安美建，2007；李秋生等，2001；王有學等，2004；楊主恩等，2005；楊少敏等，2008；周仕勇，1999）等.

區(qū)域內的部分一維速度模型是通過接收函數(shù)與面波頻散曲線聯(lián)合反演方法獲得的，其分布如圖1所示.

3 模型構建

首先根據(jù)GT事件走時數(shù)據(jù)得到各個臺站的實測走時偏差圖，然后將新疆地區(qū)劃分為1°×1°的網格，根據(jù)已有的部分一維速度模型生成初始三維介質模型.對于存在多個一維模型的網格，取其所有一維模型的均值；對于無一維模型與之對應的網格，則采用克里金空間插值方法獲得.之后，計算現(xiàn)有模型走時偏差圖，計算并記錄體波走時反演法對現(xiàn)有模型的修正量；結合現(xiàn)有模型走時偏差圖與實測走時偏差圖的主要差異，找出相應的反投影法對現(xiàn)有模型進行修正的修正量，結合其它資料判斷是否需要對這些修正數(shù)據(jù)再進行必要的校正；最后以校正后的修正數(shù)據(jù)對現(xiàn)有模型進行調整，得出一個新模型.反復進行上述步驟，直至所得三維模型的走時偏差圖與實測走時偏差圖之間不再有明顯的不同.

圖2給出了所構建三維模型在不同深度上的S波速度分布.

圖2 所構建的三維速度模型切面圖Fig.2 Profiles of modified 3Dvelocity model

4 定位實驗

為檢驗所構建三維模型是否能夠有效地減小區(qū)域震相走時殘差并提高區(qū)域地震定位精度，我們選取了4種不同的介質模型以便進行對比：IASP91一維地球模型，Crustal 2.0模型，目標區(qū)域初始三維模型，修正后的三維介質模型.

重新定位后各事件震中相對于參考事件原震中位置的偏差如圖3所示.圖中各參考事件原震中位置均位于極坐標原點（未顯示）.重新定位后的震中位置則如圖3中小圓圈所示，即圖中小圓圈的空間坐標代表各地震事件重新定位后的震中與參考事件震中的偏差.上述4種模型的平均定位偏差分別為5.0，5.8，5.1和3.8km，修正模型中的定位偏差明顯小于其它模型.另外，從圖3中不難發(fā)現(xiàn)，修正模型中定位偏差更加集中于原點且空間分布比較均勻，不像在Crustal 2.0模型中系統(tǒng)偏南，而IASP91模型中則系統(tǒng)偏東.多少有些讓人意外的是，對于以局部一維模型為基礎得到的初始三維模型，其定位偏差與一維的IASP91模型相仿，并未像預期的那樣對減小定位偏差起到作用.

對于重新定位后各事件的震源深度偏差（相對于參考事件），對IASP91模型定位結果而言，幾乎所有的定位深度都明顯小于參考事件深度，而修正模型的定位深度與參考事件深度最為接近，且基本上不再有系統(tǒng)偏差，如圖4所示.就重新定位后各事件的發(fā)震時間偏差而言，IASP91模型、Crustal 2.0模型與初始模型的時間偏差相當；而修正模型則小得多，與震源深度偏差類似.修正模型消除了其它幾個模型中發(fā)震時間系統(tǒng)偏晚的現(xiàn)象.對定位實驗中走時殘差的統(tǒng)計表明，修正模型的走時殘差分布明顯優(yōu)于其它幾個模型.圖5給出了初始模型與修正后模型的走時殘差統(tǒng)計圖.由圖5可以看出，修正后模型的走時殘差明顯小于初始模型，它們更加對稱地集中于零殘差附近，消除了初始模型中震相走時系統(tǒng)偏小的情形.

圖5 初始模型與修正后模型的走時殘差統(tǒng)計圖Fig.5 Residuals in initial 3Dvelocity model and modified 3Dvelocity model

為檢驗所構建模型在稀疏臺網條件下對小震級事件的定位能力，我們進行了相應的模擬實驗.所謂稀疏臺網條件下的小震級事件定位，形式上基本可以理解為：記錄到的地震事件的臺站比較少，因此參與定位的臺站就比較少，臺站對事件的張角可能比較大，可能沒有距震源很近的臺站等.為此，我們依然采用前面實驗中使用的參考事件，只是隨機剔除其一部分記錄臺站，使其臺站最大張角與最小震中距擴大.模擬實驗表明，此時各模型的定位偏差均較前面實驗有所擴大.相比而言，修正模型中擴大得最小，且其定位結果沒有系統(tǒng)偏差，其走時殘差的分布也較其它模型更為合理.

5 討論與結論

綜上所述，通過綜合利用目標區(qū)域內的局部一維模型與相關資料，并結合走時標定方法，我們構建了新疆地區(qū)三維介質模型.用于地震定位時，該模型的走時殘差明顯小于初始模型，且所得結果與參考事件各參數(shù)的相似程度明顯優(yōu)于初始模型及全球平均模型.其中初始模型中對各參考事件的定位偏差平均為5.1km，而修正后模型的定位偏差平均為3.8km，即從統(tǒng)計意義上講，修正后模型的定位偏差較初始模型可減少約20%.這說明本文的構建方法是行之有效的，確能達到預期目的——通過構建區(qū)域三維介質模型，減小區(qū)域震相走時殘差，提高區(qū)域地震定位精度.

令人有些意外的是，初始三維模型中的定位偏差與一維IASP91模型相仿，并未像預期的那樣對減小定位偏差起到明顯作用.但結合定位深度偏差與發(fā)震時刻偏差來看，初始三維模型仍有明顯改善.如圖4所示，其定位深度偏差明顯要比IASP91模型小很多，而且消除了系統(tǒng)偏差，對于發(fā)震時刻偏差也有類似結果.從圖5可以看出，初始三維模型中的走時有明顯的系統(tǒng)偏差，這應該是造成其震中位置偏差較大的主要原因.

另外，從原理上講，三維模型地震定位只需將普通定位方法中的走時計算部分的一維模型更換為三維模型即可，但實踐中直接計算三維模型中的走時十分耗時，實用性很差，應以查表法代替.

在模型構建的過程中，構建人員綜合其它相關資料，增加了所構建模型的合理性.但這一過程帶給構建人員的工作量相當大，盡管我們引入的對模型修正的數(shù)值參考已經大大降低了這一工作量.這也是這一方法仍有待繼續(xù)改進的方面.

馮梅，安美建.2007.中國大陸中上地殼剪切波速結構［J］.地震學報，29（4）：337－－347.

高銳，肖序常，劉訓，管燁，李秋生，盧德源，李朋武.2001.新疆地學斷面深地震反射剖面揭示的西昆侖—塔里木結合帶巖石圈細結構［J］.地球學報，22（6）：547－－552.

高銳，肖序常，高弘.2002.西昆侖—塔里木—天山巖石圈深地震探測綜述［J］.地質通報，21（1）：11－－18.

賀日政，高銳，李秋生，管燁，李朋武.2001.新疆天山（獨山子）—西昆侖（泉水溝）地學斷面地震與重力聯(lián)合反演地殼構造特征［J］.地球學報，22（6）：553－－558.

姜枚，許志琴，薛光琦，史大年.1999.青海茫崖—新疆若羌地震探測剖面及其深部構造的研究［J］.地質學報，73（2）：153－－161.

李海鷗，姜枚，王亞軍，張立樹，于更新.2006.新疆富蘊—庫爾勒剖面接收函數(shù)方法獲得的地殼上地幔結構成像［J］.地質學報，80（1）：135－－141.

李秋生，盧德源，高銳.2001.新疆地學斷面（泉水溝—獨山子）深地震探測成果綜合研究［J］.地球學報，22（6）：534－－540.

李順成，劉啟元，陳九輝，郭飆，賴院根，王繼.2005.橫跨天山的寬頻帶流動地震臺陣觀測［J］.地球物理學進展，20（4）：955－－960.

米寧，王良書，李華，徐鳴潔，李成，張勇，陳運平，于大勇.2005.天山和塔里木盆地結合部地殼上地幔速度結構［J］.科學通報，50（4）：363－－368.

錢輝，許志琴，姜枚，宿和平.2006.西昆侖接收函數(shù)反演與構造解析［J］.中國地質，33（2）：309－－316.

王有學，韓果花，姜枚.2004.阿爾泰—阿爾金地學斷面地殼結構［J］.地球物理學報，47（2）：46－－54.

胥頤，劉建華，劉福田.2006.天山—帕米爾結合帶的地殼速度結構及地震活動研究［J］.地球物理學報，49（6）：1693－－1700.

楊少敏，李杰，王琪.2008.GPS研究天山現(xiàn)今變形與斷層活動［J］.中國科學：D輯，38（7）：872－－880.

楊主恩，張先康，趙瑞武，周偉新.2005.天山中段的深淺構造［J］.地震地質，27（1）：11－－19.

趙俊猛，李植純，馬宗晉.2003.天山分段性的地球物理分析［J］.地學前沿，10（S1）：125－－131.

趙俊猛，程宏崗，裴順平，劉宏兵，張健獅，劉寶豐.2008.塔里木盆地北緣的深部結構［J］.科學通報，53（8）：946－－955.

周仕勇.1999.主地震定位法分析以及1997年新疆伽師強震群高精度定位［J］.地震學報，21（3）：43－－51.

Cormier V F.2001.Construction of 3-D earth models for station specific path corrections by dynamic ray tracing［C］∥Proceedings of the 23rdSeismic Research Review.Tucson，Arizona，USA，NNSA，1：30－－36.

Kennett B L N，Engdahl E R.1991.Travel times for global earthquake location and phase identification［J］.Geophys J Int，105（3）：429－－465.

Kennett B L N，Engdahl E R，Buland R.1995.Constraints on seismic velocities in the Earth from travel times［J］.Geophys J Int，122（1）：108－－124.

Russian Federation／United States Calibration Working Group.2002.Refinement of regional seismic event location in northern Eurasia using 3-D crustal and upper mantle velocity model［C］∥Proceedings of the 24thSeismic Research Review.Tucson，Arizona，USA，NNSA，2：176－－185.

Ryaboy V，Baumgardt D R，F(xiàn)irbas P，Dainty A M.2001.Application of 3-D crustal and upper mantle velocity model of North America for location of regional seismic events［J］.Pure Appl Geophys，158（1）：79－－103.

Schneider W A，Ranzinger K A，Balch A H，Kruse C.1992.A dynamic programming approach to first arrival traveltime computation in media with arbitrarily distributed velocities［J］.Geophysics，57（1）：39－－50.

Steck L K，Rowe C A，Begnaud M L，Phillips W S，Gee V L，Velasco A A.2004.Advancing seismic event location through difference constraints and three-dimensional models［C］∥Proceedings of the 26thSeismic Research Review.Tucson，Arizona，USA，NNSA，3：346－－355.