基于多尺度球面小波解算GPS應(yīng)變場的方法及應(yīng)用

蘇小寧， 孟國杰, 王振

1 中國地震局地震預(yù)測研究所， 北京　100036 2 同濟(jì)大學(xué)測繪與地理信息學(xué)院， 上?！?00092 3 Earthquake Research Institute, University of Tokyo， Tokyo Japan

基于多尺度球面小波解算GPS應(yīng)變場的方法及應(yīng)用

蘇小寧1,2， 孟國杰1*, 王振3

1 中國地震局地震預(yù)測研究所， 北京100036 2 同濟(jì)大學(xué)測繪與地理信息學(xué)院， 上海200092 3 Earthquake Research Institute, University of Tokyo， Tokyo Japan

摘要提出了利用多尺度球面小波解算GPS應(yīng)變場的方法，該方法通過建立不規(guī)則分布的小波基函數(shù)以凸顯測站非均勻分布的特征，獲取不同空間尺度的應(yīng)變場.通過對模擬數(shù)據(jù)進(jìn)行對比分析驗(yàn)證該方法的有效性，模擬結(jié)果顯示，多尺度球面小波應(yīng)變計算方法穩(wěn)健性較好.在位移場中加入高斯隨機(jī)誤差(均值為0.0 mm，標(biāo)準(zhǔn)差為0.5 mm)對應(yīng)變計算結(jié)果的影響可忽略不計.從原有位移場隨機(jī)抽取的60%的數(shù)據(jù)解算應(yīng)變場時，仍然能夠獲取比較可靠的結(jié)果，但如果出現(xiàn)大范圍的數(shù)據(jù)空缺，則會對應(yīng)變結(jié)果產(chǎn)生明顯影響.基于華北地區(qū)GPS測站原始觀測數(shù)據(jù)，通過高精度數(shù)據(jù)處理方法解算了該區(qū)域GPS速度場，根據(jù)多尺度球面小波方法，解算了該區(qū)域的應(yīng)變率場及其誤差，并分析了其空間分布特征.

關(guān)鍵詞球面小波； 多尺度； GPS應(yīng)變場； 跨斷裂GPS剖面； 模擬分析

1引言

隨著全球范圍大量GPS測站的布設(shè)與加密，GPS觀測資料在研究現(xiàn)今地殼形變及其動力學(xué)機(jī)制等方面得到廣泛應(yīng)用(Gan et al., 2007; Meng et al., 2008; Liang et al., 2013；Wang et al., 2011).由于GPS速度場建立于某個參考基準(zhǔn)，例如，全球參考基準(zhǔn)，區(qū)域參考基準(zhǔn)，選擇不同的參考基準(zhǔn)會導(dǎo)致速度場存在較大差異，而GPS應(yīng)變場與基準(zhǔn)選擇無關(guān).因此，地殼應(yīng)變場的解算在地殼形變和斷層活動性研究方面具有重要地位(Poyraz, 2015; Pearson and Snay, 2014; Craig and Calais, 2014).

GPS應(yīng)變場解算主要有兩種模型：一種是根據(jù)斷層幾何和運(yùn)動特征建立塊體運(yùn)動學(xué)模型，計算塊體內(nèi)部應(yīng)變場，稱為物理模型；另一種是直接利用GPS速度場進(jìn)行計算，稱為數(shù)學(xué)模型(Tape et al., 2009; 葛偉鵬等，2013；Wu et al.,2011)，目前已經(jīng)提出的數(shù)學(xué)模型主要有最小二乘配置、距離加權(quán)、多面函數(shù)、張力樣條、Kriging差值等(武艷強(qiáng)等，2009；江在森和劉經(jīng)南，2010；Shen et al.，1996；Gan et al.，2007；劉曉霞等，2014).由于解算方法不同，即使利用同一數(shù)據(jù)也難以得到相同的應(yīng)變結(jié)果(武艷強(qiáng)等，2009；江在森和劉經(jīng)南，2010)，因此如何把實(shí)際資料反映的不同空間尺度、不同時間尺度的地殼運(yùn)動與變形信息客觀地表達(dá)出來，仍然是值得探究的問題.不同于其他數(shù)學(xué)模型，通過建立不規(guī)則分布的小波基函數(shù)，多尺度球面小波方法能夠充分體現(xiàn)GPS測站不均勻分布的特征；同時也可通過選擇不同的尺度因子，研究不同空間尺度的GPS應(yīng)變場.本文分析了多尺度球面小波原理，探討了利用多尺度球面小波解算GPS應(yīng)變場的方法，通過模擬斷層閉鎖積累的應(yīng)變，驗(yàn)證了該方法的有效性，并利用本文提出的應(yīng)變計算方法研究了華北地區(qū)不同空間尺度的GPS應(yīng)變場分布特征.

