孫 悅,薛樹強(qiáng),韓保民,肖 圳

1.山東理工大學(xué)建筑工程與空間信息學(xué)院,山東 淄博 255049; 2.中國測繪科學(xué)研究院,北京 100036

GNSS的廣泛應(yīng)用,使得陸域地殼運(yùn)動(dòng)觀測網(wǎng)絡(luò)已經(jīng)非常密集,且實(shí)現(xiàn)了長期連續(xù)觀測[1-3]。1980年提出的GPS與聲學(xué)測距相結(jié)合的聯(lián)合觀測系統(tǒng)讓建設(shè)高精度海底大地測量的海底參考點(diǎn)成為可能[4]。經(jīng)過多年的發(fā)展,該方法理論與技術(shù)得到了不斷的完善[5-7]。水下定位精度主要受復(fù)雜海洋環(huán)境的影響,特別是海洋聲速場時(shí)空變化誤差的影響,例如,季節(jié)性溫度和流場變化會對聲速剖面整體產(chǎn)生影響等[8-9]。近年來,海底高精度定位方法也不斷完善[10-11],實(shí)現(xiàn)了聲速誤差的有效補(bǔ)償[12-13]。需要指出的是,海底基準(zhǔn)網(wǎng)采用對稱設(shè)計(jì),有利于減小聲速誤差影響[1,14-15]。

文獻(xiàn)[16]提出在海底布設(shè)多個(gè)基準(zhǔn)點(diǎn)構(gòu)成的基準(zhǔn)網(wǎng),其基準(zhǔn)點(diǎn)均勻地分布在半徑近似等于水深的圓上,可實(shí)現(xiàn)更為準(zhǔn)確的海底地殼運(yùn)動(dòng)監(jiān)測。海底多站聯(lián)合觀測技術(shù)一直沿用至今,已實(shí)現(xiàn)海底位移監(jiān)測和海底板塊擴(kuò)展監(jiān)測[17-18]。由于供電不足等原因,海底基準(zhǔn)站難以全部正常工作,更難以保障長時(shí)間連續(xù)觀測。倘若海底一組基準(zhǔn)站即將供電不足,則需要重新布站或更換電池,考慮到水下作業(yè)的固有成本和技術(shù)難度,通常在每個(gè)海底基準(zhǔn)站附近安裝一個(gè)新的海底基準(zhǔn)站[19]。中心點(diǎn)法是目前國際上開展海底地殼運(yùn)動(dòng)監(jiān)測普遍使用的方法,而換站補(bǔ)償是該方法必須采用的策略。當(dāng)然,原位觀測是海底基準(zhǔn)維護(hù)的理想途徑,但實(shí)現(xiàn)原位觀測具有較高的深海工程作業(yè)難度。因此,聯(lián)合多個(gè)間斷基準(zhǔn)站坐標(biāo)時(shí)序獲取海底構(gòu)造運(yùn)動(dòng)信息就成為實(shí)現(xiàn)海底大地測量監(jiān)測的重要課題,文獻(xiàn)[15]提出采用海底基準(zhǔn)站網(wǎng)的中心點(diǎn)作為虛擬觀測,通過估計(jì)虛擬觀測中心點(diǎn)的偏移量進(jìn)行形變監(jiān)測。近20年太平洋西岸附近建立的世界上最為密集的地殼運(yùn)動(dòng)觀測網(wǎng)絡(luò)清晰記錄了一系列大型地震引起的形變等,實(shí)現(xiàn)了災(zāi)害過程模擬[20-22]。文獻(xiàn)[23]利用日本東北大地震前的數(shù)據(jù)對日本海底地殼運(yùn)動(dòng)觀測系統(tǒng)進(jìn)行了評估,結(jié)果表明,海底基準(zhǔn)站網(wǎng)平面定位精度已達(dá)到厘米級,高程方向定位精度也可達(dá)到10 cm。

