GPS連續(xù)觀測站低頻擬合殘差序列空間分布的分形特征

閆 偉 牛安福

1 中國地震臺網(wǎng)中心，北京市南橫街5號，100045

隨著GPS連續(xù)站時間序列資料的積累，時間序列中所含噪聲的性質(zhì)及其定量估計問題也逐漸引起人們的重視，且發(fā)表了許多研究成果。其中一個重要的結(jié)論是，GPS時間序列中包含的絕不僅僅是白噪聲，而且白噪聲可能還不占主要地位。至于其中究竟包含哪些種類噪聲，以及它們的大小及其分布，則認(rèn)識尚有不同。Nikolaidis［1］和Agnew［2］認(rèn)為，白噪聲＋閃爍噪聲模型可能是描述GPS時間序列噪聲比較恰當(dāng)?shù)哪Ｐ停⒒谶@一噪聲模型用最大似然估計方法定量估計了分布在全球100多個GPS連續(xù)站的噪聲分量。而朱文耀等［3］則采用白噪聲＋閃爍噪聲＋隨機(jī)漫步噪聲的模型，同樣用最大似然估計法定量估計了分布在全球的178個GPS連續(xù)站的噪聲分量，但是也發(fā)現(xiàn)，并非所有站都可以解得全部的噪聲分量。黃立人［4］用功率譜密度方法定量估計了模型中各噪聲分量，并表明GPS站坐標(biāo)時間序列中白噪聲甚至不是噪聲的主要成分；同時也指出，在確定觀測結(jié)果中閃爍噪聲和白噪聲的定量估計問題時，除了對噪聲特性的判別正確之外，對時間序列建模的正確性也是一個十分重要的問題。

本文從分形的角度討論GPS連續(xù)站時間序列低頻擬合殘差的空間分布特征，求得各GPS連續(xù)站低頻殘差的分?jǐn)?shù)維及無標(biāo)度區(qū)間，并對GPS噪聲提取模型進(jìn)行討論。

1 資料概況

中國地殼運動網(wǎng)絡(luò)建立的25個GPS連續(xù)站試運行以來，除了個別站（如南海永興島的YONG 站）因為自然條件惡劣導(dǎo)致資料連續(xù)性較差之外，大部分已獲得10a以上的連續(xù)觀測資料。2003年以后又有哈爾濱（HARB）和鄭州（ZHEN）兩站先后納入網(wǎng)絡(luò)工程的基準(zhǔn)站序列。對這些連續(xù)站資料由網(wǎng)絡(luò)的數(shù)據(jù)中心進(jìn)行統(tǒng)一處理，得到單日坐標(biāo)解的時間序列。

時間序列觀測值中難免帶有粗差，考慮到單日坐標(biāo)解的誤差約在0～3mm，本文將觀測誤差大于10mm 的對應(yīng)觀測值進(jìn)行剔除。

時間序列低頻信息的提取方法主要包括多項式最小二乘擬合、小波分析提取趨勢項、移動平均值法等［5－7］，本文采用db5小波函數(shù)8次擬合，求得原始觀測值與小波低頻擬合殘差作為研究對象。

以鹽池（YANC）GPS 基準(zhǔn)站北向坐標(biāo)時間序列為例，提取結(jié)果如圖1示。

圖1 鹽池（YANC）GPS基準(zhǔn)站北分量坐標(biāo)時間序列低頻擬合殘差信息提取Fig.1 The low frequency information of the north component of YANC GPS station

為了研究低頻殘差信息的空間統(tǒng)計特征，本文將北向坐標(biāo)（N）與東向坐標(biāo)（E）在站心坐標(biāo)系內(nèi)描點，得到空間分布圖（圖2）。為了直觀地表述其空間分布特征，本文計算以數(shù)據(jù)幾何中心為圓心，不同半徑圓環(huán)內(nèi)出現(xiàn)的觀測點單位面積的頻率：

式中，N為總點數(shù)，N（r1）、N（r2）分別表示落入內(nèi)外徑r1、r2的圓環(huán)的點數(shù)，fs為單位面積的頻率。

