蘆山余震震源參數及震源區品質因子反演

溫瑞智，王宏偉，任葉飛，冀 昆

(中國地震局工程力學研究所，150080哈爾濱)

2013年4月20日四川省蘆山縣發生MS7.0級地震，蘆山主震震中位于龍門山斷裂帶西南端的彭縣-灌縣斷裂上，主震過后余震活動十分頻繁，余震沿發震斷層向主震兩側延伸，主要分布在長約32 km、寬約15～20 km、深度為5～24 km的范圍內［1］，余震破裂類型以逆沖型為主［2］.蘆山地震發震區域位于龍門山斷裂及鮮水河-小江斷裂的“Y”形交界處，在過去的40年里，距離蘆山地震震中200 km范圍內發生過多次MS6.0級以上地震，其中包括2008年MS8.0級汶川地震及其余震，這一區域地震多發人口密集，地震危險性很高.

我國測震資料相對豐富，國內研究人員大多基于Atkinson方法與遺傳算法采用測震記錄反演我國地震多發地區的品質因子及震源參數，比如甘肅、云南和四川地區［3-5］.隨著我國高質量強震動記錄的增加，采用強震動記錄反演中小地震震源參數及品質因子已經逐漸得到推廣［6-8］.我國數字強震動臺網在龍門山地區密度相對較大，蘆山余震序列中臺站多次觸發，記錄到2.0～5.4級余震176次，共收集到超過1 000組加速度記錄，這為反演蘆山余震震源參數及品質因子提供了基礎數據.本文基于參考事件的廣義反演方法估計了發震區域剪切波品質因子及震源參數，包括地震矩、拐角頻率、應力降、破裂半徑及地震波能量，進一步確定了震源參數的定標關系.

1 廣義反演方法

自由表面水平地震動剪切波加速度傅里葉幅值譜可表示為

式中:Oij(f)表示第j個臺站觀測到第i個地震的記錄剪切波水平方向加速度傅里葉幅值譜;Si(f)表示第i個地震的加速度震源譜;GS(Rij)表示幾何擴散;Rij表示第i個地震到第j個臺站的震源距;exp(-πfRij/Q(f)Vs)表示地震波非彈性衰減;Q(f)表示與頻率f相關的品質因子;Vs表示震源處剪切波速，取為3.6 km/s;Gj(f)表示第j個臺站的場地反應.

本文幾何擴散采用Atkinson等［9］給出的三段線性衰減曲線，表示為

式中:R01=1.5D，R02=2.5D，b1=1.0，b2=0.5，D表示地殼厚度.蘆山地震序列震中附近區域地殼結構復雜，地殼厚度從西北側的52.5 km減小至東南側的 41.5 km［1］，近似以平均地殼厚度47 km作為地殼厚度，則R01=70.5 km，R02=117.5 km.

式(1)兩邊同時取對數，得到線性疊加形式:

式(3)矩陣形式為

式中:A是每一行包含3個非零項(兩個1和-πfRij/Vs)的稀疏矩陣;X表示式(3)右邊所有未知量的向量;b表示式(3)左邊所有已知量的向量.該矩陣方程可詳細表示為式(5)所示形式.采用奇異值分解方法求解式(4)，在每個頻率上確定Ⅰ(地震個數)+J(臺站個數)+1(Q值)個未知數，由于存在一個未加約束的自由度，需要考慮震源與場地間的權衡問題，通常采用參考場地法或參考事件法來解決.蘆山余震序列中6個觸發基巖臺站水平/垂直譜比法估計的場地反應均不能滿足參考場地應在整個頻段上不存在場地放大的要求，因此本文采用參考事件法.

2 數據及數據處理

采用以下原則選取合適的強震記錄用于廣義反演:

1)記錄震源距范圍為25～100 km.蘆山余震震源深度主要分布于5～24 km，為避免震源深度不確定性對震源距計算的影響和減少截取的剪切波中面波引起的干擾分別設定最小和最大震源距限制［10］.

2)對原始數據進行零線校正及0.1～30 Hz巴特沃斯帶通濾波處理后，為減小噪聲對記錄的影響和避免場地出現非線性反應選取三分量峰值地面加速度均滿足5～100 cm/s2的記錄［11］.

3)為減小反演結果的離散性，在滿足1)和2)的基礎上，選擇同時滿足不少于4次地震中觸發的臺站以及不少于4個臺站記錄到的地震所對應的記錄.

