海底震后形變多基站聯合估計模型

薛樹強, 肖圳,*, 董杰, 韓保民, 孫悅,, 楊文龍

1 中國測繪科學研究院地理信息工程國家重點實驗室, 北京 100830

2 山東理工大學建筑工程與空間信息學院, 山東淄博 255049

0 引言

隨著以全球導航衛星系統(Global Navigation Satellite System, GNSS)為代表的空間大地測量觀測技術的發展和成熟,大地測量監測網絡已實現全球陸域有效覆蓋,且實現了長期連續觀測(楊元喜等,2011,2020).然而,占地球表面積71%的海洋還存在大量海底基準空白和導航定位服務盲區(楊元喜等,2017).為此,國際大地測量協會專門成立了海洋大地測量交叉委員會(Poutanen and Rózsa,2020);同時,我國在“十三五”期間也開展了海底大地測量基準試驗網建設(中國科學院,2022).新一代大地測量觀測系統(Global Geodetic Observing System, GGOS)致力于地球變化監測并為地震等災害過程研究提供高精度觀測,但由于絕大部分地震發生在板塊交接帶,且分布于海底(Bürgmann and Chadwell,2014),導致海底構造探測和地震等災害過程研究在廣闊海域缺少有效約束(Liu等,2023).因此,海底大地測量對于全球海底板塊監測、地球動力學和地震等災害過程研究具有重要意義.

最近20多年來,震后形變理論和算法得到了較大的發展(Sun and Okubo,1993;喬學軍等,2021;劉泰等,2022),國內外許多學者基于大地測量觀測對震后形變開展了大量研究工作,例如有學者利用震后黏彈性模型研究地震后對周邊區域地殼形變的影響(余建勝等,2018;劉泰等,2017;陳飛等,2020);還有學者從Sentinel-1和ALOS-2圖像中獲得同震形變場并進行一系列分析(姜衛平等,2022);也有學者利用水準監測數據建模推斷震后形變源和震后形變機制(Hao et al.,2012;Iinuma et al.,2012);Tobita結合日本東北部GNSS時間序列和指數、對數模型構建了震后形變模型(Tobita,2016).在GNSS坐標時間序列震后形變估計方面,不少學者研究采用分步估計法,即先利用震前數據估計站速度,然后通過去除震后數據中站速度引起的位移和同震形變得到震后形變(Gonzalez-Ortega et al.,2014;Tobita,2016).對GNSS時間序列進行整體建模的模型參數估計具有更為嚴密的理論基礎(Altamimi et al.,2016).此外,GNSS測站震后形變具有相同的驅動機制,所以相近測站具有相似的震后形變,從而導致GNSS時間序列中存在共模形變,因此已有學者采用主成分分析法(Principal component analysis, PCA)提取共模形變,且采用兩步法估計或迭代估計得到更為合理的震后形變參數(Savage and Svarc,1997;蘇利娜等,2019).

利用海底大地測量開展地震過程研究,已成為最近十多年來的熱點研究課題(Fujiwara et al.,2022;Iinuma et al.,2012;Ozawa et al.,2012;Watanabe et al.,2014;孫悅等,2023).目前,海底大地測量、海底板塊運動和擴張監測等通常基于全球導航衛星系統-聲學(Global Navigation Satellite System-Acoustics, GNSS-A)組合觀測方法實現,該方法最早由美國斯克利普斯海洋學研究所提出(Spiess,1980).水下定位受復雜海洋溫度、鹽度等環境參數影響,特別是海洋聲速場隨時間和空間都存在顯著的差異性,給水下精密定位帶來了極大挑戰(趙建虎和梁文彪,2019;薛樹強等,2023).為了提高水下定位精度和可靠性,通常需要在等水深的半徑上均勻布設多個海底基準站構成區域海底基準網(Yokota et al.,2016;Matsumoto and Fujita,2008),例如,美國學者提出的在海底布設3～4個海底基準站的方法一直被日本等沿用至今(Spiess,1985a,1985b).該方法經過最近20多年的發展和完善,其平面定位精度已經達到2～3 cm(Fujita et al.,2006;Watanabe et al.,2020;Yang and Qin,2021),而Sato等利用日本東北大地震前的數據對日本海底地殼運動觀測系統進行了評估,結果表明,海底基準站高程定位精度也可達到10 cm(Sato et al.,2013b).海底基準站觀測設備由于供電等問題需要新舊站更替或替換電池,由于水下作業的固有成本和技術難度,通常采用新舊站更替方式保證時間序列的連續性,且每年一般只能復測3～5次,導致海底基準站坐標時序樣本相對較小,無法與GNSS等實時連續大地測量觀測相比(Yokota et al.,2018).除上述原因之外,水下定位還會受到測量船相對于海底控制點的航跡、聲速剖面誤差的影響(馬越原等,2022;肖圳等,2023;Xue et al., 2023).

