田 曉，鄭洪艷，蘇廣利，王家慶

(中國地震局第一監(jiān)測中心，天津 300180)

地殼垂直運(yùn)動是地殼運(yùn)動的一個重要分量，由于空間大地測量技術(shù)測定的高程分量精度低于水平分量，因此利用重復(fù)高精度水準(zhǔn)測量仍然是目前研究地殼垂直運(yùn)動的主要方法。但由于水準(zhǔn)重合點(diǎn)是離散的，且水準(zhǔn)網(wǎng)內(nèi)部沒有數(shù)據(jù)，因此動態(tài)平差得到的垂直運(yùn)動速率在時間域和空間域都是不連續(xù)的，實(shí)際上地殼垂直運(yùn)動是一種連續(xù)的形態(tài)變化。為了實(shí)現(xiàn)在時空域內(nèi)客觀準(zhǔn)確地描述形變分布特征，往往需要基于離散的測點(diǎn)數(shù)據(jù)構(gòu)建合適的數(shù)學(xué)模型擬合(最佳逼近)整個地區(qū)的地殼垂直形變速率面[1]。

在我國，從20世紀(jì)80年代末期開始，黃立人、趙承坤、楊國華和陶本藻等研究了多面函數(shù)擬合法在地殼垂直運(yùn)動中的應(yīng)用，幾位學(xué)者針對多面函數(shù)擬合中核函數(shù)的選擇、平滑因子的選擇及核函數(shù)中心點(diǎn)的選擇作了細(xì)致的研究，并提出了一系列獨(dú)特的見解[2-6]。劉經(jīng)南等[7]研究了多面函數(shù)擬合法在建立我國地殼平面運(yùn)動速度場模型中的應(yīng)用。王文利等[8]利用多面函數(shù)擬合了我國陸地垂直運(yùn)動速率圖。劉同文等[9]應(yīng)用多面函數(shù)模型對西秦嶺地區(qū)的垂直形變速率作了研究和探討。秦姍蘭等[10]利用水準(zhǔn)資料研究了西秦嶺地區(qū)的垂直形變，應(yīng)用多面函數(shù)對區(qū)域整體形變特征進(jìn)行擬合。

綜上所述，研究地殼垂直運(yùn)動中，多面函數(shù)法是一種非常傳統(tǒng)并得到廣泛應(yīng)用的方法，而其他方法應(yīng)用并不多。本文應(yīng)用多面函數(shù)和移動法綜合擬合模型，并以首都圈和晉冀蒙兩個地區(qū)的垂直運(yùn)動為例，與其他幾種模型作對比研究。

1 模型方法數(shù)學(xué)原理

1.1 多面函數(shù)模型

多面函數(shù)模型是由美國Hardy教授于1977年提出的，并建議用于地殼垂直形變分析。該方法的基本思想為：任何數(shù)學(xué)表面和任何不規(guī)則的圓滑表面，總可以用一系列有規(guī)則的數(shù)學(xué)表面的和以任意精度逼近[9]。根據(jù)這一思想，假設(shè)地球表面上有一點(diǎn)(x,y)，其垂直速率值ξ(x,y)可以用下式來表示

(1)

式中，aj為待定系數(shù)；Qx,y,xj,yj是多面函數(shù)的核函數(shù)。核函數(shù)有多種類型，常用的有距離型核函數(shù)和錐面函數(shù)，其中距離型核函數(shù)表達(dá)式為

(2)

式中，δ為平滑因子，通?？扇∫恍≌龜?shù)或0；μ一般取1/2、-1/2、3/2。當(dāng)μ=1/2時，核函數(shù)稱為正雙曲面函數(shù)；當(dāng)μ=-1/2時，核函數(shù)稱為倒雙曲面函數(shù)；當(dāng)μ=3/2時，核函數(shù)稱為三次曲面函數(shù)。

錐面函數(shù)表達(dá)式為

(3)

式中，c為待定參數(shù)。

設(shè)有m個已知水準(zhǔn)點(diǎn)xi,yi，選取其中n(n≤m)個點(diǎn)作為核函數(shù)的中心點(diǎn)xj,yj，令Qij=Q(xi,yi,xj,yj)，則各數(shù)據(jù)點(diǎn)應(yīng)滿足

(4)

由此可列出誤差方程

V=QX-ξ

(5)

根據(jù)最小二乘原理可求得待定系數(shù)X，即

(6)

待定系數(shù)求出后，根據(jù)式(4)可計算測區(qū)各未知點(diǎn)的垂直形變速率值。

