陳 皓 陶庭葉 李水平 李江洋 時夢杰 高 飛

1 合肥工業(yè)大學(xué)土木與水利工程學(xué)院，合肥市屯溪路193號，230009

郯廬斷裂帶是一條位于我國東部地區(qū)的最大規(guī)模的斷裂帶，其主要表現(xiàn)為右旋走滑[1]。其內(nèi)部構(gòu)造復(fù)雜，歷史上多次發(fā)生5.0級以上地震，其中1668年發(fā)生的郯城地震達到8.5級。為探究2011年日本宮城地震是否對郯廬斷裂帶南段的運動性質(zhì)產(chǎn)生影響，對該地區(qū)的地殼形變與斷層活動特征展開研究。李彥川等[2]采用塊體負位錯模型分析郯廬斷裂帶中南段的地殼形變特征；李臘月等[1]對不同時期的速度場利用負位錯模型反演郯廬斷裂帶中南段的閉鎖程度和滑動虧損速率，認(rèn)為日本宮城地震緩解了郯廬斷裂帶中南段的應(yīng)變積累。以往的研究中，采用的站點基本來源于中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)項目(CMONOC)，然而郯廬斷裂中南段站點密度有限，無法揭示區(qū)域斷層活動的細節(jié)特征。隨著安徽連續(xù)運行參考站(Anhui continuously operating reference station，AHCORS)的建設(shè)與發(fā)展，增加了郯廬斷裂帶南段地區(qū)的站點密度，可為揭示區(qū)域活動斷層的變形特征和地震危險性提供更詳細的數(shù)據(jù)約束。

利用全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)的觀測資料建立適當(dāng)?shù)臄?shù)學(xué)模型，解算某一地區(qū)的地殼應(yīng)變特征是研究地震危險性的重要手段。多尺度球面小波函數(shù)可以根據(jù)GNSS測站的密集程度和分布特征選擇不同的尺度因子，揭示不同空間尺度下的應(yīng)變特征[3]。蘇小寧等[4]將多尺度球面小波解算GPS應(yīng)變場的結(jié)果與位錯模型結(jié)果對比，驗證多尺度球面小波解算應(yīng)變的有效性和穩(wěn)健性；徐克科等[5]提出多尺度球面小波的尺度和位置的選取、模型正則化因子以及參數(shù)估計的方法；李承濤等[6]采用多尺度球面小波解算不同空間尺度下巴顏喀拉塊體東北緣的應(yīng)變率場，并分析其與2017年九寨溝MS7.0地震的聯(lián)系。本文基于2013-01～2018-06 AHCORS和CMONOC的同期共170個測站的觀測數(shù)據(jù)，利用GAMIT/GLOBK解算其速度場，通過多尺度球面小波函數(shù)估計郯廬斷裂帶南段在不同空間尺度下的應(yīng)變參數(shù)，揭示郯廬斷裂帶南段及鄰近地區(qū)的地殼變形特征。研究結(jié)果可為該區(qū)域未來地震危險性評估提供參考。

1 GNSS數(shù)據(jù)處理

本文使用的數(shù)據(jù)主要是2013-01～2018-06期間37個AHCORS參考站的觀測文件，利用GAMIT/GLOBK(Ver.10.7)軟件對每隔7 d的觀測數(shù)據(jù)進行處理。解算時引入研究區(qū)域以外的9個歷史同期國際全球衛(wèi)星導(dǎo)航定位服務(wù)(IGS)站(BJFS、DAEJ、IRKT、KIT3、LHAZ、POL2、TCMS、URUM、WUHN)，從而與AHCORS參考站構(gòu)成一個區(qū)域網(wǎng)。

數(shù)據(jù)處理主要分為3個步驟：首先利用GAMIT軟件解算AHCORS參考站以及周邊IGS站的觀測數(shù)據(jù)，獲取單日松弛解，GAMIT單日解的解算策略如表1所示；然后利用GLOBK聯(lián)合SOPAC(Scripps Orbit and Permanent Array Center)機構(gòu)發(fā)布的全球h文件進行網(wǎng)平差處理；最后通過坐標(biāo)轉(zhuǎn)換獲取AHCORS參考站在ITRF2008框架下的坐標(biāo)和速度。

表1 數(shù)據(jù)處理策略

采用雙差相位解算方法，利用IGS的精密星歷和地球旋轉(zhuǎn)參數(shù)對其進行強約束，同時對測站坐標(biāo)和衛(wèi)星軌道進行解算。

