吳富梅，魏子卿，劉光明

1. 地理信息工程國家重點實驗室，陜西 西安 710054； 2. 西安測繪研究所，陜西 西安 710054

?

利用GNSS數(shù)據(jù)分析大港驗潮站地殼沉降

吳富梅1,2，魏子卿1,2，劉光明1,2

1. 地理信息工程國家重點實驗室，陜西 西安 710054； 2. 西安測繪研究所，陜西 西安 710054

青島大港驗潮站的地殼沉降關(guān)系到該站平均海平面的絕對變化，因而也就關(guān)系到我國高程基準面的變化。本文利用青島GNSS基準站約10年的觀測數(shù)據(jù)對該站的地殼沉降變化進行分析。首先將青島GNSS基準站納入由50個國際IGS站和43個國內(nèi)陸態(tài)網(wǎng)絡基準站組成的全球網(wǎng)中，進行單日松弛解和單日約束解解算，獲得該站坐標時間序列。然后對該站垂向坐標時間序列進行分析，利用粗差探測、偏差探測、趨勢項分析、頻譜分析等方法對粗差、偏差、趨勢項和周期項進行探測、分析，并通過時間序列模型估計獲得時間序列中的周期項振幅和偏差估值。分析表明青島GNSS基準站垂直方向近一段時間未發(fā)現(xiàn)存在顯著性的地殼沉降變化，但受到比較明顯的周年和半周年周期變化影響。結(jié)合青島大港驗潮站驗潮數(shù)據(jù)分析結(jié)果得出結(jié)論：青島大港驗潮站平均海平面的絕對上升速率是1.62 mm/a。

高程基準; GNSS基準站; 時間序列; 地殼沉降變化

1957年，我國確定將青島大港驗潮站作為國家高程基準基本驗潮站，將大港驗潮站的平均海平面定義為我國高程基準的起算面。我國目前采用的“1985國家高程基準”，是由大港驗潮站1952至1979年驗潮數(shù)據(jù)確定的。計算時以19年為周期，滑動步長為1年，得到10組平均海面值，取其平均值242.89 cm作為當?shù)攸S海平均海面高，從而求出國家水準原點的高程值為72.260 m[1-3]。

大港驗潮站平均海平面的變化直接關(guān)系到我國高程基準面的變化。引起驗潮站平均海平面發(fā)生沉降變化的因素很多，例如溫室效應造成的冰川消融、局部地區(qū)地下水的開采以及海底地形的變化等[4-5]，另外一個不可忽視的變化是地殼的沉降變化。

為了監(jiān)測我國高程基準面的變化，并研究地殼沉降變化對大港驗潮站平均海平面的影響，2005年在距離青島驗潮站約1.5 km處建立了青島GNSS連續(xù)運行觀測站，于2005年9月開始試運行。至今，青島GNSS基準站已經(jīng)累積觀測數(shù)據(jù)10年左右，對這些數(shù)據(jù)進行處理，研究其垂向坐標時間序列，分析其趨勢項，可以監(jiān)測我國高程基準面升降變化。基于此，本文首先對青島GNSS基準站約10年的觀測數(shù)據(jù)進行組網(wǎng)解算，并對準于ITRF2014框架，獲得基準站垂向坐標時間序列；然后利用多種方法對垂向坐標時間序列中的粗差、偏差、趨勢項和周期分量等分別進行處理和分析，得出青島GNSS基準站垂向時間序列不存在顯著趨勢性變化，但受到比較明顯的周年和半周年周期變化影響。最后結(jié)合青島大港驗潮站數(shù)據(jù)分析結(jié)果得出：青島大港驗潮站平均海平面的絕對上升速率是1.62 mm/a。

1 大港驗潮站平均海平面升降原理

圖1為平均海平面絕對變化原理示意圖。監(jiān)測系統(tǒng)由GNSS基準站和近鄰的驗潮站組成，它們共建在海邊的同一處穩(wěn)定的基巖上，可以認為其地殼運動的影響是相同的。

圖1 平均海平面絕對變化原理圖Fig.1 Absolute change of the mean sea level

通過分析驗潮站驗潮數(shù)據(jù)可以獲得平均海面相對于水尺零點的高度ζ，通過GNSS基準站可以獲得天線處相對于參考橢球面的大地高H，假設(shè)GNSS基準站的天線相位中心與驗潮站水尺零點高差為h，平均海平面到參考框架橢球面的高度為H′，則

