文章編號(hào)1000-5269（2024）05-0026-06 DOI：10.15958/j.cnki.gdxbzrb.2024.05.04

摘要：針對(duì)均值漂移模型在異常測(cè)站剔除過程中出現(xiàn)的“回漂”現(xiàn)象導(dǎo)致粗差未能盡數(shù)剔除、使其擬合結(jié)果有偏于實(shí)際的問題，提出基于懲罰回歸優(yōu)化均值漂移模型對(duì)異常測(cè)站進(jìn)行剔除的方法。利用優(yōu)化后的模型方法對(duì)“中國大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)”華南塊體中東部觀測(cè)的水平速度場(chǎng)數(shù)據(jù)進(jìn)行擬合對(duì)比分析，并利用剔除后的有效測(cè)站進(jìn)行形變特征反演。結(jié)果表明：利用懲罰回歸優(yōu)化的均值漂移模型能更有效地識(shí)別異常測(cè)站點(diǎn)，篩選后的測(cè)站數(shù)據(jù)能顧及塊體整體運(yùn)動(dòng)的平滑性，同時(shí)比均值漂移模型剔除具有更高的擬合精度，利用有效數(shù)據(jù)反演的形變特征也符合該區(qū)域?qū)嶋H運(yùn)動(dòng)趨勢(shì)。

關(guān)鍵詞：地殼運(yùn)動(dòng)；異常測(cè)站；均值漂移模型；懲罰函數(shù)；應(yīng)變特征 中圖分類號(hào)：P227；P313 文獻(xiàn)標(biāo)志碼：A

任何觀測(cè)結(jié)果都會(huì)包含噪聲，全球衛(wèi)星定位導(dǎo)航系統(tǒng)（GlobalNavigationSatelliteSystem，GNSS）監(jiān)測(cè)站在地殼形變速度場(chǎng)連續(xù)觀測(cè)中也8d826b8fce9328f3fba0a6b89d931f199c06d2a64435f34a68f90616a1d1d818不例外。監(jiān)測(cè)站點(diǎn)可能會(huì)受到局部干擾或局部構(gòu)造形變的影響，使同一塊體上的監(jiān)測(cè)站點(diǎn)有時(shí)會(huì)表現(xiàn)出不協(xié)調(diào)的位移［1］。人們把這類不協(xié)調(diào)位移的測(cè)站定義為異常測(cè)站。這類測(cè)站的存在會(huì)影響地殼運(yùn)動(dòng)模型擬合的精度和可靠性。所以，必須對(duì)監(jiān)測(cè)站速度場(chǎng)數(shù)據(jù)進(jìn)行合理有效的篩選，并對(duì)異常測(cè)站進(jìn)行剔除。異常測(cè)站剔除的過程究其原理，可視為粗差的探測(cè)及定位，常用的方法是以最小二乘法為基礎(chǔ)將粗差歸為函數(shù)模型，利用兩倍中誤差法、擬準(zhǔn)檢定法、均值漂移模型等方法進(jìn)行探測(cè)。文獻(xiàn)［2］利用兩倍中誤差法、擬準(zhǔn)檢定法、均值漂移模型進(jìn)行測(cè)站篩選對(duì)比，得出均值漂移抗差估計(jì)模型剔除效果明顯好于擬準(zhǔn)檢定法和兩倍中誤差法。而均值漂移量的搜索方向無法進(jìn)行限制，容易出現(xiàn)“回漂”現(xiàn)象［3］，導(dǎo)致粗差剔除不完整。針對(duì)此問題，本文提出利用平滑剪切絕對(duì)偏差（smoothlyclippedabsolutedeviation，SCAD）懲罰函數(shù)對(duì)均值漂移模型異常測(cè)站剔除方法進(jìn)行優(yōu)化，對(duì)“中國大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)”在華南塊體的GNSS觀測(cè)的速度場(chǎng)數(shù)據(jù)進(jìn)行篩選，并對(duì)篩選出的有效速度場(chǎng)數(shù)據(jù)進(jìn)行擬合分析。結(jié)果表明：此方法能有效識(shí)別異常測(cè)站點(diǎn)，一定程度上克服了“回漂”影響，利用優(yōu)化的均值漂移模型篩選的測(cè)站數(shù)據(jù)進(jìn)行擬合具有更高的擬合精度。

1測(cè)站剔除的方法

1.1基于剛性運(yùn)動(dòng)的地殼形變模型

