上海CORS站高程數(shù)據(jù)非線性建模方法研究

姜 濤

(廣東省核工業(yè)地質(zhì)局二九二大隊(duì)，廣東 河源 517001)

0 引 言

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

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

目前GPS測(cè)量技術(shù)已經(jīng)應(yīng)用于交通、建筑、導(dǎo)航、定位等領(lǐng)域，起到了越來(lái)越重要的作用[10]。為了對(duì)CORS站高程坐標(biāo)建立符合其實(shí)際運(yùn)動(dòng)規(guī)律的數(shù)學(xué)模型，本文選取了國(guó)內(nèi)上海CORS站10余年高程數(shù)據(jù)作為研究對(duì)象，對(duì)其建立非線性速度場(chǎng)模型[11]。建立非線性速度場(chǎng)模型的前提是先給CORS站高程數(shù)據(jù)建立線性最小二乘模型，其次求出線性最小二乘模型的解，最后將其作為非線性最小二乘模型的迭代初值進(jìn)行未知參數(shù)的求解，得到反映CORS站高程坐標(biāo)運(yùn)動(dòng)規(guī)律的非線性速度場(chǎng)模型。

1 CORS站高程數(shù)據(jù)建模方法

1.1 線性最小二乘建模原理

通過(guò)對(duì)CORS站高程坐標(biāo)時(shí)間序列的時(shí)序圖分析發(fā)現(xiàn)，CORS站高程方向運(yùn)動(dòng)規(guī)律存在明顯的年周期運(yùn)動(dòng)。由于不同因素的影響，CORS站高程坐標(biāo)運(yùn)動(dòng)周期性變化中可能存在的月、半年、兩年變化在時(shí)序圖中不能被準(zhǔn)確發(fā)現(xiàn)，并且不同的CORS站其周期性變化規(guī)律也不盡相同[12]。為了得到符合各個(gè)CORS站高程運(yùn)動(dòng)規(guī)律的數(shù)學(xué)模型，需進(jìn)行多次的數(shù)據(jù)擬合。包含年、半年周期項(xiàng)的線性最小二乘擬合模型如式(1)所示。

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

(1)

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

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

V=AN-Z

(2)

式(2)中，

(3)

(4)

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

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

非線性最小二乘擬合模型是在線性模型(1)的基礎(chǔ)上加入相應(yīng)的初相值，且周期項(xiàng)不是固定的數(shù)值[13]。具體模型表達(dá)式如式(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為未知參數(shù)，記為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為單位矩陣，求解未知參數(shù)S的值。

2 CORS站高程數(shù)據(jù)建模實(shí)驗(yàn)分析

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

圖1 上海站高程坐標(biāo)時(shí)序圖

2.1 線性最小二乘建模

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

圖2 上海站高程坐標(biāo)線性最小二乘擬合圖

表1 線性最小二乘建模參數(shù)值

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

2.2 非線性最小二乘建模預(yù)測(cè)

2.2.1 非線性最小二乘建模

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

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

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

圖3 上海CORS站高程坐標(biāo)非線性速度場(chǎng)擬合模型圖

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

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

2.2.2 非線性最小二乘預(yù)測(cè)

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

圖4 上海CORS站高程坐標(biāo)非線性速度場(chǎng)預(yù)測(cè)結(jié)果圖

表4 非線性最小二乘預(yù)測(cè)結(jié)果參數(shù)值

3 結(jié) 語(yǔ)

本文對(duì)上海CORS站2005～2017年高程數(shù)據(jù)進(jìn)行了建模方法的研究，首先分析對(duì)比了線性最小二乘建模和非線性最小二乘建模方法的原理，然后采用CORS站高程數(shù)據(jù)進(jìn)行了建模預(yù)測(cè)實(shí)驗(yàn)的分析，得出以下結(jié)論。

(1)線性最小二乘建模采用固定周期項(xiàng)的方法與上海CORS站高程數(shù)據(jù)實(shí)際周期變化值有一定的誤差，導(dǎo)致了建立的數(shù)學(xué)模型不能恰當(dāng)?shù)胤从掣叱虜?shù)據(jù)實(shí)際的運(yùn)動(dòng)變化規(guī)律。

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

(3)以殘差均方根誤差為建模精度評(píng)價(jià)標(biāo)準(zhǔn)得到，非線性最小二乘擬合建模精度較線性最小二乘擬合建模精度提高了32%，且非線性模型預(yù)測(cè)誤差在1.7 cm以內(nèi)，均方根誤差為0.008 7 m，預(yù)測(cè)精度高。