山東地區(qū)地震波一維速度模型研究

李 鉑，崔 鑫，葉慶東，于 澄

（1.中國科學技術大學，合肥 230026；2.山東省地震局，濟南 250014；3.中國地震局地球物理研究所，北京 100080）

0 引言

地震波速度結構的反演是地震學的經(jīng)典問題，精準的速度模型是許多研究的基礎性資料，如地震定位、大尺度構造應力場研究、地球動力學等。速度結構研究主要可分為一維速度結構和三維速度結構研究，一維速度結構經(jīng)典的研究方法主要有Herglotz－Wiechert方 法、Gutenberg 反 演 方 法、τ（p）法［1］等，主要依據(jù)大量地震資料的走時曲線計算速度結構，目前應用較少。地震波速度結構研究有2個主要的難點，一是地震波速度結構和地震定位的耦合問題，兩者會互相影響其誤差和準確性；二是地震數(shù)據(jù)的結構問題，即地震臺站深度普遍位于近地表，難于深入地下，從而難以在深度上進行有效約束。為了解決地震波速度結構和地震定位的耦合性問題，在Geiger定位方法的基礎上，1976年，Crosson［2］首次提出了震源參數(shù)和速度結構聯(lián)合反演方法，將速度結構作為未知參數(shù)與震源參數(shù)同時反演，由此彌補選取速度模型引起的誤差，得到了廣泛的應用。Aki［3－4］等人將該方法推廣到三維速度結構，Pavlis［5］，Spencer［6］，劉福田［7－8］等人相繼提出了參數(shù)分離法，大大提高了運算效率。Kissling［9］等于1994年提出可以基于走時殘差均方根最小作為目標函數(shù)進行一維速度結構和震源位置聯(lián)合反演，定位精度更高。

基于現(xiàn)在全球數(shù)字地震觀測技術的發(fā)展和觀測臺網(wǎng)的大量建設，積累了大量高質量的地震波走時數(shù)據(jù)和數(shù)字波形記錄，使得三維速度結構的反演方法取得了長足的進步。目前層析成像方法在數(shù)據(jù)處理時分為2組：波速分析和地震偏移。波速分析的目的是反演地球內部介質的波速結構，其主要方法有：體波走時層析成像［10］、面波走時層析成像［11］及接收函數(shù)法［12］。

1 區(qū)域地質構造特征

山東省位于我國地勢劃分中的第三大階梯中［13］，海拔高度除一小部分山地超過千米以外，大部分山地丘陵都在500m 左右，地勢起伏較小。山東省地層在全國地層區(qū)劃中，屬華北地層大區(qū)的晉冀魯豫地層區(qū)和秦祁昆地層區(qū)，前者可分為華北平原、魯西和魯東3個地層分區(qū)，后者為祁連－北秦嶺地層分區(qū)（魯東南地層分區(qū)）。魯西北－魯西南平原區(qū)屬華北平原地層分區(qū)，魯中南山地丘陵區(qū)屬魯西地層分區(qū)，魯東丘陵區(qū)的西北部（魯東北地區(qū)）屬魯東地層分區(qū)，魯東丘陵區(qū)的東南部（魯東南地區(qū)）屬祁連－北秦嶺地層分區(qū)（魯東南地層分區(qū)）。山東省各斷代地層發(fā)育比較齊全，自中太古代至新生代地層都有分布，地表出露以中、新生代地層為主，其次為古生代地層，元古宙地層分布局限，太古宙地層零星出露。華北平原地層分區(qū)以發(fā)育大量新生代地層區(qū)別于魯西地層分區(qū)，魯東地層分區(qū)和祁連－北秦嶺地層分區(qū)以沒有古生代地層區(qū)別于華北平原地層分區(qū)和魯西地層分區(qū)，魯東地層分區(qū)、祁連－北秦嶺地層分區(qū)與魯西地層分區(qū)（包括華北平原地層分區(qū)）三者的前寒武紀基底巖系組成特征明顯不同。