H′=H-h+ζ

(1)

(2)

式(2)描述了平均海面絕對變化。由該式可見，平均海平面的絕對變化包括兩方面的影響，即GNSS基準站求得的當?shù)仃懙卮怪边\動和驗潮資料求得的海平面相對地殼的變化的共同作用，絕對平均海面變化速度等于驗潮站垂直地殼形變與相對海平面變化速度的代數(shù)和。

根據(jù)參考文獻[6]，獲得青島大港驗潮站平均海平面相對地殼的運動速度是1.62 mm/a。需要求得GNSS基準站地殼垂直運動速度，才能獲得青島驗潮站平均海面的絕對變化。

2 青島GNSS基準站數(shù)據(jù)預處理

在對青島GNSS基準站垂向坐標時間序列分析之前，首先對原始觀測數(shù)據(jù)進行預處理，主要包括以下兩步：

2.1 單日松弛解解算

將全球50個穩(wěn)定、可靠、具有連續(xù)觀測數(shù)據(jù)的IGS站與國內(nèi)包含SDQD在內(nèi)的43個陸態(tài)網(wǎng)絡基準站組成一個全球網(wǎng)，采用GAMIT軟件對該網(wǎng)2005年第295 d至2015年第365 d的觀測數(shù)據(jù)進行單日松弛解解算，在解算中國際IGS站坐標給予1 m約束，國內(nèi)基準站坐標給予10 m約束。除去個別天因觀測數(shù)據(jù)異常無法解算或者無觀測數(shù)據(jù)，共計獲得3095個單日松弛解。

2.2 單日約束解解算

獲得單日松弛解后，首先去約束[7-9]

(Σunc)-1(Σest)-1-(Σconst)-1

(3)

式中，Σunc代表去約束之后的協(xié)方差陣；Σest代表松弛解獲得的協(xié)方差陣；Σconst代表GAMIT解算時加入的約束。

然后以50個IGS站的ITRF2014坐標為參考基準，采用如下公式加入最小約束[7-9]，將整網(wǎng)坐標對準于ITRF2014

(4)

式中，N=(Σunc)-1；ΔX=XC-Xapr，Xapr是坐標近似值(GAMIT解算時采用的概略坐標)，XC是解算獲得的站坐標；K是單日松弛解獲得的自由項；B=(ATA)-1AT

解算中，2006年第298 d和2007年第86 d因觀測數(shù)據(jù)問題，解算結(jié)果異常而剔除，這樣共計獲得3293 d垂向坐標時間序列，如圖2所示。

圖2 青島GNSS基準站垂向坐標時間序列Fig.2 Height time series of Qingdao GNSS station

3 青島GNSS基準站垂向坐標時間序列分析

在進行地殼沉降分析之前，需要對GNSS垂向坐標時間序列進行“清理”，獲得比較“純凈”的時間序列。GNSS垂向坐標時間序列中包含各種信號和誤差影響，主要包括地殼運動引起的長期變化(地殼沉降)、與地球物理相關(guān)的季節(jié)性變化、與技術(shù)數(shù)據(jù)處理相關(guān)的一些誤差以及由于地震、基墩、接收機故障等引起的偏差或者粗差[10]。

3.1 粗差探測

垂向坐標時間序列標準化殘差[11-12]

(5)

式中，vi為殘差；σi為均方差因子。

3.2 偏差探測

采用STARS算法進行偏差探測，STARS算法最早用于氣象數(shù)據(jù)偏差探測[10,13-14]。STARS算法的基本原理是基于t分布

(6)

如果RSII+l-1,i>0，就判斷vi發(fā)生偏差跳變；如果RSIi+1-1,i<0，就判斷vi沒有發(fā)生偏差跳變。

3.3 時間序列頻譜分析

采用離散傅里葉變換進行時間序列頻譜分析

(8)

式中，W(j)為頻譜值；i為虛數(shù)。

3.4 坐標時間序列模型

垂向坐標時間序列一般用如下模型描述

(9)

通過以上4步不斷迭代，探測粗差和偏差，分析周期項。具體步驟如下：