2多尺度球面小波解算GPS應(yīng)變場

2.1多尺度球面小波基本原理

球面小波是小波理論從無限的平面空間到有限的球面空間的延伸，是小波分析理論的前沿研究分支，已在地磁場和重力場研究中得到了廣泛應(yīng)用(馬濤等，2003；Holschneider et al., 2003; Chambodut et al., 2005)，同時也可用于建立GPS地殼運(yùn)動速度場和應(yīng)變場模型.首先定義在單位球上的小波母函數(shù)，然后通過作用在母函數(shù)上的仿射變換(平移、伸縮)得到小波基函數(shù).任何二維平面上的小波母函數(shù)通過逆映射變換都可以得到其球面上的表達(dá)式(Bogdanova et al., 2005)，所以平面上的小波母函數(shù)都有其球面的逆映射形式.本文選擇用高斯函數(shù)的差(Difference of Gaussians，DOG) (包伯成等，2009)形成的球面DOG小波作為基函數(shù)，其表達(dá)式如式(1)所示：

(1)

小波母函數(shù)的伸縮使得小波方法具有不同的分辨率.對于單個測站而言，不同時刻的數(shù)據(jù)代表時間分辨率，而同一時刻不同位置的測站則代表空間分辨率.將待處理信號用球面小波的方法在不同尺度

上進(jìn)行分解，進(jìn)而得到對應(yīng)的平滑信息和細(xì)節(jié)信息.

2.2分解尺度的選擇

通常GPS測站在空間上分布不均勻，直接用GPS速度場解算應(yīng)變率場時不同測站的權(quán)重不同，導(dǎo)致在觀測點(diǎn)密集的地方權(quán)重較大，故權(quán)重的選擇對解算結(jié)果產(chǎn)生很大的影響(江在森和劉經(jīng)南，2010).分布不均勻的GPS測站也可視為不同間隔的空間序列，測站的疏密程度代表了不同的空間尺度.通過選擇建立不規(guī)則分布的小波基函數(shù)，多尺度球面小波方法充分考慮了測站不均勻分布的特征.選擇分解尺度的原則是：在測站較密集的區(qū)域分解尺度較大，在測站較稀疏的區(qū)域分解尺度較小.如果選擇等邊三角形作為圖形單元表達(dá)整個球面，那么不同尺度下三角形的平均邊長如表1所示(Wang and Dahlen, 1995)，設(shè)立基函數(shù)需要滿足的條件是在圖形單元覆蓋的區(qū)域至少包含3個GPS測站.當(dāng)只包含3個GPS測站時，最小圖形單元對應(yīng)的尺度因子為最大分解尺度，而最小分解尺度可根據(jù)觀測網(wǎng)的范圍確定.

圖1為實(shí)際分布的GPS測站所確定的最大分解尺度.從該圖可以看出，只有兩個小區(qū)域能夠達(dá)到測站間距為6.89 km的10級分解，測站間距為13.78 km的9級分解較為零散地分布于整個區(qū)域，測站間距為27.55 km的8級分解覆蓋了大部分區(qū)域，整個區(qū)域最小分解尺度可達(dá)到6級，所以小波分解的最大分解尺度可以選擇6～9級，從而在不同區(qū)域獲取不同空間分辨率的應(yīng)變信息.

表1　球面三角形網(wǎng)格的幾何特征

2.3模型參數(shù)的解算

在球面上GPS速度場可表示為

(2)

已知任何平方可積空間的標(biāo)量函數(shù)f(f∈L2(S2))可以寫成

圖1　由GPS站分布確定的球面小波最大分解尺度Fig.1　Maximum resolved scale of spherical wavelet by the locations of GPS sites

(3)

因此，利用式(3)將式(2)的分解問題轉(zhuǎn)換為3個標(biāo)量函數(shù)的求解問題.把式(3)帶入式(2)可得

(4)

2.4應(yīng)變場及其誤差解算

Savage等(2001)、孟國杰等(2008)均給出了球坐標(biāo)系下三維位移與應(yīng)變的微分表達(dá)式.考慮到地球自由表面的應(yīng)變計算實(shí)質(zhì)仍然是二維問題，而GPS測站不在同一球面，因此在利用三維應(yīng)變計算公式解算二維球面應(yīng)變時需要將不在同一個球面的測站速度投影到同一球面，但由于參與應(yīng)變計算的測站很難滿足同量值的垂向運(yùn)動，因此如何投影是非常復(fù)雜的問題；此外，相對水平形變產(chǎn)生的應(yīng)變，同量級的垂向形變對水平應(yīng)變的影響可忽略不計.基于此兩點(diǎn)，本文暫不考慮垂向速度對水平向應(yīng)變率的影響，此時應(yīng)變張量計算如式(5)所示：

