顧及抗差性的GNSS多尺度地殼運動場構建方法

蘇小寧，朱慶，孟國杰，孫海麗，閆浩文

1 蘭州交通大學測繪與地理信息學院，蘭州 730070 2 西南交通大學地球科學與環境工程學院，成都 611756 3 中國地震局地震預測研究所，北京 100036 4 首都師范大學資源環境與旅游學院，北京 100048

0 引言

GNSS(Global Navigation Satellite System)觀測技術在研究全球地殼運動和構造活動等方面發揮了重要的作用，“中國地殼運動觀測網絡”、“中國大陸構造環境監測網絡”、首都圈GNSS觀測網和省屬地震氣象測繪部門的CORS(Continuously Operating Reference Stations)網等覆蓋全國和區域的多源觀測數據，為獲取中國大陸較高空間分辨率的變形場提供了豐富的資料(Wang and Shen，2020).由于GNSS速度場相對于某個參考基準，因此選擇不同參考基準導致速度場的表現形式差異巨大.而應變率場與參考基準選擇無關，通過對比應變積累和地震矩釋放的關系，在地震危險性分析和長期預報方面發揮著重要作用(Ojo et al., 2021)，這其中的關鍵是解算出能反映測站空間不規則分布特征且精準的應變率場.

目前解算應變場的方法可總結為兩類：一類是建立塊體-斷層模型，計算塊體內部均勻應變，稱為物理模型(申重陽等，2002；Tape et al., 2009)；另一種是利用GNSS速度場，構建空間連續的速度場模型，進而解算應變率場，稱為數學模型(Tape et al., 2009).可見，在數學模型中關鍵是如何精準構建空間連續的速度場模型，常用的方法有距離加權(Shen et al，1996)、多面函數(武艷強等，2009)、Kriging插值(劉曉霞等，2014)、最小二乘配置(江在森和劉經南，2010)和多尺度球面小波(Tape et al., 2009；程鵬飛等，2015；蘇小寧等，2016)等，這些方法在確定觀測值權重時或將GNSS速度場看作為等精度的觀測值，或利用擬合的速度場形式誤差表征其不等精度的特征.但實際上經常存在由于觀測數據質量較差導致的誤差較大或者土層觀測點標墩不穩定產生的局部速度差異，也就是在觀測數據中存在少量粗差.此時利用現有方法構建速度場模型時，少量粗差將會對解算的應變場產生重要的影響，即上述方法均不具備抗差功能.

本文基于能夠體現GNSS測站不規則分布特征的多尺度球面小波方法構建GNSS速度場，獲取不同空間尺度的地殼運動場，在此基礎上通過抗差迭代算法顧及速度場中存在的少量粗差對建立速度場模型的影響，通過在實測數據中人為加入粗差的模擬計算驗證方法的抗差性，結合川滇地區最新的GNSS速度場數據給出了區域的速度場和應變率場模型.

1 基于球面小波算法的多尺度地殼運動場建模

1.1 速度場函數模型

球面上觀測點的三維速度可以以南北向、東西向和垂直向3個分量來表示，本文只考慮水平向速度場，因此只涉及到南北向和東西向2個分量，考慮垂直向時的方法與本文所述相同.地面觀測站(λ，φ)的地殼運動速度ν(λ，φ)可表示為

(1)

基于球面小波算法構建多尺度地殼運動場模型的本質是采用小波基函數作為估計速度場的框架，利用對小波母函數進行平移和伸縮等仿射變換后得到的小波基函數來表示不同分量的速度場，因此公式(1)可以變換為

(2)

(3)

1.2 空間球面小波基函數

以高斯函數差建立的球面小波基函數(Difference of Gaussians，DOG)在構建多尺度地殼運動場模型時得到了廣泛應用(蘇小寧等，2016)，DOG小波的函數模型為(包伯成等，2009)

(4)

1.3 空間球面小波分解

