基于改進多面函數(shù)和球面整體應(yīng)變法計算華北中部地區(qū)地殼水平應(yīng)變

彭釗 王寶鎖 邵永謙 李恩建

1)天津市地震局，天津 300201 2)上海市地震局，上海 200062

0 引言

華北地區(qū)處于歐亞板塊東緣，地質(zhì)構(gòu)造環(huán)境復(fù)雜，是我國大陸強震活動較頻繁的地區(qū)之一。自GPS技術(shù)引入地震監(jiān)測預(yù)測領(lǐng)域以來，華北地區(qū)一直是GPS重點觀測區(qū)域(甘衛(wèi)軍等，2007)，很多學(xué)者利用GPS資料對該區(qū)的水平運動與應(yīng)變特征進行了研究(喬學(xué)軍等，2008；劉峽等，2010；黨學(xué)會等，2015；姚宜斌等，2015)。然而，華北地區(qū)形變信息較弱、干擾較強，計算地區(qū)應(yīng)變場時，除收集采用高精度、高密度的輸入數(shù)據(jù)外，應(yīng)使用合理的計算模型，在避免引入系統(tǒng)誤差的同時進行精度評定。

目前GPS應(yīng)變場計算主要采用三角形法或插值法等方法，這類方法模型簡明，但存在以下問題：①在數(shù)學(xué)上均存在一定的計算偏差；②對原始數(shù)據(jù)粗差抵抗性較差；③無法得到應(yīng)變率計算精度以評價可靠性。

采用球面整體應(yīng)變計算法計算應(yīng)變率在數(shù)學(xué)上是無偏的(石耀霖等，2006；楊國華等，2010)，但是，采用該方法需要位移或速度的數(shù)學(xué)解析式，而GPS站點分布不連續(xù)，其站點速率離散，采用多面函數(shù)法擬合GPS速度場可以獲得區(qū)域運動速度解析式(劉經(jīng)南等，2001)，進而應(yīng)用球面整體應(yīng)變計算法。同時，由于得到了速度和應(yīng)變率的數(shù)學(xué)解析式，可以根據(jù)誤差傳播率計算應(yīng)變結(jié)果的誤差，以評價結(jié)果的精度和可信度。

本文綜合采用改進多面函數(shù)法和球面整體應(yīng)變法構(gòu)建無偏應(yīng)變計算模型。利用華北中部地區(qū)CMONOC的GPS連續(xù)站數(shù)據(jù)，基于連續(xù)形變假設(shè)，采用改進多面函數(shù)法得到該區(qū)域GPS水平速率解析式，再用球面整體應(yīng)變計算法計算華北中部地區(qū)地殼EW向線應(yīng)變率、SN向線應(yīng)變率、NE向剪應(yīng)變率、面膨脹率，并初步評估其精度。

1 應(yīng)變場計算模型及其精度評定

1.1 改進多面函數(shù)法擬合水平速度場

采用多面函數(shù)擬合GPS水平速度場，可以由離散的GPS站點的位移或速率計算得到區(qū)域水平速度解析式，其數(shù)學(xué)模型為(Hardy，1978)

(1)

式中，B、L為緯度和經(jīng)度；n為節(jié)點數(shù)；α為待定系數(shù)；θ(B，L；Bi，Li)為核函數(shù)。用矩陣的形式描述多面函數(shù)法，即

S=θα

(2)

根據(jù)正規(guī)方程可得

α=(θTθ)-1θTS

(3)

本文采用基于Tikhonov正則化的改進多面函數(shù)法(彭釗等，2019)，該方法在經(jīng)典最小二乘準(zhǔn)則的基礎(chǔ)上，引入正則化矩陣，刪去核函數(shù)中的平滑系數(shù)，與原方法相比，其優(yōu)點有：①規(guī)避了原方法平滑系數(shù)的不確定性；②去除了節(jié)點數(shù)不得大于已測點數(shù)的約束條件，采用更多的擬合節(jié)點，既可以有效提高擬合精度，又可以避免因站點分布不均造成的偏差，對于二維水平速度場擬合有重要價值。

