Lg波衰減模型中建模誤差的定量評估＊

胡巖松陳友麟劉瑞豐，劉偉

1） 中國北京 100081 中國地震局地球物理研究所

2 ） 美國馬里蘭州格林貝爾特 20770 Array信息技術公司先進技術事業(yè)部

3） 中國西安 710054 中國地震局第二監(jiān)測中心

引言

地震Lg波對地殼的橫向變化非常敏感，在穩(wěn)定的地殼區(qū)域，如克拉通和地盾等，Lg波可以傳播至3 000 km以外.但當Lg波穿過大陸地殼較薄區(qū)、地殼厚度突變區(qū)、巨厚的沉積盆地或存在部分熔融物質的地區(qū)時，Lg波會呈現出明顯的衰減差異（Nuttli，1980；Kennett，1986；Mitchell，Hwang，1987；Bostock，Kennett，1990；Mitchellet al，1997；Rodgerset al，1997；Schurret al，2003；Xieet al，2004）.例如在青藏高原北部的強烈衰減區(qū)，Lg 波的傳播不超過 700 km （McNamaraet al，1996；Xie，2002）.地震波振幅衰減的多少可以用地下介質品質因子Q值來衡量，它是了解地下裂隙的數量、孔隙密度與分布，孔隙中流體含量以及地殼溫度的重要參數.例如青藏高原的Lg波Q值極低（約為100——200），Chen和Xie （2017）認為這可能是青藏高原地殼內的高溫造成的.由于構造因素導致的Q值變化比速度更加敏感，因此借助Lg波對Q值進行研究可以為了解地下介質提供更多的約束和解釋（Akinciet al，1995；Rietbrock，2001；Gung，Romanowicz，2004）.

Q值層析成像是一個非常復雜的問題，因為Q值反演所使用的地震波振幅數據不僅受有限的Q值的影響，而且還受三維速度結構中波前擴展所導致的未知的三維幾何擴散的影響（Nolet，1987）.在反演Q值時，一般先用簡化過的一維幾何擴散修正觀測到的地面運動振幅，然后再取對數，將其結果作為Q值反演的輸入數據，所得到的輸入數據稱為“折衷振幅數據”（Menkeet al，2006；Chen，Xie，2017）.即使通過一維幾何擴散修正，目前未知的三維幾何擴散影響依舊保留在輸入數據中.Chen和Xie （2017）將這些殘留的三維幾何擴散的影響當作數據誤差來處理并稱其為建模誤差，相較于振幅測量不準導致的誤差，建模誤差才是Q值數據的主要誤差來源（Xie，1998）.本研究使用雙臺法測量Q值，并將兩個臺站相對地震的方位角差也引入了建模誤差.

在許多地球物理反演研究中，直接運用中心極限定理隱含地假設輸入數據的誤差呈正態(tài)分布，反演問題就能簡化至用最小二乘算法求解.中心極限定理指出：如果觀測數據的誤差來源很多，無論某一種來源的誤差呈何種形態(tài)分布，所有來源誤差的綜合效應可以假設呈正態(tài)分布（盛驟等，2008）.如果一組數據在對數域呈正態(tài)分布，那么數據本身就可能不是正態(tài)分布，所以在使用振幅數據的Q值反演中，對于數據的正態(tài)性檢驗不可忽略.Chen和Xie（2017）首次開展了對“折衷振幅數據”建模誤差的統(tǒng)計分析工作，使用中國東南地區(qū)的Lg波，得出“折衷振幅數據”的建模誤差可以用正態(tài)分布描述的結論，此結論是僅適用于當前研究區(qū)或所使用的數據，還是一個具有普適性的結論，尚待商榷.本文將使用在川滇及其周邊地區(qū)收集到的更為密集的地震數據，對Lg波的降幅數據誤差進行更詳盡的統(tǒng)計分析工作.