首先假定垂向坐標時間列中存在趨勢項，利用式(9)去常數(shù)項、趨勢項、周期項(假定有周年和半周年周期項，不考慮偏差)，對殘差進行第1次粗差探測，設(shè)定臨界值為c=8，探測出4個比較大的粗差，如圖3所示。

在此基礎(chǔ)上，進行第1次偏差探測，設(shè)定t0=15(根據(jù)實際數(shù)據(jù)檢驗獲得)，共計探測出3個偏差(發(fā)生時刻在2009.249、2010.909、2010.994)，如圖4所示。

圖3 垂向坐標時間序列及粗差(第1次探測)Fig.3 Height time series and outliers(after the first detection)

圖4 垂向坐標時間序列及偏差(第1次探測)Fig.4 Height time series and offsets(after the first detection)

再次利用式(9)去常數(shù)項、趨勢項、周期項和偏差，對殘差進行第2次粗差探測，設(shè)定臨界值是c=3，探測出30個粗差，如圖5所示。

圖5 垂向坐標時間序列及粗差(第2次探測)Fig.5 Height time series and outliers(after the second detection)

進行第2次偏差探測，設(shè)定t0=10(根據(jù)實際數(shù)據(jù)檢驗獲得)，又探測出3個偏差，如圖6所示(其中有3個是之前探測出的)。

圖6 垂向坐標時間序列及偏差(第2次探測)Fig.6 Height time series and offsets(after the second detection)

將該時間序列剔除粗差、調(diào)整偏差、去常數(shù)項和趨勢項，對缺失的數(shù)據(jù)進行修補之后進行頻譜分析，如圖7所示。從圖中可以看出，垂向坐標時間序列具有較明顯的周年(0.998 2年)、半年周期(0.477年)影響，另外還有一個不明原因、周期接近3個月(0.239年)的噪聲影響。

圖7 垂向坐標時間序列頻譜分析Fig.7 Spectrum analysis of height time series

至此，在垂向坐標時間序列中共計探測出粗差30個，偏差6個，并且明確存在周年和半周年周期影響，為下一步地殼沉降分析奠定基礎(chǔ)。

4 大港驗潮站地殼沉降分析

坐標水平分量因板塊地殼水平運動一定存在趨勢項，但垂向則未必，因此在估計垂向趨勢項之前有必要進行趨勢項檢驗。

4.1 趨勢項檢驗

采用非參數(shù)Mann-Kendall檢驗[15]，提出如下假設(shè)

(10)

統(tǒng)計量ZS服從正態(tài)分布

(11)

Mann-Kendall檢驗統(tǒng)計量S為

(12)

式中，n是時間序列長度；xi、xj分別是第i、j個垂向坐標值；sgn是符號函數(shù)。

(13)

式中，m是重復點的集群數(shù)；ti是第i個集群的重復點數(shù)。

存在顯著性趨勢項變化。

4.2 趨勢項估計

為了進一步驗證垂向坐標時間序列中不存在顯著性趨勢項，采用式(12)，考慮常數(shù)項、趨勢項、周年和半周年(接近3個月的周期，因原因不明未考慮)、粗差以及偏差等各項因素，計算獲得趨勢項速度是-0.1 mm/a，標準差是0.005 mm/a，這進一步驗證了垂向坐標時間序列不存在顯著性趨勢項變化。

同時估計出周年和半周年周期振幅分別為4.684 mm、1.112 mm。圖8給出垂向坐標時間序列及其擬合曲線圖，表1給出6個偏差的估值及發(fā)生時間。

圖8 青島GNSS基準站垂向坐標時間序列及其擬合曲線Fig.8 Height time series and the fitting curve of Qingdao GNSS station

偏差1偏差2偏差3偏差4偏差5偏差6發(fā)生歷元/年2008.8542009.2492010.9092010.9232010.9942015.739偏差估值/mm-2.35020.15417.62113.088-40.13111.483

結(jié)合圖8和表1可以看出，在2009.249、2010.909、2010.923、2010.994發(fā)生較大偏差跳變。據(jù)了解青島站在2009年曾經(jīng)因天線故障更換過接收機天線，但是在天線信息中沒有記錄，所以會在2009.249發(fā)生跳變。2010年陸態(tài)網(wǎng)絡部分站點進行了整改，但是整改是在2010.495～2010.632歷元之間，與2010年發(fā)生的3處跳變沒有關(guān)系，從圖8中也可看出，這3處跳變更像比較短時間內(nèi)的異常，因此判定是由觀測數(shù)據(jù)異常引起的。另外，2008.249處比較小的跳變可能是由觀測數(shù)據(jù)異常引起的，2015.739處的跳變也許因該站點2015年的拆卸有關(guān)，還需更多的數(shù)據(jù)研究。