引入Tikhonov正則化，原理是在偏差函數(shù)后增加一個正則化項，在正規(guī)方程里，通過在待定參數(shù)中引入正則化矩陣實現(xiàn)，即

α=(θTθ+ΓTΓ)-1θTS

(4)

取Γ=aIn，即單位矩陣的a倍，稱為L2正則化。

此時式(4)變?yōu)?/p>

α=(θTθ+λIn)-1θTS

(5)

式中λ=a2為正則化系數(shù)，In為n階單位矩陣。應(yīng)用時，根據(jù)泛化誤差最小化原則確定λ的取值。

1.2 球面整體應(yīng)變計算法計算應(yīng)變率

球面整體應(yīng)變計算法基于球面坐標(biāo)而非直角坐標(biāo)計算應(yīng)變，在球面上計算應(yīng)變，直接套用直角坐標(biāo)系下的應(yīng)變算式是有偏的。在已知平面運動解析式的情況下，基于球面坐標(biāo)通過微分方程求解應(yīng)變的計算式為(楊博等，2010；楊國華等，2010)

(6)

式中，εn、εe、εne、εΔ分別為NS向線應(yīng)變率、EW向線應(yīng)變率、NE向剪應(yīng)變率、面膨脹率；Vn、Ve分別為N、E向速率；R為地球平均曲率半徑。式(6)為應(yīng)變計算的無偏算式，它適用于球面非均勻介質(zhì)連續(xù)應(yīng)變場的計算與分析。

采用改進多面函數(shù)法，可由離散點的N、E向的位移計算得到北、東向運動的解析式，進而計算研究區(qū)域內(nèi)任一點N、E向速率Vn、Ve及其偏導(dǎo)。

1.3 精度評定

根據(jù)式(1)、式(6)，由離散GPS站點速率值計算得到區(qū)域應(yīng)變率的數(shù)學(xué)解析式，進而可以根據(jù)誤差傳播率來計算應(yīng)變率結(jié)果的誤差。

由式(6)和誤差傳播律可求得各應(yīng)變率分量方差為

(7)

根據(jù)式(2)，各方向速率及其偏導(dǎo)數(shù)的方差為

(8)

參數(shù)α的單位權(quán)方差為

(9)

式中，V為殘差；P為權(quán)矩陣；r為自由度。

其方差陣為

(10)

2 計算結(jié)果及分

2.1 計算過程

采用華北中部地區(qū)93個CMONOC站(45連續(xù)站+48區(qū)域站)2011～2017年的連續(xù)時序數(shù)據(jù)計算各站點NS、EW兩方向WGS84坐標(biāo)系下的水平運動速率，根據(jù)穩(wěn)定歐亞大陸的歐拉矢量將其換算到相對于穩(wěn)定歐亞大陸的運動速率(圖1)。數(shù)據(jù)來自中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(www.cgps.ac.cn)。

圖 1 站點分布及其水平速率

采用改進多面函數(shù)法，根據(jù)式(2)、(5)得到區(qū)域內(nèi)速率解析式；根據(jù)球面整體應(yīng)變計算法代入式(6)，計算0.1°×0.1°間隔點的EW向線應(yīng)變率、SN向線應(yīng)變率、NE向剪應(yīng)變率和面膨脹率(劉經(jīng)南等，2002)；再利用式(7)～(10)計算各點應(yīng)變率的精度，因所選取的站點都屬于CMONOC站，觀測條件和精度相仿，觀測時間一致，故采用等權(quán)重方式進行計算。

2.2 改進多面函數(shù)法的擬合精度