另一方面，地球物理層析成像方法主要包括同時迭代重構法（simultaneous iterative reconstruction technique）、反投影法（back-projection）和最小平方 QR 分解法（least-squares QR）等.當Q值層析成像所使用的振幅數據符合最小二乘算法的要求時，這些算法就能夠被廣泛使用在Q值層析成像中（Peiet al，2006；Xieet al，2006；Hearnet al，2008； Pasyanoset al，2009；Zhaoet al，2010，2013；Zhao，Xie，2016）.這些方法都是利用一個簡單的初始Q值模型，使用交互式的方式改進模型以使數據得到更好的擬合，從而使反演快速收斂，得到令數據殘差最小化的Q值最終模型，然后再通過“棋盤格”測試（G?ket al，2003；Zhaoet al，2010，2013；Lianget al，2014）、合成“點擴散函數”（Daiet al，2020）或者反映射線路徑的采樣密度的區(qū)域（Zhao ，Xie，2016）等各種分辨率測試方法來找出這些模型的近似分辨率和近似誤差.需要指出的是，這些方法并不能定量地給出模型的分辨率和誤差值，但是在實際應用這些模型時，定量的模型分辨率和誤差是非常重要的.例如，對模型的分辨率和誤差定量分析的結果有助于推斷出所得到的Q模型的可靠部分，從而通過可靠的Q模型恢復特征來推斷地殼屬性.盡管理論上可以量化Q模型的分辨率和誤差，并檢查它們之間的影響（Backus，Gilbert，1970），但實際應用中它們會受到計算能力的限制.大規(guī)模矩陣奇異值（singular value decomposition，縮寫為SVD）分解的數值計算需要使用迭代方法，但由于有限的計算精度，所得到的奇異矢量在迭代過程中會逐漸失去正交性，從而使最終的SVD分解是不完整的，即并非所有的奇異矢量都還保持正交.如果在迭代過程中保持所有奇異矢量的正交性，則需要龐大的計算量和內存.Chen和Xie （2017）開發(fā)了一種基于PROPACK算法的新的2D層析成像方法，該算法的優(yōu)點在于它在迭代過程中只針對失去正交性的奇異矢量做調整，由于失去正交性的只是部分奇異矢量，因此計算量和內存需求大量減少.而正因為PROPACK算法得到的所有奇異矢量都保持了正交性，所以它們不僅可以用來計算模型參數，還可以計算模型的分辨率矩陣和協方差矩陣（Larsen，1998）.

本文擬利用雙臺法得到川滇及其周邊地區(qū)地震數據的”折衷振幅數據”，并對其建模誤差統(tǒng)計樣本進行了提取和分析，采用 K-S （kolmogorov-Smirnov），Q-Q （Quartile-Quartile）圖和正態(tài)分布圖等三種統(tǒng)計檢驗方法得到研究區(qū)建模誤差的統(tǒng)計特征，以期驗證Q值反演數據中建模誤差統(tǒng)計特征的普適性.在證明建模誤差近似服從正態(tài)分布的基礎上，使用最小二乘法求解Q值層析成像問題，并采用基于PROPACK 算法的 2D 層析成像方法得到模型的參數、分辨率矩陣與協方差矩陣，以期得到定量的模型分辨率和誤差.最后，采用不同的層析成像方法來求解Q值模型，通過所得到模型結果的異同來檢驗新層析成像方法的有效性和可靠性.

1 數據和方法

1.1 數據篩選

川滇地區(qū)位于華南地塊、印度板塊和青藏高原的交會處（95°E——108°E，22°N——35°N），區(qū)域內深大斷裂縱橫交錯，地震活動強烈.川滇地區(qū)也是印度板塊和歐亞板塊碰撞帶東南緣的一部分及特提斯——喜馬拉雅造山系的轉折點（鄧起東等，2002），川滇地區(qū)的主要斷裂帶和地塊分布如圖1所示.本研究收集了來自國家數字測震臺網數據備份中心（鄭秀芬等，2009）提供的2007年7月至2013年5月期間262個固定臺站（圖1b）記錄到的地震波形數據.臺站分布在川滇及其周邊地區(qū)（圖1a），挑選了震中距范圍在200——3 000 km、震級在M3.5——6.5之間且具有清晰的P波初至和Lg波波包發(fā)育的地震事件資料用于Lg波Q值的計算.通過手動拾取震相后，再進行去均值、去趨勢、濾波等數據處理，選取133次地震的24 001條地震波垂直分量波形數據用來計算Lg波頻譜.Wei和Zhao （2019）在川滇地區(qū)分別對頻率為0.5，1.0和1.5 Hz的Lg波進行反演，得出在該區(qū)域內QLg沒有很強的頻率依賴性.所以在本研究中，我們只分析了1 Hz附近的Lg譜幅值數據，使用一個中心頻率為1.0 Hz的窄帶寬（0.5——1.5 Hz）來獲得平均譜.