傳統(tǒng)的塊體運(yùn)動(dòng)模型不考慮塊體內(nèi)部形變，將塊體水平運(yùn)動(dòng)視為繞歐拉軸的有限旋轉(zhuǎn)運(yùn)動(dòng)，其模型為［4-5］

VEVN=r-cosλsinφ-sinλsinφcosφsinλ-cosφ0ωxωyωz（1）

式中：VE表示板塊任意一點(diǎn)在東西方向的運(yùn)動(dòng)速率，VN表示南北方向的運(yùn)動(dòng)速率；r是地球平均半徑；λ為大地經(jīng)度，φ為大地緯度；ωx，ωy，ωz定義為歐拉旋轉(zhuǎn)參數(shù)。相應(yīng)歐拉極參數(shù)計(jì)算公式為

ω=ω2x+ω2y+ω2zl=atan（ωyωx）φ=atan（ωzω2x+ω2y）（2）

式中：ω為歐拉經(jīng)度；l為歐拉緯度；φ為旋轉(zhuǎn)角。利用式（1）對(duì)速度場(chǎng)數(shù)據(jù)進(jìn)行擬合處理，利用觀測(cè)值與擬合值求差，得出速率殘差值，其服從正態(tài)分布，符合隨機(jī)誤差特性，因此利用兩倍中誤差作為閾值進(jìn)行測(cè)站篩選是合理的方法，但其模型本身未能顧及塊體整體位移的平滑性，從而導(dǎo)致剔除效果減弱。

1.2基于均值漂移模型的測(cè)站篩選［6-7］

將式（1）簡(jiǎn)化為測(cè)量平差問題的Gauss-Markov模型

L=AX+ΔΔ～Nn（0，σ20P-1）（3）

設(shè)未知參數(shù)向量估值為X^，L為觀測(cè)值向量，V=AX^－L為誤差方程，則最小平差解為X^=（ATPA）-1ATPL。設(shè)可靠性矩陣為R=I－A（ATPA）－1ATP，且RL=－V。假設(shè)第i個(gè)觀測(cè)值存在異常，則含有粗差的均值漂移模型為

L=AX+r+ΔΔ～Nn（0，σ20P-1）r=（0，…，ri，…，0）（4）

式中：ri為漂移量。根據(jù)最小二乘法得出均值漂移模型的抗差估計(jì)的解：

X^=（ATPA）-1ATP（5）

式中：為等價(jià)觀測(cè)值。

利用迭代法求出未知參數(shù)解，即

X^k+1=（ATPA）-1ATPk+1（6）

觀測(cè)值迭代抗差解計(jì)算表達(dá)式可表示為

k+1=k－rk，rk=（rk1，…，rkn）

rki=0，k－1i≤c

rk-1i，k－1i>c（7）

式中：k－1i為第k－1步的標(biāo)準(zhǔn)化殘差，c值視實(shí)際情況而定，一般取值范圍在0～3。rk－1i為rk－1=-（diag（R））－1k－1的第i個(gè)分量，diag（R）表示R的對(duì)角矩陣，k－1=k－1－AX^k－1。若觀測(cè)值無粗差影響，則RL=－V且V很??；若觀測(cè)值中含有粗差，設(shè)R（L+r）=－，r為異常值，則有Rr－V=-，又因?yàn)閂相比于很小且R－1≈0，故r=-（diag（R））－1，此時(shí)粗差值便可定位。利用均值漂移模型剔除異常測(cè)站具體步驟見文獻(xiàn)［2］。

1.3基于SCAD懲罰回歸的測(cè)站篩選［8-10］

由式（4）確定的異常值識(shí)別量r是稀疏的，因此在利用上述方法檢測(cè)異常測(cè)站時(shí)容易出現(xiàn)“回漂”現(xiàn)象，故利用SCAD懲罰函數(shù)對(duì)其漂移量施加懲罰，得到的目標(biāo)函數(shù)如式（8）：

FnX，r;λn=‖Y－AX－r‖22+n∑ni=1Pλnri（8）

式中：‖·‖2為l2范數(shù)；∑ni=1Pλn（ri）為懲罰函數(shù)項(xiàng);n為樣本數(shù)量，即觀測(cè)站數(shù);λn>0是懲罰函數(shù)中的閾值參數(shù)。