需要注意的是，應(yīng)用多面函數(shù)擬合垂直形變速率時，要根據(jù)測區(qū)大小和地形情況確定合適的核函數(shù)、平滑因子和中心點(diǎn)個數(shù)。

1.2 移動曲面模型

移動曲面模型的基本原理是以每一個待定點(diǎn)為中心，選取待定點(diǎn)周圍適量的已知點(diǎn)，擬合出一個多項(xiàng)式曲面，進(jìn)而內(nèi)插出待定點(diǎn)上的函數(shù)值，是一種特殊的局部函數(shù)逼近法。

為了選取鄰近的數(shù)據(jù)點(diǎn)，以待定點(diǎn)Pxp,yp為圓心，以R為半徑作圓，凡落在圓內(nèi)的數(shù)據(jù)點(diǎn)即被選用。應(yīng)用時要根據(jù)實(shí)際情況選擇合適的半徑R。

定權(quán)時要根據(jù)已知點(diǎn)與待定點(diǎn)間的距離進(jìn)行定權(quán)，常采用的權(quán)有以下幾種形式

在進(jìn)行擬合前，可先將參與擬合的各點(diǎn)坐標(biāo)進(jìn)行中心化處理，即將各點(diǎn)坐標(biāo)減去待定點(diǎn)P的坐標(biāo)，這樣轉(zhuǎn)換后就變成以P為原點(diǎn)的坐標(biāo)

(7)

(8)

最后求待定點(diǎn)P的垂直形變速率時，直接代入坐標(biāo)0,0即可，從而在很大程度上方便了計算，而且有利于改善方程的數(shù)值穩(wěn)定性和提高解算精度[11]。

1.3 多面函數(shù)和移動法加權(quán)綜合模型

該方法是基于綜合預(yù)報理論提出的，它并不是一種全新的方法，而是首先采用多面函數(shù)和移動法進(jìn)行單獨(dú)擬合，然后對其擬合值進(jìn)行加權(quán)綜合而得到最終擬合值[12]。用數(shù)學(xué)模型可表示為

(9)

文獻(xiàn)[12]中并沒有給出該模型解的具體解法和表達(dá)式，本文給出W解的具體解法如下：

先構(gòu)造拉格朗日極值函數(shù)

令

(10)

K=WTFTY-FW

(11)

另外，由式(10)可得

(12)

由式(11)和式(12)可采用迭代法解得W

(2) 將W0代入式(11)計算K。

(3) 將K代入式(12)計算W。

(4) 判斷W-W0<10-6,若是，退出循環(huán)；若否，令W0=W繼續(xù)重復(fù)步驟(2)—(4)。

經(jīng)過以上迭代過程即可解得W為W*，則相應(yīng)綜合模型擬合結(jié)果為

(13)

1.4 多面函數(shù)和移動法綜合擬合模型

多面函數(shù)法對整體性變化的擬合效果較好，但需要確定合適的模型參數(shù)，如果選取不當(dāng)，容易出現(xiàn)擬合結(jié)果波動較大的現(xiàn)象，對局部變化擬合效果較差；而移動擬合法是以待定點(diǎn)為中心，以一定距離為半徑，用一個多項(xiàng)式曲面擬合待定點(diǎn)函數(shù)值的局部擬合法，具有較好的靈活性[13]。為此，可以結(jié)合兩者優(yōu)點(diǎn)，采用先進(jìn)行多面函數(shù)擬合，再對殘差進(jìn)行移動法擬合的綜合模型擬合垂直形變速率。這很類似于組合法確定似大地水準(zhǔn)面中的“移去—恢復(fù)”法，具體步驟如下：

(2) 求出m個擬合點(diǎn)的速率擬合殘差v

v=ξN-ξ

(14)

(15)

2 算例與分析

2.1 擬合方案介紹

本文選用首都圈地區(qū)2001年和2005年，以及晉冀蒙地區(qū)2002年和2006年各兩期水準(zhǔn)數(shù)據(jù)進(jìn)行研究，并對兩個地區(qū)的數(shù)據(jù)采用自由網(wǎng)平差法進(jìn)行動態(tài)平差，計算垂直運(yùn)動速率。兩段動態(tài)平差的單位權(quán)中誤差分別為1.03和1.15 mm/a，說明這兩個區(qū)域的數(shù)據(jù)整體實(shí)測精度比較好。剔除平差后速率中個別量值過大的突變點(diǎn)，首都圈地區(qū)共選用719個具有平差后速率的數(shù)據(jù)點(diǎn)(如圖1所示)，晉冀蒙地區(qū)共有627個具有平差后速率的數(shù)據(jù)點(diǎn)(如圖2所示)。由于對一個區(qū)域垂直形變場的擬合研究，主要為了得到未測點(diǎn)的垂直形變速率(即水準(zhǔn)環(huán)內(nèi)部區(qū)域的垂直形變速率)，因此為了更充分地得出結(jié)論，采用以下4種方案分別進(jìn)行擬合計算：