圖1 研究所用地震事件震中（a）和研究區(qū)臺站（b）分布圖圖中紅色矩形框代表研究區(qū)域和臺站所在區(qū)域；黑線為主要地塊分界線，灰線為主要斷裂帶，下同Fig.1 Distribution of epicenters of events （a） and seismic stations （b） in this studyRed rectangular represents the extent of the study area；the black solid lines indicate the major tectonic block boundaries and the grey dashed lines indicate faults，the same below

1.2 Q值的測量方法

在特定頻率f處觀測到的Lg波振幅可以通過正演建模得到，即：

式中，A（f）表示振幅項，S（f）為震源譜（Serenoet al，1988），G（ ?）是幾何擴散項， ? 表示震中距，v表示Lg波的群速度，Q（f）為衰減值，R（f）為臺站的場地響應項（Streetet al，1975；Xie，1993；Yang，2002；Chen，Xie，2017）.真實的幾何擴散與3D地殼結構有關，但一維幾何擴散可以簡化為震中距 ? 的函數（Streetet al，1975），表示為

對于Lg波來說，一維幾何擴散的冪指數取0.5，臨界距離 ?0設為100 km.Q值是通過指數函數影響振幅的，經過1D幾何擴散校正并取對數后得到

式（3）等號左邊是反演數據，稱為“折衷振幅數據”.本研究使用雙臺法 （two station based method，縮寫為TS），將式（1）應用于與震源位于同一大圓弧的兩個臺站，因為這兩個臺站記錄了同一次地震，而路徑衰減的差別主要體現在兩個臺站之間，所以兩個臺站的振幅比直接消除了震源譜，并求出兩個臺站之間的Q值（Xie，Mitchell，1990；Xieet al，2006；Zoret al，2007；Gallegoset al，2014；Chen，Xie，2017）.但傳統(tǒng)雙臺法不能消除場地效應的影響，只是假設兩個臺站的場地效應相同.為了避免必須進行Q值和震源參數之間的權衡所導致的Q 模型的非唯一性 （Menkeet al，2006），研究中又使用了逆雙臺法（reversed two-station method，縮寫為 RTS）（Chunet al，1987）和逆雙事件法（reversed two-event method，縮寫為 RTE）（Chen，Xie，2017）.逆雙臺法是在雙臺法的基礎上，從雙臺路徑的對拓方向再引進一次地震，兩個臺站記錄的兩個地震的交叉振幅比，可以將兩個臺站的場地效應也直接消除，從而得到兩個臺站之間的Q值.逆雙事件法是逆雙臺法的拓展，即將逆雙臺法中臺站和地震的位置互換就是逆雙事件法，它可以測量兩個地震事件之間的Q值.

2 數據誤差的統(tǒng)計分析

2.1 建模誤差的統(tǒng)計方法

雙臺法的Q值計算表示為

式中，Q（ef）為進行一維幾何擴散校正后包含誤差的Q（f）的估計值，為兩個臺站之間的有效距離，為雙臺法的振幅項與幾何擴散項乘積的比值.本研究中將建模誤差的定義范圍進行了拓展，主要包括：① 殘留的三維幾何擴散導致的誤差；② 不能完全去除臺站響應部分造成的誤差；③ 由Sn尾波引起的“噪聲”和地震事件前無法模擬的背景噪聲造成的誤差；④ 兩個臺站和地震的方位角差別導致的誤差.

Q（ef）中包含的誤差可以表述為（Chen，Xie，2017）：

2.2 建模誤差樣本采樣

圖2 基于雙臺法的各個路徑上Q值重復測量次數分布（a） TS 方法；（b） RTS 方法；（C） RTE 方法Fig.2 Distribution of the number of repeated two-stations based method over individual paths（a） TS method；（b） RTS method；（c） RTE method

2.3 統(tǒng)計學方法檢驗結果

原始的 εTS， εRTS和 εRTE的樣本分布具有中心峰值的特征，是顯著的正態(tài)分布的特征，如圖3左側所示，然而遠離峰值的分布在正負方向上都有較長的尾部，這些分布與Nolet （1987）討論的遠震P波走時殘差分布非常相似； ε的特征也符合Chen和Xie （2017）的統(tǒng)計特征結果.