空間球面小波分解的原則是根據觀測站的分布建立一系列球面小波基函數.通過對球面進行三角形分解，可以獲得在球面上均勻分布的球面三角形網格，分解的層次與球面小波基函數的尺度相對應，也對應于不同邊長的球面三角形，具體形式可參考 Wang和Dahlen(1995).首先根據分解尺度建立空間規則分布的小波基函數，再根據地表觀測點的空間密度對規則分布的小波基函數進行篩選，使得最終的小波基函數集在空間分布上與地表觀測站的空間密度密切相關，進而構建的模型能夠凸顯測站不規則分布的特征.分解的具體步驟如下：

(1)選取研究區經緯度的覆蓋范圍形成球面矩形，計算球面矩形對角線的長度，搜索小于該長度且距離最近的三角形邊長，確定最小分解尺度；

(2)根據不同分解尺度，建立多層次且規則分布的空間小波基函數，標注對應的球面三角形邊長；

(3)根據地表觀測點分布，對所有小波基函數逐個搜索，以小波基函數所在位置為圓心和分解尺度對應的球面三角形邊長為半徑，確定小波基函數所對應的“空間支撐”范圍，統計“空間支撐”內地面觀測站的個數，篩選出“空間支撐”內部多于3個站點的小波基函數，獲得模型可用的最大分解尺度、空間不均勻分布的小波基函數及其數量；

(4)統計分析不同分解尺度小波基函數的有效覆蓋范圍，選擇有效覆蓋范圍超過70%的最大分解尺度為最終的球面小波分解尺度，形成在空間不規則分布且不同分解尺度的小波基函數集.

1.4 觀測方程病態問題及其最優化解

根據地表觀測站的空間分布建立的球面小波基函數個數為N個，因此與其對應的系數，也就是待求參數的個數也是N個，而觀測點的個數為M個，也就是觀測量個數為M.在實際建立速度場模型時，觀測點個數往往小于待求參數個數，導致觀測方程的參數求解成為對欠定方程的求解問題.同時，由母函數經過仿射變換得到的DOG小波基函數為非正交基，導致法方程秩虧，涉及到不適定方程的求解問題.

在建立速度場模型時，可以通過最優化算法中的正則化方法，求解出待求參數向量，目標函數和參數的最優解分別為

(5)

m=(GTPDG+λ2I)-1GTPDd,

(6)

式中S(m)為目標函數，通過以目標函數為最小準則獲取參數的最優解.G為由小波基函數組成的系數矩陣，m為待求參數向量，d為地表GNSS觀測點觀測的速度場向量，CD為觀測數據的方差-協方差矩陣，其逆矩陣為權矩陣PD，考慮到不同位置測站計算的地殼運動速度為獨立觀測量，則CD和PD均為對角矩陣.CM為待求參數的方差-協方差矩陣，本文采用Tikhonov正則化方法(徐天河和楊元喜，2003)，CM的具體形式為λ2I，I為單位對角陣，λ為正則化參數.利用公式(6)，通過多次迭代之后，則可以求解出參數m的平差值.可以看出，由于對法方程矩陣的對角元素進行了修正，顯然該平差值為有偏估計值.

從公式(6)可以看出，在參數最優解求解的過程中，如果已知基函數分布和最優化解算方法，則參數的解只依賴于如何確定觀測數據的權矩陣.以往普遍采用擬合的地表測站GNSS速度場的形式誤差給定，但如果在觀測數據中存在異常值或者粗差值，利用形式誤差確定權矩陣的方法并不能精確地表達其隨機特征，此時公式(6)不具備抗差功能，則無法獲得精確的參數解，最終建立的速度場模型值得商榷，因此，有必要對公式(6)進行修正，使其具備抗差功能.

2 顧及抗差迭代算法的多尺度地殼運動場模型構建

一個好的估計方法不僅要具有“穩健性”，還得具有“抗干擾性”.“穩健性”是指當假定模型與實際模型有微小差異時，其估計方法的性能只受到微弱的影響；“抗干擾性”是指當觀測樣本中混入少量粗差時，估計值受到的影響不大(楊元喜等，2002).地表不規則分布的地殼運動速度場屬于獨立觀測，且不等精度，適合采用按等價權原理作抗差估計，目的是解決不等精度獨立觀測值權函數的確定問題(周江文，1989；楊元喜，1994).

