上海CORS站高程數據非線性建模方法研究

姜 濤

(廣東省核工業地質局二九二大隊，廣東 河源 517001)

0 引 言

國家連續運行參考站(continuously operating referencestations，CORS)目前已在國內大部分省市建立，且積累了大量的連續觀測數據[1-2]。這些數據為各省市的測量、導航、地震研究等提供了重要的數據信息，其中CORS站高程坐標數據的信息最為重要[3-5]。為此，研究CORS站高程坐標方向的運動規律具有一定的實際價值和意義。

從國內不同的CORS站高程坐標時序圖可以看出，CORS站高程坐標運動規律存在周期性變化，且不同的CORS站其高程坐標運動規律是不同的[6]。造成不同CORS站高程坐標運動周期性變化規律的因素有很多種，其中包括地球潮汐、地殼運動、大氣層變化、GPS測量技術等[7]。先前對CORS站高程數據建立數學模型時，常采用線性最小二乘的方法，該方法事先將固定的整年周期值賦給線性最小二乘模型周期項，導致建立的模型周期值與實際周期值存在一定的偏差[8]。這種建立CORS站高程數據數學模型的方法明顯地認為CORS站高程坐標運動只存在線性運動，而忽視了其非線性運動特征[9]。

目前GPS測量技術已經應用于交通、建筑、導航、定位等領域，起到了越來越重要的作用[10]。為了對CORS站高程坐標建立符合其實際運動規律的數學模型，本文選取了國內上海CORS站10余年高程數據作為研究對象，對其建立非線性速度場模型[11]。建立非線性速度場模型的前提是先給CORS站高程數據建立線性最小二乘模型，其次求出線性最小二乘模型的解，最后將其作為非線性最小二乘模型的迭代初值進行未知參數的求解，得到反映CORS站高程坐標運動規律的非線性速度場模型。

1 CORS站高程數據建模方法

1.1 線性最小二乘建模原理

通過對CORS站高程坐標時間序列的時序圖分析發現，CORS站高程方向運動規律存在明顯的年周期運動。由于不同因素的影響，CORS站高程坐標運動周期性變化中可能存在的月、半年、兩年變化在時序圖中不能被準確發現，并且不同的CORS站其周期性變化規律也不盡相同[12]。為了得到符合各個CORS站高程運動規律的數學模型，需進行多次的數據擬合。包含年、半年周期項的線性最小二乘擬合模型如式(1)所示。

Yi=A+B×ti+C1×sin(2π×ti)+C2×sin(4π×ti)

(1)

式(1)中，Yi為模型擬合值，A為常數項，B為線性速度值，ti代表時間，C1表示一年周期項對應的振幅值，C2表示半年周期項對應的振幅值，周期值分別為T1=1,T2=0.5，該周期值表示線性模型給定的固定周期值為一年和半年。

A,B,C1,C2為4個未知數，設CORS站高程坐標序列觀測值為Z。利用間接平差原理得到的誤差方程如式(2)所示。

V=AN-Z

(2)

式(2)中，

(3)

(4)

按照最小二乘原理VTPV=min，P為單位矩陣，求解未知參數N的值。

1.2 非線性最小二乘建模原理

非線性最小二乘擬合模型是在線性模型(1)的基礎上加入相應的初相值，且周期項不是固定的數值[13]。具體模型表達式如式(5)所示。

Yi=A+B×ti+C1×sin(2π×f1×ti+φ1)+C2×sin(4π×f2×ti+φ2)

(5)

式(5)中A,B,C1,f1,φ1,C2,f2,φ2為未知參數，記為S。

將非線性模型(5)在近似值N0處線性化，線性化后的誤差方程如式(6)所示。

V=BS-Z

(6)

式(6)中，S=[ABC1f1φ1C2f2φ2]T，

qi=sin(2×π×ti+φ1)，wi=2×π×ti×cos(2×π×ti+φ1)，ei=cos(2×π×ti+φ1)，

ri=sin(4×π×ti+φ2)，di=4×π×ti×cos(4×π×ti+φ2)，fi=cos(4×π×ti+φ2)，

最后按照最小二乘原理VTPV=min，P為單位矩陣，求解未知參數S的值。

2 CORS站高程數據建模實驗分析

本文從SOPAC網站下載上海CORS站2005～2017年高程坐標數據，并繪制上海CORS站高程坐標時序如圖1所示。其中，縱軸表示的是上海站高程數據值，單位為米，該高程數據由原始高程數據減去平均值得到。橫軸表示的是高程數據對應的觀測時間為2005～2017年。從圖1上海站高程坐標時序圖可以看出：高程坐標數據隨著時間的變化具有周期性，且不是線性變化，即不存在直線上升或下降的趨勢[14]。這就為高程數據建立數據模型提供了相應的思路。

圖1 上海站高程坐標時序圖

2.1 線性最小二乘建模