(5)

式中εθ=?uθ/? θ，εφ=?uφ/? φ，εθ φ=?uθ/? φ和εφ θ=?uφ/? θ為速度梯度張量，r是地球半徑，θ是余緯，φ是經(jīng)度.值得注意的是，式(5)定義的坐標(biāo)系以南向?yàn)檎?，而GPS速度場通常以北向?yàn)檎?，在計算時需要對不同坐標(biāo)系的速度和應(yīng)變張量進(jìn)行坐標(biāo)旋轉(zhuǎn).

(6)

3應(yīng)變計算方法有效性驗(yàn)證

利用位錯模型分別正演走滑、逆沖和傾滑三類斷層在地表產(chǎn)生的位移場和應(yīng)變場，然后采用球面小波方法，由位移場數(shù)據(jù)計算應(yīng)變張量.Baxter等(2011)認(rèn)為走滑斷層在斷層附近能夠引起與計算方法無關(guān)的應(yīng)變失真現(xiàn)象，為檢驗(yàn)本文應(yīng)變計算方法的有效性，僅對逆沖和傾滑兩類斷層進(jìn)行分析.限于篇幅，本文僅展示對傾滑斷層的分析結(jié)果.

由位錯模型正演的應(yīng)變稱為“理論應(yīng)變”.由正演的位移場，利用本文的應(yīng)變計算方法解算的應(yīng)變?yōu)椤坝嬎銘?yīng)變”.通過“理論應(yīng)變”和“計算應(yīng)變”對比，分析本文應(yīng)變計算方法的有效性.對于逆沖和傾滑斷層，均考慮位移場加入隨機(jī)誤差和不加入隨機(jī)誤差、測站空間均勻分布和不均勻分布等不同情況，其目的是研究本文應(yīng)變計算方法對誤差的敏感程度，以及數(shù)據(jù)分布不均勻性對應(yīng)變計算結(jié)果的影響.其中增加誤差的方法是在正演的地表位移場中加入均值為0.0mm，標(biāo)準(zhǔn)差為0.5mm的高斯隨機(jī)白噪聲.產(chǎn)生測站空間不均勻分布的位移場的方法有兩種，一是從正演的均勻分布的位移場中隨機(jī)抽取其中的60%，二是在該區(qū)域有GPS測站的位置正演位移場(圖2).

3.1基于位錯模型正演的位移場和應(yīng)變場

設(shè)定一條傾角為45°的傾滑斷層，斷層起點(diǎn)為(115°E，41°N)，終點(diǎn)為(120°E，40°N)，寬度為40 km，假定地表至地下15 km處完全鎖定，15 km以下存在20 mm的逆沖和20 mm左旋的滑動量.利用均勻各向同性半無限空間的位錯模型(Okada, 1985)，計算由于斷層閉鎖在地表產(chǎn)生的采樣間隔為 0.25°×0.25°位移場和應(yīng)變場.圖3為傾滑斷層位錯在地表產(chǎn)生的主應(yīng)變/面應(yīng)變和最大剪應(yīng)變.為顯示清晰，對主應(yīng)變進(jìn)行了重采樣，從圖中可以看出梯度較大的應(yīng)變場主要分布在斷層附近.

3.2利用多尺度球面小波方法解算的應(yīng)變場

在構(gòu)建DOG小波基函數(shù)時，需要通過對尺度因子增加比例系數(shù)，以調(diào)節(jié)基函數(shù)的強(qiáng)度和覆蓋范圍，使得應(yīng)變場既能夠分布光滑，也能夠突出小區(qū)域的細(xì)節(jié)，該系數(shù)類似于高斯距離加權(quán)法的平滑半徑和張力樣條的張力參數(shù).Tape等(2009)認(rèn)為該系數(shù)的選擇范圍應(yīng)為1～2，通過多次重復(fù)試算，本文最終選擇1.25.

圖2　由位錯模型正演的位移場(a) 原始位移場； (b) 加入誤差的位移場； (c) 加入誤差并隨機(jī)提取其中60%的位移場； (d) 實(shí)際GPS測站的位移場.Fig.2　Displacement fields from dislocation model(a) Original displacements； (b) Displacements with stochastic errors added； (c) Randomly extractd 40 percent of displacements, with stochastic errors introduced； (d) Displacements at real GPS sites.

利用正演的位移場，由本文的應(yīng)變計算方法解算應(yīng)變場，結(jié)果如圖4所示，小波尺度因子選擇3～9.與圖3對比可以看出，該結(jié)果與理論應(yīng)變在應(yīng)變梯度比較大的區(qū)域幾乎一致，在離斷層較遠(yuǎn)的區(qū)域存在微小差異，因此，可以認(rèn)為二者在空間分布特征和量值上都非常接近.