2.1 顧及抗差性的最優化算法原理

在公式(5)所示最優化算法的目標函數中，第一項為觀測數據殘差的加權平方和，顯然個別異常大殘差的出現會導致平方和顯著增大.為了達到平方和最小的目標函數，平差值必然會遷就那些異常觀測值或者粗差值，使得平差值偏離最優值，因此個別異常值或粗差會對整個估值產生較大的影響.抗差估計的本質是利用增長較慢的函數代替平方和函數，得到具有較好抗粗差性的估值.

(7)

2.2 顧及抗差性的迭代最優化算法實現步驟

(1)確定待估參數的初值，確定權矩陣的初值，給出觀測數據的方差-協方差矩陣；

(2)按照公式(5)和公式(6)所示的傳統最優化算法估計參數的平差值，計算觀測數據與模型值的殘差；

(5)對參數平差值與初值的差值進行評價，如果其符合限差要求，則結束求解；否則重復步驟(2)—(4)進行迭代求解，直至待求參數平差值的差值符合要求為止.其中第k+1次的迭代解為

(8)

3 川滇地區地殼運動場模型

3.1 多源數據融合的川滇地區GNSS震間速度場

中國大陸觀測地殼運動的GNSS測站主要來自于中國大陸構造環境監測網絡，該網絡包括260個連續觀測站和2000個流動觀測站，其中32個連續觀測站和1056個流動觀測站觀測起始于1998年，其余測站觀測起始于2009年，所有的數據均已經獲得了超過10年的觀測資料.除此之外，在一些地震活動性較強的區域，例如首都圈地區、山西地塹帶和川滇地區等，也進行了GNSS流動觀測網的加密觀測；省屬的地震、氣象和測繪部門為了滿足區域性CORS網的需求，在所屬范圍也布設了連續觀測站，這些測站中的一部分能夠滿足地殼運動觀測的需求.Wang和Shen(2020)對以上數據進行了收集和整理，同時也收集了中國大陸周邊的GNSS觀測資料，通過高精度數據處理與震間形變場提取，給出了中國大陸高空間分辨率的GNSS震間速度場.數據處理采用GAMIT/GLOBK軟件(Herring et al., 2015, 2018)，所有數據解算采用統一的策略，確保解算的坐標參數的高精度和強自洽性.在GNSS原始觀測數據解算階段，將地面測站坐標、地球定向參數(極移和協調時間時)和衛星軌道一并估算，有效地消除了由于地球物理模型或者軌道參數變化所帶來的系統性偏差(Wang and Shen, 2020).在時間序列建模階段，同時考慮了大地震同震和震后位移的影響；考慮到不同天線類型和可能的天線設置誤差，對GAMIT給出的測站坐標誤差進行了重新分配，并利用有色噪聲隨機模型構造觀測數據的方差-協方差矩陣，體現出時間序列中的時間相關性特征(Tao et al., 2021).

在Wang和Shen(2020)給定的GNSS速度場中，中國大陸構造環境監測網絡的數據截止于2015年，本文采用與其一致的觀測數據解算策略、時間序列建模策略和地殼運動解算策略，補充處理了該觀測網2017、2018、2019和2020年的觀測數據，增加了觀測數據的長度，也提高了速度場的穩定性，更新的川滇地區GNSS速度場如圖1所示.

圖1 川滇地區相對于歐亞參考框架的GNSS速度場，誤差橢圓置信區間為95%Fig.1 GNSS velocity field with respect to Eurasian reference frame in the Sichuan and Yunnan regions. The confidence interval of error ellipse is 95%

圖2 川滇地區GNSS測站分布的球面小波分解尺度Fig.2 Decomposition scale of spherical wavelet based on the distribution of GNSS stations in the Sichuan and Yunnan regions

3.2 基于球面小波算法的多尺度地殼運動場模型