在對 ε的樣本進行了幾次篩選試驗后，確定了任何 ε偏離中位數超過兩個 ε的標準差的樣本都被剔除的篩選準則，即這也符合統(tǒng)計學上的“68-95-99.7”準則（Pukelsheim，1994），68.27% （n＝1），95.45% （n＝2）和 99.73% （n＝3）的值是分別距離樣本數據的平均值小于一個、兩個和三個標準差內的百分比，超過三個標準差的數據可能是異常值，并且如果有許多數據超過三個標準差，則有理由質疑假定的正態(tài)性（Wheelan，2013）.兩個標準差范圍代表著95%的置信水平，可以同時平衡樣本選擇的精度和準確度.如表1所示，經過數據篩選后，對于所有的樣本 ε ，剔除的數據量不超過總體的12%.方差Va（rεTS），Va（rεRTS）和Va（rεRTE）分 別 為0.025 1，0.006 7和0.011 0， εTS的 方 差 遠 大 于 εRTS和εRTE的方差，說明 εTS的 分布比 εRTS和 εRTE更廣， V a（rεRTS）和 V a（rεRTE）之間的微小差異很可能是由每個方法收集的重復測量樣本數不同造成的， V a（rεTS）的值最大可以歸因于TS方法在測量中不能消除場地效應的影響.

表1 使用 TS，RTS 和 RTE 方法估計 ε的詳細信息Table 1 Estimated ε information using TS， RTS and RTE methods

K-S （Kolmogorov-Smirnov）檢驗是一種非參數測試，測試分布的正態(tài)性時，將已知的正態(tài)概率分布與數據生成的分布（經驗分布函數）進行量化比對（Massey，1951；Vrbik，2020），研究中的零假設為數據是正態(tài)分布.本研究采取單樣本K-S檢驗，p值為否定零假設的最低顯著性水平（Vrbik，2020）.K-S檢驗通過p值與顯著性水平α進行比較，從而來決定是否拒絕零假設.如表1所示，p值都大于顯著性水平，所以在95%的置信水平下接受零假設，建模誤差是正態(tài)分布.剔除異常數據中新的ε樣本的K-S檢驗圖如圖3，右側圖所示.使用篩選后的εTS樣本計算累積分布函數（cumulation distribution function，縮寫為CDF）如圖3a右圖所示，并與正態(tài)分布的CDF進行比較，觀測和篩選的經驗累積分布函數（empirical cumulation distribution function，縮寫為 ECDF）與理論最佳正態(tài)累積分布函數（normal cumulation distribution function，縮寫為NCDF）之間的相關系數是99.99% （見表1），通過K-S檢驗不僅證明了篩選準則的可靠性，還證明了篩選后的 εTS樣 本非常接近于正態(tài)分布函數.εRTS和 εRTE的樣本的檢測結果如圖3b，c右圖所示，也得出了相同的結論.

圖3 建模誤差 ε 樣本分布（左）和 K-S 檢測圖（右）（a） TS 方法；（b） RTS 方法；（c）RTE 方法.左圖中兩條黑色直線標記了數據篩選準則|ε-median（ε）|＞2std（ε）的區(qū)間，紅色線表示最佳的正態(tài)概率密度函數值；右圖中灰色線代表著95%置信區(qū)間的CDF的上下界限，綠色小線段代表著K-S檢驗的經驗CDF和理論CDF之間的最大垂向差值Fig.3 Distribution of samples of modeling error ε （left）and results of K-S test （right）（a） TS method；（b） RTS method ；（c） RTS method.The black lines in the left figure mark the interval of |ε-median（ε）|＞2std（ε），which we use as criteria for data screening and the best-fit theoretical normal PDF of the screened samples is plotted over the distribution map with the red line.The grey lines in the right figure indicate lower and upper 95 percent confidence bounds for the CDF and the green segment indicates the maximum distance from the empirical CDF to the theoretical normal CDF which is used in the K-S test