3.3數(shù)據(jù)誤差對解算應(yīng)變的影響

考慮到GPS測站速度是含有誤差的，因此在正演的位移場中加入均值為0.0 mm、標(biāo)準(zhǔn)差為0.5 mm的高斯隨機(jī)白噪聲，以研究本文應(yīng)變計算方法對于誤差的敏感程度.加入誤差的位移場如圖2b所示，從該圖可以看出，相對于原始位移場，加入誤差之后的位移場一致性稍差，但分布形態(tài)并未改變.以加入誤差的位移場作為輸入數(shù)據(jù)，由本文應(yīng)變計算方法計算應(yīng)變場，結(jié)果如圖5所示.對比圖3和圖5可以看出，二者的主應(yīng)變/面應(yīng)變和最大剪應(yīng)變分布特征均十分相似，量值也基本一致，可以認(rèn)為二者非常接近.說明本文應(yīng)變計算方法對于高斯隨機(jī)分布的偶然誤差不敏感，方法具有較好的穩(wěn)健性.

圖3　由位錯模型正演的應(yīng)變場(單位：10-8)(a) 主應(yīng)變和面應(yīng)變； (b) 最大剪應(yīng)變.Fig.3　Displacement and strain fields from dislocation model (unit: 10-8)(a) Principal strain and dilatation strain fields； (b) Maximum shear strain fields.

圖4　利用多尺度球面小波方法解算的應(yīng)變場(原始位移場數(shù)據(jù)，單位：10-8)(a) 主應(yīng)變和面應(yīng)變； (b)最大剪應(yīng)變.Fig.4　Strain distribution using multi-scale spherical wavelet (pure displacement fields from dislocation model, unit: 10-8)(a) Principal and dilatation strain； (b) Maximum shear strain.

圖5　基于多尺度球面小波方法解算的應(yīng)變場(位移場加入誤差,單位：10-8)(a) 主應(yīng)變和面應(yīng)變； (b)最大剪應(yīng)變.Fig.5　Strain and displacement distribution using multi-scale spherical wavelet (displacement fields with stochastic errors added, unit: 10-8)(a) Principal and dilatation strains； (b) Maximum shear strain.

3.4測站空間分布對解算應(yīng)變的影響

GPS觀測網(wǎng)的測站分布通常是不規(guī)則的，因此有必要分析測站空間分布的不均勻性對本文應(yīng)變計算方法有效性的影響.從圖2b位移場中隨機(jī)剔除其中40%的數(shù)據(jù)，剔除后的位移場如圖2c所示.可以看出剔除之后位移場并未出現(xiàn)區(qū)域性的空缺.由圖2c所示的位移場，利用本文應(yīng)變計算方法計算的應(yīng)變場，結(jié)果如圖6所示.對比圖3和圖6可以看出，二者的主應(yīng)變/面應(yīng)變和最大剪應(yīng)變分布形態(tài)仍然十分相似，雖然圖6中應(yīng)變的量值僅為圖3的85%，但仍然可以認(rèn)為二者非常接近.說明利用60%的位移場數(shù)據(jù)能夠恢復(fù)真實(shí)的應(yīng)變場，表明本文應(yīng)變計算方法在數(shù)據(jù)比較稀疏時仍然能夠獲得好的應(yīng)變解算結(jié)果.

為進(jìn)一步分析數(shù)據(jù)空間分布對應(yīng)變計算結(jié)果的影響，在該條斷層周邊區(qū)域存在GPS測站的位置正演了位移場，并加入均值為0.0 mm、標(biāo)準(zhǔn)差為0.5 mm的高斯隨機(jī)白噪聲，該位移場如圖2d所示，由此圖可以看出，在斷層北側(cè)測站比較稀疏，平均間距約為120 km；斷層?xùn)|東南側(cè)存在一個400 km×400 km的空區(qū)，其他區(qū)域分布較為均勻.由圖2d所示的位移場，利用本文應(yīng)變計算方法計算的應(yīng)變場，結(jié)果如圖7所示.對比圖7與圖3可以看出，二者的應(yīng)變分

圖6　基于多尺度球面小波方法解算的應(yīng)變場(隨機(jī)提取60%的位移場并加入誤差,單位：10-8)(a) 主應(yīng)變和面應(yīng)變； (b)最大剪應(yīng)變.Fig.6　Displacement distribution by multi-scale spherical wavelet approach (displacement fields added by stochastic errors and deleted to 60%, randomly, unit: 10-8)(a) Principal and dilatation strains； (b) Maximum shear strain.