采用改進多面函數(shù)法，計算區(qū)域內(nèi)速率解析式，核函數(shù)取基于大地長的正雙曲型核函數(shù)(楊國華等，1990；黃立人等，1993)，以全部93個站點作為節(jié)點。經(jīng)過試算，NS向正則化系數(shù)取3.5，EW向正則化系數(shù)取0.9時，模型泛化誤差最小，計算此時各站點的建模速率，其與實際速率的絕對值差值見表1。

采用原多面函數(shù)法對站點速率進行擬合，取平滑系數(shù)為1，其余參數(shù)不變，其結(jié)果與改進多面函數(shù)的擬合結(jié)果的對比見表2。

表1 建模速率與實際速率的絕對差值(單位：mm/a)

方向最小值最大值平均值標(biāo)準(zhǔn)差NS向0.0302.0200.5250.714EW向0.0043.6620.6480.977

表 2 兩種方法的擬合精度比較(單位：mm/a)

擬合方法NS向EW向殘差均值標(biāo)準(zhǔn)差殘差均值標(biāo)準(zhǔn)差多面函數(shù)法0.6240.7930.7761.142改進多面函數(shù)法0.5250.7140.6480.977誤差降低比例/%18.811.119.714.2

由表 1、2可知，改進多面函數(shù)法各站點建模速率與實際速率相差不大，其平均殘差為0.5mm/a左右，擬合精度相對于原方法有10%以上的提高，說明采用改進多面函數(shù)法建立的區(qū)域速度解析式符合整體速度場的趨勢且精度高于原方法。計算1°×1°間隔點模型速率，區(qū)域模型推估速率分布如圖2 所示。

圖 2 華北中部地區(qū)模型推估速率分布

2.3 華北中部地區(qū)NS、EW向線應(yīng)變率及其精度

計算得到的華北中部地區(qū)NS、EW向線應(yīng)變率及其中誤差如圖3、4 所示。

圖 3 華北中部地區(qū)NS向線應(yīng)變率(a)及其中誤差(b)

圖 4 華北中部地區(qū)EW向線應(yīng)變率(a)及其中誤差(b)

在NS向，區(qū)域整體以張性為主，應(yīng)變率相對較小，黃河下游帶及向南延伸的區(qū)域呈相對張性高值，晉冀蒙交界及其以西地區(qū)呈較小的壓性狀態(tài)，區(qū)域內(nèi)部呈近NW向的階梯分布；在EW向，張渤帶兩側(cè)及以西地區(qū)呈壓性狀態(tài)，山西帶、黃河下游帶及其以南地區(qū)呈張性狀態(tài)，郯廬帶山東段以東為張性相對高值地區(qū)，應(yīng)變高值分布與大型斷裂帶分布基本一致，可能表明EW向線應(yīng)變受構(gòu)造運動控制較強。

2.4 華北中部地區(qū)NE向剪應(yīng)變率和面膨脹率及其精度

計算得到的華北中部地區(qū)NE向剪應(yīng)變率及其中誤差如圖5 所示，面膨脹率及其中誤差如圖6 所示。

圖 5 華北中部地區(qū)NE向剪應(yīng)變率(a)及其中誤差(b)

圖 6 華北中部地區(qū)面膨脹率(a)及其中誤差(b)

NE向剪應(yīng)變率區(qū)間為-10.5×10-9～14.0×10-9/a，中誤差的區(qū)間為0.2×10-9～2.3×10-9/a，面膨脹率區(qū)間為-21.0×10-9～20.7×10-9/a，中誤差的區(qū)間為0.3×10-9～2.5×10-9/a，計算結(jié)果精度較高。

區(qū)域內(nèi)剪應(yīng)變分布大體呈東負西正的特征，西部山西帶中部和邢臺地區(qū)，東部郯廬帶山東段數(shù)值稍大，區(qū)域極值出現(xiàn)在唐山-秦皇島地區(qū)，反映出該區(qū)域相對于首都圈其他地區(qū)較強烈的旋剪形變；面膨脹率分布整體表現(xiàn)為北壓南張的格局，除晉冀蒙交界及其以西地區(qū)呈面收縮外，華北中部大部分地區(qū)呈面膨脹，山西帶南部和郯廬帶山東段數(shù)值較大。