研究中還使用Q-Q （Quartile-Quartile）分位數圖來檢測 ε是否服從正態(tài)分布.Q-Q分位數圖是一種圖形統(tǒng)計檢驗方法，用此方法可以找出兩組數據是否來自同一分布.如果在Q-Q分位數圖中縱軸表示理論正態(tài)分布分位數，橫軸上輸入驗證數據，擬合的數據如果聚類在y＝x的線上，這表明樣本數據呈正態(tài)分布（Marden，2004；Wheelan，2013）.如圖4a所示， εTS的Q-Q分位數圖顯示 εTS與正態(tài)分布基本一致，僅左下角和右上角存在一定偏差，代表著樣本中剔除了長尾異常，即尾部未采樣，但總體上可以證明 ε服從正態(tài)分布.εRTS和 εRTE的樣本的檢測結果如圖4b，c所示，也得出了相同的結論.

圖4 樣本 ε的 Q-Q 圖（a） TS 方法；（b） RTS 方法；（c） RTE 方法Fig.4 Q-Q plots of samples of ε（a） TS method；（b） RTS method；（c） RTE method

研究中還使用了其它圖形檢驗方法如正態(tài)概率圖以及數值檢驗方法如Shapiro-Wilk測試應用于檢驗篩選的 ε樣本，所有的這些方法都是為了檢驗隨機樣本在一定置信水平下符合正態(tài)分布的假設而設計的，并互相驗證測試結果.本研究中的檢驗結果與Chen和Xie （2017）開展的Q值模型建模誤差的統(tǒng)計特征一致，但本次研究中使用K-S檢驗、Q-Q分位圖并沒有像Chen和Xie （2017）的結果一樣對TS和RTS方法的正態(tài)性檢測失效，因為這些方法是為相對較?。ú怀^幾千）的樣本量設計的，川滇研究區(qū)的Q值模型的建模誤差樣本都在幾千的范圍內，可以很好地呈現出正態(tài)性.所以通過多種統(tǒng)計檢驗方法，都證明了本文所采用的數據篩選準則是合適的，可以剔除建模誤差樣本的異常值，產生一個接近于完美正態(tài)分布的新的數據集，從而使用最小二乘算法求解Q值模型.

3 川滇地區(qū)Q值模型及評價

3.1 Q值層析成像方法

研究區(qū)Lg波Q值的橫向變化可以用二維衰減層析成像解決，將研究區(qū)域劃分為M個小單元來參數化Q值的橫向變化，假設每個單元中的Q值恒定.在收集N條唯一路徑上的??/Q（f）測量值之后，對每條路徑都可以構建線性方程（Chen，Xie，2017）

式中，n＝1，2，···，N代表著第n條射線路徑，m＝1，2，···，M代表著第m個模型的小單元.Q（nf）是第n條射線路徑上的Q（f）的估計值，第n條射線路徑長度為 ?n，當它穿過研究區(qū)域的第m個小單元的時候，相交的小段距離長度為 ?nm.式（6）的矩陣形式為

式中，數據矢量d中的元素為 ?n/Q（nf），轉換矩陣中的元素為 ?nm，而模型參數矢量m中的元素為1 /Qm（f）.轉換矩陣A是一個稀疏矩陣，因為任意第n條射線路徑穿過單元格數量比總的單元格數量M少得多.另外，實際反演的模型值是1/Q，只是按照習慣在圖像中轉換為Q.

一旦證明了Q值反演的數據集受正態(tài)分布的誤差影響，則可以使用最小二乘法來反演1/Q值模型.數據的方差可以由 ε的方差得出：

為了更好地反映通過改變 ?n/Q（nf）測量值之間的方差來量化模型的不確定性，通過式（8）可以計算基于雙臺方法的任何特定路徑上的?/Q（ef）］值.在研究中引入一個N×N的數據協方差矩陣C，它的第n個對角線元素由 V a［r?n/Q（nf）］給出，然后歸一化得到下式：

最小二乘意義下求解模型m需要求解矩陣A的逆矩陣，但是在層析成像中這通常需要很大的計算量和內存，取而代之的是各種迭代技術，如反投影法（Xieet al，2004，2006）.反投影方法主要著力于迭代更新一個測試解m，直到預測的數據不匹配而停止改善.研究中使用的新的2D層析成像方法不是迭代技術，而是通過PROPACK算法（Larsen，1998）對矩陣A進行精確的SVD分解，從而求得矩陣A的逆矩陣.

3.2 川滇地區(qū)Q值模型和分辨率