為了直觀地表述其空間分布特征，本文對頻率結(jié)果進(jìn)行歸一化（最大值為1，最小值為0），并進(jìn)行填色顯示，如圖2所示（黑色圓點為低頻殘差序列空間分布；黑色圓環(huán)為統(tǒng)計區(qū)域；顏色條表示單位面積內(nèi)觀測值出現(xiàn)頻率的相對大小，頻率最高值取1，最低值取0）。

從圖2可以看出，YANC 基準(zhǔn)站低頻殘差的空間分布具有極其不均勻性，其單位面積的頻率以幾何中心為密集圓心向外發(fā)散性衰減。本文引入分形方法定量描述其分布特點。

圖2 1999～2011年YANC低頻殘差信息的空間分布及其相對的頻率分布Fig.2 The distribution of the YANC station’s low frequency residual during the period of 1999－2011 and the relative frequency in differentrings area

2 分形參數(shù)計算方法

2.1 根據(jù)測度關(guān)系求分維

實際測定分形維數(shù)的方法大致可以分成如下5類［8］：改變觀察尺度求維數(shù)，根據(jù)測度關(guān)系求維數(shù)，根據(jù)相關(guān)函數(shù)求維數(shù)，根據(jù)分布函數(shù)求維數(shù)，根據(jù)頻譜求維數(shù)。本文采用根據(jù)測度關(guān)系求維數(shù)的方法。

分布于空間的點的集合為以某點為中心、以r為半徑的圓。設(shè)包含于此圓內(nèi)部的點的總數(shù)為M（r），如果點的分布是均勻直線，即有M（r）∝r1；如果點的分布是均勻平面，即有M（r）∝r2。因此，將其一般化，如滿足：

就可以說點的分布的分形維數(shù)就是D。一般情況下，選取空間集合幾何中心或重心作為圓心［9］。

在實際計算中，本文采取以高頻點位幾何中心為圓心，取不同的半徑進(jìn)行個數(shù)統(tǒng)計并求解其空間分形維數(shù)D。

2.2 無標(biāo)度區(qū)間

在分形維數(shù)的數(shù)學(xué)定義中，要求碼尺趨于零時的極限存在。但是對于不同學(xué)科中的分形以及自然界存在的分形，一般來說并不存在無窮嵌套結(jié)構(gòu)，而只存在有限的嵌套層次，所以，碼尺趨于零的這個要求在測量中很難實現(xiàn)。近幾年來國內(nèi)外研究［10－16］表明，對實體分形而言，測量的分形維數(shù)值隨碼尺而變化，也就是說，對同一分形體由于測量精度采樣數(shù)目的原因可能存在分維結(jié)果的不確定性。在實際計算過程中，需要截取線性段較好的區(qū)間作為其無標(biāo)度區(qū)間。本文采用截取線性段的方法是在保證盡量多點的情況下人工拾取線性段。線性段上下限對應(yīng)的圓的半徑即為無標(biāo)度區(qū)范圍。拾取結(jié)果如圖3（a）所示。

代入本文中YANC連續(xù)站低頻殘差空間分布并拾取線性段數(shù)據(jù)（圖3（b）），用最小二乘線性擬合理論計算得到，直線斜率k＝1.807±0.007，若人工拾取的“線性段”符合線性特征，那么此時k即為分維值D及其解算誤差。無標(biāo)度區(qū)間即為該線性段對應(yīng)圓環(huán)的半徑0～2.9mm，落入此區(qū)間的擬合殘差點共2 706個，占總點數(shù)的61.781%。

圖3 YANC基準(zhǔn)站低頻擬合殘差空間分布分維計算及其無標(biāo)度區(qū)拾取結(jié)果Fig.3 The fractal result of YANC site

3 顯著性檢驗

用最小二乘法求線性回歸方程并不需要事先假設(shè)自變量與因變量之間一定存在線性相關(guān)關(guān)系。就最小二乘法本身而言，對于任意試驗數(shù)據(jù)都可確定相應(yīng)的線性方程，只有y與x之間確實存在線性相關(guān)關(guān)系時，用最小二乘法求出的線性回歸方程才能近似地表示它們之間的線性相關(guān)關(guān)系。因此，必須檢驗y與x之間是否符合線性關(guān)系，即進(jìn)行線性回歸的顯著性檢驗。本文采用直線擬合相關(guān)系數(shù)與F統(tǒng)計檢驗零假設(shè)進(jìn)行顯著性檢驗。