山東省表層構造首先反映在地貌特征上，是由魯中和半島地區(qū)的低山丘陵及環(huán)繞的堆積平原、陸架海域［14］構成的現(xiàn)代地貌格局。沂沭斷裂帶縱貫山東中部，如“刀劈斧砍”將山東一分為二。蘇魯造山帶則奠定了魯東地區(qū)基底構造線的總體格局。齊河－廣饒斷裂和聊城－蘭考斷裂則是分劃魯西地塊和華北坳陷平原的構造帶。因此，山東的地質塊體所反映的地表構造格局具有一坳（濟陽坳陷）、兩塊（魯西地塊、膠北地塊）、兩帶（沂沭斷裂帶、蘇魯造山帶）及一域（黃、渤海陸架海域）6大構造塊體格局。

圖1 山東地區(qū)斷層分布圖

2 理論與算法

地震波的到時是臺站坐標（s）、震源參數(shù)（h，包括發(fā)震時刻和地理坐標）和速度場（m）的非線性函數(shù)，

一般而言，真正的震源參數(shù)和速度場都是未知的。僅有到時和臺站坐標可以測量得到，我們不能直接解出（1）式。所以我們必須對未知參數(shù)做有根據(jù)的推測。使用初始速度模型，我們可以對初始震中到臺站進行射線追蹤，并計算理論到達時間（tcalc）.觀測到時與理論到時的差，即為走時殘差（tres），可以將其視為當前震中和速度模型與真實震中和速度模型的差（Δ）。為了計算震源和模型參數(shù)之間的偏差，我們需要知道觀測走時對所有參數(shù)的相關性。對于（1）式進行泰勒展開，即可獲得基于震中和速度模型的走時殘差。

根據(jù)矩陣理論，震源參數(shù)和速度模型之間的耦合關系寫為如下形式：

式中：t為走時殘差向量；

H為走時對震源參數(shù)的偏導數(shù)矩陣；

h為震源參數(shù)的擾動量向量；

M為走時對速度模型的偏導數(shù)矩陣；

m為速度模型的擾動量向量；

e為走時誤差向量，包括觀測走時拾取誤差以及由于臺站坐標造成的理論走時誤差、不適當?shù)乃俣饶Ｐ秃驼鹪磪?shù)的誤差、方程線性化等的誤差；

A為所有偏導數(shù)矩陣。A的結構為

d為震源參數(shù)和速度模型的擾動量，d＝（Δh1，Δh2，Δh3，Δh4，Δm1，…，Δmn）。

在定位地震時忽略公式（3）中Mm 的作用可能引入震源位置的系統(tǒng)誤差［15－16］。與此類似，在等式中忽略Hh可能導致速度參數(shù)的誤差［17－18］。因此，同時使用震源和模型參數(shù)并不一定能獲得真正的震中和速度模型。

除非我們已經(jīng)“猜”出了正確的震源坐標。地震數(shù)據(jù)和層析成像會更新震源和速度參數(shù)。Thurber［15］在1992 指出，聯(lián)合反演震源－速度模型，要比輪換使用獨立的震中和速度模型擾動可信度更高，即先用速度模型進行地震定位，再采用震源位置修正震源模型，反復迭代得到最終的震中和速度模型。Kissling［19］1988年在加利福尼亞長谷地區(qū)數(shù)據(jù)測試中展示了這種分離變量法在進行地震層析成像中的高效率性。

通常，公式（3）的解是個最小二乘解，該解使得誤差（eTe）的加權組合最小，eTe即為整個算法的目標函數(shù)［20］?；谡龖B(tài)分布誤差和擾動模型假設，最小二乘解能給出符合初始模型的可信度最好的解。通過搜索初始參考模型的鄰域，初始模型的缺陷會導致三維層析成像的人工缺陷，特別是解不通過模型更新和反復迭代優(yōu)化的時候。因為參數(shù)空間常常包含幾千個未知數(shù)，簡單的迭代過程可能不能避免局部最小，從而不能給出全局最小。而且A 的維數(shù)很大，分離變量法也可能會有困難。若有可能，計算模型、數(shù)據(jù)分辨率矩陣和模型協(xié)方差矩陣，能夠對解了解更詳細。否則震中參數(shù)和速度模型的權衡就很難。