可視化的棋盤格測試幾乎是所有的層析成像方法均采用的一種評估反演可靠性的方法.為了測試研究使用的Lg波頻譜振幅數據集能否很好地恢復Lg波衰減特征，使用棋盤格測試作為檢測模型的可靠性方法之一，如圖5所示.測試了不同網格尺寸的模型參數化范圍，最終選擇1°×1°作為網格反演大小.周連慶等（2008）研究表明，川滇地區(qū)最高Q值不超過800，最低Q值不小于100，因此設置了Q值交錯為100和800的棋盤格來生成合成數據.合成數據的路徑位置和數量與實際數據相同，然后通過合成數據來恢復原始的棋盤圖案來評估反演的質量.基于SVD的新的層析成像方法可以完整地恢復出頻率為1 Hz的QLg模型的棋盤格，并且恢復的棋盤格圖案與射線覆蓋范圍（圖6a）一致.

圖5 （a） 1°×1°棋盤格結果；（b）研究區(qū)棋盤格恢復Fig.5 （a） The checkboard with 1°×1° grid；（b） The recovery of the original checkboard

盡管棋盤格測試說明反演過程是穩(wěn)定可靠的，但它對模型分辨率大小的說明是相對模糊的，僅能說明射線覆蓋的好壞，但是密集的射線覆蓋并不一定意味著模型的高分辨率.一種極端情況是，如果模型地區(qū)僅僅被同一方向的射線覆蓋，而缺乏交叉方向的射線，那么模型則不會被很好地解析.

新的層析成像方法同時給出了模型的分辨率矩陣，可以通過構建Backus-Gilbert擴散函數（spread function，縮寫為 SF）來定量評價模型的分辨率（Chen，Xie，2017），通過分辨率矩陣的元素所控制的均方根距離，可以對模型的有限分辨率進行近似度量，從而給出每個反演格點的可分辨范圍.圖6b顯示出本研究中通過Backus-Gilbert擴散函數對川滇地區(qū)Q模型的分辨率進行評估，在地震射線覆蓋良好的研究區(qū)域中，分辨精度可以達到100 km.通過棋盤格測試和Backus-Gilbert擴散函數都證明了研究中得到的川滇地區(qū)Lg波的1/Q模型是可靠的.

圖6 QL g模型的 Backus-Gilbert擴散函數（a） 基于雙臺方法的射線路徑覆蓋圖；（b） 用Backus-Gilbert擴散函數來表示研究區(qū)域的模型分辨率Fig.6 Backus-Gilbert spread functions of the tomographic Lg Q model（a） The path coverage of two-station based method over the studied region；（b） Spread function of the tomographic Q model

基于SVD的新的層析成像方法和反投影方法獲得的川滇地區(qū)的地殼Q值模型如圖7所示.由于兩種方法原理不同，但都需要做光滑處理，為方便比較，兩種方法采用了相同的光滑參數.在相同的光滑處理下，不同方法的QLg的區(qū)域變化似乎相當一致，最大值大于700，最小值小于80，橫向非均勻性較強，模式具有與該地區(qū)速度結構相似的特征（如Tian 等（2021））.由圖7可以看出，川滇地區(qū)西部及其鄰近地區(qū)Q值普遍低于東部，四川盆地出現明顯的東西向變化與Chen和Xie （2017）的結果相似，同時盆地東部的高Q值與Zhao等（2013）估算的整個四川盆地相似.揚子克拉通的Q值很高（低衰減），松潘甘孜褶皺帶Q值很低（高衰減）.本研究中得到的QLg結果與Wei和Zhao （2019）的川滇地區(qū)1 Hz的QLg的區(qū)域分布基本一致.在滇中地塊西側和松潘——甘孜褶皺帶西側，反投影法得到的結果似乎有更精細的結果，但是從分辨率結果可以看出，滇中地塊的Backus-Gilbert 擴散函數值（SF）約在100 km，而松潘——甘孜褶皺帶西側以及三江斷層西側的SF值在200至300 km.所以在相同的分辨率下進行構造討論，新方法相比于反投影方法能夠體現出更精細的地殼橫向變化結果.不同的反演方法涉及不同的參數選擇，在大多數情況下參數的設置是經驗決定的，兩種方法的比較過程中盡可能使用類似的參數，所以細節(jié)方面不盡相同.影響層析成像結果的因素很多，例如反投影方法需要設置一個初始模型，不同的初始條件，反演的結果也可能不同.新方法可以通過PROPACK算法直接對轉換矩陣進行SVD分解，所以不需要初始模型，其優(yōu)越性表現在模型分辨率的定量分析，從而決定了模型所反映的結構的可靠性.

