傅彥博，孫付平，朱新慧，劉 婧

信息工程大學(xué)，河南 鄭州 450001

隨著各類空間大地測量技術(shù)的迅猛發(fā)展，毫米級地球參考框架正逐步地建立并應(yīng)用于大地測量領(lǐng)域中。溫度變化被發(fā)現(xiàn)為GPS測站垂向周年變化作出了一定比例的毫米量級貢獻(xiàn)[1-3]；季節(jié)性溫度變化引發(fā)的測站水泥墩與基巖的熱脹冷縮被發(fā)現(xiàn)與測站周年變化存在很強的正相關(guān)性，會造成測站高程方向的季節(jié)性變化[2-8]；大氣負(fù)荷、水文負(fù)載、非潮汐海洋負(fù)載等負(fù)載變化是引起測站垂向周年變化的主要因素[9-17]；高階電離層延遲是造成中緯度測站垂向周年，東西方向半周年以及低緯測站垂向半周年運動的主要原因之一[18-20]；文獻(xiàn)[21]統(tǒng)計了全球492個GPS測站坐標(biāo)時間序列的周期項，發(fā)現(xiàn)周年項在全球具有普遍性；文獻(xiàn)[22]發(fā)現(xiàn)了南北半球周年項初始相位存在約180°的系統(tǒng)差。由于GPS測站的非線性變化成因復(fù)雜，目前國際上還未建立包含多種機制影響的非線性變化理論改正模型。若基于現(xiàn)有全球GPS實測數(shù)據(jù)，分析測站垂向坐標(biāo)時間序列中周期項的全球分布規(guī)律，并根據(jù)此規(guī)律構(gòu)建全球測站垂向非線性變化統(tǒng)計改正模型，將會從整體上削弱全球測站垂向非線性運動，并為從全球角度研究測站垂向非線性變化提供參考。

本文在全球范圍內(nèi)選取了461個連續(xù)觀測時段在2年以上的GPS測站(其中北半球357個，南半球104個)，提取了測站垂向周年項的振幅與初始相位，分析得出了周年項振幅與初始相位的全球分布規(guī)律函數(shù)，基于此函數(shù)構(gòu)建了全球GPS測站垂向周年變化統(tǒng)計改正模型，并對模型的改正效果進(jìn)行了測試與評估。

1 數(shù)據(jù)源與提取周期項

1.1 數(shù)據(jù)源

本文采用的測站坐標(biāo)殘差數(shù)據(jù)可通過ITRF2008官方網(wǎng)站(http:∥itrf.ensg.ign.fr/ITRF_solutions/2008/)下載。為了便于后續(xù)數(shù)據(jù)分析，本文對原始數(shù)據(jù)進(jìn)行了粗差剔除與零均值化處理。

1.2 周期項提取

根據(jù)諧波分析的基本思想，含有若干周期項的時間序列可以看成由一系列不同頻率的周期波疊加而成，因此只需要確定某一周期波的準(zhǔn)確頻率以及參考?xì)v元，就能計算出其對應(yīng)的振幅和初始相位[3]。本文采用周期圖法對測站垂向坐標(biāo)殘差時間序列進(jìn)行了頻譜分析，探測并統(tǒng)計了周期項。鑒于篇幅，只列出了BJFS站垂向坐標(biāo)殘差時間序列的頻譜圖，如圖1所示。

圖1 BJFS站垂直方向坐標(biāo)時間序列的頻率譜Fig.1 Frequency spectrum of vertical coordinates time series of BJFS station

圖1中能量最高頻點的頻率為0.002 79 Hz，可以得出所對應(yīng)的周期為358.4 d，接近于一周年。確定了周期項所對應(yīng)的頻率后，可以利用最小二乘擬合的方法得出周期項的振幅和初始相位[23]。

1.3 全球GPS測站垂向時間序列周期項統(tǒng)計

本文對全球461個GPS測站垂向時間序列所包含的各類周期項進(jìn)行了統(tǒng)計，如表1所示。

表1 全球GPS測站各類周期項統(tǒng)計

在實際探測中，發(fā)現(xiàn)測站時間序列中很少存在嚴(yán)格的周年項，主要為類周年項，因此本文在提取周年振幅與初始相位時，將周期介于350～380 d 之間的周期項都視為周年項[21-24]。

從表1可以看出，周年項在全球具有普遍性，因此研究周年項振幅與初始相位的全球分布規(guī)律，是建立全球測站垂向非線性變化統(tǒng)計改正模型的重要手段。

2 全球測站周年項分布規(guī)律

建立基于周年項的統(tǒng)計改正模型，需要分析周年項振幅與初始相位的全球分布規(guī)律?？紤]到南北半球存在以下差異：①南北半球存在非對稱構(gòu)造變化[25]。②南北半球的海洋與陸地面積比例不同，導(dǎo)致南北半球的氣候變化、氣象條件、GPS測站部署有所差異[26]。因此，分別對南北半球測站垂向周年項的全球分布規(guī)律進(jìn)行研究。

