凌聰聰 廖超明 楊翼飛 周聰林 陳香萍

1 廣西壯族自治區(qū)自然資源信息中心,南寧市中新路2號,530029

2 南寧師范大學自然資源與測繪學院,南寧市明秀東路175號,530001

廣西位于華南塊體西南端,是云貴高原至東南沿海山地丘陵的過渡區(qū),區(qū)域西北以高山地為主,東南多為丘陵和平原,呈自西北向東南掀斜的地貌特征[1]。區(qū)域內活動主要受NE向和NW向斷裂控制,這2個方向的斷裂構成“X”型格架,為廣西“山”型構造體系提供了構造環(huán)境[2]。開展廣西區(qū)域現(xiàn)今地殼構造活動研究,對認識中國東南沿海地區(qū)地殼構造活動及探討廣西區(qū)域構造格局的形成機制等具有重要科學意義。

針對廣西區(qū)域速度場與應變場研究成果較為陳舊的問題,本文利用廣西地區(qū)長時序連續(xù)運行參考站(GXCORS)觀測數(shù)據構建區(qū)域三維速度場模型,并聯(lián)合廣西及周邊地區(qū)136個GNSS站速度場數(shù)據分析廣西地區(qū)現(xiàn)今應變特征。

1 數(shù)據處理

選取2009-01～2019-06廣西地區(qū)83個GXCORS站觀測數(shù)據,使用GAMIT/GLOBK軟件進行解算,獲得ITRF2014參考框架下站點坐標時間序列。利用CATS軟件[3],基于白噪聲(WN)+閃爍噪聲(FN)組合模型估計站點的三維速度,同時收集相同參考框架下廣西及周邊地區(qū)中國大陸構造環(huán)境監(jiān)測網絡工程53個連續(xù)GNSS站的三維速度。

1.1 CORS站數(shù)據處理

利用中國及周邊地區(qū)10個IGS站作為框架站進行GAMIT基線處理,得到單日松弛基線解H文件,主要參數(shù)見表1。GLOBK坐標時間序列聯(lián)合了全球IGS1、IGS2、IGS3子網解,選用igb14核心站作為框架站,并以ITRF2014參考框架坐標作為先驗值。對于因地震或天線更換引起的站點坐標變化,使用地震文件(igb14_comb.eq)進行修正。采用3倍中誤差法和滑動四分位距法剔除坐標時間序列的粗差點。

1.2 CORS站速度場

GNSS站坐標時間序列是給定時間范圍內的歷元坐標解組合,GNSS站坐標與時間的關系式一般為:

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

esin(4πti)+fcos(4πti)+

(1)

式中,ti為站點坐標對應的觀測歷元,a、b為站點平均位置和速度,c、d為站點周年運動特征系數(shù),e、f為測站半周年運動特征系數(shù),gj為各種原因導致的階躍式坐標突跳,H為階躍函數(shù),Tgj為發(fā)生坐標突跳的時間歷元,vi為站點坐標時間序列的殘差值。

基于ITRF2014框架下的GXCORS站坐標時間序列結果,使用CATS軟件估計站點三維速度。為合理估計站點的速度及其不確定性,在計算站點速度時考慮WN及FN的影響。圖1、2分別為GXCORS站水平向和垂直向速度場分布,三維分量速率估值V及相應RMSE統(tǒng)計見表2(單位mm/a)。

色帶為站點觀測數(shù)據的時間跨度,下同

表2 GXCORS站N、E、U方向速率估值及相應RMSE

由表2可知,估計得到的GXCORS站三維速度場整體精度較高,水平方向精度明顯優(yōu)于垂直方向,其中觀測數(shù)據跨度為3～4 a的站點的3個方向速率估值精度較低。結合圖1可以看出,GXCORS站水平速度場整體為SEE向,優(yōu)勢方向為N106.3°E,平均速率約為34.95 mm/a,與已有研究結果基本一致[2]。另外,隨著經度的增大,N方向速率逐漸增大,E方向速率基本保持不變。由圖2可以看出,區(qū)域垂直方向構造運動整體相對較弱,除JZ14站下沉速率異常外,GXCORS站垂直方向速率在-2.28～1.96 mm/a之間,總體呈區(qū)域性緩慢隆升和沉降的運動態(tài)勢,與He等[4]的研究結果基本一致。