根據GNSS測站分布，依據空間球面小波分解原則，獲得川滇地區的球面小波分解尺度，如圖2所示.研究區內最大分解尺度的最小值為6，最大分解尺度的最大值為9.除西南角和東南角的邊緣區域，其他區域的最大分解尺度均在7以上.最大分解尺度能達到9的區域主要沿著主要構造帶分布，例如小江斷裂帶南端、紅河斷裂與麗江—小金河斷裂交匯的區域和鮮水河斷裂與龍門山斷裂南端交匯的三岔口區域.川滇菱形塊體覆蓋的主要區域最大分解尺度都能達到8，因此本文最終選擇的最大分解尺度為8，分解尺度的選擇范圍為3～8.

正則化參數的求解采用留一交叉驗證法(Leave-one-out Cross Validation，LOOCV)，首先去除一個觀測點數據，計算待求參數，然后計算該點觀測值與模型估值之間的差值，通過全局搜索，選擇出最優的正則化參數，進而求解出待求參數的最優解,正則化參數的搜索過程如圖3所示.選擇1～1010作為搜索范圍，可以看出該過程搜索到了函數的全局最小值，進而提取出正則化參數λ的取值，圖3中橫坐標取正則化參數λ的常用對數，其在東西分量和南北分量分別為533和849.

圖3 基于LOOCV算法的最優化算法正則化參數Fig.3 Selection of regulation parameters based on optimized LOOCV method

確定球面小波分解尺度和正則化參數之后，通過多次迭代即可求解出待估參數的最優化值，進而正演出研究范圍內空間任意一點的地殼運動速度，實現地殼形變場建模，給出空間連續的地殼運動圖像，如圖4所示，圖中冷色(藍色)代表小速率值，暖色(紅色)代表大速率值.可以看出，安寧河—則木河—小江斷裂組成的構造帶是川滇菱形塊體地殼運動速率的分界帶.

圖4 川滇地區地殼運動場模型(a) 東西分量; (b) 南北分量.Fig.4 Crustal deformation model in the Sichuan and Yunnan regions(a) The EW component; (b) The NS component.

3.3 顧及抗差迭代算法的多尺度地殼運動場模型

抗差迭代算法可有效削弱異常觀測值和粗差對地殼形變場建模的影響，進而獲得準確的應變率場.通過多次抗差迭代計算，在所有參與建模的367個GNSS測站中，檢測出74個降權觀測值，占比為20%，檢測出7個除權觀測值，占比為2%，圖5中藍色三角形站點為正常觀測站點，黑色三角形站點為降權觀測站點，紅色三角形站點為除權觀測站點.

圖5 抗差迭代算法給出的川滇地區三類GNSS測站分布圖Fig.5 The distribution of three types of GNSS stations in Sichuan and Yunnan regions based on robust iterative algorithm

圖6給出了利用抗差迭代算法對部分測站進行降權和除權之后，基于多尺度球面小波算法構建的川滇地區地殼運動場模型，可以看出考慮抗差性前后計算的速度場空間分布形態基本一致，安寧河—則木河—小江斷裂是川滇菱形塊體的東邊界這一基本特征未發生明顯改變.為了進一步分析抗差迭代算法在構建地殼形變場時的作用，對本文顧及抗差性的方法與傳統多尺度球面小波方法獲得的地殼運動場模型進行對比，結果如圖7所示，圖中也給出了降權站點和除權站點的分布.可以看出，二者差異的量值在1.0 mm·a-1以內，說明雖然在構建速度場模型時考慮了抗差性，但由于原有速度場本身數據質量較好，因此得到的結果是穩定的，二者的差異較小，而差異性較大的區域總是覆蓋著被降權或者除權的測站.

4 討論

4.1 GNSS速度場中存在粗差測站的原因

利用目視法對基于GNSS觀測資料解算的中國大陸地殼運動場進行分析，發現不可避免地存在一些與周圍測站的運動特征明顯不協調的測站.其產生的原因主要包括：(1)因觀測時間跨度太短而無法擬合出真正的震間運動；(2)觀測標墩建于土層之上，由于季節變化等因素導致標墩本身存在非線性的變化；(3)測站觀測環境的改變導致觀測噪聲增加或運動特征偏離真值.