方案1：均勻選取首都圈地區(qū)143個點(diǎn)作為外部檢核點(diǎn)(約占總點(diǎn)數(shù)的1/5)，其余576個點(diǎn)作為擬合點(diǎn)參與函數(shù)模型的參數(shù)擬合(針對路線上數(shù)據(jù))。

方案2：選取首都圈地區(qū)5條水準(zhǔn)環(huán)內(nèi)部路線全部點(diǎn)作為外部檢核點(diǎn)(共80個點(diǎn))，其余639個點(diǎn)作為擬合點(diǎn)參與函數(shù)模型的參數(shù)擬合(針對環(huán)內(nèi)部數(shù)據(jù))。

方案3：均勻選取晉冀蒙地區(qū)125個點(diǎn)作為外部檢核點(diǎn)(約占總點(diǎn)數(shù)的1/5)，其余502個點(diǎn)作為擬合點(diǎn)參與函數(shù)模型的參數(shù)擬合(針對路線上數(shù)據(jù))。

方案4：選取晉冀蒙地區(qū)5條水準(zhǔn)環(huán)內(nèi)部路線全部點(diǎn)作為外部檢核點(diǎn)(共60個點(diǎn))，其余567個點(diǎn)作為擬合點(diǎn)參與函數(shù)模型的參數(shù)擬合(針對環(huán)內(nèi)部數(shù)據(jù))。

本文采用外符合精度及檢核點(diǎn)速率擬合值與實(shí)際觀測值之差(殘差)的絕對值大于2 mm/a的點(diǎn)的個數(shù)來衡量擬合的精度。外符合精度σ的表達(dá)式為

(16)

式中，Δi表示檢核點(diǎn)速率擬合值與實(shí)際觀測值之差；n為檢核點(diǎn)個數(shù)，并稱 Δi>2 mm/a的點(diǎn)為病態(tài)點(diǎn)。

圖1 首都圈地區(qū)垂直形變速率(相對于北京)

圖2 晉冀蒙地區(qū)垂直形變速率(相對于大同)

2.2 各種模型擬合結(jié)果比較

對首都圈地區(qū)的擬合計算有方案1和方案2，各種模型擬合結(jié)果見表1。

表1 首都圈地區(qū)各種模型擬合結(jié)果比較

為了直觀地看出各種模型的擬合效果，給出不同模型在方案1和方案2中擬合檢核點(diǎn)速率殘差情況，如圖3和圖4所示。

圖3 方案1不同模型檢核點(diǎn)擬合殘差

圖4 方案2不同模型檢核點(diǎn)擬合殘差

對晉冀蒙地區(qū)的擬合計算有方案3和方案4，各種模型擬合結(jié)果見表2。

表2 晉冀蒙地區(qū)各種模型擬合結(jié)果比較

為了直觀地看出各種模型的擬合效果，給出不同模型在方案3和方案4中擬合檢核點(diǎn)速率殘差情況，如圖5和圖6所示。

圖5 方案3不同模型檢核點(diǎn)擬合殘差

圖6 方案4不同模型檢核點(diǎn)擬合殘差

2.3 算例結(jié)果分析

(1) 從表1和表2中可以發(fā)現(xiàn)，4種模型中兩種綜合模型的擬合精度優(yōu)于兩種單一模型，但是特性不同。多面函數(shù)和移動法加權(quán)綜合模型雖然精度上優(yōu)于多面函數(shù)模型，但是病態(tài)點(diǎn)個數(shù)并不具有優(yōu)勢，這是因?yàn)榧訖?quán)綜合會使得單一模型的各種誤差傳遞到最終結(jié)果中。而多面函數(shù)和移動法綜合擬合模型不管在精度上還是病態(tài)點(diǎn)個數(shù)上均優(yōu)于多面函數(shù)模型。

(2) 多面函數(shù)和移動法綜合擬合模型在兩個地區(qū)4種不同情況下都是4種模型中擬合精度最高的，病態(tài)點(diǎn)數(shù)也是最少的，并且從圖4—圖7可以看出多面函數(shù)和移動法綜合擬合模型的擬合殘差是最穩(wěn)定的。