綜上所述,目前國際上主要采用中心點(diǎn)法和新舊基準(zhǔn)站換站補(bǔ)償方法,實(shí)現(xiàn)海底高精度、長時(shí)序、連續(xù)觀測,并以此為基礎(chǔ)開展海底地殼運(yùn)動(dòng)監(jiān)測。雖然虛擬中心點(diǎn)方法可用于海底區(qū)域地殼運(yùn)動(dòng)監(jiān)測和形變分析,但無法獲取海底基準(zhǔn)站網(wǎng)中各基準(zhǔn)站協(xié)議坐標(biāo),無法為海洋大地測量和海洋導(dǎo)航提供參考框架,也難以實(shí)施精細(xì)化質(zhì)量控制,因此筆者認(rèn)為這不是海底基準(zhǔn)建立與維持的最佳途徑。為此,筆者借鑒國際地球參考框架(ITRF)構(gòu)建方式[24],提出了基于協(xié)議參考?xì)v元的站坐標(biāo)與站速度聯(lián)合估計(jì)方法,可為水下導(dǎo)航定位服務(wù)提供實(shí)時(shí)、動(dòng)態(tài)、高精度參考基準(zhǔn)。

1 海底坐標(biāo)時(shí)序多站聯(lián)合處理模型

假設(shè)整體海底基準(zhǔn)網(wǎng)與地殼存在一個(gè)共同的運(yùn)動(dòng)速度,設(shè)定一個(gè)參考?xì)v元(建議采用觀測周期內(nèi)的中間時(shí)刻),并將每個(gè)站點(diǎn)在該參考?xì)v元處的坐標(biāo)作為未知參數(shù),與海底基準(zhǔn)網(wǎng)的整體運(yùn)動(dòng)速度進(jìn)行聯(lián)合估計(jì),簡稱聯(lián)合估計(jì)模型。設(shè)有I個(gè)海底基準(zhǔn)站,第i個(gè)海底基準(zhǔn)站在t0時(shí)刻的坐標(biāo)xi0,假設(shè)所有海底基準(zhǔn)站在該時(shí)刻具有相同的速度v0,則t時(shí)刻的坐標(biāo)為

(1)

式中,xi0(i=1,2,…,I)為第i個(gè)站在參考?xì)v元t0的協(xié)議坐標(biāo),為未知參數(shù);速度v0為在參考?xì)v元t0的協(xié)議速度;δPSD為地質(zhì)災(zāi)害引起的非線性部分,如地震的震后形變等。一旦獲得各測站的協(xié)議坐標(biāo)和協(xié)議速度,無論該測站是否處于工作狀態(tài),都可實(shí)現(xiàn)該基準(zhǔn)站坐標(biāo)維持?？紤]到數(shù)值計(jì)算需要和大地測量習(xí)慣,上述時(shí)間單位建議采用年。

為便于討論,下面僅考慮平面方向時(shí)間序列,并暫不考慮測站非線性變化(以E方向?yàn)槔?。對于第i個(gè)海底基準(zhǔn)站,假設(shè)存在J個(gè)歷元觀測,則可構(gòu)造以下觀測方程

(2)

式中,ε為觀測與模型不符值,主要包括觀測隨機(jī)誤差、區(qū)域構(gòu)造形變及地震震后形變等影響。等號左側(cè)為觀測量,記為L,右側(cè)xi0和ve0為待估量,則式(2)的最小二乘解為

(3)

式中,P為觀測權(quán)陣;A為觀測設(shè)計(jì)矩陣;1J×1為由1構(gòu)成的(J×1)維列向量;0J×1為由0構(gòu)成的(J×1)維列向量。

海底觀測時(shí)間序列偏離線性運(yùn)動(dòng)模式,一般是由地震等地質(zhì)災(zāi)害事件引起的。此時(shí),可考慮采用地震的震后形變模型來構(gòu)建測站非線性運(yùn)動(dòng)。對于高程方向,還應(yīng)考慮年周期、半年周期等非線性周期性信號[25]。然而,受限于海底定位精度及復(fù)測周期,目前海底基準(zhǔn)站網(wǎng)還很難實(shí)現(xiàn)高程方向時(shí)序分析。因此,本文主要考慮海底水平構(gòu)造運(yùn)動(dòng),并將海底地震等引起的異常當(dāng)成粗差觀測處理。無論是觀測異常,還是各種地質(zhì)災(zāi)害引起的時(shí)間序列異常,都可使用抗差估計(jì)策略予以處理。本文采用以下IGGⅢ抗差估計(jì)[26-28]