由于山東地區(qū)地質構造比較復雜、地殼上的地幔速度結構區(qū)分明顯［21］，所以如果采用簡單的一維速度模型，實際速度模型差異對定位結果肯定會有較大的影響。利用Kissling方法計算得到一維速度模型，包括臺站校正值，應用于地震定位中，很大程度上可消除初始模型與真實速度模型的差異，進而提高地震定位的精度。本文采用Kissling 計算方法，利用Kissling等人編寫的程序Velest來計算山東地區(qū)P波的一維速度模型。

3 觀測資料

2007年山東數(shù)字地震臺網(wǎng)開始運行，截止目前，接入山東地震臺網(wǎng)的地震臺站為63個，其中寬頻帶地震儀39臺，短周期地震儀24臺，臺站分布較均勻，平均臺間距50km。

為了獲得盡可能好的反演效果，我們通過以下標準進行地震事件的選擇：

（1）震級大于ML2.0級；

（2）震中位于山東及其鄰區(qū)，且震中位于陸上；

（3）地震定位精度為1。

同時，為了取得山東局部地質地區(qū)的地震波速度模型，我們針對3 個地區(qū)選取了部分地震事件。主要是：河南范縣和山東鄄城交界地區(qū)，郯廬帶兩側地區(qū)。河南范縣和山東鄄城交界地區(qū)主要是以聊城－蘭考斷裂作為研究區(qū)域，郯廬帶則以昌邑－大店斷裂作為中心選取地震事件。

經(jīng)過篩選 總計有45次地震事件，919條地震射線。震中分布圖如圖2所示，其中紅色實心圓為地震震中，綠色三角形為地震臺站。圖中用顏色標示了山東及其鄰區(qū)的高程信息，對照顏色與高程色塊圖，可以看出山東及其鄰區(qū)大部分都是平原（紫色區(qū)域），中部有山脈分布（綠色區(qū)域）：

圖2 震中與臺站分布圖

4 速度模型確定

Kissling方法基本計算過程為：在參考已有先驗信息的基礎上，選取不同的速度模型作為初始模型，同時反演震源參數(shù)和速度模型，包括臺站校正值，再在初始模型的基礎上對層厚和速度進行微調，經(jīng)過反復計算得到走時均方根殘差最小的一維速度模型作為我們最后求得的一維速度模型［22］。選取不同的速度模型反復計算旨在求得多值解空間中的全局最小速度模型，避免陷入局部最小。

山東及其鄰區(qū)相關的速度模型主要有：趙仲和［23］1983年給出 的MDBJ模型，陳立華［24］1990年給出的華北地區(qū)速度模型（下述簡稱“陳立華模型”），以及IASP91全球速度模型，如圖3所示：

圖3 3種不同的P波速度模型

由于陳立華模型相對其他模型較新，而且使用了更多的地震數(shù)據(jù)，其他模型則或者是全球模型，較為粗糙，又或者研究對象為首都圈地區(qū)，因而山東地區(qū)的地震數(shù)據(jù)結構較差?；谝陨峡紤]，我們最終使用陳立華模型作為初始模型進行反演。由于山東的地震事件深度普遍較淺，我們選取的事件中，最深的地震位于地下10km，這自然影響了地震速度模型的反演，因為我們只能得到地震射線最低點以上的速度結構，對于更深的地方，對速度模型沒有分辨。為此，我們需要求得地震射線的最低點，以便參考確定速度模型的深度下限［25］。經(jīng)計算，地震射線的平均深度為17.1km，我們以此作為參考，取略深于這一深度的速度模型，開始進行迭代計算。斷層層數(shù)為8層，最下面一層為均勻半無限空間，反演過程不改變層厚度，只進行層中波速的計算。計算到同時滿足下列3項條件下迭代中止：①震源位置，臺站校正值和速度值已沒有很大的變化；②所有地震均方根殘差相對于第一次重新定位后的結果明顯地減??；③計算的一維速度模型和臺站校正值有實際的地質意義且沒有違反先驗信息。據(jù)此，我們最終計算出一維P波速度模型［26］，如圖4中的黑色加粗實線，該速度模型從地表到30km 深度處分為5層。從圖4中可以發(fā)現(xiàn)，該模型與陳立華模型在淺層近似程度較高，但是陳立華模型在21～26km 處有低速層，本文模型無法反映該低速層的存在，這可能源于山東地區(qū)的特殊地球物理背景，也有可能是因為該深度范圍的地震樣本數(shù)較少，無法給出較好的速度約束。初始的地震到時數(shù)據(jù)的殘差均方根為1.08s，我們給出的震中和速度模型使殘差均方根降低至0.62s。