采用GLOBK解算時間序列和速度場的處理策略如下：1)引入SOPAC的7個全球h文件(igs1～7)，將其與GAMIT基線解算得到的AHCORS參考站的h文件按照1.0的權(quán)重因子合并；2)利用GLOBK進行網(wǎng)平差，以IGS站在ITRF2008框架下的坐標(biāo)和速度為基準(zhǔn)，求得AHCORS參考站的速度。

由于ITRF2008框架以地球質(zhì)心為基準(zhǔn)，該框架下的點位運動速度包含了板塊的運動速度，難以準(zhǔn)確表達AHCORS參考站間的相對運動和形變，因此需要將ITRF2008框架下的運動速度轉(zhuǎn)換到歐亞框架下。本文采用Wang等[7]提出的相對于ITRF2008框架下的歐拉矢量(-0.087 mas/a，-0.514 mas/a，0.741 mas/a)，將AHCORS參考站在ITRF2008框架下的速度減去歐拉矢量所得到的速度，獲得其在歐亞框架下的速度。AHCORS參考站在ITRF2008框架下的速度場如圖1(a)所示，在歐亞框架下的速度場如圖1(b)所示。從圖1(b)可以看出，安徽區(qū)域在歐亞框架下的速度主要往東偏南方向運動，與袁鵬等[8]的結(jié)論一致。表2(單位mm/a)為歐亞框架下AHCORS參考站在N、E、U方向的速度估值和中誤差統(tǒng)計。

圖1 郯廬斷裂帶南段及周邊區(qū)域水平運動速度場Fig.1 Horizontal velocity field on the southern of theTanlu fault and its surrounding areas

表2 歐亞框架下AHCORS參考站在N、E、U方向速度估值及中誤差統(tǒng)計

2 聯(lián)合AHCORS參考站的多尺度球面小波應(yīng)變場解算

2.1 多尺度球面小波基本原理

球面小波具有將某一給定空間進行局部化的特點，可以用來建立速度場和應(yīng)變場的模型。其主要原理為，在半徑為1的單位球上定義一個小波母函數(shù)，然后將母函數(shù)通過平移和伸縮得到小波基函數(shù)[4]。本文選擇高斯函數(shù)的差(difference of gaussians，DOG)形成的球面DOG小波作為基函數(shù)，其表達式為：

(1)

式中，尺度a=1/2q，q為尺度因子，q的值越大，尺度越?。沪脼橛^測點與球面小波中心間的夾角，一般在0°～180°之間；α是一個大于1的常數(shù)，本文根據(jù)實驗結(jié)果將α取值為1.25。

2.2 多尺度球面小波估計速度

在使用多尺度球面小波解算郯廬斷裂帶南段應(yīng)變場前，需要根據(jù)GNSS測站分布的疏密情況來選擇分解尺度。多尺度球面小波方法具有可將給定空間局部化的特點，可有效消除測站分布不均勻帶來的影響。選擇分解尺度的原則為：在測站分布密集的區(qū)域選擇較大的分解尺度，在測站分布相對稀疏的區(qū)域選擇較小的分解尺度。

本文根據(jù)AHCORS參考站和郯廬斷裂帶南段附近的CMONOC站點的疏密程度來選擇分解尺度，如圖2所示，圖中黑色實心小圓點為測站。從圖2可以看出，小波分解的最大分解尺度為6～9級，測站間距為21.8 km的9級分解尺度區(qū)域共有8個，其余大部分區(qū)域是測站間距為43.6 km的8級分解尺度。

圖2 由AHCORS參考站分布確定的球面小波分解尺度Fig.2 Spherical wavelet decomposition scale based on AHCORS distribution

在球面上二維GPS速度場的表達式為：

(2)

(3)

式中，ν(θ,φ)為測站速度，gk(θ,φ)為DOG小波基函數(shù)，ak和bk為待估參數(shù)。由于gk(θ,φ)是經(jīng)過平移和伸縮得到的，屬于非正交基，存在著法方程系數(shù)矩陣秩虧的問題，因此解算出的待估參數(shù)不唯一。本文可采用廣義交叉驗證法(generalized cross-validation，GCV)確定正則化參數(shù)，結(jié)果如圖3所示。

圖3 正則化參數(shù)Fig.3 Regularization parameter