2.4 等值線圖繪制

利用多面函數(shù)和移動法綜合擬合模型對首都圈和晉冀蒙地區(qū)的垂直形變速率進(jìn)行擬合，并進(jìn)行格網(wǎng)化處理，格網(wǎng)間距設(shè)為1′×1′，得到兩個地區(qū)的垂直形變速率等值線圖，如圖7和圖8所示。

從圖7和圖8可以發(fā)現(xiàn)首都圈和晉冀蒙兩個地區(qū)用多面函數(shù)和移動法綜合擬合模型得到的垂直形變速率等值線圖與圖1和圖2的速率圖反映的升降情況吻合度較好，說明該模型可以用于這兩個地區(qū)的垂直形變場擬合。

圖7 首都圈地區(qū)垂直形變速率等值線

圖8 晉冀蒙地區(qū)垂直形變速率等值線

3 結(jié) 論

(1) 兩種綜合模型的擬合精度優(yōu)于兩種單一模型，但多面函數(shù)和移動法加權(quán)綜合模型會使得單一模型的各種誤差傳遞到最終結(jié)果中，從而導(dǎo)致病態(tài)點(diǎn)個數(shù)較多。而多面函數(shù)和移動法綜合擬合模型不管在精度上還是病態(tài)點(diǎn)個數(shù)上都優(yōu)于多面函數(shù)模型。

(2) 多面函數(shù)和移動法綜合擬合模型在首都圈和晉冀蒙兩個地區(qū)4種情況下都是擬合精度最高、病態(tài)點(diǎn)數(shù)最少、擬合殘差最穩(wěn)定的，并且應(yīng)用該模型擬合這兩個地區(qū)的垂直形變場得到的等值線圖與速率圖吻合度較好，具有一定的實(shí)際應(yīng)用價值。

[1] 郝明.基于精密水準(zhǔn)數(shù)據(jù)的青藏高原東緣現(xiàn)今地殼垂直運(yùn)動與典型地震同震及震后垂直形變研究[D].北京:中國地震局地質(zhì)研究所，2012.

[2] 黃立人,楊國華,劉天奎.用速率面擬合法研究華北部分地區(qū)的現(xiàn)今地殼垂直運(yùn)動[J].地殼形變與地震，1989,9(3):26-34.

[3] 趙承坤，黃立人.速率面擬合法中核函數(shù)中心點(diǎn)的選擇[J].地殼形變與地震，1991,11(2):48-54.

[4] 陶本藻,王新洲,于正林,等.用于垂直形變模型的多面函數(shù)擬合法的試驗(yàn)研究[J].地殼形變與地震，1992，12(1)：1-12.

[5] 楊國華,黃立人.速率面擬合法中多面函數(shù)幾個特性的初步數(shù)值研究[J].地殼形變與地震，1990,10(4)：70-82.

[6] 黃立人,陶本藻,趙承坤.多面函數(shù)擬合在地殼垂直運(yùn)動研究中的應(yīng)用[J].測繪學(xué)報，1993,22(1):25-32.

[7] 劉經(jīng)南,施闖,姚宜斌,等.多面函數(shù)擬合法及其在建立中國地殼平面運(yùn)動速度場模型中的應(yīng)用研究[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2001，26(6):500-503.

[8] 王文利,陳士銀,董鴻聞,等.利用多面函數(shù)擬合中國陸地垂直運(yùn)動速率圖[J].測繪通報,2002(8):6-8.

[9] 劉同文,楊志強(qiáng),王慧敏.西秦嶺地區(qū)垂直形變的分析研究[J].測繪科學(xué)，2014，39(4)：61-63.

[10] 秦姍蘭,王慶良,季靈運(yùn),等.利用水準(zhǔn)資料研究西秦嶺地區(qū)的垂直形變[J].大地測量與地球動力學(xué)，2012，32(2):16-19.

[11] 徐紹銓,張華海,楊志強(qiáng),等.GPS測量原理及應(yīng)用[M].武漢：武漢大學(xué)出版社，2008.

[12] 趙超英,張勤.地殼垂直形變場的綜合逼近模型[J].大地測量與地球動力學(xué)，2003，23(2)：42-44.

[13] 田曉，鄭洪艷，許明元,等.一種改進(jìn)的適用于不同地形的GPS高程擬合模型[J].測繪通報，2017(1):35-38.