(4)

(5)

IGGⅢ方案擁有正常權(quán)段、可疑降權(quán)段,以及淘汰權(quán)段,可充分利用觀測數(shù)據(jù),具有較強(qiáng)抗差性,其中,k1=1.5,k2=2.5為常用推薦值。

(6)

如圖1所示,在舊海底基準(zhǔn)網(wǎng)將要供電不足時(shí)安裝新海底基準(zhǔn)網(wǎng),二者都是局部基準(zhǔn)網(wǎng),共同構(gòu)成整體基準(zhǔn)網(wǎng)。局部基準(zhǔn)網(wǎng)中心點(diǎn)之間存在一定的差值。為了確保上述中心點(diǎn)法監(jiān)測海底地殼運(yùn)動(dòng)的連續(xù)性,需要通過新舊海底基準(zhǔn)網(wǎng)同步觀測,得到中心點(diǎn)間的差值Δx,以后歷元啟動(dòng)新海底基準(zhǔn)站網(wǎng),并將Δx補(bǔ)償?shù)叫潞５谆鶞?zhǔn)網(wǎng)中心坐標(biāo),即可保證中心點(diǎn)坐標(biāo)時(shí)序的連續(xù)性[19]。

圖1 更換海底基準(zhǔn)站

換站補(bǔ)償后的中心點(diǎn)法是海底構(gòu)造監(jiān)測的一種有效方法,采用重心基準(zhǔn)有利于區(qū)域構(gòu)造形變分析,但由于缺少如ITRF控制下的長期基準(zhǔn)約束,不利于長期、大尺度海底構(gòu)造分析。中心點(diǎn)法可以平滑觀測噪聲,但某個(gè)站點(diǎn)異?？杀恢行狞c(diǎn)平均化處理削弱而無法對其進(jìn)行探測和分離,或者因其中少數(shù)站點(diǎn)異常,導(dǎo)致整個(gè)基準(zhǔn)站組視為異常。此外,中心點(diǎn)坐標(biāo)時(shí)序?yàn)樘摂M觀測時(shí)序,從而無法獲取各站點(diǎn)在協(xié)議參考?xì)v元處的站坐標(biāo)估計(jì)和站速度。

2 試驗(yàn)結(jié)果分析

文獻(xiàn)[19]公布了7個(gè)觀測站在2011—2020年的GNSS-A實(shí)測數(shù)據(jù),其空間分布如圖2所示。本文基于這7個(gè)站的觀測數(shù)據(jù)驗(yàn)證本文方法的有效性。

圖2 試驗(yàn)區(qū)海底基準(zhǔn)站分布

2.1 中心點(diǎn)法及換站補(bǔ)償

中心點(diǎn)法及換站補(bǔ)償對GNSS-A數(shù)據(jù)的處理策略見表1,大致分成了3個(gè)環(huán)節(jié)。

表1 數(shù)據(jù)處理策略

以MYGI站E方向?yàn)槔f明換站引起的時(shí)間序列間斷問題。MYGI站在海底有8個(gè)基準(zhǔn)站,編號分別是:M01、M03、M04、M05、M12、M13、M14、M15,在2011—2020年觀測期間發(fā)生了3次換站,本文通過開源的GARPOS軟件解算得到10年的虛擬中心點(diǎn)數(shù)據(jù)。中心點(diǎn)坐標(biāo)如圖3所示,不同海底基準(zhǔn)站網(wǎng)的中心點(diǎn)之間相差較大。因此換站導(dǎo)致原始海底觀測時(shí)序無法直接用于海底地殼運(yùn)動(dòng)監(jiān)測。