根據(jù)多尺度球面小波法可以模擬郯廬斷裂帶南段區(qū)域的速度場。圖4為多尺度球面小波估計速度場的矢量圖，圖4(a)為實際速度場與估計速度場的對比，圖4(b)為殘差矢量圖。圖5為多尺度球面小波估計速度場的標(biāo)量圖，圖5(a)為實際的速度場，圖5(b)為利用球面小波理論估計的速度場，圖5(c)為二者之間的殘差。從圖4和圖5可以看出，利用球面小波法估計得到的速度場與實際速度場的誤差在2 mm/a以內(nèi)，整體運動情況也大致相同。因此，多尺度球面小波法可以有效估計速度場。

圖4 多尺度球面小波估計速度場與實際速度場的矢量對比Fig.4 Vector comparison of velocity field estimated by multi-scale spherical wavelet and actual velocity field

圖5 多尺度球面小波估計速度場與實際速度場的標(biāo)量對比Fig.5 Scalar comparison of velocity field estimated by multi-scale spherical wavelet and actual velocity field

2.3 應(yīng)變場解算

在利用GCV法確定模型的正則化參數(shù)后，可利用Savage等[9]、孟國杰等[10]提出的在球面坐標(biāo)系下應(yīng)變與位移的微分公式，來求得目標(biāo)點的應(yīng)變張量。由于垂直方向上產(chǎn)生的形變對于水平方向的應(yīng)變影響較小，且將不處于同一球面的測站投影到同一球面來進行應(yīng)變解算的過程較為復(fù)雜，本文不予考慮。應(yīng)變張量的計算表達式為：

(4)

3 郯廬斷裂帶南段地殼形變與斷層活動分析

3.1 應(yīng)變場分析

利用多尺度球面小波法計算郯廬斷裂帶南段的面膨脹率結(jié)果如圖6所示，(a)～(d)4幅圖分別表示空間分解尺度6～9下的面膨脹率,圖中面膨脹率為正值表示膨脹，負值表示壓縮。當(dāng)qmax=6時，要求最小站點間的距離不超過174.4 km。該圖可以明顯反映出整個郯廬斷裂帶南段大區(qū)域的構(gòu)造特點，面膨脹率最大值為9.995×10-9/a。當(dāng)qmax=9時，要求最小站點間的距離不超過21.8 km，此時面膨脹率最大值為1.976×10-8/a?？傮w上看，郯廬斷裂帶南段以西除阜陽、徐州、棗莊等地，大部分地區(qū)呈面膨脹狀態(tài)，以東地區(qū)大多表現(xiàn)為面壓縮狀態(tài)。比較不同空間尺度下的面膨脹率可以發(fā)現(xiàn)，某些局部區(qū)域存在差異。因此，在小的分解尺度下可以反映出大尺度的應(yīng)變積累特征，在大的分解尺度下能夠反映出小區(qū)域范圍內(nèi)的應(yīng)變積累特征。

圖6 基于多尺度球面小波解算的面膨脹率Fig.6 The dilatation rate of multiscale spherical wavelet solution

比較最大分解尺度qmax=8和qmax=9時的應(yīng)變特征結(jié)果可以發(fā)現(xiàn)，兩者間的差異較小，故本文將針對最大分解尺度qmax=8時的應(yīng)變特征進行分析。

圖7為qmax=8時郯廬斷裂帶南段及周邊地區(qū)的主應(yīng)變率/面膨脹率、最大剪應(yīng)變率的分布。從圖7(a)可以看出，郯廬斷裂帶泗洪至莒南段西側(cè)的徐州附近主要表現(xiàn)為E-W向的拉張，商丘北部表現(xiàn)為NE-SW方向的拉張、NW-SE方向的壓縮；商丘附近則相反，轉(zhuǎn)變?yōu)镹E-SW方向的壓縮、NW-SE方向的拉張。宿松、岳西等地表現(xiàn)為S-N向壓縮，黃山、婺源等地表現(xiàn)為近E-W向的拉張。郯廬斷裂帶的泗洪至莒南段東側(cè)大多表現(xiàn)為S-N向的壓縮。整體上看，郯廬斷裂帶南段西側(cè)的主應(yīng)變比東側(cè)的大，其中最大張應(yīng)變位于郯廬斷裂帶西側(cè)的徐州附近，主應(yīng)變的最大值為2.327×10-8/a。商丘、岳西、婺源等應(yīng)變方向變化異常地區(qū)，均處于面膨脹與面壓縮的交界地帶，地質(zhì)構(gòu)造可能較為復(fù)雜，存在著較大的應(yīng)變積累，后期有發(fā)生地震的可能性，應(yīng)予以關(guān)注。從圖7(b)可以看出，最大剪應(yīng)變在徐州附近的表現(xiàn)最為異常，達到4.044×10-8/a。其次，商丘、棗莊、阜陽、黃山等地最大剪應(yīng)變率也較大?？傮w上看，最大剪應(yīng)變率沿郯廬斷裂帶南段兩側(cè)呈對稱分布，其西北和東南地區(qū)為高值區(qū)域，存在著應(yīng)變積累，應(yīng)予以關(guān)注。