圖2 ITRF2014框架下GXCORS站垂直方向速度場

2 三維速度場模型構建

2.1 Kriging插值法

Kriging插值法主要依據協(xié)方差函數(shù)對隨機過程構建數(shù)學模型,進而實現(xiàn)對空間位置值的預測[5],其考慮空間位置對空間屬性的影響,利用插值點影響范圍的已知數(shù)據估計未知值,實現(xiàn)空間不均勻離散點數(shù)據的無偏最優(yōu)線性估計。Kriging插值法的基本公式為:

(2)

常用變異函數(shù)模型有指數(shù)函數(shù)模型、球狀模型和高斯模型等,其中速度場插值常采用球狀模型作為變異函數(shù),其表達式為:

(3)

式中,h為滯后距,c為基臺值,a為變程。利用式(3)擬合計算變異函數(shù)的球狀模型參數(shù),然后利用最小二乘法求解觀測站權系數(shù),最后將權系數(shù)代入方程組求得待定觀測站的估值。

2.2 水平速度場模型構建

將83個GXCORS站水平方向速率作為輸入值,利用Kriging插值法構建30′×30′分辨率的廣西區(qū)域水平速度場模型。如圖3所示,廣西區(qū)域整體表現(xiàn)為SEE向(N106.3°E)水平運動,其中E方向平均速率為33.55 mm/a,N方向平均速率為-9.76 mm/a,整體速率為34.95 mm/a,表明構建的廣西區(qū)域水平速度場模型與GXCORS站實測速度場一致。

圖3 廣西區(qū)域30′×30′水平方向速度場模型

為驗證構建的水平速度場模型的精度,分別對模型進行內符合精度和外符合精度檢驗。計算83個GXCORS站模型插值與實測值的殘差,并將殘差RMSE作為模型內符合精度的評價指標;收集廣西區(qū)域6個中國大陸構造環(huán)境監(jiān)測網絡工程連續(xù)GNSS站的實測速度值,計算其與模型插值的殘差及殘差RMSE,作為模型外符合精度的評價指標。廣西區(qū)域速度場模型內符合精度檢驗結果見圖4,外符合精度檢驗結果見表3。

圖4 內符合精度檢驗中模型插值與實測值的速率殘差分布

表3 外符合精度檢驗中模型插值與實測值的速率殘差統(tǒng)計

由圖4可知,水平速度場模型插值與實測速率的殘差在±0.35 mm/a之間,殘差RMSE為±0.11 mm/a,說明利用Kriging插值法構建的廣西區(qū)域水平速度場模型的擬合精度較高。由表3可知,除GXWZ、GXBS站N方向速率殘差絕對值大于1 mm/a外,其他測站水平方向速率殘差絕對值均小于1 mm/a,外符合精度為±0.98 mm/a,進一步說明本文構建的廣西區(qū)域水平速度場模型可靠性較高。

2.3 垂向速度場模型構建

剔除JZ14站的異常垂向速度后,將82個GXCORS站垂向速度作為輸入,利用Kriging插值法構建廣西區(qū)域現(xiàn)今地殼運動垂向速度場模型,結果見圖5。可以發(fā)現(xiàn),桂西局部(百色西北)地區(qū)沉降相對明顯,平均沉降速率約為-1.29 mm/a;桂東南(南寧至貴港、欽州、玉林)和桂中北部(柳州)平原及盆地區(qū)域沉降相對緩慢,平均沉降速率約為-0.92 mm/a;防城港西北至崇左東部及來賓中部等局部區(qū)域隆升,平均隆升速率約為0.95 mm/a。

圖5 廣西區(qū)域垂向速度場模型