根據上述3個原則選取了25個臺站在34次地震中記錄到的262組記錄，臺站記錄震源距分布見圖1，圖中陰影區域表示5 km間隔范圍內記錄數量分布直方圖圍成的區域，64表示最大頻數.由圖1可見，縱坐標51DXY及以上臺站(A組)，記錄震源距大多不超過70 km，主要分布于25～45 km范圍內;縱坐標 51DXY以下臺站(B組)，震源距主要分布于80～90 km范圍內.A、B兩組臺站記錄震源距差異決定了傳播路徑不同，由于品質因子對地殼介質橫向和垂直方向的不均勻性特別敏感，A、B兩組臺站記錄品質因子差異較大［12］.

圖1 臺站記錄震源距分布

本文在A組臺站數據中選取了30次地震中15個臺站記錄到的179組記錄用于廣義反演.臺站及余震震中分布見圖2.所有傳播路徑分布均勻且交織在一起，方位角對震源影響及傳播介質不均勻性對品質因子影響可忽略.A組記錄震源距不超過R01(即70.5 km)，幾何擴散GS(R)取為R-1.

圖2 臺站、地震震中位置及震源與臺站間傳播路徑分布

廣義反演方法需截取記錄的剪切波部分，剪切波抵達時間和結束時間分別定義為Husid函數中地震波能量開始急劇增加的點和累加均方根函數開始下降的點對應的時間［13］.在截取的剪切波前后各乘以10%的剪切波持時的余弦邊瓣窗口以消除截斷誤差，計算剪切波兩個水平方向的傅里葉譜，采用b=20的Konno和 Ohmachi窗口函數平滑傅氏譜［14］，并以兩個水平方向譜的矢量和作為記錄水平傅氏譜.

3 參考事件震源譜

Ren等［11］以62WIX臺站作為參考場地基于廣義反演方法估計了96個汶川余震加速度震源譜.本文選取2013年4月20日09:37:29發生的4.9級地震作為參考事件，從96個汶川余震中選取與參考事件震級相差不超過 0.1(4.9±0.1)的余震13個.采用下文4.1節提到的網格搜索方法確定選取的13個余震的平均位移震源譜的理論震源參數，地震矩M0=2.427×1023dyne·cm、拐角頻率fc=0.89 Hz、高頻滾降系數 γ=2.0，以該理論震源譜作為參考事件的實際震源譜用于廣義反演.

4 結果分析及討論

4.1 震源參數

基于參考事件的廣義反演方法分離了上文所選15個臺站的場地反應、30個地震的加速度震源譜以及區域品質因子.反演得到的余震位移震源譜見圖3，圖中可見，位移震源譜較符合Aki［15］提出的ω2震源譜模型.采用震源參數(M0、fc)及γ表示的實際位移震源譜的理論形式:

圖3 反演得到的30次余震位移震源譜

式中:RΘΦ為點源輻射圖型因子，隨方位角呈對稱性變化，在此取平均值0.55;V表示地震波水平分量含剪切波能量的比例，取為0.707;F為半空間表面放大，本文取為1;ρs為震源處介質密度，取為2 700 kg/m3;R0為參考距離，通常取為1 km.

其中矩震級MW與地震矩有如下關系［16］:

采用網格搜索方法在0.5～20 Hz頻段內搜索滿足理論位移震源譜與實際位移震源譜相對面積差最小的MW、fc及γ，即

式中:n表示0.5～20 Hz內所有頻率點數;fi表示第i個點的頻率;SS(fi)與SO(fi)分別表示搜索的理論位移震源譜及基于廣義反演方法確定的實際位移震源譜.理論震源譜的搜索過程中MW變化范圍為余震震級±0.5，變化步長0.01;fc變化范圍為0.01～5.0 Hz，變化步長為 0.01 Hz;γ 變化范圍為2±0.3，變化步長 0.1.

最小二乘擬合地震矩與拐角頻率可得:

Aki［15］認為M0fc3是一個與應力降有關的常數，故將式(9)斜率固定為-3.00后擬合得:

擬合結果見圖4，M0fc3=1.73×1023dyne·cm·s-3，略低于Dutta等［17］研究阿拉斯加中南部地區中小地震得到的M0fc3=2.09×1023dyne·cm·s-3和Hassani等［13］研究伊朗中東部地區中小地震得到的M0fc3=2.48×1023dyne·cm·s-3.蘆山余震平均應力降為 3.14 MPa，略低于喻畑等［8］研究13次汶川余震 (MW＞5.0)震源參數得到的3.80 MPa平均等效應力降.