圖7　多尺度球面小波方法解算的應(yīng)變(真實(shí)測站位移場并加入誤差,單位：10-8)(a) 主應(yīng)變和面應(yīng)變； (b) 最大剪應(yīng)變.Fig.7　Strain and displacement distribution using multi-scale spherical wavelet (displacements with stochastic errors introduced in real GPS sites, unit: 10-8)(a) Principal and dilatation strains； (b) Maximum shear strain.

布特征有相似性，基本反映了斷層的活動特征，但在量值上圖7中的應(yīng)變僅為圖3的63%，說明二者存在較大差異.該結(jié)果似乎與圖7的結(jié)論存在矛盾，但考慮到應(yīng)變本質(zhì)反映的是形變梯度，因此當(dāng)測站間距跨越了變形區(qū)域，形變就可能會被忽略，其計算的應(yīng)變也會與真實(shí)結(jié)果存在較大差異.

4華北地區(qū)應(yīng)變率場分布特征

將多尺度球面小波應(yīng)變計算方法應(yīng)用于華北地區(qū)(其范圍如圖8所示)實(shí)測GPS站速度數(shù)據(jù)分析.測站由兩部分組成: (1) 跨張-渤斷裂和郯-廬斷裂布設(shè)的6個GPS連續(xù)觀測站及2條區(qū)域站剖面，其中連續(xù)站如表2所示，區(qū)域站每年觀測兩期，在2008-2012年期間進(jìn)行了9期觀測； (2) 中國構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(簡稱“陸態(tài)網(wǎng)”)GPS基準(zhǔn)站和復(fù)測區(qū)域站.以上共計28個連續(xù)站，308個復(fù)測區(qū)域站.考慮到華北地區(qū)受到2011年日本東北大地震同震(王敏等，2011；Cheng et al., 2014)，以及震后余滑和松弛的持續(xù)性影響，為保證速度場的精度，在分析時未選擇日本東北大地震之后的觀測資料，同時剔除觀測周期不足兩年的測站，最終得到220個站的速度，扣除區(qū)域整體運(yùn)動后，速度場結(jié)果如圖8所示，從該圖可以看出，華北地區(qū)的運(yùn)動特征在不同區(qū)域存在差異性.在張-渤斷裂帶以北，GPS速度場一致性比較好，在其以南比較雜亂.

圖8　華北地區(qū)GPS速度場(無整體旋轉(zhuǎn)基準(zhǔn))Fig.8　GPS velocity fields in North China with respect to No-Net-Rotation reference frame

圖9　基于多尺度球面小波方法解算的不同尺度GPS面應(yīng)變率場(單位：10-8/a)Fig.9　GPS strain rate fields in different scale calculated by multi-scale spherical wavelet (units：10-8/a)

測站經(jīng)度(°E)緯度(°N)開始時間結(jié)束時間CHCH115.8440.912008-2342013-122SHUN116.5840.282008-2012013-122JIUL115.0241.542008-2402013-122ZUOL115.1440.162009-1722013-122FEHU119.1039.762010-290持續(xù)觀測TUOJ120.7538.162011-015持續(xù)觀測

圖9給出了不同空間尺度的面應(yīng)變率，約定壓縮為正，膨脹為負(fù).從圖9a至9d空間解析度依次增加，對應(yīng)的尺度因子分別為6～9，其中圖9a要求最小站點(diǎn)距離小于110.2 km，最大面應(yīng)變率為1.7×10-8/a；圖9d要求最小站點(diǎn)距離小于12.8 km，最大面應(yīng)變率為5.3×10-8/a.可以看出不同空間尺度應(yīng)變率在一些區(qū)域存在差異，小的分解尺度得到的應(yīng)變率場反映出大尺度應(yīng)變積累特征，大的分解尺度得到的應(yīng)變率場能夠反映出小區(qū)域的應(yīng)變積累特征，說明本文應(yīng)變計算方法具有依據(jù)站點(diǎn)密度計算不同空間尺度應(yīng)變的優(yōu)點(diǎn).

圖10給出了最大分解尺度為8時的主應(yīng)變/面應(yīng)變率和最大剪應(yīng)變率，從該圖可以看出，主應(yīng)變率顯示該區(qū)域的形變場主要受控于近東西向的擠壓；面應(yīng)變率呈現(xiàn)壓縮和拉張交替出現(xiàn)的特征，其中廊坊—天津一帶呈面壓縮，張-渤斷裂帶面壓縮和面膨脹交替出現(xiàn).最大剪應(yīng)變率的分布顯示，該區(qū)域存在兩個最大剪應(yīng)變率較大的區(qū)域，分別分布于沿著張-渤斷裂帶，保定與滄州之間.該結(jié)果與Wu等(2006)和楊國華等(2013)的應(yīng)變率場結(jié)果在大空間尺度上比較一致，但在局部區(qū)域存在差異.可以看出，利用本文應(yīng)變計算方法可以獲取小空間尺度的應(yīng)變積累特征.