2.1 周年項振幅全球分布規(guī)律

圖2和圖3分別展示了北半球、南半球測站垂向周年項振幅隨緯度的分布規(guī)律。從中可看出，南北半球垂向周年項振幅隨緯度的總體變化趨勢都呈現(xiàn)了單一的凹凸性，且總體趨勢較為平緩。由于二次項函數(shù)的凹凸性單一，不存在拐點，且函數(shù)值隨著橫坐標(biāo)的變化不會出現(xiàn)劇烈的震蕩，結(jié)合振幅隨緯度的分布特征，對垂向周年項振幅全球分布規(guī)律進(jìn)行建模時，用二次項擬合的方法最合適；考慮到函數(shù)在赤道處的連續(xù)性，使用分段擬合的方法，最終得出了南北半球振幅的分布規(guī)律函數(shù)。北半球垂向周年項振幅分布規(guī)律函數(shù)為

AN(B)=-0.002 417B2+0.206B+3.2

(1)

式中，B為緯度，下同。

南半球垂向周年項振幅分布規(guī)律函數(shù)為

AS(B)=-0.001 378B2+0.091B+3.2

(2)

對南北半球的振幅分布規(guī)律函數(shù)進(jìn)行了t檢驗(顯著性水平為0.05，零假設(shè)為擬合函數(shù)值與振幅值來自同一分布)，檢驗結(jié)果如表2所示。

表2 南北半球振幅分布規(guī)律函數(shù)顯著性檢驗結(jié)果

經(jīng)過顯著性檢驗，在顯著性水平為0.05的條件下，南北半球振幅擬合函數(shù)的H值都為0，則接受零假設(shè)，即模型擬合值與真實振幅值為同一分布，說明二次項擬合函數(shù)能準(zhǔn)確刻畫振幅的分布特征。

從圖2和圖3中可以看出兩個明顯的現(xiàn)象：①隨著緯度的增加，南北半球測站垂向周年項振幅呈現(xiàn)先增大后減小的規(guī)律；②北半球的垂向周年項振幅普遍高于南半球。

針對第1個現(xiàn)象，由于環(huán)境負(fù)載與測站垂向坐標(biāo)季節(jié)性溫度變化存在強相關(guān)性[15]，對于季節(jié)性變化劇烈的中緯度地區(qū)的GPS測站，其環(huán)境負(fù)載導(dǎo)致的垂向位移普遍高于其他地區(qū)；同時，環(huán)境負(fù)載是引起測站垂向周年變化的主要因素[9]。因此，中緯度地區(qū)的周年振幅高于低、高緯度測站。

針對第2個現(xiàn)象，由于地表環(huán)境負(fù)載是引起測站垂向周年變化的主要因素，同時北半球測站環(huán)境負(fù)載的垂向絕對平均坐標(biāo)值普遍高于南半球測站[9]，因此導(dǎo)致北半球測站的周年項振幅普遍高于南半球測站。

2.2 周年項初始相位全球分布規(guī)律研究

提取周年項的初始相位時，需要統(tǒng)一各測站垂向時間序列的參考?xì)v元[3]。在兼顧時間序列長度與測站數(shù)目的基礎(chǔ)上，本文選用2002年1月1日為參考?xì)v元。由于部分測站時間序列的觀測時段不包含參考?xì)v元，因此參與初始相位全球分布規(guī)律分析的北半球、南半球測站分別為298、89個，測站數(shù)目總和占全球測站數(shù)目的84%。

圖4和圖5分別展示了北半球、南半球測站垂向周年項初始相位隨緯度的分布規(guī)律。

根據(jù)南北半球垂向周年項初始相位的分布特征，同2.1節(jié)，使用分段二項式擬合的方法得出了垂向周年項初始相位的全球分布規(guī)律函數(shù)。北半球垂向周年項初始相位分布規(guī)律函數(shù)為

φN(B)=-0.000 624B2+0.051 2B-0.25

(3)

南半球垂向周年項初始相位分布規(guī)律函數(shù)為

φS(B)=0.000 615 8B2-0.045 57B-0.25

(4)

對南北半球的初始相位分布規(guī)律函數(shù)進(jìn)行了t檢驗(顯著性水平為0.05，零假設(shè)為擬合函數(shù)值與振幅值來自同一分布)，檢驗結(jié)果如表3所示。

表3 南北半球初始相位分布規(guī)律函數(shù)顯著性檢驗結(jié)果

根據(jù)表3可知，模型擬合值與真實初始相位值為同一分布，說明二次項擬合函數(shù)能準(zhǔn)確刻畫初始相位的分布特征。