圖7 川滇地區(qū)及其鄰近地區(qū)的 Q 值模型（a） 新的層析成像方法結果；（b） 反投影方法結果Fig.7 The Lg Q model of Sichuan-Yunnan area and its neighborhood（a） Result from the new tomographic method；（b） Result from the back projection method

3.3 模型的不確定性

通過SVD奇異值分解，得到了Lg波1/Q模型的協方差矩陣.因為證實了數據誤差服從正態(tài)分布，誤差的大小完全可以由誤差的方差表示（式（8））.反演過程中，數據包含的誤差會傳遞到模型中，反映在模型的協方差矩陣中，從而由數據誤差導致的1/Q模型的誤差可以由模型的協方差矩陣完全量化.Q值層析成像模型的不確定度由模型的分辨率和誤差共同決定.一階近似下，用δQ/Q來標定模型的相對誤差，這個值可以使用協方差矩陣對角線元素的平方根來近似估計.在川滇地區(qū)，1 Hz的Q模型相對誤差如圖8所示，一個標準差下的模型的相對誤差小于9%，在射線路徑覆蓋良好的區(qū)域，相對誤差小于3%，這得益于研究中嚴格的數據篩選準則和基于SVD的新的層析成像方法.完整的Q值層析成像模型的協方差矩陣可以用來定量估計特定路徑預測Q值的誤差，當使用路徑Q值校正振幅時，這些誤差估計可用于建立地震事件識別和震源大小估計的置信水平.

圖8 每個模型單元的相對不確定性Fig.8 The relative uncertainty plot at each model cell

4 討論與結論

使用川滇地區(qū)及其鄰近地區(qū)2007年7月至2013年5月的地震數據，對Lg波的“折衷振幅數據”進行了嚴格的正態(tài)性檢驗工作，這項工作在應用最小二乘算法之前十分重要，但經常被忽略，尤其在應用類似“折衷振幅數據”這樣被處理過的數據時.用一個隨機變量 ε來表示建模誤差，將“折衷振幅數據”誤差統(tǒng)計分析的復雜任務轉化為更容易處理的對 ε樣本的統(tǒng)計分析工作，獲得了如下結論：

1）Q層析成像使用“折衷振幅數據”作為輸入數據導致的建模誤差主要由幾何擴散函數、地震事件前和相關相位前尾波的噪聲以及場地效應的方位變化影響，通過多種統(tǒng)計檢驗方法證明，以前Q值成像研究中僅僅被隱含假設為正態(tài)分布的建模誤差是真實服從正態(tài)分布的，盡管那些研究可能使用了不同的數據篩選方法.ε的分布以一個擬正態(tài)分布的峰為主，輔以向峰兩側擴散的弱長尾異常值.

2） 通過恰當的數據篩選準則，可以剔除建模誤差樣本的異常值，從而產生一個接近于完美正態(tài)分布的新的數據集，用于在最小二乘意義下求解Q值模型.

3） 研究中所使用的川滇地區(qū)地震數據表明，方差 V a（rεTS），Va（rεRTS）和 V a（rεRTE） 分別為0.025 1，0.006 7 和 0.011 0， V a（rεTS） 值最大，這可以歸因于 TS 方法在測量中不能消除場地效應的影響.

4） 研究中使用了基于SVD分解的新的Q層析成像方法來反演Q值模型，相比于反投影等迭代方法，它最大優(yōu)點是可以得到模型分辨率矩陣以及協方差矩陣，能夠把模型的分辨率和誤差定量化.

5） 新方法比反投影法可以獲得更精細的地殼Q值橫向結構，在射線覆蓋良好的區(qū)域，新方法的QLg模型可分辨達到100 km左右，模型相對誤差小于3%.研究區(qū)定量的不確定度結果表明，新的Q層析成像方法是穩(wěn)定可靠的，可以用來解釋地下介質的橫向變化，同時通過使用模型的協方差矩陣，可以定量估計特定路徑上Q預測的誤差，這樣的定量誤差估計可以在未來繼續(xù)開展在一定置信水平下的地震動預測、震源譜估計、震源類型識別和震源大小的估計.

國家測震臺網數據備份中心為本文提供了數據支持，各位審稿專家為本文提出了寶貴的修改意見，作者在此一并表示感謝.