相關(guān)系數(shù)計算公式為：

其中，n為樣本數(shù)，x和y分別為兩變量的變量值。代入本文中YANC 擬合殘差分人工截取的線性段數(shù)據(jù)，得到互相關(guān)系數(shù)r＝0.997 25＞0.99，表明兩變量具有較強(qiáng)的相關(guān)性，但并不能說明一定符合線性關(guān)系，還需要進(jìn)行線性回歸方程的顯著性檢驗。

對 于 樣 本 值（x1，y1），（x2，y2），…，（xn，yn），線性回歸方程可寫成：

若假設(shè)式（4）符合實際，則b不應(yīng)為零。因此我們需要檢驗假設(shè)：

本文使用F檢驗方法進(jìn)行檢驗：

4 討 論

基于以上方法，本文計算了全國27 個GPS連續(xù)站的分形參數(shù)及顯著性檢驗參數(shù)，結(jié)果見表1。

從表1 可以看出，27 個GPS 連續(xù)站低頻擬合殘差時間序列空間分布在無標(biāo)度區(qū)內(nèi)均符合分形特征（分維數(shù)為非整數(shù)），且分維值介于1.71和1.97之間。表明GPS殘差序列中包含一部分自然分形規(guī)律作用下的觀測結(jié)果，從分形的角度可以不將其作為觀測噪聲，如果籠統(tǒng)地將其代入噪聲模型進(jìn)行計算，則有可能引入系統(tǒng)誤差。這對于我們更加精細(xì)地認(rèn)識和分析GPS觀測結(jié)果提供了佐證。

為了更加細(xì)致地分析其中的計算結(jié)果，將其中部分?jǐn)?shù)據(jù)以統(tǒng)計圖的形式展示（圖4）。

圖4（a）表示了各個臺站的分維值及其無標(biāo)度區(qū)間。可以看出，各個臺站計算得到的分維值不盡相同，且存在一個波動范圍（1.71～1.97），可能與各臺站所處的位置及環(huán)境不同有關(guān)。黃立人指出，建立模型的完善性可能關(guān)系到結(jié)果的適用性。如果在涉及指數(shù)的模型時我們選取其分維值作為指數(shù)，可能可以減少建模時因模型的不適用而產(chǎn)生系統(tǒng)誤差。

表1 各臺站分形參數(shù)及其顯著性檢驗結(jié)果Tab.1 The results of fractal parameters and thesignificant test of each sites

圖4 各臺站無標(biāo)度區(qū)分形參數(shù)及其顯著性檢驗結(jié)果Fig.4 The fractal parameters in the scaleless range and the result of its significant test

圖4（b）表示各個臺站落入無標(biāo)度區(qū)間的點數(shù)占總點數(shù)的百分比。可以看出，除BJFS、LUZH 和XIAM 外，其余24 個臺站的百分比都大于50%，表明多數(shù)點符合自然分形規(guī)律，即符合冪律為D的分形規(guī)律。

由表1中顯著性檢驗及圖4（c）可以看出，除DXIN、QION 和XIAA 三個連續(xù)站殘差序列無標(biāo)度區(qū)內(nèi)相關(guān)系數(shù)在0.99～0.995之間外，其他24個臺的相關(guān)系數(shù)都達(dá)到了0.995以上；27個連續(xù)站無標(biāo)度區(qū)內(nèi)的F檢驗在α＝0.005時全部拒絕零假設(shè)。說明線性程度較好，能夠利用最小二乘方法對提取的無標(biāo)度區(qū)進(jìn)行線性擬合，也說明在雙對數(shù)坐標(biāo)系下選取的圓的半徑與圓內(nèi)的點的總數(shù)存在較好的線性關(guān)系，即GPS殘差序列存在分形特征，從顯著性檢驗的角度支持以上結(jié)論。

由于分形理論中無標(biāo)度區(qū)（雙對數(shù)坐標(biāo)系下線性段）非常易于拾取，為我們提取GPS信息成分提供了另一種工具，這將有助于我們更加精細(xì)地認(rèn)識和分析GPS觀測結(jié)果。