圖4 本文得到的P波速度模型

為了研究山東各地區(qū)的速度模型，我們選取了3個地區(qū)分別計算其對應的速度模型，分別為：魯西模型（震中位于河南范縣及其鄰區(qū)），郯廬帶西側模型，郯廬帶東側模型。經(jīng)計算可得到圖5所示的速度模型。

圖5 魯西、郯廬帶P波速度模型對比

如圖5 所示，在地殼淺層，郯廬帶西側速度稍低，郯廬帶東側則相對具有稍高的速度，這可能與區(qū)域的地質背景有關，如郯廬帶西側地區(qū)的沉積層覆蓋較厚，因而速度較低，東側地區(qū)則相對較薄，基巖出露的情況也較多，因而速度較高。

魯西地區(qū)相對于郯廬帶兩側地區(qū)，其淺層速度較低，而在10km 深度左右速度則較高。淺層的低速可能源于地表較厚的沉積層，而10km 深度的速度差異則反映了魯西和魯東地區(qū)花崗巖層巖石的差異。

圖6 原始震中與修訂震中圖

我們再將原始的地震震中位置和與最終速度模型聯(lián)合給出的震中位置進行做圖比較，如圖6所示，紅色圓圈為原始地震震中，橙色圓圈為最終速度模型對應的震中位置，綠色三角形為地震臺站?？梢园l(fā)現(xiàn)，在郯廬斷裂帶兩側區(qū)域，新的震中位置相對于原始震中位置距離斷層更近。

5 討論與結論

（1）通過Kissling提出的計算最小一維速度模型方法，利用山東地震臺網(wǎng)的觀測走時資料，確定了山東及其鄰區(qū)的P 波一維速度模型以及臺站校正值。

（2）采用一維速度模型，對選取的地震重新定位。定位后的地震走時殘差有一定的降低，地震定位精度在經(jīng)度、緯度、深度方向上有了提高。表明了一維速度模型要優(yōu)于其他速度模型。同時走時殘差降低的幅度有限，說明用一維速度結構近似真實的地殼結構仍然過于簡化，應該使用更復雜更合理的速度模型進行聯(lián)合定位，才能期望得到更高精度的結果。

（3）地震深度和震相的選擇限制了速度結構的深度。由于山東地區(qū)地震事件普遍位于近地表，深度普遍小于10km，而我們選取的是直達P震相，其射線的最低點往往不超過20km，因而無法獲取更深處的地震波速度信息。在這個意義上，PmP 與Pn震相具有重要意義，因為這2個震相攜帶了地殼深部、莫霍面以及上地幔頂部的相關信息。未來的研究工作中應考慮加入這2個震相的到時數(shù)據(jù)。

［1］ Geiger，L.（1912）.Probability method for the determination of earthquake epicenters from arrival time only.Bull.St.Louis.Univ.，8，60－71.

［2］ Crosson，R.S.（1976）.Crustal structure modeling of earthquake data，1，Simultaneous least squares estimation of hypocenter and velocity parameters.Geophys.Res.，81（17）：3036－3046.

［3］ Aki，K.，Lee，W.H.K.（1976）.Determination of three－dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes，part 1：A homogeneous initial model.Geophys.Res.，81（23）：4381－4399.

［4］ Aki，K.，et al.（1977）.Determination of the three－dimensional seismic structure of the lithosphere.Geophys.Res.，82（2）：277－296.