日本近20年建立了世界上最為密集的海底地殼運動觀測網絡(Ozawa et al.,2002;Ozawa et al.,2012;Yokota et al.,2016),該網絡清晰記錄了一些大型地震引起的形變等,許多學者在此基礎上對2011年日本東北部地震(MW9.0)進行了研究,例如,有學者實現了災害過程模擬(Iinuma et al.,2012; Sun et al.,2014),研究了日本東北部地震導致的同震滑動分布和板塊間耦合(Iinuma et al.,2012;Sato et al.,2013a).日本學者利用海底大地測量也開展了大量地震過程研究,并在2020年公開了7個海底站的十年尺度監測資料,這為本文研究提供了數據條件.本文借鑒國際地球參考框架2020(ITRF 2020)震后形變研究思想(Altamimi et al.,2016),結合海底基準站坐標時序樣本小、存在中斷等實際問題,構建了以站速度和震后形變作為公共參數的多站聯合時序分析模型;后又提出了將N-E方向震后形變弛豫因子作為公共參數的震后形變參數聯合估計模型,并提出了求解站坐標、站速率等線性參數和震后形變非線性參數混合模型的“一步法”;最后采用日本海底基準站網觀測時序驗證了本文的主要研究結果.

1 海底震后形變模型

海底震后形變會影響海底基準站坐標時序分析,對海底基準維持產生不利影響.為此,本文嘗試引入海底基準參考歷元,并構建海底基準站坐標、速度、震后形變等數聯合估計模型.由于海底基準站局部網通常布設在地殼構造穩定區域,且分布以水深為半徑的圓上,可假設有K個海底基站,第i個海底基站在參考歷元t0時刻的平面坐標為xi0(i=1,2,…,K),所有海底基準站具有相同的速度v0和相同的震后形變,則t時刻的海底基準站網坐標可表示為

(1)

其中,v0稱為海底基準站在參考歷元t0的速度,δPSD為震后形變,εi(t)為觀測誤差.上述模型將鄰近站點的海底震后形變作為公共參數,因此我們稱其為海底震后形變網解模型.考慮海底基準站定位精度相對較低,且海底復測周期長導致觀測樣本有限,所以本文建議采用上述多站聯合的網解模型.相對于GNSS測站相距幾十千米,海底基準站間相距僅有幾百米到幾千米,故而將站速度和震后形變參數作為公共參數是合理的.相比于基于主成分分析方法的震后共模形變提取方法(蘇利娜等,2019),本文認為上述公共參數法對于平滑觀測噪聲可具有更強的約束性.

海底基準站隨板塊存在近似線性運動的同時,還會受地震等地質事件影響,存在非線性運動.震后形變通常采用指數或對數函數模擬,如指數函數模型可表示為(Altamimi et al.,2016,2018)

(2)

其中,y=[KEKNτ]T,KE,KN分別為E方向和N方向的振幅因子,τ分別為震后形變弛豫因子(我們假設N方向和E方向具有相同的弛豫因子),T0為地震發生時刻.已有學者利用指數和對數組合模型開展震后形變分析,并給出了相應參數的物理意義(Tobita,2016).在實踐中,分別對測站坐標時序分量構建震后形變模型,會得到不同的弛豫時間因子.然而,對于同一個地震事件,不同方向的震后弛豫因子應該是相同的.因此,本文采用N-E時序聯合估計地震形變弛豫時間因子.需要指出,雖然理論上地震會影響全球測站坐標時序,但無論是從計算成本還是從實際應用方便角度,對震后形變參數實施假設檢驗,當震后形變參數顯著時方可視為模型誤差予以參數化估計,即上述模型應建立在震后形變識別基礎上.當采用假設檢驗方法對附加參數進行識別時,需要充分考慮檢驗參數的相關性問題.對于長時序分析問題,參數間的相關性一般很低,此時建議采用t檢驗對對逐個參數施加假設檢驗(陶本藻,1986).此外,對于多維模型參數選取問題,亦可采用更為高效的AIC或BIC信息熵準則(Akaike,1973;薛樹強和楊元喜,2013;薛樹強等,2022).

