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

閆 偉 牛安福

1 中國地震臺(tái)網(wǎng)中心，北京市南橫街5號(hào)，100045

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

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

1 資料概況

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

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

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

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

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

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

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

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

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

圖2 1999～2011年YANC低頻殘差信息的空間分布及其相對(duì)的頻率分布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ù)計(jì)算方法

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

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

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

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

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

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

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

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

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

3 顯著性檢驗(yàn)

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

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

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

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

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

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

4 討 論

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

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

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

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

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

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

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

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

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

致謝：感謝阿爾及爾大學(xué)Sid－Ali Ouadfeul教授的建設(shè)性意見及其與布魯門大學(xué)Lena A Brinkhoff博士有益的討論，感謝中國地殼觀測(cè)網(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)出的全球高程震蕩運(yùn)動(dòng)及季節(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ù)觀測(cè)站的噪聲分析［J］.地震學(xué)報(bào)，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測(cè)量與數(shù)據(jù)處理［M］.武漢：武漢大學(xué)出 版 社，2005（Li Zhenghang，Huang Jinsong.GPS Measurement and Data Processing［M］.Wuhan：Wuhan University Press，2005）

［6］張燕，吳云，施順英，等.GPS時(shí)間序列揭示地震前兆的初步探索［J］.大地測(cè)量與地球動(dòng)力學(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測(cè)量的精度及其應(yīng)用［J］.測(cè)繪通報(bào)，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é)報(bào)，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é)報(bào)，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é)報(bào)，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