圖3 MYGI站原始中心點(diǎn)時(shí)間序列

對虛擬中心點(diǎn)的偏移量時(shí)序進(jìn)行最小二乘擬合,即可得到海底基準(zhǔn)網(wǎng)的整體運(yùn)動(dòng)趨勢和線性速度信息。非線性部分主要代表地殼形變或地震等異常。如圖4所示(以E方向?yàn)槔?,可以看到MYGI站整體符合線性運(yùn)動(dòng)趨勢和地殼整體構(gòu)造運(yùn)動(dòng)特征,其中擬合直線公式為y=vt+dx,v代表局部海底基準(zhǔn)網(wǎng)中心點(diǎn)的運(yùn)動(dòng)速度,dx是最小二乘擬合得到的坐標(biāo)改正數(shù),文中所有中心點(diǎn)速度擬合直線均由上述方法繪制。

圖4 MYGI站中心位移與最小二乘擬合

將新舊海底基準(zhǔn)站網(wǎng)中心之間的差值補(bǔ)償?shù)叫潞５谆鶞?zhǔn)網(wǎng)的中心點(diǎn)上(圖5),可以看到補(bǔ)償之后海底基準(zhǔn)站網(wǎng)的中心點(diǎn)時(shí)序已具備了海底地殼運(yùn)動(dòng)信號提取能力。

圖5 MYGI站補(bǔ)償中心點(diǎn)時(shí)間序列

在使用重心基準(zhǔn)前,先基于初始中心點(diǎn)計(jì)算海底基準(zhǔn)站網(wǎng)的運(yùn)動(dòng)速度,正常情況下的運(yùn)動(dòng)速度為厘米級。由表1可知,MYGW觀測站利用初始中心點(diǎn)計(jì)算出來的運(yùn)動(dòng)速度過大,與實(shí)際不符,即基于初始中心點(diǎn)來研究海底運(yùn)動(dòng)一般是無效的。表2給出了重心基準(zhǔn)下求得的虛擬中心點(diǎn)運(yùn)動(dòng)速度,相對于初始中心點(diǎn)結(jié)果更符合地殼運(yùn)動(dòng)實(shí)際情況。

表2 水平方向初始中心解重心基準(zhǔn)解的中心速度

2.2 聯(lián)合估計(jì)模型解

下文繼續(xù)以MYGI站為例,使用本文模型計(jì)算整體基準(zhǔn)網(wǎng)速度。通過GARPOS軟件確定海底坐標(biāo),圖6中擬合直線公式為y=vt+dx,v代表整體海底基準(zhǔn)網(wǎng)運(yùn)動(dòng)速度,dx表示最小二乘擬合得到的協(xié)議坐標(biāo)改正數(shù)的平均值,文中所有動(dòng)態(tài)平差速度擬合直線均由上述方法繪制?？梢钥闯?受觀測異常影響,絕大部分海底基準(zhǔn)站位移時(shí)序都處于擬合直線的左下側(cè),這是不合理的,筆者認(rèn)為這主要是由觀測異常引起的。

圖6 MYGI站海底基準(zhǔn)站位移與模型擬合

接下來比較本文聯(lián)合估計(jì)模型解中的站速度與重心基準(zhǔn)下中心點(diǎn)速度的差異,比較結(jié)果見表3,兩種結(jié)果之間的差值大多在毫米級。因此,利用聯(lián)合估計(jì)模型計(jì)算出來的速度也是有效的。

表3 中心速度與聯(lián)合估計(jì)模型速度

