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

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

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

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

引言

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

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

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

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

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

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)篩選

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

圖1 研究所用地震事件震中（a）和研究區(qū)臺站（b）分布圖圖中紅色矩形框代表研究區(qū)域和臺站所在區(qū)域；黑線為主要地塊分界線，灰線為主要斷裂帶，下同F(xiàn)ig.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）表示振幅項(xiàng)，S（f）為震源譜（Serenoet al，1988），G（ ?）是幾何擴(kuò)散項(xiàng)， ? 表示震中距，v表示Lg波的群速度，Q（f）為衰減值，R（f）為臺站的場地響應(yīng)項(xiàng)（Streetet al，1975；Xie，1993；Yang，2002；Chen，Xie，2017）.真實(shí)的幾何擴(kuò)散與3D地殼結(jié)構(gòu)有關(guān)，但一維幾何擴(kuò)散可以簡化為震中距 ? 的函數(shù)（Streetet al，1975），表示為

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

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

2 數(shù)據(jù)誤差的統(tǒng)計(jì)分析

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

雙臺法的Q值計(jì)算表示為

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

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

2.2 建模誤差樣本采樣

圖2 基于雙臺法的各個路徑上Q值重復(fù)測量次數(shù)分布（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)計(jì)學(xué)方法檢驗(yàn)結(jié)果

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

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

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

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

圖3 建模誤差 ε 樣本分布（左）和 K-S 檢測圖（右）（a） TS 方法；（b） RTS 方法；（c）RTE 方法.左圖中兩條黑色直線標(biāo)記了數(shù)據(jù)篩選準(zhǔn)則|ε-median（ε）|＞2std（ε）的區(qū)間，紅色線表示最佳的正態(tài)概率密度函數(shù)值；右圖中灰色線代表著95%置信區(qū)間的CDF的上下界限，綠色小線段代表著K-S檢驗(yàn)的經(jīng)驗(yàn)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）分位數(shù)圖來檢測 ε是否服從正態(tài)分布.Q-Q分位數(shù)圖是一種圖形統(tǒng)計(jì)檢驗(yàn)方法，用此方法可以找出兩組數(shù)據(jù)是否來自同一分布.如果在Q-Q分位數(shù)圖中縱軸表示理論正態(tài)分布分位數(shù)，橫軸上輸入驗(yàn)證數(shù)據(jù)，擬合的數(shù)據(jù)如果聚類在y＝x的線上，這表明樣本數(shù)據(jù)呈正態(tài)分布（Marden，2004；Wheelan，2013）.如圖4a所示， εTS的Q-Q分位數(shù)圖顯示 εTS與正態(tài)分布基本一致，僅左下角和右上角存在一定偏差，代表著樣本中剔除了長尾異常，即尾部未采樣，但總體上可以證明 ε服從正態(tài)分布.εRTS和 εRTE的樣本的檢測結(jié)果如圖4b，c所示，也得出了相同的結(jié)論.

圖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

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

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

3.1 Q值層析成像方法

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

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

式中，數(shù)據(jù)矢量d中的元素為 ?n/Q（nf），轉(zhuǎn)換矩陣中的元素為 ?nm，而模型參數(shù)矢量m中的元素為1 /Qm（f）.轉(zhuǎn)換矩陣A是一個稀疏矩陣，因?yàn)槿我獾趎條射線路徑穿過單元格數(shù)量比總的單元格數(shù)量M少得多.另外，實(shí)際反演的模型值是1/Q，只是按照習(xí)慣在圖像中轉(zhuǎn)換為Q.

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

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

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

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

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

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

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

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

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

圖7 川滇地區(qū)及其鄰近地區(qū)的 Q 值模型（a） 新的層析成像方法結(jié)果；（b） 反投影方法結(jié)果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模型的協(xié)方差矩陣.因?yàn)樽C實(shí)了數(shù)據(jù)誤差服從正態(tài)分布，誤差的大小完全可以由誤差的方差表示（式（8））.反演過程中，數(shù)據(jù)包含的誤差會傳遞到模型中，反映在模型的協(xié)方差矩陣中，從而由數(shù)據(jù)誤差導(dǎo)致的1/Q模型的誤差可以由模型的協(xié)方差矩陣完全量化.Q值層析成像模型的不確定度由模型的分辨率和誤差共同決定.一階近似下，用δQ/Q來標(biāo)定模型的相對誤差，這個值可以使用協(xié)方差矩陣對角線元素的平方根來近似估計(jì).在川滇地區(qū)，1 Hz的Q模型相對誤差如圖8所示，一個標(biāo)準(zhǔn)差下的模型的相對誤差小于9%，在射線路徑覆蓋良好的區(qū)域，相對誤差小于3%，這得益于研究中嚴(yán)格的數(shù)據(jù)篩選準(zhǔn)則和基于SVD的新的層析成像方法.完整的Q值層析成像模型的協(xié)方差矩陣可以用來定量估計(jì)特定路徑預(yù)測Q值的誤差，當(dāng)使用路徑Q值校正振幅時(shí)，這些誤差估計(jì)可用于建立地震事件識別和震源大小估計(jì)的置信水平.

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

4 討論與結(jié)論

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

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

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

3） 研究中所使用的川滇地區(qū)地震數(shù)據(jù)表明，方差 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 方法在測量中不能消除場地效應(yīng)的影響.

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

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

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