圖5、圖6之間最大的差異在于擬合曲線的凹凸。在中緯度地區(qū)，南北半球測站垂向周年項初始相位之間存在約180°的系統(tǒng)差，即垂向周年位移趨勢相反。由于南北半球測站的季節(jié)性溫度變化相反，測站環(huán)境負(fù)載與測站季節(jié)性溫度變化存在強相關(guān)性[15]，且環(huán)境負(fù)載是導(dǎo)致測站垂向周年變化的主要因素[9]，因此在同一時間段內(nèi)，季節(jié)性溫度變化最劇烈的中緯度地區(qū)會呈現(xiàn)南北半球測站垂向周年位移趨勢相反的現(xiàn)象，低、高緯度地區(qū)測站的季節(jié)性變化相對較弱，因此其垂向周年位移不會呈現(xiàn)明顯的相反趨勢。

3 建模及改正效果分析

3.1 建立統(tǒng)計改正模型

在只考慮周年項的條件下，統(tǒng)計改正模型時間序列M(t)為

M(t)=A(B)·cos[ω1at+φ(B)]

(5)

GPS測站時間序列采樣周期為7 d/次，可以得出周年角頻率ω1a

(6)

根據(jù)全球GPS測站垂向周年項振幅與初始相位的全球分布規(guī)律函數(shù)，建立了全球GPS測站垂向周年變化統(tǒng)計改正模型。北半球的GPS測站垂向周年變化統(tǒng)計改正模型為

(7)

式中，B為測站緯度；t0為測站時間序列的起始?xì)v元；t為擬合值所對應(yīng)的歷元，下同。

南半球的GPS測站垂向周年變化統(tǒng)計改正模型為

(8)

3.2 模型改正效果分析

利用建立的統(tǒng)計改正模型對全球461個測站的垂向時間序列進(jìn)行了改正，并且根據(jù)時間序列的RMS對改正效果進(jìn)行了評價。

設(shè)原始時間序列為O(t)，模型序列為M(t)，改正后的剩余殘差序列為R(t)，原始時間序列所包含歷元數(shù)為n。則原始時間序列O(t)的RMS為

(9)

改正模型序列M(t)的RMS為

(10)

剩余殘差序列R(t)的RMS為

(11)

本文利用了模型改正比例km及剩余殘差比例kr，評估模型改正效果

(12)

本文在全球低、中、高緯度各選取了具代表性的4個測站，分別統(tǒng)計了模型改正比例及剩余殘差比例，如表4所示。

表4 所選各測站信息及改正效果統(tǒng)計

圖6是低緯度4個測站的改正效果(其中藍(lán)色線表示原始時間序列，黑色線表示模型改正序列， 紅色線表示剩余殘差序列，下同)。圖7是中

緯度4個測站的改正效果。圖8是高緯度4個測站的改正效果。

圖2 北半球垂向周年項振幅隨緯度分布規(guī)律Fig.2 Distribution law of vertical annual terms’ amplitudes of northern hemisphere with latitude

圖3 南半球垂向周年項振幅隨緯度分布規(guī)律Fig.3 Distribution law of vertical annual terms’ amplitudes of southern hemisphere with latitude

圖4 北半球垂向周年初始相位隨緯度分布規(guī)律Fig.4 Distribution law of vertical annual terms’ initial phases of northern hemisphere with latitude

圖5 南半球垂向周年項初始相位隨緯度分布規(guī)律Fig.5 Distribution law of vertical annual terms’ initial phases of southern hemisphere with latitude

圖6 低緯度測站改正效果Fig.6 Correction effect of stations located in low latitudes

圖7 中緯度測站改正效果Fig.7 Correction effect of stations located in middle latitudes

圖8 高緯度測站改正效果Fig.8 Correction effect of stations located in high latitudes

結(jié)合表4與圖6—圖8，可以看出中緯度測站的平均模型改正比例最大，平均剩余殘差比例最小，總體改正效果最好；低、高緯度測站的改正效果均不如中緯度測站?？赡苡幸韵聝牲c原因：①中緯度地區(qū)由于季節(jié)性溫度效應(yīng)較為明顯[6-7]，環(huán)境負(fù)載、熱膨脹效應(yīng)等因素與季節(jié)性溫度變化存在強相關(guān)性[7-9]，且環(huán)境負(fù)載是導(dǎo)致測站垂向周年運動的主要因素[9]，因此垂向周年變化是中緯度地區(qū)測站垂向位移的主要周期性變化；低、高緯度測站由于受季節(jié)性溫度效應(yīng)影響較小，其他周期性變化在測站垂向位移中的占比高，導(dǎo)致中緯度測站垂向周年振幅隨緯度的分布比低、高緯度測站更為集中，改正效果更好。②與原因①同理，中緯度測站垂向周年初始相位值分布比低、高緯度測站分布更為集中，因此與低、高緯度測站相比，中緯度測站的改正效果更好。