［5］ Pavlis，G.，Booker，J.R.（1980）.The mixed discrete－continuous inverse problem：Application of the simultaneous determination of earthquake hypocenters and velocity structure.Geophys.Res.，85（B9）：4801－4810.

［6］ Spencer，C.，Gubbins，D.（1980）.Travel－time inversion for simultaneous earthquake location and velocity structure determination in laterally varying media.Geophys.J.Roy.Astr.Soc.，63（1）：95－116.

［7］ 劉福田.震源位置和速度結構的聯(lián)合反演（I）——理論和方法［J］.地球物理學報，1984，27（2）：167－175.

［8］ 李強，劉福田.一種橫向不均勻介質中地震基本參數(shù)的測定方法［J］.中國地震，1991，7（3）：54－63.

［9］ Kissling E，et al（1994）.Initial reference models in local earthquake tomography.J.Geophys.Res.，99，19 635－19 646.

［10］ 雷建設.淺析中國地區(qū)體波走時層析成像研究進展［J］.西北地震學報，2000，22（4）：471－475.

［11］ 黃忠賢.中國東部海域巖石圈結構面波層析成像［J］.地球物理學報，2009，52（3）：653－662.

［12］ 黃海波.利用遠震接收函數(shù)方法研究南海西沙群島下方地殼結構［J］.地球物理學報，2011，54（11）：2788－2798.

［13］ 蔣科植，朱介壽，程先瓊.利用面波方法研究歐亞大陸及周邊海域巖石圈結構［J］.華北地震科學，2012，30（2）：1－7.

［14］ 王帥軍，王夫運，張成科，等.北京及附近地區(qū)地殼速度結構與構造［J］.華北地震科學，2011，29（2）：6－12.

［15］ Thurber，C.H.（1992）.Hypocenter－velocity structure coupling in local earthquake tomography.Phys.EarthPlanet.Inter.，75：55－62.

［16］ Eberhart－Phillips，D.，and A.J.Michael（1993）.Three－dimensional velocity structure，seismicity，and fault structure in the Parkfield region，central California.Geophys.Res.，98，737－758.

［17］ Michael，A.J.（1988）.Effects of three－dimensional velocity structure on the seismicity of the 1984 Morgan Hill，California，aftershock sequence.Bull.Seismol.Soc.Am.，78：1199－1221.

［18］ Vander Hilst，R.D.，and W.Spakman（1989）.Importance of the reference model in linearized tomography and images of subduction below the Caribbean plate.Geophys.Res.Lett.，16，1093－1096.

［19］ Kissling，E.，et a1（1995）.Improved seismic velocity reference model from local earthquake data in Northwestern Italy.Terra.Nova.，7，528－534.

［20］ Spakman，W.，and G.Nolet（1988）.Imaging algorithms，accuracy and resolution in delay time tomography in Mathematical Geophysics.In Vlaar N.J.，et al.（Eds.），Mathematicalgeophysics：asurveyofrecentdevelopmentsinseismologyandgeodynamics（pp.155－187）.Dordrecht：Reidel Publ.Co..

［21］ 于利民，刁桂苓，李欽祖.由深源遠震體波記錄反演華北北部地殼上地幔速度結構［J］.華北地震科學，1983，（3）：11－20.

［22］ 傅淑芳，劉寶誠.地震學教程［M］.北京：地震出版社，1991：447－480.

［23］ 趙仲和.北京地區(qū)地震參數(shù)與速度結構的聯(lián)合測定［J］.地球物理學報，1983，26（2）：131－139.

［24］ 陳立華，宋仲和.華北地區(qū)地殼上地幔P波速度結構［J］.地球物理學報，1990，33（5）：540－547.

［25］ 萬永革，李鴻吉.京津唐張地區(qū)速度結構和震源位置聯(lián)合反演的遺傳算法［J］.地震地磁觀測與研究，1996，（3）：57－66.

［26］ 王小龍，劉淵源，余國政，等.利用遠震接收函數(shù)反演烏江彭水電站地震臺下方地殼厚度［J］.華北地震科學，2011，29（2）：13－18.