2.5 應(yīng)變率計算結(jié)果可靠性分析

以應(yīng)變率與應(yīng)變率中誤差的比值作為應(yīng)變結(jié)果可靠性的評價指標(biāo)，以2倍中誤差作為應(yīng)變計算結(jié)果真實可信的閾值，即

(11)

式中，ε為應(yīng)變率；σε為應(yīng)變率中誤差。

表3應(yīng)變結(jié)果可靠性的統(tǒng)計

評價指標(biāo)NS線應(yīng)變率/%EW線應(yīng)變率/%NE剪應(yīng)變率/%E≥2908693

分別計算各格網(wǎng)點線應(yīng)變率和剪應(yīng)變率及其中誤差的比值，其統(tǒng)計結(jié)果見表3。

由表3 可知，NS、EW向線應(yīng)變率和NE向剪應(yīng)變率的計算結(jié)果評價為可信的比例都在85%以上，結(jié)合圖1～3，除了邊緣地區(qū)誤差稍高以外，大多數(shù)地區(qū)的應(yīng)變率計算結(jié)果都真實可靠。

3 討論與總結(jié)

對弱形變地區(qū)的應(yīng)變場采用無偏的計算方法是必要的，因為無偏模型能避免引入難以估計和消除的系統(tǒng)誤差，而這樣的系統(tǒng)誤差對弱形變地區(qū)來說不可忽視，同時精度和可靠性評定也有助于評價計算結(jié)果的質(zhì)量，并方便與其他研究者的結(jié)果進行對比。

本文綜合采用改進多面函數(shù)法和球面整體應(yīng)變法構(gòu)建無偏應(yīng)變計算模型，推導(dǎo)計算及精度評定公式，并給出可靠性的評價指標(biāo)?；谌A北中部93個CMONOC站2011～2017年的連續(xù)時序數(shù)據(jù)，計算華北中部地區(qū)線應(yīng)變率及剪應(yīng)變率，并計算相應(yīng)的中誤差及其比值，誤差分析表明結(jié)果具有較高的精度和可靠性，證明采用基于改進多面函數(shù)和整體球面應(yīng)變法的無偏應(yīng)變計算模型計算弱應(yīng)變應(yīng)變率是可行的，且精度和可信度都較高，應(yīng)變計算結(jié)果表明，華北中部地區(qū)應(yīng)變率整體較小，高值區(qū)主要集中在大型斷裂帶及其附近區(qū)域，應(yīng)變率梯度較高的區(qū)域主要有唐山地區(qū)、郯廬帶山東段、山西帶、張渤帶等；區(qū)域整體以張性為主，晉冀蒙交界及其以西地區(qū)呈壓性狀態(tài)；唐山-秦皇島地區(qū)具有區(qū)域內(nèi)最高的剪應(yīng)變率，表現(xiàn)出較強的旋剪形變。

本文采用的是同時間跨度CMONOC連續(xù)站和部分區(qū)域站的數(shù)據(jù)，計算使用等權(quán)方式，后續(xù)研究可通過補充更多區(qū)域站數(shù)據(jù)和多期流動GPS觀測數(shù)據(jù)來提高應(yīng)變場空間分辨率，反映更多細節(jié)信息，但是，不同區(qū)域站數(shù)據(jù)和流動觀測數(shù)據(jù)的時間跨度和數(shù)據(jù)精度可能不同，需要采用不等權(quán)的方式進行計算。另外，本文僅就應(yīng)變率場的分布特征進行了分析，在地質(zhì)及地球物理構(gòu)造解釋方面存在不足，后續(xù)可通過結(jié)合地質(zhì)、測震和其他地球物理信息對區(qū)域地殼構(gòu)造變形及其動力學(xué)解釋進行進一步研究。