對構建的廣西區(qū)域垂向速度場模型進行內符合精度和外符合精度檢驗,結果見圖4和表3?？梢钥闯?垂向速度場模型插值與實測值的殘差在±0.40 mm/a之間,內符合精度檢驗結果為±0.09 mm/a;除GXBS站和GXNN站垂向速率殘差絕對值大于1 mm/a外,其他測站垂向速率殘差絕對值均小于1 mm/a,外符合精度檢驗結果為±0.66 mm/a。

3 區(qū)域地殼應變參數(shù)反演

3.1 格網距離加權法

格網距離加權法[6-7]根據各站點的位置和速度,結合站點至格網點的距離信息,計算格網點的應變率。將區(qū)域分為連續(xù)均勻的格網,利用站點水平速度計算各格網點的速度梯度:

(4)

式中,ui為站點速度,ti為格網點的擬合速度,eij為速度梯度張量,xj為站點j與格網點的距離。將各站點按其到單元格網中心的距離進行格網加權[6],可表示為:

(5)

式中,W為加權因子,d為站點到單元格網中心的距離,α為距離平滑因子,表示站點影響隨距離衰減的程度。計算區(qū)域應變參數(shù)的主要步驟為:1)將站點速度代入式(4),得到間接平差函數(shù)模型;2)根據式(5)計算加權因子,并與站點速度誤差一起構成平差隨機模型;3)將間接平差函數(shù)模型與平差隨機模型結合,利用最小二乘法計算各格網點的水平速度梯度;4)根據速度梯度張量和應變參數(shù)計算區(qū)域應變參數(shù):

(6)

δ=ε1+ε2

(7)

γ1=ε11-ε22,γ2=ε12

(8)

式中,ε1、ε2分別為最大、最小主應變率,εij為應變張量,δ為面應變率,γ1、γ2分別為第1、第2剪切應變。

3.2 應變場特征分析

基于廣西及周邊區(qū)域136個GNSS站(83個GXCORS站和53個陸態(tài)站)監(jiān)測到的水平速度場,利用格網距離加權法計算50 km×50 km分辨率區(qū)域現(xiàn)今地殼應變參數(shù),并繪制剪切應變率、面應變率和主應變分布,見圖6、7。

圖6 廣西區(qū)域剪切應變和主應變場

3.2.1 主應變特征分析

從圖6、7可以看出,桂北地區(qū)以近SN向主壓應變?yōu)橹?最大主壓應變位于貓兒山地區(qū),約為-4.21×10-9/a;桂中及東西鄰近區(qū)域(23.5°N附近)主壓應變逐步變小,并過渡為近SN向的主張應變;桂南地區(qū)往南NS向主張應變逐漸增大,且出現(xiàn)微小的近EW向主壓應變;桂東地區(qū)以NW向主張應變?yōu)橹?局部(梧州市、賀州市)兼NE向主壓應變;桂西地區(qū)以NE向主張應變?yōu)橹?兼微小NW向主壓應變。桂北地區(qū)自北向南主壓應變軸由近NS向逐漸轉為NW向(西部)和NE向(東部)。桂中地區(qū)的廣西“山”字型構造弧頂北部主壓應變軸逐漸轉為近EW向,且主壓應變及主張應變率接近最小,分別約為0.13×10-9/a和-0.26×10-9/a。桂中往南主壓應變逐漸轉為主張應變,到桂南地區(qū)主張應變率逐漸增大,最大應變位于桂南地區(qū)合浦-北流斷裂中部附近,約為5.07×10-9/a。主張應變軸在108.5°E附近以NS向為主,并向東西兩側呈扇形展開,在桂東南地區(qū)呈NNW向主張應變,桂西南地區(qū)呈NNE向主張應變。

3.2.2 剪切應變特征分析