當存在n個觀測歷元時,可構建觀測方程：

j=1,2,…,n

(3)

其中,xij=xi(tj)為海底基準站在觀測歷元tj的坐標,可由GNSS-A觀測技術獲取.上述多站聯合網解模型及N-E時序聯合估計模型,也是解決目前海底基準站小樣本坐標時序觀測和觀測精度相對較低問題的有效途徑.

2 海底基準站坐標與速度聯合估計模型

Eij=Ei0+vE0(tj-t0)+δEPSD(tj)+εEij,Nij=Ni0+vN0(tj-t0)+δNPSD(tj)+εNij,

(4)

線性化可得

dLi,j=ai,jdxi+bi,jdy,

(5)

其中,dxi=[dEidNidvEdvN]T,dy=[dKEdKNdτ]T,

(6)

(7)

(8)

當存在K個海底基準站和n個觀測歷元時,即可構成如下矩陣形式的觀測模型：

dL=CdX+ε,

(9)

(10)

其中,

(11)

為K個測站在n個觀測歷元的累計設計矩陣.考慮當非線性迭代廣泛存在的不適定問題,建議考慮附加一定先驗約束.

基于高斯-牛頓法可迭代計算線性化觀測模型(9)的非線性最小二乘解,即

(12)

w(dLi)=

(13)

為等價權因子,σ0為觀測單位權方差.IGGIII方案擁有正常權段、可疑降權段、以及淘汰權段,可以充分利用觀測數據,具有較強的抗差性,其中,k1=1.5,k2=2.5為常用推薦值(楊元喜,1996).

采用類似的方法,可以分別獲取E方向和N方向的震后形變振幅和弛豫時間因子,下文不再贅述.此外,擬合震后形變模型也可采用以下“兩步法”,即

(1)首先,忽略模型(5)中的震后形變項,獲取觀測時序殘差V,即

V=[I-A(ATPA)ATP]dL;

(14)

(2)將上述殘差向量作為虛擬觀測量,采用高斯-牛頓法求解震后形變參數,即

(15)

其中,dSk=V-δPSD(y,t),δPSD(y,t)為全部觀測的非線性函數模型項.

需要指出,上述兩步法需要考慮殘差之間的相關性,并需要迭代計算以消除分步估計的影響.對于小樣本觀測,忽略殘差的相關性可能會對震后形變參數估計產生較大影響.考慮到海底基準站一般每年復測3～5次,觀測樣本較小會導致殘差相關性較大,本文建議采用理論更為嚴密的“一步法”,即采用模型(12)一次獲取站坐標及震后形變參數.

3 試驗結果分析

2011年3月11日,日本海底發生了9.0級大地震,國內外學者圍繞日本震后地球物理反演開展了大量研究工作(Ozawa et al.,2012;Sun et al.,2014;陳飛等,2020).本文利用日本公布的7個觀測站在2011—2020年的GNSS-A實測數據(空間分布如圖1所示)驗證本文方法的有效性.這7個站均受到2011年日本大地震的影響,以FUKU站為例,E方向的中心點站坐標時序如圖2所示,可以看到在2011年3月日本大地震發生后坐標受到地震影響,隨著時間的推移逐漸恢復線性趨勢.

圖1 日本海底基準站分布

圖2 受震后影響的FUKU站中心點坐標時序

使用本文提出的海底震后形變網解模型以及兩步法分別對日本公開的7個測站坐標時序求解各測站的站速度和震后形變等參數,水平方向的震后形變參數 (振幅、弛豫因子)如表1所示.

表1 兩步法結果與本文模型解對比Table 1 Comparison between two-step method and the proposed model solution