本文所采用的觀測數據均超過了10年，超過解算震間變形時對時間跨度4.5年的要求(Blewitt et al., 2013).觀測環境的改變在測站的日志文件中都會有明確記載，在數據預處理階段已經被考慮到.因此，為了分析測站被降權處理和被除權處理的原因，本文重點對這些測站的觀測標墩和數據質量進行了進一步分析.其中測站標墩分為基巖點和土層點，通常認為基巖點能夠精準代表地殼的運動特征，而土層點會受到土壤隨著季節變化的熱脹冷縮和土層本身不穩定等因素的影響(譚偉杰等，2017)，產生局部效應而可能無法精準代表地殼的運動特征.被降權和除權處理的81個站點全部為土層站點，可見觀測結果中存在明顯的局部影響.通過觀測數據形式誤差的大小體現出數據質量的差異性，分析權不變、降權和除權三類測站形式誤差的差異性.計算三類數據東西向和南北向形式誤差的均值，權不變、降權和除權三類數據東西向形式誤差的均值分別為0.61、0.70和0.93 mm·a-1，南北向形式誤差的均值分別為0.52、0.60和0.83 mm·a-1，可以看出三類數據的形式誤差存在差異，表明其時間序列的離散度逐漸增大，數據質量逐漸降低.

圖6 顧及抗差迭代算法的川滇地區地殼運動場模型(a) 東西分量; (b) 南北分量.Fig.6 Crustal deformation model in Sichuan and Yunnan regions considering robust iterative algorithm(a) The EW component; (b) The NS component.

圖7 顧及抗差迭代算法與傳統方法構建的地殼運動場在東西分量的差值Fig.7 The difference of crustal deformation model in the EW component between our and traditional method

4.2 抗差迭代算法的抗差性分析

為了進一步檢驗抗差迭代算法的抗差性能，在速度場空間變化較大的川滇地區人為地引入粗差，分析粗差引入前后速度場模型的差異性特征.引入粗差的方式有在單個測站引入粗差和在相鄰的幾個測站引入粗差兩種，具體如圖8所示，其中測站H335、H387和H380為3個在區域上獨立的含有粗差的測站，測站H346、H111、H112、SCYY和H103為在區域上具有相關性的5個含有粗差的測站.用含有粗差的測站構建運動場模型，并與不含粗差的數據構建的結果進行對比，分析本文算法的抗差性.

圖8 速度場中單個測站和區域測站加入人為誤差的兩種方案Fig.8 Two schemes of adding gross errors artificially of individual station and regional stations in velocity field

圖9 加入人為粗差前后東西向速度場模型的差異分布Fig.9 The difference of crustal deformation model in EW component before and after adding gross errors artificially

圖9給出了對二者速度場求差的結果，圖中-0.5～0.5 mm·a-1以內的值用同一個色標表示，可以看出，其中測站H335和H387兩個單獨的存在粗差的測站經過抗差迭代算法之后，對建模的結果并未產生影響，而另一個獨立的H380測站的粗差卻對建模結果產生了明顯的影響，進一步對測站分布分析可以看出，在H335和H387兩個測站10 km以內的范圍還有其他測站，而H380測站在10 km以內的范圍內沒有其他測站，也就是說抗差迭代算法從測站H335和H387檢測出了粗差，而從H380測站沒有檢測到粗差，把含有粗差的觀測值當作好的觀測值對待，這是由于該測站周邊沒有其他觀測站作為對比，即便其本身含有粗差也被當做了有效值.由測站SCYY、H346、H103、H111和H112等5個測站組成的含有粗差的區域對建模結果也產生了明顯的影響，表明區域性存在的粗差會對建模結果產生影響.

4.3 抗差迭代算法對解算應變率場的貢獻

