聯(lián)合安徽CORS分析郯廬斷裂帶南段地殼形變與斷層活動(dòng)特征

陳 皓 陶庭葉 李水平 李江洋 時(shí)夢(mèng)杰 高 飛

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

2.2 多尺度球面小波估計(jì)速度

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

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

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

在球面上二維GPS速度場(chǎng)的表達(dá)式為：

(2)

(3)

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

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

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

圖4 多尺度球面小波估計(jì)速度場(chǎng)與實(shí)際速度場(chǎng)的矢量對(duì)比Fig.4 Vector comparison of velocity field estimated by multi-scale spherical wavelet and actual velocity field

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

2.3 應(yīng)變場(chǎng)解算

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

(4)

3 郯廬斷裂帶南段地殼形變與斷層活動(dòng)分析

3.1 應(yīng)變場(chǎng)分析

利用多尺度球面小波法計(jì)算郯廬斷裂帶南段的面膨脹率結(jié)果如圖6所示，(a)～(d)4幅圖分別表示空間分解尺度6～9下的面膨脹率,圖中面膨脹率為正值表示膨脹，負(fù)值表示壓縮。當(dāng)qmax=6時(shí)，要求最小站點(diǎn)間的距離不超過(guò)174.4 km。該圖可以明顯反映出整個(gè)郯廬斷裂帶南段大區(qū)域的構(gòu)造特點(diǎn)，面膨脹率最大值為9.995×10-9/a。當(dāng)qmax=9時(shí)，要求最小站點(diǎn)間的距離不超過(guò)21.8 km，此時(shí)面膨脹率最大值為1.976×10-8/a。總體上看，郯廬斷裂帶南段以西除阜陽(yáng)、徐州、棗莊等地，大部分地區(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時(shí)的應(yīng)變特征結(jié)果可以發(fā)現(xiàn)，兩者間的差異較小，故本文將針對(duì)最大分解尺度qmax=8時(shí)的應(yīng)變特征進(jìn)行分析。

圖7為qmax=8時(shí)郯廬斷裂帶南段及周邊地區(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)變?cè)谛熘莞浇谋憩F(xiàn)最為異常，達(dá)到4.044×10-8/a。其次，商丘、棗莊、阜陽(yáng)、黃山等地最大剪應(yīng)變率也較大?？傮w上看，最大剪應(yīng)變率沿郯廬斷裂帶南段兩側(cè)呈對(duì)稱分布，其西北和東南地區(qū)為高值區(qū)域，存在著應(yīng)變積累，應(yīng)予以關(guān)注。

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

3.2 速度剖面分析

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

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

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

4 結(jié) 語(yǔ)

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

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

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

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