通過對上海站高程坐標時序圖分析發現，高程坐標數據隨著時間的變化其周期性可能存在一年、半年變化。為此，本文在建立線性最小二乘模型時，首先假定模型周期值為固定一年、半年數值，即T1=1,T2=0.5。然后按照線性最小二乘模型(1)進行高程數據的擬合，得到的上海站高程數據線性最小二乘擬合結果如圖2所示，其中縱軸表示高程坐標數據(包括觀測高程數據和線性最小二乘擬合建模所求解的高程擬合值)黑色的點表示的是原始高程數據，紅色曲線表示的是線性最小二乘擬合曲線圖。從圖2可以看出，上海站高程坐標運動存在周期性的變化規律，且周期性變化主要是一年變化為主、半年變化占比較少。表1給出了線性最小二乘建模初值、線性最小二乘解及建模精度評價指標均方根誤差的值[15]。

圖2 上海站高程坐標線性最小二乘擬合圖

表1 線性最小二乘建模參數值

由表1線性最小二乘建模參數值可以看出：線性最小二乘擬合建模所需的周期項為給定的周期值，分別是一年、半年數值，由計算出的振幅值C1和C2數值可得，上海CORS站高程數據周期性變化一年周期項對應的振幅值明顯大于半年周期項對應的振幅值，這就說明上海CORS站高程數據周期性變化主要為一年變化，半年變化占比非常少。

2.2 非線性最小二乘建模預測

2.2.1 非線性最小二乘建模

前文中已經提到線性最小二乘建模方法中給定固定周期值的方法與實際CORS站高程數據周期性變化值存在一定的偏差，為此，本文建立了非線性最小二乘擬合模型(5)。該模型以線性最小二乘模型的解作為迭代初值，其中非線性最小二乘模型迭代初值如表2所示。

表2 非線性最小二乘建模迭代初值

將表2中給出的非線性最小二乘模型迭代初值代入非線性最小二乘模型(5)，進行不斷地迭代解算得到上海站高程數據非線性最小二乘擬合建模結果如圖3所示：縱軸代表高程數據值，橫軸表示的是高程值對應的時間。

圖3 上海CORS站高程坐標非線性速度場擬合模型圖

從圖3上海CORS站高程坐標非線性速度場擬合模型圖可以看出：擬合曲線能夠較好地描述CORS站高程數據的變化情況，即擬合曲線能夠反映出CORS站高程數據的周期性變化的規律。并且從圖3中也可以看出，擬合模型計算出來的值與實際觀測的高程值誤差在2 cm以內。為了說明非線性速度場模型可以反映CORS站高程數據周期性變化的實際情況，本文計算出非線性最小二乘擬合模型中的未知參數值，如表3所示。從表3可以看出，非線性最小二乘模型中的未知參數f1，f2的數值并不是整數值，即周期性變化的數值并不是固定的整年或整半年數值，這就說明了非線性最小二乘模型將周期項作為未知數進行迭代是符合CORS站高程數據周期性變化規律的。另外，通過對比非線性最小二乘擬合模型殘差均方根誤差與線性最小二乘擬合模型誤差發現，非線性最小二乘擬合模型相較于線性最小二乘擬合模型建模精度提高了32%。

表3 非線性最小二乘建模解算未知參數值

2.2.2 非線性最小二乘預測

為了驗證本文提出的非線性最小二乘速度場建模方法的可行性，作者在本小節以上海CORS站2005～2014年的高程數據利用非線性速度場模型進行擬合建模，然后求解出非線性模型的未知參數，得到非線性速度場模型的數學表達式。將上海CORS站2015～2017年的高程數據值作為實際觀測值，并與非線性預測模型求解出的預測值進行比較，得到的預測結果如圖4所示，前一段綠色的點代表的為2005～2014年的CORS站高程數據值，黑色的曲線代表的為利用非線性最小二乘方法對該段高程數據的擬合。后一段黑色的點代表的是待預測的2015～2017年的CORS站高程數據觀測值，紅色的曲線為非線性速度場模型預測值。相應的非線性模型參數值和預測結果精度值見表4，其中RMSE/m為預測值于實際值誤差的均方根誤差統計值。從圖4的預測曲線可以看出，預測曲線數據值的變化與CORS站高程數據實際值的變化基本一致，且預測時間越短，預測的精度越高，但從表4可看出總體預測誤差在1.7 cm以內，說明本文建立的非線性速度場模型具有一定可行性與正確性。

圖4 上海CORS站高程坐標非線性速度場預測結果圖

表4 非線性最小二乘預測結果參數值

3 結 語

本文對上海CORS站2005～2017年高程數據進行了建模方法的研究，首先分析對比了線性最小二乘建模和非線性最小二乘建模方法的原理，然后采用CORS站高程數據進行了建模預測實驗的分析，得出以下結論。

(1)線性最小二乘建模采用固定周期項的方法與上海CORS站高程數據實際周期變化值有一定的誤差，導致了建立的數學模型不能恰當地反映高程數據實際的運動變化規律。

(2)上海CORS站高程數據線性最小二乘和非線性最小二乘擬合建模結果圖都反映出高程運動存在的周期性變化以一年變化為主，半年變化占比較少。

(3)以殘差均方根誤差為建模精度評價標準得到，非線性最小二乘擬合建模精度較線性最小二乘擬合建模精度提高了32%，且非線性模型預測誤差在1.7 cm以內，均方根誤差為0.008 7 m，預測精度高。