由表3可知,此時(shí)水平方向的中心速度與聯(lián)合估計(jì)模型速度之間平均相差5 mm/a左右,相差較大。下文分析速度差值較大的站點(diǎn)。例如,CHOS觀測站在E方向上的中心速度與聯(lián)合估計(jì)模型速度相差1.94 cm/a,誤差較大。因此,本文著重分析CHOS觀測站E方向時(shí)序。CHOS站兩種方法各自的殘差分布情況如圖7所示。由圖7可知,2011年前幾次觀測數(shù)據(jù)的模型殘差和中心殘差均較大,因此,可認(rèn)為2011年的觀測應(yīng)視為異常觀測?？紤]到2011年該區(qū)域發(fā)生了Mw 9.0級地震,該異?？赡苁堑卣鸬恼鸷笮巫儭椥运沙诘仍蛞鸬?。此外,在2018年,聯(lián)合估計(jì)模型的殘差較大,而中心殘差沒有明顯異常?？梢哉J(rèn)為聯(lián)合估計(jì)模型對觀測異常更為敏感。

圖7 CHOS站模型擬合殘差與中心擬合殘差比較

為消除異常觀測影響,本文采用IGGⅢ抗差估計(jì),即對聯(lián)合估計(jì)模型與中心點(diǎn)法采用相同策略進(jìn)行抗差估計(jì),結(jié)果見表4。

表4 CHOS站E方向抗差前后速度

由表4可知,抗差之后聯(lián)合估計(jì)模型速度與中心速度之間的差值已經(jīng)減小到了毫米級,且本文方法抗差前后的結(jié)果變化也相對較小,說明本文方法更穩(wěn)健?？共詈髢煞N方法的殘差比較如圖8所示,說明抗差對聯(lián)合估計(jì)模型與中心點(diǎn)法都有明顯的改善作用,經(jīng)過抗差估計(jì)可以把絕大部分殘差約束在±0.1 m以內(nèi)。

圖8 CHOS站抗差后模型擬合殘差與中心擬合殘差比較

結(jié)合圖7、圖8可以看出,聯(lián)合估計(jì)模型和中心點(diǎn)法在抗差過程中都對2011年的觀測數(shù)據(jù)進(jìn)行了處理,但本文提出的聯(lián)合估計(jì)模型,可對數(shù)據(jù)處理進(jìn)行精細(xì)抗差,即可只對單個(gè)有問題的海底基準(zhǔn)站觀測數(shù)據(jù)進(jìn)行刪除或降權(quán);而中心點(diǎn)法中的抗差估計(jì),則是直接對中心點(diǎn)進(jìn)行刪除或降權(quán),這樣處理比較粗略,其對粗差的敏感度不如聯(lián)合估計(jì)模型。具體表現(xiàn)為聯(lián)合估計(jì)模型在抗差過程中可檢測到多個(gè)歷元觀測異常,而中心點(diǎn)法則只檢測到2011年附近的觀測異常。

下文給出所有觀測站的抗差估計(jì)結(jié)果,兩種方法的比較見表5。此時(shí)聯(lián)合估計(jì)模型的速度估計(jì)與中心點(diǎn)法的速度估計(jì)之間的絕對差值平均接近3.1 mm/a,有了較大改善。

表5 抗差后水平方向兩種方法運(yùn)動(dòng)速度

為評定中心點(diǎn)法與聯(lián)合估計(jì)模型的結(jié)果精度,本文計(jì)算了這兩種方法進(jìn)行抗差最小二乘時(shí)的標(biāo)準(zhǔn)差(表6)。

表6 中心點(diǎn)解與動(dòng)態(tài)模型解的標(biāo)準(zhǔn)差

由表6可知,大多數(shù)站的聯(lián)合估計(jì)模型標(biāo)準(zhǔn)差小于中心點(diǎn)法標(biāo)準(zhǔn)差,但也存在反常情況。主要是由于中心點(diǎn)法對整組觀測進(jìn)行淘汰或降權(quán),而聯(lián)合估計(jì)模型對數(shù)據(jù)的處理要更為精細(xì),僅對存在異常的基準(zhǔn)站點(diǎn)進(jìn)行異常觀測質(zhì)量控制。反常發(fā)生在MYGI站的N方向,這里中心點(diǎn)法的標(biāo)準(zhǔn)差要小于聯(lián)合估計(jì)模型的標(biāo)準(zhǔn)差,兩種方法的最小二乘和抗差最小二乘的速度擬合直線如圖9、圖10所示。圖10中MYGI站在N方向上的海底基準(zhǔn)站位移時(shí)序較為散亂,其精度相對較低也是合理的。而圖9的中心點(diǎn)時(shí)序也不穩(wěn)定,可能是在抗差最小二乘過程中刪除較多的觀測數(shù)據(jù),所以得到的標(biāo)準(zhǔn)差較小。