圖7 qmax=8時的主應(yīng)變率/面膨脹率和最大剪應(yīng)變率Fig.7 Principal strain rate/ dilatation rate and maximum shear strain rate at qmax=8

3.2 速度剖面分析

為進一步分析郯廬斷裂帶南段及其周邊地殼形變和斷層活動情況，本文橫跨郯廬斷裂帶南段設(shè)置3個速度剖面，選擇的剖面位置如圖1(b)所示，剖面長度為600 km，寬度為100 km。圖8為3個速度剖面的結(jié)果，由北向南依次排列，圖中圓點表示AHCORS參考站的速度，三角形表示CMONOC站點的速度。表3為斷層滑動速率的估值結(jié)果。從圖8和表3可以看出，郯廬斷裂帶南段主要表現(xiàn)為右旋走滑，其北部存在著拉張，拉張分量向南逐漸遞減，在最南部轉(zhuǎn)變?yōu)閴嚎s的狀態(tài)。平行斷裂的運動速率為0.10～1.21 mm/a， 垂直斷裂的運動速率為-0.52～0.76 mm/a。郭良遷等[11]利用塊體-應(yīng)變模型反演得到的郯廬斷裂帶表現(xiàn)為右旋擠壓，李彥川等[2]利用負位錯模型得到的郯廬斷裂帶南部的特征為右旋走滑為主兼拉張。上述結(jié)論與表3的a、b段結(jié)果基本吻合。由于本文還結(jié)合了AHCORS的速度場，所以b段的拉張分量較弱，更能反映該區(qū)域的斷裂帶特征。c段剖面處于郯廬斷裂帶桐城-廬江段，宋傳中等[12]通過地質(zhì)學(xué)分析認(rèn)為，郯廬斷裂帶桐城-廬江段在新第三紀(jì)以來表現(xiàn)為逆沖擠壓的特征，與本文得到的c段結(jié)論較為一致。

圖8 郯廬斷裂帶南段的3個速度剖面Fig.8 Three velocity profiles on the southernsegment of the Tanlu fault

表3 郯廬斷裂帶南段斷層滑動速率(右旋拉張為正)

4 結(jié) 語

本文采用GAMIT/GLOBK解算2013-01～2018-06期間AHCORS參考站在歐亞框架下的速度場，同時結(jié)合CMONOC速度場利用多尺度球面小波法計算郯廬斷裂帶南段及周邊地區(qū)的應(yīng)變率，并對橫跨郯廬斷裂帶南段的3個速度剖面進行分析。得出以下結(jié)論：

1)安徽地區(qū)在歐亞框架下主要朝東偏南方向運動。本文利用多尺度球面小波法模擬出的速度場與實際速度場的誤差在2 mm/a以內(nèi)，可有效反映該區(qū)域的地殼運動情況。

2)總體上看，郯廬斷裂帶南段以西除阜陽、徐州、棗莊等地外，大部分地區(qū)呈面膨脹狀態(tài)，以東地區(qū)大多表現(xiàn)為面壓縮狀態(tài)。根據(jù)最大分解尺度下的主應(yīng)變率/面膨脹率和最大剪應(yīng)變率的分布情況可以看出，應(yīng)變方向變化異常的地區(qū)，均出現(xiàn)在面膨脹與面壓縮的交界地帶。最大剪應(yīng)變率沿郯廬斷裂帶南段兩側(cè)呈對稱分布，其西北和東南地區(qū)為高值區(qū)域，構(gòu)造特征復(fù)雜，存在著較大的應(yīng)變積累，有發(fā)生地震的可能，應(yīng)予以關(guān)注。

3)郯廬斷裂帶南段主要表現(xiàn)為右旋走滑，且北部存在著拉張，拉張分量向南逐漸遞減，在最南部轉(zhuǎn)變?yōu)閴嚎s的狀態(tài)。平行斷裂的運動速率為0.10～1.21 mm/a-1，垂直斷裂的運動速率為-0.52～0.76 mm/a。