圖10　華北地區(qū)GPS應(yīng)變率場(單位：10-8/a)(a) 主應(yīng)變率和面應(yīng)變率； (b) 最大剪應(yīng)變率.Fig.10　GPS strain rate fields in North China (units：10-8/a)(a) Principal and dilatation strain rates； (b) Maximum shear strain rate.

圖11　華北地區(qū)GPS面應(yīng)變率場中誤差Fig.11　Standard deviation of GPS strain rate field in North China

圖11為多尺度球面小波方法計算的面應(yīng)變率中誤差，從該圖可以看出，面應(yīng)變率的中誤差均小于1×10-8/a，如果取2倍中誤差作為面應(yīng)變率可以分辨的閾值，那么圖10中幾個典型區(qū)域應(yīng)變結(jié)果均是可信的.

5結(jié)論

如何由地表大地測量得到的位移、基線變化等觀測量獲取地殼應(yīng)變一直是地學(xué)研究中的熱點(diǎn)問題，本文將多尺度球面小波方法應(yīng)用于GPS應(yīng)變場解算，通過對模擬數(shù)據(jù)的分析，驗(yàn)證了該方法的有效性.多尺度球面小波應(yīng)變計算方法穩(wěn)健性較好，對數(shù)據(jù)中加入的高斯隨機(jī)誤差反映不敏感，隨機(jī)誤差對計算結(jié)果的影響基本可以忽略.從規(guī)則分布的位移場數(shù)據(jù)隨機(jī)抽取其中的60%解算應(yīng)變場，結(jié)果顯示，即使測站分布比較稀疏，該方法仍然能夠獲取比較可靠的結(jié)果；但如果位移場數(shù)據(jù)出現(xiàn)大區(qū)域的空缺，則會對應(yīng)變結(jié)果產(chǎn)生顯著影響.

基于華北地區(qū)新布設(shè)的跨斷層GPS連續(xù)觀測站和區(qū)域剖面站，以及“陸態(tài)網(wǎng)”GPS基準(zhǔn)站和復(fù)測區(qū)域站數(shù)據(jù)，解算了該區(qū)域的GPS速度場，并利用本文的應(yīng)變計算方法計算了華北地區(qū)的應(yīng)變率場及其誤差.結(jié)果顯示，廊坊—天津一帶面應(yīng)變呈壓縮特征，張-渤斷裂帶呈壓縮和膨脹交替特征，該區(qū)域存在兩個最大剪應(yīng)變率較大的區(qū)域，分別分布于張-渤斷裂帶和保定與滄州之間.此外，不同空間尺度的應(yīng)變率存在差異性，利用本文提出的應(yīng)變計算方法可以獲得小空間尺度的應(yīng)變積累特征.

致謝中國構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)國家數(shù)據(jù)中心提供了部分原始GPS觀測數(shù)據(jù).部分圖形采用GMT軟件繪制.伍吉倉教授和Carl Tape博士對本研究提出了許多寶貴意見，在此一并表示感謝.

References

Bao B C, Hu W, Liu Z, et al. 2009. Dynamical analysis of DOG wavelet mapping with dilation and translation.ActaPhysicaSinica(in Chinese), 58(4): 2240-2247.

Baxter S C, Kedar S, Parker J W, et al. 2011. Limitations of strain estimation techniques from discrete deformation observations.Geophys.Res.Lett., 38(1): L01305, doi: 10.1029/2010GL046028. Bogdanova I, Vandergheynst P, Antoine J P, et al. 2005. Stereographic wavelet frames on the sphere.Appl.Comput.HarmonicAnal., 19(2): 223-252. Chambodut A, Panet I, Mandea M, et al. 2005. Wavelet frames: an alternative to spherical harmonic representation of potential fields.Geophys.J.Int., 163(3): 875-899. Cheng J, Liu M, Gan W J, et al. 2014. Seismic impact of theMW9.0 Tohoku earthquake in Eastern China.Bull.Seismol.Soc.Am., 104(3): 1258-1267, doi: 10.1785/0120130274.

Craig T J, Calais E. 2014. Strain accumulation in the New Madrid and Wabash Valley Seismic Zones from 14 years of continuous GPS observation.J.Geophys.Res., 119(12): 9110-9129, doi: 10.1002/2014JB011498.