由于弛豫因子主要與地殼黏滯系數和彈性厚度有關,結合區域構造特征可知(見表1),KAMN、KAMS、MYGI三個站所在塊體具有相同運動趨勢,而MYGW、FUKU、CHOS三個站又與上述三個站具有反向運動趨勢,而TOS2距離震中較遠,又表現出于上述三個站不同的運動趨勢,然而,兩步法得到的震后弛豫因子多數均在0.2年左右,他們之間的區分度很小,這意味著利用小樣本時序殘差震后形變參數容易低估弛豫因子的數值.相比之下,一步法得到的弛豫因子結果區分度較好.另外,因為不同海底基準站相對于地震發生時刻的復測時間不同,所以基于不同海底基準站坐標時序分析得到的震后形變振幅和弛豫因子并不在統一的時間參考下.為此,我們統計了日本海底基準站相對于地震發生時刻的復測時間信息,如圖3所示.

圖3給出的7個海底基準站復測時間分布圖表明,TOS2站捕獲的地震形變弛豫因子會相對較小,其他站也會存在微小差異,我們稱之為弛豫因子偏差.為此,我們將計算得到弛豫因子偏差(如表2所示)補償到各測站時序提取的弛豫因子上,得到統一到地震發生起始時刻的震后形變弛豫因子.

圖3 7個海底基準站復測時間分布

表2 各測站弛豫因子偏差Table 2 Relaxation factor deviation of each station

將各測站弛豫因子偏差補償到計算出來的弛豫因子參數后的兩步法與海底震后形變網解模型一步法(分為聯合估計N-E方向弛豫因子和分別估計N、E方向弛豫因子)結果對比如表3所示.

表3 補償后兩步法與海底震后形變網解模型參數結果對比Table 3 Comparison of model parameters between compensated two-step method and post earthquake deformation network

由表3可知,兩步法得到的弛豫因子依然沒有較好的區分度.需要指出,對于GNSS基準站連續坐標時序的大樣本情形,殘差的相關性可以忽略,但上述試驗結果表明,對于海底基準站復測小樣本觀測情形,殘差的相關性可對震后形變估計產生很大影響.一步法(7參數法)即N-E聯合估計得到的弛豫因子明顯可以根據數據的大小分為兩組(前三個站為一組,后四個站為一組),這主要是因為,兩組觀測站空間分布差異性決定的,即處于兩個不同的塊體上,這一結論與實際站速度觀測整體也是基本符合的(如圖4),因此本文認為研究海底震后形變參數時使用一步法(7參數法)得到的結果相比于兩步法要更準確.

圖4 海底基準站運動

4 結論

對海底局域網內的多個基準站坐標時序聯合處理,不僅可以提高海底板塊運動的監測精度,也可以提高海底震后形變信息提取的敏感性.受限于觀測成本和觀測技術所限,一般每年只能對海底基準開展3～5次復測,相對GNSS而言其復測精度也偏低,為此,本文建議將多個測站的站速度和震后形變參數作為公共參數處理,并采用“一步法”對線性參數和非線性參數進行聯合解算.利用海底基準站復測坐標時序研究震后弛豫因子時,需要充分考慮地震發生時刻與海底基準站復測時刻的偏差影響.

當采用“兩步法”提取震后形變信息時,受海底基準坐標觀測時序的樣本大小限制,必然導致其在第一步時間序列線性擬合階段獲取的殘差時間序列存在較大的相關性,若忽略這種相關性影響,則會影響其在第二步提取震后形變信息.因此,本文建議采用“一步法”同時對測站坐標、速度和震后形變參數做出聯合估計.我們認為,地震引起地殼形變的弛豫時間不應因坐標系的定義不同而不同,因此,為了更進一步提高震后形變信息提取的可靠性,我們建議將N、E坐標時序中的震后形變弛豫因子作為公共參數進行估計.

研究發現,日本2011年發生了MW9.0大地震后,位于不同海底板塊上的海底測站感應到的海底震后形變振幅和弛豫時間存在較為明顯的差異性,同一板塊上的震后弛豫時間具有較好的一致性,但當采用“兩步法”提取震后弛豫因子時,則很難獲取與之一致的分析結果.因此本文認為,研究海底震后形變參數時選擇一步法會得到比“兩步法”更準確的結果.