從圖6可以看出,區(qū)域剪切應變總體呈四周高、中間低的特征,其中剪切應變率高值區(qū)對應主應變率高值區(qū),分別為桂北三江、龍勝、資源一帶、桂東梧州至玉林(蒼梧、龍圩、岑溪、北流)一帶、桂南北海、欽州、防城港一帶和桂西靖西、德保一帶。剪切應變率低值區(qū)位于23.5°N附近的來賓平原,該區(qū)域主應變率最小,也是主壓應變向主張應變轉化的過渡區(qū)。桂南地區(qū)剪切應變高值區(qū)位于NE向防城-靈山斷裂帶、合浦-北流斷裂帶與NW向百色-合浦斷裂帶交會處,其中欽州灣附近剪切應變率最大達6.43×10-9/a。桂東地區(qū)剪切應變高值區(qū)位于NE向合浦-北流斷裂帶及NNE向賀街-夏郢斷裂帶附近,其中梧州附近剪切應變率最大達5.56×10-9/a。桂北地區(qū)剪切應變高值區(qū)位于近NS向龍勝-永福斷裂帶北段,其中貓兒山附近剪切應變率高達6.15×10-9/a。桂西地區(qū)剪切應變高值區(qū)位于NW向靖西-崇左斷裂帶西北段,最大剪切應變率約為4.65×10-9/a。剪切應變率低值區(qū)位于NE向南寧-桂林斷裂帶與NW向巴馬-博白斷裂帶交會處,該區(qū)域恰好位于廣西“山”型構造弧頂北部,最大剪切應變率約為(0.5～2)×10-9/a。

在廣西“山”字型構造研究中,張文佑[8]利用材料力學中橫梁彎曲原理說明了該構造區(qū)域的受力特征,指出“山”字構造弧形與脊軸之間的中和帶地層受應力最小,因此常有盆地和平原發(fā)育。結合圖2和6可以看到,桂中區(qū)域的應變參數(shù)反演結果較好地監(jiān)測到了廣西“山”字型構造的現(xiàn)今應變特征。

3.2.3 面應變特征分析

由圖7可以看出,區(qū)域面應變整體呈北部區(qū)域壓縮、南部區(qū)域拉張的應力狀態(tài)。其中,面應變壓縮區(qū)主要位于桂林、柳州、河池、百色北及來賓北一帶,最大為-4.68×10-9/a;面應變拉張區(qū)域主要位于崇左、欽州、防城港、北海、玉林和梧州一帶,最大為5.09×10-9/a。結合主應變場特征可知,在印度板塊的碰撞下,青藏高原中部物質E向及SE向擠出,直接作用于華南塊體西面,造成北部區(qū)域擠壓變形;而南部區(qū)域拉張變形可能與南海北側張性活動相關[9]。結合廣西區(qū)域地形地貌來看,桂西、桂西北、桂北和桂東北地區(qū)以山地為主,桂中、桂東南、桂南和桂西南地區(qū)則主要為盆地及平原,反演得到的面應變表現(xiàn)出北部地殼壓縮增厚、南部地殼拉張變薄,總體上與區(qū)域現(xiàn)代地貌一致。

圖7 廣西區(qū)域面應變和主應變場

3.3 區(qū)域應變與震情時空分布分析

2012～2022年廣西境內發(fā)生2級以上地震29次,其中有感地震7次,小震17次,中強震5次,平均震源深度為7.4 km,具體信息見表4。結合圖6可以看出,97%的地震發(fā)生在剪切應變率(3.0～5.0)×10-9/a的中高值區(qū)域,這與已有研究認為區(qū)域地震多發(fā)生在剪切應變高值區(qū)向低值區(qū)過渡帶上的結論[10-12]一致。

表4 研究區(qū)10 a來2級以上地震信息

由圖6和表4可知,2012～2022年研究區(qū)發(fā)生的29次2級以上地震震中主要分布在NW向和NE向2組斷裂帶附近。以斷裂帶為單元,發(fā)震中心的時空遷移具有一定規(guī)律,一般自斷裂帶一端向另一端遷移,然后逐步轉折,且按時序先后發(fā)震。如NE向南寧-桂林斷裂帶的發(fā)震中心自北東端(柳州)向南西端(崇左)遷移,然后轉至北東端(江南區(qū));NW向巴馬-博白斷裂帶、靖西-那坡斷裂帶的發(fā)震中心自北西端向南東端遷移,然后轉至北西端。另外,從2013～2021年發(fā)生的5次中強震的遷移路線來看,震中先從桂西平果縣遷移至桂東蒼梧縣,然后轉向南至桂東南北流市,最后往西北遷至桂西靖西市、德?？h。