由文獻(xiàn)［9］可知，SCAD懲罰函數(shù)可通過局部線性近似擬合，即Pλn（x）≈P′λn（x^（0））x，x^（0）為x的初始估計(jì)量，P′λn（x）為Pλn（x）的一階導(dǎo)數(shù)，定義如式（9）：

P′λn（x）=λnI（x≤λn）+（aλn－x）+（a－1）λnI（x>λn）（9）

式中：x>0，a>2，I（·）為示性函數(shù)，（aλn-x）+=max（0，aλn-x）。利用式（9）對(duì)式（8）進(jìn)行擬合后可得加懲罰的目標(biāo)函數(shù)如式（10）：

Fn（X，r;λn）=‖Y-AX-r‖22+n∑ni=1P′λn（r^（0）i）ri（10）

式中：（0）i為ri的初始估計(jì)值，X和r的估計(jì)量如式（11）。

（X^，r^）=argminLnX∈Rp-1，Δ∈Rn（X，r;λn）（11）

基于SCAD懲罰函數(shù)優(yōu)化的均值漂移模型的參數(shù)解算采用坐標(biāo)下降法，解算步驟如下：

第1步：設(shè)置初始值（X^（0），r（0）），令k=1。

第2步：令r（k）=argminHn（X^（k－1），r；λn1，λn2），利用公式（12）確定懲罰函數(shù)Hn的最小值，并設(shè)定r^（k）=（r^k1，…，r^kn）T，令i=1，…，n。

r^（k）i=argminHn（X^（k－1），（r^k1，…，r^ki-1，ri，r^（k-1）i-1）T;λn1，λn2）

（12）

第3步：令X^（k）=argminHn（X，r^k；λn1，λn2），利用公式（13）確定懲罰函數(shù)Hn的最小值，并設(shè)定X^（k）=（X^k0，…，X^kn）T，令j=1，…，p。

X^（k）i=argminHn（（X^k0，…，X^kj-1，Xj，X^（k-1）j-1，…，X^（k-1）p）T，r^（k）;λn1，λn2）

（13）

第4步：令k=k+1，重復(fù)第2步和第3步，直至漂移量為0，即為正常測(cè)站，不為0的為異常測(cè)站，剔除。

2測(cè)站剔除及形變特征分析

本文采用1998—2014年間“中國大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)”觀測(cè)的華南塊體水平速度場(chǎng)數(shù)據(jù)，通過解算得到的華南塊體中東部地區(qū)在ITRF2008框架下統(tǒng)一速度場(chǎng)數(shù)據(jù)精度優(yōu)于2mm/a，該區(qū)域監(jiān)測(cè)站分布及水平速度場(chǎng)如圖1、圖2所示［11-12］。從圖2可以看出，此區(qū)域常年運(yùn)動(dòng)區(qū)域穩(wěn)定，以6.783mm/a的速度向東移動(dòng)，故以剛性運(yùn)動(dòng)模型為基礎(chǔ)進(jìn)行測(cè)站篩選。

為驗(yàn)證本文方法，用以下3種方法對(duì)異常測(cè)站進(jìn)行剔除：

方案Ⅰ：利用塊體的剛性運(yùn)動(dòng)模型對(duì)未進(jìn)行篩選的水平速度場(chǎng)數(shù)據(jù)進(jìn)行擬合分析。

方案Ⅱ：利用均值漂移模型對(duì)異常測(cè)站進(jìn)行剔除，用塊體的剛性運(yùn)動(dòng)模型對(duì)有效數(shù)據(jù)進(jìn)行擬合分析。

方案Ⅲ：利用懲罰回歸優(yōu)化的均值漂移模型對(duì)異常測(cè)站進(jìn)行剔除，用塊體的剛性運(yùn)動(dòng)模型對(duì)有效數(shù)據(jù)進(jìn)行擬合分析。

3種方案的擬合結(jié)果列于表1，同時(shí)給出方案Ⅱ和方案Ⅲ各自的異常測(cè)站分布點(diǎn)圖及速率殘差，見圖3—圖7；利用塊體線性應(yīng)變模型對(duì)方案Ⅲ擬合的結(jié)果進(jìn)行應(yīng)變特征反演，如圖8所示。

由計(jì)算結(jié)果得出如下分析：

1）利用方案Ⅱ、方案Ⅲ對(duì)異常測(cè)站進(jìn)行剔除，應(yīng)剔除的異常測(cè)站點(diǎn)數(shù)分別為23和31。從圖3、圖4可以看出：優(yōu)化后的均值漂移模型剔除的異常測(cè)站分布與未優(yōu)化前相比多剔除了8個(gè)；但在分布上，與均值漂移模型保持一致，間接證明了利用懲罰回歸剔除異常測(cè)站的可行性。