致謝：感謝阿爾及爾大學(xué)Sid－Ali Ouadfeul教授的建設(shè)性意見及其與布魯門大學(xué)Lena A Brinkhoff博士有益的討論，感謝中國地殼觀測網(wǎng)絡(luò)中心提供GPS基準(zhǔn)站單日解數(shù)據(jù)。

［1］Nikolaidis R.Observation of Geodetic and Seismic Deformation with the Global Positioning System［D］.San Diego：University of California，2002

［2］Agnew D C.The Time－Domain Behavior of Power－Law Noise［J］.Geophys Res Let，1992，19：333－336

［3］朱文耀，符養(yǎng)，李彥.GPS高程導(dǎo)出的全球高程震蕩運動及季節(jié)變化［J］.中國科學(xué)D輯，2003，33：470－481（Zhu Wenyao，F(xiàn)u Yang，Li Yan.Global Elevation Oscillatory Motion and Seasonal Variations of GPS Height Derived［J］.Science in China Series D：Earth Sciences，2003，33：470－481）

［4］黃立人，符養(yǎng).GPS連續(xù)觀測站的噪聲分析［J］.地震學(xué)報，2007，29（2）：197－202（Huang Liren，F(xiàn)u Yang.Analysis on the Noises from Continuously Monitoring GPS Sites［J］.Acta Seismologica Sinica，2007，29（2）：197－202）

［5］李征航，黃勁松.GPS測量與數(shù)據(jù)處理［M］.武漢：武漢大學(xué)出 版 社，2005（Li Zhenghang，Huang Jinsong.GPS Measurement and Data Processing［M］.Wuhan：Wuhan University Press，2005）

［6］張燕，吳云，施順英，等.GPS時間序列揭示地震前兆的初步探索［J］.大地測量與地球動力學(xué)，2005，25（3）：96－99（Zhang Yan，Wu Yun，Shi Shunying，et al.The GPS Time Series Revealed Preliminary Exploration of Earthquake Precursors［J］.Journal of Geodesy and Geodynamics，2005，25（3）：96－99）

［7］岳建平，岳東杰.工程GPS測量的精度及其應(yīng)用［J］.測繪通報，1999（11）：27－29（Yue Jianping，Yue Dongjie.Precision Engineering Measurement of GPS and Its Application［J］.Bulletin of Surveying and Mapping，1999（11）：27－29）

［8］高安秀樹.分?jǐn)?shù)維［M］.北京：地震出版社，1989（Tayayusa H.Fractals［M］.Beijing：Earthquake Press，1989）

［9］張濟(jì)中.分形［M］.北京：清華大學(xué)出版社，2011（Zhang Jizhong.Fractals［M］.Beijing：Tsinghua University Press，2011）

［10］Pande C S.Fractal Characterization of Fracture Surfaces［J］.Acta Metall，1987，35（7）：1 633－1 637

［11］Mu Z Q，Lung C W.Studies on the Fractal Dimension and Fracture Toughness of Steel［J］.Journal of Physics D：Applied Physics，1988，21：848－850

［12］謝和平，陳至達(dá).巖石的連續(xù)損傷力學(xué)模型探討［J］.煤炭學(xué)報，1998，13（1）：33－42（Xie Heping，Chen Zhida.Study of a Rock Model with Continuous Damage Mechanics［J］.Journal of China Coal Society，1998，13（1）：33－42）

［13］董連科，王曉偉，王克鋼，等.分形用于材料斷裂韌性研究的不定性問題［J］.高壓物理學(xué)報，1990（2）：118－129（Dong Lianke，Wang Xiaowei，Wang Kegang，et al.Fractal Is Used to Study the Fracture Toughness of Materials not Qualitative Problem［J］.Chinese Journal of High Pressure Physics，1990（2）：118－129）

［14］Lung C W，Zhang S Z.Fractal Dimension of the Fractured Surface of Materials［J］.Physica D，1989，38：242－245

［15］龍期威.分形圖形周界和面積的關(guān)系［J］.高壓物理學(xué)報，1990（4）：259－262（Long Qiwei.The Relationship between the Fractal Perimeter and Area［J］.Chinese Journal of High Pressure Physics，1990，4（4）：259－262）

［16］Xie H.Studies on Fractal Models of the Micro－Fracture of Marble［J］.Chinese Sci Bulletin，1989，34（15）：1 292－1 296