利用震級與地震釋放能量的關系計算得到,2012～2022年研究區(qū)29次2級以上地震共釋放1.805×1013J的能量。以發(fā)震斷裂帶為統(tǒng)計單元,NW向巴馬-博白斷裂帶、百色-合浦斷裂帶及靖西-那坡斷裂帶累積釋放1.001×1013J的能量;NNE向賀街-夏郢斷裂帶在2 a左右累積釋放能量達7.944×1012J。研究區(qū)地震能量的釋放主要分布在桂西NW向主控斷裂帶和桂東NNE向賀街-夏郢斷裂帶,約占總釋放能量的99.5%,桂西和桂東局部區(qū)域的應變能得到有效釋放。目前,桂北和桂南地區(qū)處于主應變和剪應變高值區(qū),隨著區(qū)域應變能的逐步積累,發(fā)生中強地震的可能性進一步增大。

3.4 區(qū)域構造動力機制分析

在青藏高原E向和SE向推擠力作用下,華南塊體西部經由川滇菱形塊體轉向SEE向,并施加在桂西區(qū)域[13]。在北面,華北塊體的圍限阻擋導致NS向產生擠壓,并作用于桂北區(qū)域,使廣西區(qū)域現(xiàn)今呈西部NW-SE向、北部SN向的主壓應變與面壓縮變形狀態(tài);在東面,菲律賓海板塊的W向仰沖碰撞透過南海東部作用于桂東區(qū)域;在南面,南海北側的張裂應變直接作用于桂南及北部灣區(qū)域[2,9],使廣西區(qū)域西南部、南部及東南部呈現(xiàn)主應變軸扇形展開的主張應變與面拉張變形狀態(tài)。研究區(qū)中部(23.5°N附近)的主壓應變率和主張應變率變化最小,是主壓應變向主張應變轉換的構造過渡帶。

4 結 語

1)ITRF2014框架下廣西區(qū)域水平方向平均運動速率為34.95 mm/a,優(yōu)勢方向N106.3°E;垂向運動速率為-2.28～1.96 mm/a,表現(xiàn)出區(qū)域性隆升或沉降。Kriging插值法構建的廣西區(qū)域水平向速度場模型的內、外符合精度分別為±0.11 mm/a和±0.98 mm/a,垂向速度場模型的內、外符合精度分別為±0.09 mm/a和±0.66 mm/a,說明本文構建的區(qū)域速度場模型精度可靠。

2)聯(lián)合廣西及周邊136個GNSS站水平向速度場數(shù)據反演區(qū)域現(xiàn)今地殼應變參數(shù)發(fā)現(xiàn),研究區(qū)北部呈近NS向主壓應變與壓縮面應變狀態(tài);中部地區(qū)(23.5°N附近)主壓應變、主張應變、剪應變和面應變均最小,為主壓應變向主張應變的過渡帶;南部主張應變自北向南逐漸增大,且主張應變呈扇形展開,面應變表現(xiàn)為拉張狀態(tài);剪切應變場呈四周高值、中間過渡區(qū)低值的變化特征。地震高發(fā)區(qū)主要分布在剪切應變高值向低值轉換的過渡區(qū)域。

3)研究區(qū)近10 a來的地震震中主要集中在NW向、NE向2組斷裂帶附近,發(fā)震中心呈自斷裂帶一端向另一端遷移,后又逐步轉折的特點。29次2級以上地震累積釋放能量1.805×1013J,其中桂西NW向主控斷裂帶和桂東NNE向賀街-夏郢斷裂帶共釋放99.5%的應變能。隨著區(qū)域應變能的逐步積累,桂南和桂北地區(qū)的應變高值區(qū)發(fā)生中強地震的可能性進一步增大。

致謝：感謝中國地震臺網中心及國家地震科學數(shù)據中心提供數(shù)據。