彭釗 王寶鎖 邵永謙 李恩建

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

0 引言

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

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

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

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

1 應(yīng)變場(chǎng)計(jì)算模型及其精度評(píng)定

1.1 改進(jìn)多面函數(shù)法擬合水平速度場(chǎng)

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

(1)

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

S=θα

(2)

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

α=(θTθ)-1θTS

(3)

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

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

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

(4)

取Γ=aIn，即單位矩陣的a倍，稱為L(zhǎng)2正則化。

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

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

(5)

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

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

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

(6)

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

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

1.3 精度評(píng)定

根據(jù)式(1)、式(6)，由離散GPS站點(diǎn)速率值計(jì)算得到區(qū)域應(yīng)變率的數(shù)學(xué)解析式，進(jìn)而可以根據(jù)誤差傳播率來計(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 計(jì)算結(jié)果及分

2.1 計(jì)算過程

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

圖 1 站點(diǎn)分布及其水平速率

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

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

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

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

表1 建模速率與實(shí)際速率的絕對(duì)差值(單位：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改進(jìn)多面函數(shù)法0.5250.7140.6480.977誤差降低比例/%18.811.119.714.2

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

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

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

計(jì)算得到的華北中部地區(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)變率相對(duì)較小，黃河下游帶及向南延伸的區(qū)域呈相對(duì)張性高值，晉冀蒙交界及其以西地區(qū)呈較小的壓性狀態(tài)，區(qū)域內(nèi)部呈近NW向的階梯分布；在EW向，張渤帶兩側(cè)及以西地區(qū)呈壓性狀態(tài)，山西帶、黃河下游帶及其以南地區(qū)呈張性狀態(tài)，郯廬帶山東段以東為張性相對(duì)高值地區(qū)，應(yīng)變高值分布與大型斷裂帶分布基本一致，可能表明EW向線應(yīng)變受構(gòu)造運(yùn)動(dòng)控制較強(qiáng)。

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

計(jì)算得到的華北中部地區(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，計(jì)算結(jié)果精度較高。

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

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

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

(11)

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

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

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

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

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

3 討論與總結(jié)

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

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

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