Gan W J, Zhang P Z, Shen Z K, et al. 2007. Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements.J.Geophys.Res., 112(B8): B08416, doi: 10.1029/2005JB004120.Ge W P, Wang M, Shen Z K, et al. 2013. Intersiesmic kinematics and defromation patterns on the upper crust of Qaidam-Qilianshan block.ChineseJ.Geophys. (in Chinese), 56(9): 2994-3010, doi: 10.6038/cjg20130913.

Holschneider M, Chambodut A, Mandea M. 2003. From global to regional analysis of the magnetic field on the sphere using wavelet frames.Phys.EarthPlanet.Inter., 135(2-3): 107-124.

Jiang Z S, Liu J N. 2010. The method in establishing strain field and velocity field of crustal movement using least squares collocation.ChineseJ.Geophys. (in Chinese), 53(5): 1109-1117, doi: 10.3969/j.issn.0001-5733.2010.05.011.

Liang S M, Gan W J, Shen C Z, et al. 2013. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements.J.Geophys.Res., 118(10): 5722-5732, doi: 10.1002/2013JB010503.

Liu X X, Jiang Z S, Wu Y Q. 2014. The applicability of Kriging interpolation method in GPS velocity gridding and strain calculating.GeomaticsandInformationScienceofWuhanUniversity(in Chinese), 39(4): 457-461.Ma T, Li F, Yue J L. 2003. Application of wavelet analysis in geophysics and geodesy.ProgressinGeophysics(in Chinese), 18(1): 49-52.

Meng G J, Ren J W, Wang M, et al. 2008. Crustal deformation in western Sichuan region and implications for 12 May 2008MS8.0 earthquake.Geochemistry,Geophysics,Geosystems, 9(11): Q11007, doi: 10.1029/2008GC002144. Meng G J, Ren J W, Wu J C, et al.2008. Computation of strain and rotation tensor and their uncertainty for small arrays in spherical coordinate system.ActaSeismologicaSinica，30(1):67-75.

Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space.Bull.Seismol.Soc.Am., 75(4): 1135-1154.

Pearson C F, Snay R A. 2014. Strain partitioning along the western margin of North America.JournalofStructuralGeology, 64: 67-78.

Poyraz F. 2015. Determining the strain upon the eastern section of the North Anatolian fault zone (NAFZ).Arab.J.Geosci., 8(3): 1787-1799, doi: 10.1007/s12517-014-1301-6. Savage J C, Gan W J, Svarc J L. 2001. Strain accumulation and rotation in the Eastern California Shear Zone.J.Geophys.Res., 106(B10): 21995-22007. Shen Z K, Jackson D D, Ge B Z. 1996. Crustal deformation across and beyond the Los Angeles basin from geodetic measurements.J.Geophys.Res., 101(B12): 27957-27980, doi: 10.1029/96JB02544.Tape C, Musé P, Simons M, et al. 2009. Multiscale estimation of GPS velocity fields.Geophys.J.Int., 179(2): 945-971. Wang M, Li Q, Wang F, et al. 2011. Far-field coseismic displacements associated with the 2011 Tohoku-oki earthquake in Japan observed by Global Positioning System.ChineseSci.Bull., 56(23): 2419-2424, doi: 10.1007/s11434-011-4588-7.

Wang Q, Qiao X J, Lan Q G, et al. 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan.NatureGeoscience, 4(9): 634-640, doi: 10.1038/ngeo1210.

Wang Z, Dahlen F A. 1995. Spherical-spline parameterization of three-dimensional earth models.Geophys.Res.Lett., 22(22): 3099-3102.

Wu J C, Tang H W, Chen Y Q, et al. 2006. The current strain distribution in the North China Basin of eastern China by least-squares collocation.JournalofGeodynamics, 41(5): 462-470, doi: 10.1016/j.jog.2006.01.003.

Wu Y Q, Jiang Z S, Yang G H, et al. 2009. The application and method of GPS strain calculation in whole mode using least square collocation in sphere surface.ChineseJ.Geophys. (in Chinese), 52(7): 1707-1714, doi: 10.3969/j.issn.0001-5733.2009.07.005.

Wu Y Q, Jiang Z S, Yang G H, et al. 2011. Comparison of GPS strain rate computing methods and their reliability.J.Geophys.Int., 185(2): 703-717.Yang G H, Yang B, Chen X, et al. 2013. The basic characteristics of spatial variation of 3D deformation field in current North

China.GeomaticsandInformationScienceofWuhanUniversity(in Chinese), 38(1): 31-35.

附中文參考文獻(xiàn)

包伯成, 胡文, 劉中等. 2009. DOG小波映射伸縮和平移的動力學(xué)分析. 物理學(xué)報, 58(4): 2240-2247.