圖9 MYGI站的本文方法速度擬合

圖10 MYGI站的中心點(diǎn)法速度擬合

圖11是所有觀測站分別用兩種方法計(jì)算出來的運(yùn)動(dòng)趨勢示意圖,其中,深色箭頭是聯(lián)合估計(jì)模型的運(yùn)動(dòng)方向及速度大小,淺色箭頭是中心點(diǎn)法的運(yùn)動(dòng)方向及速度大小。可以看出,觀測站按照運(yùn)動(dòng)方向大致分為了兩組,分別向兩個(gè)不同的方向運(yùn)動(dòng)。這主要是板塊構(gòu)造引起的,且研究結(jié)果和其他學(xué)者的研究結(jié)果具有很好的一致性。

3 結(jié) 論

本文構(gòu)建的海底基準(zhǔn)站網(wǎng)聯(lián)合估計(jì)模型,參考國際地球參考框架(ITRF)建立與維持策略,采用明確的參考?xì)v元站坐標(biāo)與站速度作為未知參數(shù),對海底基準(zhǔn)站網(wǎng)時(shí)序觀測進(jìn)行動(dòng)態(tài)數(shù)據(jù)處理。不但可以研究形變問題,還可以獲取各站點(diǎn)在協(xié)議參考?xì)v元的站坐標(biāo)和站速度,用于導(dǎo)航應(yīng)用。雖然兩種方法所采用的基準(zhǔn)和換站處理方法均不同,但計(jì)算出來的海底基準(zhǔn)站網(wǎng)速度非常接近。中心點(diǎn)法需要在換站情況下新舊兩個(gè)海底基準(zhǔn)站網(wǎng)同時(shí)觀測來進(jìn)行中心點(diǎn)補(bǔ)償,而聯(lián)合估計(jì)模型則不需要新舊基準(zhǔn)網(wǎng)聯(lián)合觀測,這對海底基準(zhǔn)站網(wǎng)維持的觀測要求相對較低。需要指出,在不同觀測時(shí)期,海底基準(zhǔn)站網(wǎng)定位可能采用不同的ITRF參考框架,在海底基準(zhǔn)站時(shí)序觀測數(shù)據(jù)處理時(shí)需要考慮上述參考框架間的差異。

海底地震事件或觀測異常可引起海底基準(zhǔn)站時(shí)序觀測出現(xiàn)異常,即出現(xiàn)脫離測站時(shí)序線性運(yùn)動(dòng)趨勢的現(xiàn)象。因此,若要精確估計(jì)海底基準(zhǔn)站網(wǎng)的線性運(yùn)動(dòng)速度,需要對這些異常觀測進(jìn)行質(zhì)量控制。然而,由于海底基準(zhǔn)站網(wǎng)一般每年只能復(fù)測2～4次,觀測樣本數(shù)量有限,為此,本文提出抗差最小二乘,有利于改進(jìn)中心點(diǎn)時(shí)序和各測站原始觀測時(shí)序的處理,且可提高測站速度估計(jì)的可靠性。

研究發(fā)現(xiàn),由于日本2011年發(fā)生了Mw 9.0級大地震,期間海底基準(zhǔn)站時(shí)序觀測存在明顯的觀測異常。因此,對于明確地震引起的震后形變,可采用ITRF 2014指數(shù)-對數(shù)模型予以處理,這也是本文的后續(xù)研究內(nèi)容之一。此外,MYGW站和MYGI站相鄰,但存在明顯不同的測站運(yùn)動(dòng)方向,這意味著該區(qū)域存在較為復(fù)雜的地質(zhì)構(gòu)造運(yùn)動(dòng)。