圖4 地震矩與拐角頻率關系(陰影區域表示加減一倍標準差范圍)

根據Brune［18］提出的圓盤形應力脈沖震源模式，計算了余震震源半徑r及應力降Δσ.地震矩與震源半徑r的關系見圖5，圖中同時給出了余震的應力降分布范圍(1.0～10.0 MPa).應力降與震級沒有明顯的相關關系，與Moya等［19］研究阪神地震余震得到的結論一致.

圖5 地震矩與震源半徑散點分布圖

基于Izutani-kanamori理論譜方法同時考慮地震波能量補償將積分上限提高至最大拐角頻率(拐角頻率上限為 5 Hz)的 10 倍，即 50 Hz［20］.本文分別計算了地震波能量Es和視應力σA.

地震波能量與地震矩的關系見圖6，雙對數坐標下線性擬合得到 logEs=0.82logM0-7.7，相關系數R=0.91，擬合直線斜率0.82 接近1.0，近似有Es∝M0，固定斜率為1擬合得

Es/M0平均值為2.19×10-12，與 Hassani等［13］得到的Es/M0=2.50×10-12較為接近，但明顯大于Dutta 等［17］得到的Es/M0=1.20×10-12.

圖6 地震波能量與地震矩關系

本文確定的視應力為 0.17～2.21 MPa，低于程萬正等［21］研究2000—2004年四川地區中小地震得到的視應力變化范圍 0.1～10.0 MPa，推斷可能是發生于主震破裂面的余震初始應力降低所致.根據應力降與視應力的計算結果，線性擬合得Δσ=6.78σA-1.05，R=0.91，表明應力降越高視應力越大，單位地震矩輻射的地震波能量越大，本文余震σA/Δσ的平均值為0.258.

不考慮面波震級MS和地方震級ML之間的換算，通過最小二乘法擬合可得到地震面波震級與矩震級的關系:MW=0.876 2MS+0.563 9(R=0.83)，如圖7所示，本文確定的矩震級與面波震級關系與Ren等［11］反演汶川余震得到的結果較為一致，整體上矩震級擬合結果與面波震級一致，對于小震級余震，矩震級略高于面波震級;對于中等震級余震，矩震級略低于面波震級.

圖7 矩震級與面波震級關系

4.2 品質因子

品質因子Q(f)通常表示為Q0fn的形式，本文反演確定的蘆山余震區域剪切波品質因子見圖8，擬合得到0.5～20 Hz頻段品質因子為Q(f)=31.867f1.0375，蘆山地震近場區地震波衰減表現出吸收快且與頻率依賴性強的特點(低Q0高n).蘆山發震區與汶川地區臨近，其頻率相關的品質因子對頻率的依賴性相近，即n值較一致，一定程度上可推斷本文反演的品質因子較為可靠.與四川其他地區相比［4，5，11］，本文反演確定的蘆山余震區域品質因子結果較小，主要原因有兩方面:1)廣義反演方法確定的品質因子反映了所有傳播路徑的平均非彈性衰減，本文選取的余震震源深度主要分布于10～25 km范圍內，該研究區域地殼平均厚度約為47 km，震源距集中在25～45 km范圍，地震波傳播路徑主要集中在上地殼，而其他研究成果采用較遠震源距的記錄，地震波傳播路徑主要集中在下地殼，眾所周知地殼深處Q值相對較大;2)本文研究區域集中于蘆山地震發震區，該區域地殼活動強烈，一般認為，在地殼活動相對穩定地區Q值相對較高，在地殼活動強烈的地區Q值相對較小.

圖8 頻率相關的品質因子

5 結論

1)采用蘆山強余震記錄基于廣義反演方法估了余震震源參數 (M0、fc、r、Δσ、Es、σA)及該區域剪切波品質因子Q(f)，并給出了震源參數間的定標關系.

2)蘆山余震M0fc3=1.73×1023dyne·cm·s-3，相應于3.14 MPa的平均應力降，略低于汶川地震余震平均等效應力降;應力降主要在0.1～10.0 MPa范圍內變化，與震級沒有明顯的相關性;視應力與應力降正相關，由蘆山余震視應力低于2000—2004年四川余震結果推斷發生于主震斷層面上的余震初始應力可能較低.蘆山余震震源參數M0fc

3及Es/M0與 Hassani等［13］研究伊朗中東部中小震級地震的結果十分接近，初步推斷兩地區余震震源破裂過程以及區域構造環境可能存在相似性.