葛偉鵬, 王敏, 沈正康等. 2013. 柴達(dá)木—祁連山地塊內(nèi)部震間上地殼塊體運(yùn)動特征與變形模式研究. 地球物理學(xué)報, 56(9): 2994-3010, doi: 10.6038/cjg20130913.

江在森, 劉經(jīng)南. 2010. 應(yīng)用最小二乘配置建立地殼運(yùn)動速度場與應(yīng)變場的方法. 地球物理學(xué)報, 53(5): 1109-1117, doi: 10.3969/j.issn.0001-5733.2010.05.011.

劉曉霞, 江在森, 武艷強(qiáng). 2014. Kriging方法在GPS速度場網(wǎng)格化和應(yīng)變率場計算中的適用性. 武漢大學(xué)學(xué)報·信息科學(xué)版, 39(4): 457-461.

馬濤, 李斐, 岳建利. 2003. 小波分析在地球物理及大地測量中的應(yīng)用. 地球物理學(xué)進(jìn)展, 18(1): 49-52.

孟國杰，任金衛(wèi)，伍吉倉等. 2008. 球坐標(biāo)系中圖形單元應(yīng)變與旋轉(zhuǎn)張及其誤差解算. 地震學(xué)報，30(1):67-75.

王敏, 李強(qiáng), 王凡等. 2011. 全球定位系統(tǒng)測定的2011年日本宮城MW9.0級地震遠(yuǎn)場同震位移. 科學(xué)通報, 56(20): 1593-1596.

武艷強(qiáng), 江在森, 楊國華等. 2009. 利用最小二乘配置在球面上整體解算GPS應(yīng)變場的方法及應(yīng)用. 地球物理學(xué)報, 52(7): 1707-1714, doi: 10.3969/j.issn.0001-5733.2009.07.005.

楊國華, 楊博, 陳欣等. 2013. 華北現(xiàn)今三維形變場空間變化的基本特征. 武漢大學(xué)學(xué)報, 38(1): 31-35.

(本文編輯胡素芳)

Methodology and application of GPS strain field estimation based on multi-scale spherical wavelet

SU Xiao-Ning1,2， MENG Guo-Jie1*， WANG Zhen3

1InstituteofEarthquakeScience,ChinaEarthquakeAdministration,Beijing100036,China2CollegeofSurveyingandGeo-informatics,TongjiUniversity,Shanghai200092,China3EarthquakeResearchInstitute,UniversityofTokyo，TokyoJapan

AbstractA methodology was proposed to estimate spherical GPS strain based on multi-scale spherical wavelet. By constructing an irregular distribution of wavelet basis functions to highlight non-uniform distribution of stations, this methodology can be used to calculate the strain fields at different spatial scales. The methodology is verified and evaluated using synthetic displacements. The simulated results show that the proposed methodology is robust. If Gaussian errors, with a mean value of 0.0 mm and a standard deviation of 0.5 mm are introduced into displacements, its effect to derived strain field is negligible. If 60% data are extracted from the original displacement field to calculate its strains, the obtained strain field is still reliable. Nonetheless, if large scale of data is missed, it will bring noticeable effect to the derived strain field. Based on the GPS raw data in North China, we obtain a displacement field utilizing a precise data processing strategy. Subsequently, we calculate strain rate field and its associated uncertainty, which spatial characteristics are finally analyzed.KeywordsSpherical wavelet； Multi-scale； GPS strain fields； Cross-fault GPS profile； Simulation analysis

基金項(xiàng)目國家自然科學(xué)基金(41174004，41461164004)，國家國際科技合作專項(xiàng)項(xiàng)目(2010DFB20190，2015DFR21100)，中國地震局地震預(yù)測研究所基本科研業(yè)務(wù)費(fèi)專項(xiàng)重點(diǎn)項(xiàng)目(2012IES010302)聯(lián)合資助.

作者簡介蘇小寧，男，1986年生，助理研究員，從事GPS數(shù)據(jù)處理與地殼形變、斷層活動反演等研究.E-mail:tongjisxn@163.com *通訊作者孟國杰，男，1969年生，研究員，長期從事地殼形變與地震位錯反演研究.E-mail:guojiemeng@126.com

doi:10.6038/cjg20160504 中圖分類號P228

收稿日期2015-04-17，2016-01-22收修定稿

蘇小寧， 孟國杰, 王振. 2016. 基于多尺度球面小波解算GPS應(yīng)變場的方法及應(yīng)用.地球物理學(xué)報，59(5):1585-1595,doi:10.6038/cjg20160504.

Su X N, Meng G J, Wang Z. 2016. Methodology and application of GPS strain field estimation based on multi-scale spherical wavelet.ChineseJ.Geophys. (in Chinese),59(5):1585-1595,doi:10.6038/cjg20160504.