從圖10可以看出，川滇地區重要應變積累的區域沿著主要的活動構造帶分布，其中主要的一個條帶是安寧河—則木河—小江斷裂帶，這也是川滇菱形塊體的東邊界，紅河斷裂與麗江—小金河斷裂交匯的區域是另外一個重要的應變積累區域，該分布特征與已有的研究成果保持一致(Qu et al., 2018)，但在應變集中區的量值上普遍高于已發表的結果(Pan and Shen, 2017；李長軍等，2019)，這主要是由于本文采用的多尺度球面小波的應變計算方法能夠凸顯出測站不均勻分布的特征，獲取不同空間分辨率的應變率場，進而提取出應變積累的集中區域.圖10也給出本文方法與傳統多尺度球面小波方法解算的應變率場的差異，從第二不變量的量值來看，應變差異值的量值遠低于應變的量值，表明兩個結果均是穩健的，對比圖10和圖7可以看出，應變率差異較大的區域也是速度場差異較大的區域，同時也是存在被降權或除權測站的區域.

4.4 川滇地區多尺度應變積累特征

基于多尺度球面小波算法解算應變率場的優點之一是能夠獲取不同空間尺度的應變率場，從而體現出不同規模構造運動的應變積累特征.小的分解尺度能夠獲取大規模構造運動的應變積累，例如塊體邊界的應變積累；大的分解尺度能夠獲取小規模構造運動的應變積累，例如塊體邊界應變積累的差異性分布或塊體內部次級斷層的應變積累.圖11分別給出了分解尺度為3～7和8的應變率第二不變量，二者之和為區域的總應變，總體特征與已有結果保持一致(Jin et al, 2019; Zhang et al., 2019).從分解尺度為3～7的應變率場可以看出，應變積累的主要區域為塊體運動的邊界區域，作為川滇菱形塊體東邊界的安寧河—則木河—小江斷裂帶，其是川滇地區主要的應變積累區域，其中以安寧河斷裂和則木河斷裂為東邊界的大涼山次級塊體內部存在應變積累，沿著小江斷裂的應變積累穿過紅河斷裂之后一直向西南方向延伸.其他主要的應變積累區域是位于紅河斷裂與麗江—小金河交匯的區域和中緬邊界的幾條西南向斷裂分布區域，主要包括畹町斷裂和南汀河斷裂等.分解尺度為8的應變積累分布顯示區域內除了沿著主要斷裂的應變積累外，還存在局部小區域的應變積累，這些應變積累主要體現了沿著斷裂應變積累的差異性特征.

圖10 川滇地區應變率第二不變量(a) 本文方法解算的應變率場; (b) 本文方法與傳統多尺度球面小波方法解算結果的差異.Fig.10 The second invariant of strain rates in Sichuan and Yunnan regions(a) The strain rate field calculated by our method; (b) The difference between our and traditional method.

圖11 川滇地區應變率第二不變量(a) 分解尺度為3～7; (b) 分解尺度為8.Fig.11 The second invariant of strain rates in Sichuan and Yunnan regions(a) The decomposition scale is from 3 to 7; (b) The decomposition scale is 8.

5 結論

利用長期積累的GNSS觀測資料經過高精度統一數據處理提取地殼運動圖像，然后解算GNSS應變率場，進而對應變積累與地震矩釋放進行對比分析，這些是開展區域地震危險性分析的主要工作，并在實際地震潛勢研究中得到了廣泛應用和驗證.由于GNSS時間序列中存在與長期構造運動無關的信息致使擬合的速度場無法真實表達實際的震間地殼運動特征，進而影響到應變率場的結果.

本文提出了顧及抗差性的速度場模型構建方法，考慮到了測站因局部運動和數據質量較差在速度場中產生的粗差.以解算的川滇地區最新GNSS觀測資料為例，在觀測的速度場中檢測出了20%的降權測站和2%的除權測站，對這些測站進一步分析發現，粗差產生的原因是由于觀測數據質量較差和測站為土層觀測點而存在局部非構造變形，進而獲得穩健性較好的速度場模型.通過在觀測數據人為加入粗差檢驗本文方法的抗差性，結果表明單個存在的粗差對速度場模型不會產生影響，而區域性存在的粗差會對結果產生影響，表明本文的方法具有好的抗差性.本文進一步分析了顧及抗差性算法對解算應變率場的貢獻，進而給出了川滇地區精準的應變率場.

致謝感謝中國地震臺網中心提供的中國大陸構造環境監測網絡GNSS觀測數據.