2）從圖5、圖6、圖7可以看出：利用未進(jìn)行異常測(cè)站剔除的速度場(chǎng)的速率殘差在東西方向和南北方向最大值均在4mm/a左右，未經(jīng)優(yōu)化的均值漂移模型剔除異常測(cè)站后的擬合結(jié)果顯示在東西方向最大為2.7mm/a，南北方向上的速率殘差最大值為3mm/a，而經(jīng)懲罰函數(shù)優(yōu)化后的均值漂移模型剔除異常測(cè)站后的擬合結(jié)果顯示在東西方向和南北方向上的速率殘差最大值均在2mm/a左右；對(duì)比方案Ⅱ、方案Ⅲ，無論從東西方向還是從南北方向的擬合，方案Ⅲ殘差最大值均小于方案Ⅱ；從表1中的速率殘差標(biāo)準(zhǔn)差可以看出方案Ⅰ未經(jīng)測(cè)站剔除的擬合精度最差，方案Ⅲ剔除后擬合精度最優(yōu)；另外，利用方案Ⅲ解算的歐拉參數(shù)值和歐拉矢量值最趨近于方案Ⅰ，即趨近于未進(jìn)行異常測(cè)站剔除的擬合結(jié)果，更好地說明方案Ⅲ更能顧及塊體整體運(yùn)動(dòng)的平滑性。綜上所述，無論是從速率殘差殘差圖對(duì)比，還是從殘差標(biāo)準(zhǔn)差及解算參數(shù)的對(duì)比，方案Ⅲ對(duì)異常測(cè)站剔除的效果最優(yōu)，即利用懲罰回歸優(yōu)化的均值漂移模型對(duì)速度場(chǎng)數(shù)據(jù)處理的效果更好。

3）從圖8可以看出，華南塊體中東部地區(qū)存在不同程度、不同方向的壓縮和擴(kuò)張。其主壓應(yīng)變從西至東逐漸減弱，由南北方向按逆時(shí)針轉(zhuǎn)動(dòng)；其主張應(yīng)變從西至東、從北至南逐漸擴(kuò)大，從北至南呈順時(shí)針方向轉(zhuǎn)動(dòng)，此結(jié)果與張靜華、李延興等的研究結(jié)果大體一致［13］。從研究區(qū)域地理位置分析，該區(qū)域所處歐亞板塊，受印度板塊擠壓，整體向東移動(dòng)；另外，該區(qū)域北邊界還受到華北塊體擠壓，南邊界受到南海塊體自南向北的推擠力；在東南方向受到菲律賓海板塊向西北方向的仰沖作用力。因此，內(nèi)部形變變化復(fù)雜。由此可見，圖8反演結(jié)果接近于實(shí)際。

3結(jié)論

GNSS測(cè)站的篩選與剔除是研究地殼形變時(shí)不可缺少的重要環(huán)節(jié)。本文針對(duì)GNSS的對(duì)地觀測(cè)速度場(chǎng)數(shù)據(jù)出現(xiàn)的異常測(cè)站，提出利用SCAD懲罰函數(shù)改進(jìn)均值漂移模型的抗差估計(jì)對(duì)異常測(cè)站進(jìn)行剔除的方法。均值漂移模型雖然能有效地檢測(cè)粗差位置，但不可避免地會(huì)出現(xiàn)“回漂”現(xiàn)象，而SCAD懲罰函數(shù)改進(jìn)的均值漂移模型能更準(zhǔn)確地解決線性模型數(shù)據(jù)中異常值的檢測(cè)問題，在一定程度上克服了一般方法未能克服的隱差現(xiàn)象及均值漂移模型的“回漂”現(xiàn)象，對(duì)后續(xù)地殼形變模型的擬合、塊體應(yīng)變反演的精度大有提升，也為地殼形變異常測(cè)站的剔除提供了一種新的思路。

參考文獻(xiàn)：

[1]

李亞萍，孫付平，朱新慧，等.抗差主成分估計(jì)在板塊運(yùn)動(dòng)參數(shù)解算中的應(yīng)用［J］.測(cè)繪工程，2016，25（8）：38-41.