3)0.5～20 Hz頻段內本文估計的剪切波品質因子Q(f)=31.867f1.0375，分析表明地震波衰減吸收較快且與頻率的依賴性較強;本文獲得的蘆山地震區域品質因子低于四川其他地區的結果，主要是由于所用記錄震源距相對較小且該區域地殼活動較為活躍.另外，需要說明的是，本文研究區域位于青藏高原與四川盆地過渡地帶，存在地殼介質橫向不均勻性.廣義反演確定的品質因子僅體現該區域的平均水平，該方法無法考慮單一區域內品質因子的橫向變化性.考慮到本文研究區域相對較小(半徑約70 km)，暫不考慮地殼介質橫向不均勻性對本文反演結果的影響.

［1］陳晨，胥頤.蘆山MS7.0級地震余震序列重新定位及構造意義［J］.地球物理學報，2013，56(12):4028-4036.

［2］林向東，葛洪魁，徐平，等.近場全波形反演:蘆山7.0級地震及余震矩張量解［J］.地球物理學報，2013，56(12):4037-4047.

［3］華衛，陳章立，鄭斯華，等.三峽水庫地區震源參數特征研究［J］.地震地質，2010，32(4):533-542.

［4］喬慧珍，張永久，程萬正.川西北地區介質衰減特性研究［J］.地震地磁觀測與研究，2006，27(4):1-7.

［5］張永久，喬慧珍，程萬正.四川盆地地區介質衰減特性研究［J］.地震研究，2007，30(1):43-48.

［6］章文波，謝禮立，郭明珠.利用強震記錄分析場地的地震反應［J］.地震學報，2001，23(6):604-614.

［7］劉杰，鄭斯華，黃玉龍.利用遺傳算法反演非彈性衰減系數、震源參數和場地響應［J］.地震學報，2003，25(2):211-218.

［8］喻畑，李小軍.汶川地震余震震源參數及地震動衰減與場地影響參數反演分析［J］.地震學報，2012，34(5):621-632.

［9］ATKINSON G M，MEREU R F.The shape of ground motion attenuation curves in southeastern Canada［J］.Bulletin of the Seismological Society of America，1992，82(5):2014-2031.

［10］MCNAMARA D，MEREMONTE M，MAHARREY J Z，et al.Frequency-dependent seismic attenuation within the Hispaniola island region of the Caribbean sea ［J］.Bulletin of the Seismological Society of America，2012，102(2):773-782.

［11］REN Yefei，WEN Ruizhi，HIROAKI Yamanaka，et al.Site effects by generalized inversion technique using strong motion recordings ofthe 2008 Wenchuan earthuake［J］.Earthquake Engineering and Engineering Vibration，2013，12(2):165-184.

［12］PETUKHIN A，IRIKURA K.A method for the separation of source and site effects and the apparent Q structure from strong motion data ［J］.Geophysical Research Letters，2000，27(20):3429-3432.

［13］HASSANI B，ZAFARANI H，FARJOODI J，et al.Estimation of site amplification，attenuation and source spectra of S-waves in East-Central Iran ［J］.Soil Dynamic and Earthquake Engineering，2011，31(2011):1397-1413.

［14］KONNO K，OHMACHI T.Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components ofmicrotremor ［J］.Bulletin of the Seismological Society of America，1998，88(1):228-241.

［15］AKI K.Scaling law of seismic spectrum［J］.Journal of Geophysical Research，1967，72(4):1217-1231.

［16］HANKS T C，KANAMORI H.A moment magnitude scale［J］.Journal of Geophysical Research，1979，84(B5):2348-2350.

［17］DUTTA U，BISWAS N，MARTIROSYAN A，et al.Estimation of earthquake source parameters and site response in Anchorage，Alaska from strong motion network data using generalized inversion method ［J］.Physics of the Earth and Planetary Interiors，2003，137:13-29.

［18］BRUNE J N.Tectonic stress and the spectra of seismic shear-waves from earthquake ［J］. Journal of Geophysical Research，1970，75(26):4997-5009.

［19］MOYA A，IRIKURA K.Estimation of site effects and Q factor using a reference event［J］.Bulletin of the Seismological Society of America，2003，93(4):1730-1745.

［20］華衛，陳章立，鄭斯華.2010年4月14日青海玉樹7.1級地震序列中小地震輻射能量的估計［J］.地球物理學進展，2012，27(1):8-17.

［21］程萬正，陳學忠，喬慧珍.四川地震輻射能量和視應力的研究［J］.地球物理學進展，2006，21(3):692-699.