5 結(jié) 論

青島大港驗潮站的地殼沉降直接關(guān)系到驗潮站平均海平面變化的絕對值，研究青島大港驗潮站地殼垂向運動對我國高程基準基準的維持有重要意義。本文將青島GNSS基準站約10年的觀測數(shù)據(jù)納入全球網(wǎng)進行綜合處理，并進行垂向坐標時間序列分析，得出如下結(jié)論：

(1) 青島GNSS基準站垂直方向受到比較明顯的周年和半周年周期項影響，振幅分別為4.684 mm、1.112 mm。另外還受到不明原因的接近3個月的周期因素影響，需要進一步分析其成因。

(2) 青島GNSS基準站垂直方向發(fā)生6次偏差或者跳變，分別在2008.854、2009.249、2010.909、2010.923、2010.994、2015.739歷元，偏差估值分別為-2.350 mm、20.154 mm、17.621 mm、13.088 mm、-40.131 mm、11.483 mm，其中，2009.249跳變是由更換天線引起的，2015.739跳變也許與拆卸有關(guān)，其余跳變由觀測異常引起。

(3) 在現(xiàn)有觀測時間序列中，未發(fā)現(xiàn)青島GNSS基準站垂向存在顯著趨勢性變化，結(jié)合文獻[6]，可以得出青島大港驗潮站平均海平面的絕對上升速率是1.62 mm/a。

[1] 郭海榮, 焦文海, 楊元喜, 等. 1985國家高程基準的系統(tǒng)差[J]. 武漢大學學報(信息科學版), 2004, 29(8): 715-719. GUO Hairong, JIAO Wenhai, YANG Yuanxi, et al. Systematic Error of the 1985 National Height Datum[J]. Geomatics and Information Science of Wuhan University, 2004, 29(8): 715-719.

[2] 陳宗鏞. 計算平均海面的方法[J]. 山東海洋學院學報, 1960(1): 65-73. CHEN Zongyong. The Method of Calculating the Mean Sea Level[J]. Journal of Ocean Institute of Shandong, 1960(1): 65-73.

[3] 陳宗鏞. 潮汐學[M]. 北京: 科學出版社, 1980. CHEN Zongyong. Tide[M]. Beijing: Science Press, 1980.

[4] 陳俊勇. 對我國建立現(xiàn)代大地坐標系統(tǒng)和高程系統(tǒng)的建議[J]. 測繪通報, 2002(8): 1-5. CHEN Junyong. On the Chinese Modern Geodetic Coordinate System and Height System[J]. Bulletin of Surveying and Mapping, 2002(8): 1-5.

[5] 文援蘭, 楊元喜. 我國近海平均海面及其變化的研究[J]. 武漢大學學報(信息科學版), 2001, 26(2): 127-132. WEN Yuanlan, YANG Yuanxi. Research on Sea Level and Rising Trend of Coastal Waters in China[J]. Geomatics and Information Science of Wuhan University, 2001, 26(2): 127-132.