［2］范成成，張俊，雷前坤，等.均值漂移模型在地殼運(yùn)動(dòng)異常測(cè)站剔除中的應(yīng)用［J］.測(cè)繪通報(bào)，2020（5）：80-84.

［3］高向東，黎揚(yáng)進(jìn)，劉秀航，等.改進(jìn)均值漂移算法的焊縫特征點(diǎn)識(shí)別分析［J］.華南理工大學(xué)學(xué)報(bào)（自然科學(xué)版），2019，47（4）：132-137.

［4］李延興，黃珹，胡新康，等.板內(nèi)塊體的剛性彈塑性運(yùn)動(dòng)模型與中國大陸主要塊體的應(yīng)變狀態(tài)［J］.地震學(xué)報(bào)，2001，23（6）：565-572.

［5］王東東.基于GPS速度場(chǎng)的地殼彈塑性形變模型及應(yīng)用［D］.貴陽：貴州大學(xué)，2018：6-7.

［6］張?zhí)教剑畞喞?，鐘先?基于均值漂移模型的異常值檢測(cè)方法［J］.上海理工大學(xué)學(xué)報(bào)，2018，40（2）：116-120.

［7］范成成.顧及先驗(yàn)信息的地殼形變分析模型研究［D］.貴陽：貴州大學(xué)，2020：12-13.

［8］ZOUH，LIRZ.One-stepsparesestimatesinnonconcavepenalizedlikelihoodmodels［J］.AnnalsofStatistics，2008，36（4）：1509-1533.

［9］潘瑩麗，劉展，宋廣雨.基于SCAD懲罰回歸的異常值檢測(cè)方法［J］.統(tǒng)計(jì)與決策，2022，38（4）：38-42.

［10］萬超.基于SCAD稀疏約束的字典學(xué)習(xí)算法及其在信號(hào)處理中的應(yīng)用［D］.廣州：廣東工業(yè)大學(xué)，2021.

［11］ZHAOB，HUANGY，ZHANGCH，etal.CrustaldeformationontheChinesemainlandduring1998—2004basedonGPSdata［J］.GeodesyandGeodynamics，2015，6（1）：7-15.

［12］石玉濤，高原.華南塊體中部上地殼剪切波分裂特征［J］.地球物理學(xué)報(bào)，2022，65（9）：3268-3279.

［13］張靜華，李延興，郭良遷，等.華南塊體的現(xiàn)今構(gòu)造運(yùn)動(dòng)與內(nèi)部形變［J］.大地測(cè)量與地球動(dòng)力學(xué)，2005（3）：57-62.

（責(zé)任編輯：曾晶）

Abstract：

Inviewoftheproblemthatthe“backdrift”phenomenonofthemeandriftmodelintheprocessofanomalousstationexclusionleadstothefailuretoeliminateallthegrosserrors，andthefittingresultsarebiasedtoreality，amethodofeliminatingabnormalstationsbasedonpenaltyregressionoptimizationmeandriftmodelisproposed.TheoptimizedmodelmethodwasusedtofitandcomparethehorizontalvelocityfielddataobservedinthecentralandeasternpartsoftheSouthChinablockbythe“Chinesemainlandtectonicenvironmentmonitoringnetwork”，andthedeformationcharacteristicswereinvertedbyusingtheeffectivestationsafterelimination.Theresultsshowthatthemeandriftmodeloptimizedbypenaltyregressioncanidentifyabnormalstationsmoreeffectively，andthefilteredstationdatacantakeintoaccounttheoverallmotionsmoothnessoftheblockandhaveahigherfittingaccuracythanthemeandriftmodel.Thedeformationcharacteristicsinvertedbyusingeffectivedataalsoconformtotheactualmovementtrendoftheregion.

Keywords：

crustalmovement;anomalystation;meanshiftmodel;penaltyfunction;straincharacteristics

收稿日期：2023-10-08

基金項(xiàng)目：貴州省青年科技人才成長(zhǎng)項(xiàng)目（黔教合KY字［2022］132號(hào)，黔教合KY字［2022］135號(hào)）；貴州工程應(yīng)用技術(shù)學(xué)院一流專業(yè)建設(shè)項(xiàng)目（ZY202104）

作者簡(jiǎn)介：范成成（1995—），男，講師，碩士，研究方向：板塊運(yùn)動(dòng)及地殼形變分析理論，E-mail：970864235@qq.com.

*通訊作者：范成成，E-mail：970864235@qq.com.