從圖6—圖8可看出，測站原始時間序列的主周期變化與模型序列中的周年變化并不完全契合，因此原始序列中非線性變化的不是嚴(yán)格的周年變化，而是類周年變化，印證了文獻(xiàn)[27]的結(jié)論。

本文統(tǒng)計了全球461個GPS測站的改正效果，如圖9和圖10所示。

圖9 全球GPS測站模型改正比例統(tǒng)計Fig.9 Statistics of model correction ratio of global GPS stations

圖10 全球GPS測站剩余殘差比例統(tǒng)計Fig.10 Statistics of residuals ratio of global GPS stations

根據(jù)統(tǒng)計結(jié)果，在所有測站中，剩余殘差比例達(dá)50%～70%的測站占多數(shù)，可知統(tǒng)計改正模型能削弱多數(shù)測站30%～50%的垂向坐標(biāo)殘差。在模型改正比例達(dá)30%～40%的測站中，存在一部分模型改正比例達(dá)到130%～140%的測站，是因為此類測站的垂向周年項振幅小于模型擬合值。綜合分析統(tǒng)計結(jié)果，改正模型能大幅度削弱GPS測站的垂向坐標(biāo)殘差并提升大規(guī)模測站數(shù)據(jù)處理效率。

3.3 統(tǒng)計改正模型評價

3.3.1 模型優(yōu)點

(1) 本文建立的統(tǒng)計改正模型能削弱全球多數(shù)測站30%～50%的垂向坐標(biāo)殘差，精確分離測站的線性與非線性運動。

(2) 模型能統(tǒng)一對全球測站垂向坐標(biāo)進(jìn)行改正，提高了數(shù)據(jù)處理效率。

(3) 模型能直觀地顯示測站垂向周年項振幅以及初始相位的全球分布規(guī)律，為研究全球范圍內(nèi)周年性地球物理運動提供參考。

3.3.2 模型缺點

(1) 由于全球測站所處地區(qū)的地球物理條件各異，引起周年項的因素復(fù)雜，部分測站的周年項振幅與模型擬合值相差較大，導(dǎo)致其改正效果不理想。

(2) 二次曲線受局部異常值影響較大，異常值的存在會影響二項式擬合效果。

(3) 由于數(shù)據(jù)觀測時段限制，有16%的測站不包含用于提取初始相位的參考?xì)v元，導(dǎo)致建模時數(shù)據(jù)利用率降低，影響此類測站的改正效果。

3.3.3 模型改進(jìn)展望

(1) 半周年項也較為普遍地存在于全球GPS測站垂向坐標(biāo)殘差時間序列中?？偨Y(jié)半周年項的全球分布規(guī)律，優(yōu)化現(xiàn)有統(tǒng)計改正模型，對一些含有半周年項的測站將有更好的改正效果。

(2) 研究測站水平方向非線性變化項的全球分布規(guī)律，建立GPS測站三維非線性變化統(tǒng)計改正模型。

(3) 由于二項式擬合對異常值較為敏感，今后改進(jìn)模型時，可以考慮使用抗差性更好的方法進(jìn)行更精確的擬合。

(4) 對全球GPS測站數(shù)據(jù)使用最新的精密數(shù)據(jù)處理策略進(jìn)行重新處理，使所有測站的坐標(biāo)時間序列中包含相同的參考?xì)v元，得出更精確的周年項初始相位分布規(guī)律函數(shù)，精化全球GPS測站垂向周年變化統(tǒng)計改正模型。

4 結(jié) 論

本文基于實測數(shù)據(jù)，研究發(fā)現(xiàn)了全球GPS測站垂向坐標(biāo)殘差時間序列中周年項的全球分布規(guī)律，建立了全球GPS測站垂向周年變化統(tǒng)計改正模型，分析了模型的改正效果，得出如下結(jié)論。

(1) 周年項是全球GPS測站中普遍存在的周期項，南北半球周年項振幅、北半球周年項初始相位隨緯度的增加呈現(xiàn)先增后減的規(guī)律，南半球周年項初始相位的變化規(guī)律與北半球相反，南北半球中緯度地區(qū)測站周年項初始相位之間存在約180°的系統(tǒng)差，低、高緯度地區(qū)測站的周年項初始相位之間則不存在明顯系統(tǒng)差。

(2) 本文建立的統(tǒng)計改正模型能削弱全球多數(shù)測站30%～50%的垂向坐標(biāo)殘差，可以精確分離測站的線性與非線性運動，說明建立統(tǒng)計改正模型是必要且可行的。

(3) 測站周年項的周期并不嚴(yán)格等于一周年，各個測站周年項的真實周期存在微小偏差。