[6] 吳富梅, 魏子卿, 李迎春. 大港驗潮站潮汐分析與國家高程基準面變化[J]. 測繪學報, 2015, 44(7): 709-716. DOI: 10.11947/j.AGCS.2015.20140110. WU Fumei, WEI Ziqing, LI Yingchun. Analysis of Tidal Data for Dagang Tidal Gauge and Study of the Changes for the National Height Datum[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(7): 709-716. DOI: 10.11947/j.AGCS.2015.20140110.

[7] ALTAMIMI Z, SILLARD P, BOUCHER C. ITRF2000: A New Release of the International Terrestrial Reference Frame for Earth Science Applications[J]. Journal of Geophysical Research, 2002, 107(B10): 2214. DOI: 10.1029/2001JB000561.

[8] ALTAMIMI Z, COLLILIEUX X, LEGRAND J, et al. ITRF2005: A New Release of the International Terrestrial Reference Frame Based on Time Series of Station Positions and Earth Orientation Parameters[J]. Journal of Geophysical Research, 2007, 112(B9): B09401. DOI: 10.1029/2007JB00 4949.

[9] ALTAMIMI Z, COLLILIEUX X, MéTIVIER L. ITRF2008: An Improved Solution of the International Terrestrial Reference Frame[J]. Journal of Geodesy, 2011, 85(8): 457-473. DOI: 10.1007/s00190-011-0444-4.

[10] 明鋒, 楊元喜, 曾安敏, 等. 顧及有色噪聲的GPS位置時間序列中斷探測法[J]. 武漢大學學報(信息科學版), 2016, 41(6): 745-751. MING Feng, YANG Yuanxi, ZENG Anmin, et al. Offset Detection in GPS Position Time Series with Colored Noise[J]. Geomatics and Information Science of Wuhan University, 2016, 41(6): 745-751.

[11] KOCH K R, YANG Y. Konfidenzbereiche und Hypothesentests für Robuste Parametersch?tzungen[J]. Zeitschrift für Vermessungswesen, 1998, 123(1): 20-26.

[12] KOCH K R. Parameter Estimation and Hypothesis Testing in Linear Models[M]. Berlin: Springer -Verlag, 1988: 271-307.

[13] SERGEI N, RODIONOV. A Sequential Algorithm for Testing Climate Regime Shifts[J]. Geophysical Research Letters, 2004, 31(9): L09204.

[14] MING Feng, YANG Yuanxi, ZENG Anmin, et al. Analysis of Seasonal Signals and Long-term Trends in the Height Time Series of IGS Sites in China[J]. Science China Earth Sciences, 2016, 59(6): 1283-1291. DOI: 10.1007/s11430-016-5285-9.

[15] KENDALL M G. Rank Correlation Measures[M]. London: Charles Griffin, 1975: 202.

(責任編輯：叢樹平)

Crustal Subsidence Analysis from GNSS Data for Dagang Tidal Station in Qingdao

WU Fumei1,2，WEI Ziqing1,2，LIU Guangming1,2

1. State Key Laboratory of Geo-information Engineering, Xi’an 710054, China; 2. Xi’an Research Institute of Surveying and Mapping, Xi’an 710054, China

Crustal subsidence at Dagang tidal station affects the Huanghai mean sea level which defines the national height datum. The crustal subsidence at Dagang tidal station located 1.5 km away from Qingdao GNSS fiducial station is analyzed, by using approximately 10 years of GNSS observations. Firstly, the daily relax solutions and daily constrained solutions for coordinates of this station are computed in a combined adjustment involving a global GNSS net of 50 IGS stations and a regional net of 43 CMONOC stations. Then, a time series spanning 10 years for station coordinates aligned to ITRF2014 is obtained. And next, the time series of height components is analyzed to detect outliers, offsets, trends and frequencies and the amplitudes of periodic parts and the values of offsets are obtained. It is shown that no trend is found in the series. But it is found that the annual and semi-annual signals influence the changes of the series obviously. Finally, it is concluded in combination with the analysis of tidal data for Dagang tidal station that the mean sea level at Dagang tidal station is ascending with a speed of 1.62 mm/a.

height datum; GNSS station; time series; crust subside

The National Natural Science Foundation of China (Nos. 41374003; 41474015; 41674025); The National Key Research and Development Program of China (No. 2016YFB0501701)

WU Fumei(1981—)，female, PhD, majors in geodesy datum research and kinematic geodetic data processing.

吳富梅，魏子卿，劉光明.利用GNSS數(shù)據(jù)分析大港驗潮站地殼沉降[J].測繪學報，2017,46(4)：430-435.

10.11947/j.AGCS.2017.20160506. WU Fumei，WEI Ziqing，LIU Guangming.Crustal Subsidence Analysis from GNSS Data for Dagang Tidal Station in Qingdao[J]. Acta Geodaetica et Cartographica Sinica,2017,46(4):430-435. DOI:10.11947/j.AGCS.2017.20160506.

P224

A

1001-1595(2017)04-0430-06

國家自然科學基金(41374003; 41474015; 41674025); 國家重點研發(fā)計劃(2016YFB0501701)

2016-10-13

吳富梅(1981—)，女，博士，研究方向為大地基準研究和動態(tài)大地測量數(shù)據(jù)處理。

E-mail： wfm8431812@163.com

修回日期： 2016-12-15