應(yīng)變計算與分析的若干問題及有關(guān)偏差的修正＊

楊國華武艷強楊 博張風(fēng)霜韓月萍周海濤董運鴻

1)中國地震局第一監(jiān)測中心,天津 300180

2)中國地震局地震預(yù)測研究所,北京 100036

應(yīng)變計算與分析的若干問題及有關(guān)偏差的修正＊

楊國華1)武艷強2)楊 博1)張風(fēng)霜1)韓月萍1)周海濤1)董運鴻1)

1)中國地震局第一監(jiān)測中心,天津 300180

2)中國地震局地震預(yù)測研究所,北京 100036

對多種平面與球面應(yīng)變場計算模型進行分析,指出何種模型屬于無偏差應(yīng)變計算模型,何種模型為有偏差應(yīng)變計算模型,何種模型偏差可以修正,何種模型只適用于均勻介質(zhì)的應(yīng)變計算而不能外延,何種模型可以用于作連續(xù)應(yīng)變場的求解與分析;給出利用應(yīng)變結(jié)果計算球面差異運動 (位移)的修正關(guān)系,并論證了球面應(yīng)變模型中非微分項的含義。

GNSS;水平應(yīng)變計算;模型適應(yīng)性;偏差修正;應(yīng)變?yōu)V波

1 引言

在應(yīng)變場描述的早期研究中,由于研究的空間尺度較小,通常是將局部球面運動投影到高斯平面,然后在平面上利用運動 (位移)與應(yīng)變的關(guān)系求解研究區(qū)的應(yīng)變場。隨著 GNSS的引入,研究區(qū)的空間尺度得以擴大,在平面上的研究已顯得不適宜,因此球面應(yīng)變計算模型相繼出現(xiàn)[1-11]。然而,這些模型在理論上并非一致,適用范圍也不盡相同,有的存在偏差,這在一定程度上影響了計算結(jié)果。此外,由于對應(yīng)變的廣譜性沒有得到充分的肯定與了解,也使得應(yīng)變分析變得混亂。為此,本文試圖澄清這一問題,使得應(yīng)變計算與分析準確、有效。

2 水平應(yīng)變場計算模型或計算式

2.1 平面水平應(yīng)變場計算模型或計算式

在平面直角坐標系下,進行地形變分析所用的應(yīng)變模型或計算式可歸為 4類。一類是在已知平面運動解析式的情況下通過微分方程求解應(yīng)變,第二類是根據(jù)離散點的平面運動通過建立應(yīng)變與運動的關(guān)系方程求解應(yīng)變,第三類是利用離散邊長的相對變化量求解兩個正交方向上的正應(yīng)變和它們之間的剪應(yīng)變,第四類是基于Delaunay的三角單元法計算應(yīng)變。

若已知東、北向運動的解析式 fx(x,y)、fy(x, y),求解應(yīng)變的計算式為:

若已知離散點東、北向的運動 (位移)vx(xj,yj)、vy(xj,yj),根據(jù)運動與應(yīng)變關(guān)系方程可求解應(yīng)變:

若已知離散邊長相對變化量εj和方位角φj,根據(jù)耶格法建立方程可求解應(yīng)變:

若已知Delaunay三角單元的三個節(jié)點坐標 (xi, yi)和相應(yīng)的水平運動 vx(xi,yi)、vy(xi,yi),通過建立運動與位置的平面方程[9]:

然后利用式(1)獲得三角形單元的平均應(yīng)變,即

2.2 球面水平應(yīng)變場計算模型或計算式

已知東、北向運動的解析式 fe(λ,φ)、fn(λ,φ),求解應(yīng)變的計算式為:

式中R為地球的平均半徑。

已知離散點東、北向的運動 ve(λj,φj)、vn(λj, φj),根據(jù)運動與應(yīng)變關(guān)系方程可求解應(yīng)變:

已知離散大地線長相對變化量εj和方位角φj,根據(jù)耶格法建立方程可求解應(yīng)變:

已知Delaunay三角單元的三個節(jié)點坐標 (λj, φj)和相應(yīng)的水平運動 ve(λj,φj)、vn(λj,φj),同樣仿照平面模式獲得三角形單元的平均應(yīng)變。

2.3 不同情況下計算模型的差異與適用性

2.3.1 解析法

平面模型式(1)從理論上講僅是忽略了應(yīng)變計算的高次項,由于高次項的數(shù)值非常小,在通常應(yīng)用中將其忽略不會產(chǎn)生不可接受的問題,所以可以認為是一個嚴密的算式。它可適用于平面非均勻介質(zhì)的連續(xù)應(yīng)變計算。為了與平面模型對比,現(xiàn)將球面模型式(4)做等價變換,為:

式中,sλ和 sφ分別為經(jīng)、緯向大地線長度。從形式上則不難發(fā)現(xiàn),東向應(yīng)變多了算項,剪切應(yīng)變也多了算項,而且是非微分項。由此使人產(chǎn)生誤解,認為球面應(yīng)變計算與運動基準有關(guān);由于運動基準的相對性,感覺難以適從,此外這還與通常的應(yīng)變分析中與基準無關(guān)的認識相悖,所以蓋模型沒有得到廣泛應(yīng)用。實際上這是由于坐標系統(tǒng)的不一致所出現(xiàn)的現(xiàn)象。

為了簡要地說明這一問題,假定坐標的“北”向(“經(jīng)”線)間是平行的,類似于卯酉圈,這時質(zhì)元可向“北”作無形變運動;然而,用現(xiàn)行球面坐標的定義卻無法給出直接的描述。換句話說,現(xiàn)行的北向運動僅是上述無形變“北”向運動的投影,因此質(zhì)元向“北”向運動時,現(xiàn)行北向運動的投影描述則必然顯示經(jīng)向間的距離在縮短,即出現(xiàn)了偏差。若假定質(zhì)元的“北”向運動為 fN(λ,φ),當Δ λ→0時,則 fn(λ,φ)→fN(λ,φ)、Δ λ→sinΔ λ,那么由于北向投影所造成經(jīng)度方向運動的變化量為:

相應(yīng)地東向應(yīng)變偏差則為:

與式(9)中東向應(yīng)變的非微分項大小相等,方向相反。對于北向應(yīng)變不論運動質(zhì)元向何方向運動,其應(yīng)變量均不存在偏差,但剪切應(yīng)變卻不然,這是因為在球面上無形變質(zhì)元運動,其角速度不變,線速度在平行圈間(如緯度)是不同的。然而,實際計算剪應(yīng)變是用線速度,所以必然出現(xiàn)偏差,這種偏差只存在于東向運動,其大小為:

相應(yīng)地剪切應(yīng)變偏差則為:

剪切應(yīng)變偏差大小與剪切應(yīng)變算式中的非微分項相等,符號相反。

分析表明,在球面上計算應(yīng)變,直接套用直角坐標系下的應(yīng)變算式是有偏的。式 (4)則是應(yīng)變計算的無偏算式,它適合于球面非均勻介質(zhì)連續(xù)應(yīng)變場的計算與分析。

然而,盡管我們已經(jīng)明了非微分項源于對偏差的修正,但運動必然涉及到基準問題。由上述推導(dǎo)可知,這里的運動在球面上帶有某種絕對性,故目前無法給出嚴格精確的結(jié)果?，F(xiàn)假定地球表面在水平運動時所產(chǎn)生地球表面積變化之和為零,那么在全球范圍內(nèi)的運動若滿足

則可認為該運動在地球表面上具有絕對性質(zhì)的運動。式中D是整個地球表面的空間范圍;dm是質(zhì)量元;r、V是 dm的位置矢量和速度矢量。ITRF參考框架早期是以這一準則而獲得的,ITRF2000、ITRF2005雖不嚴格符合這一要求,但偏差很小。這就是說,當前應(yīng)變修正量計算所用的運動結(jié)果依ITRF參考框架所給出的運動結(jié)果方可。因此,式(4)的恰當使用可獲得球面連續(xù)無偏應(yīng)變場。

2.3.2 離散運動法

當無法獲得運動場的解析式但研究區(qū)內(nèi)具有離散測點及其相應(yīng)的運動時,通常通過運動與應(yīng)變的關(guān)系方程求解應(yīng)變。在平面上利用式 (2),在球面上利用式(5)。式(2)為嚴密關(guān)系方程,由此所獲得的應(yīng)變是無偏差的,但它的適用范圍只是在均勻介質(zhì)條件下。如果將求解 3個應(yīng)變參量,通過一定的方法轉(zhuǎn)變?yōu)榍蠼?3個應(yīng)變參量解析式,則可適用于連續(xù)應(yīng)變場的求解,但目前尚未獲得較為理想的結(jié)果,能獲取的主要是低頻信息或中低頻以內(nèi)的信息。就式(5)而言,實質(zhì)上與式 (2)一致,整體運動 3個參量、應(yīng)變 3個參量,都適用相對運動量計算應(yīng)變參量,所不同的是坐標系統(tǒng)。由上述分析可知,由于坐標系統(tǒng)非直角坐標系統(tǒng),故后者必然存在偏差,若偏差得以修正,則它的適用范圍與式(2)同。

2.3.3 線應(yīng)變法

式(3)與式 (6)實際是通過不同方向的應(yīng)變量求解兩個正交方向上的正應(yīng)變和它們之間的剪應(yīng)變。不論在平面還是在球面上,均不存在偏差。但如果用三角形邊長的應(yīng)變量求解兩個正交方向上的正應(yīng)變和剪應(yīng)變,并通過空間上各個三角形所求的正應(yīng)變和剪應(yīng)變進行連續(xù)應(yīng)變場模擬分析是不適宜的,這時則存在不可修正的偏差。為了解釋這一問題現(xiàn)舉例說明。假定連續(xù)形變塊體A內(nèi)有5個監(jiān)測點,其差異運動結(jié)果如圖 1所示。很顯明,這是一個存在變形的塊體。若運動恰好分別符合以兩個三角形重心為圓心、半徑相同時所給出的旋轉(zhuǎn)運動結(jié)果,那么依式 (3)、(6)計算應(yīng)變則必然為零。由此得出了無形變的結(jié)論,這顯然是錯誤的。所以,這種方法只能用于均勻介質(zhì)形變塊體的應(yīng)變求解,不適用非均勻介質(zhì)塊體的連續(xù)應(yīng)變求解。

圖 1 塊體A的相對運動示意圖Fig.1 Sket of the relative motion of blockA

2.3.4 Delaunay三角單元法

Delaunay三角單元法是當前使用較多的一種方法,并“以元代面”的方式使應(yīng)變在空間上連續(xù)起來。在平面上該方法計算結(jié)果其一與三角形的組成有關(guān)(即多解性);其二相鄰單元的應(yīng)變計算彼此間互不顧及;其三,由于三角形的頂點分布不均勻 (即非等邊),所以計算所得到的應(yīng)變結(jié)果并不在同一頻譜域內(nèi),這是該方法在使用中應(yīng)注意或存在的最大問題。因為不同頻譜域的應(yīng)變值相差可以大到以數(shù)量級計,所以在進行空間對比分析時,很容易得出不正確的結(jié)論。此外,如果使用式 (2)、(5)進行計算,由圖 1可知,不僅上述問題依舊,就連續(xù)應(yīng)變場分析來說還存在失真問題。

3 由球面應(yīng)變結(jié)果計算位移 (運動)的關(guān)系式

在某些情況下或分析研究的需求,需要利用連續(xù)應(yīng)變(率)參量求解連續(xù)位移參量。在平面坐標系統(tǒng)下根據(jù)式(2)有

由式(15)可獲取研究區(qū)的相對位移。但在球面上則不能直接用無偏連續(xù)應(yīng)變參數(shù)計算,故還需由式(16)才能獲得球面坐標系下相對位移的正確結(jié)果。

4 有偏模型獲得水平應(yīng)變場的修正

應(yīng)變計算中,由于大部分學(xué)者更多地是使用式(5)或式 (2)。通過分析可知,盡管有偏,但可修正。對于較高頻的應(yīng)變 (量級在 10-8及以上)根據(jù)實際情況可以考慮忽略這一修正,但對低頻應(yīng)應(yīng)變 (量級在 10-9)則不能忽略,因為它們基本上處在同一個量級上。修正算式如下:

該式所用的運動結(jié)果是 ITRF參考框架下的運動結(jié)果。筆者以 1999—2007年 GNSS平均運動速率為基礎(chǔ)進行計算,結(jié)果表明:中國大陸東向應(yīng)變率修正量大小分布的變化范圍為 -2.5～3.0,在空間的變化過程為西大東小,西部地區(qū)的變化梯度較大,東部較小;東、北向剪切應(yīng)變修正量的變化范圍為0.9～2.8,從變化趨勢看為北大南小。此外,以上述3個應(yīng)變量為基礎(chǔ)進一步計算,如最大主應(yīng)變、最小主應(yīng)變、最大剪切應(yīng)變和面應(yīng)變等都存在相應(yīng)的變化。由此可知,應(yīng)變修正量在低頻應(yīng)變及弱形變地區(qū)(如中國東部地區(qū))分析中是不可缺少的。

5 應(yīng)變分析中的廣譜性問題

當前在應(yīng)變分析中的廣譜性問題多被忽略或沒有得到應(yīng)有的重視,造成了利用同一運動結(jié)果而得到的應(yīng)變結(jié)果差異甚大,讓人們有無所適從之感。這除了存在模型選用不恰當或模型不完善的問題外,還與數(shù)據(jù)采用點的空間分布與取用有關(guān),尤其觀測值存在明顯的誤差干擾則更是如此。我們不能否認的一個問題是,依據(jù)當前的測站密度與觀測誤差,如果不進行進一步處理,其應(yīng)變信息與誤差干擾的量級基本相當,尤其在中國東部地區(qū)。我們不應(yīng)該回避這一問題,同時我們也必須清楚地殼介質(zhì)是非各向同性的,應(yīng)變具有廣譜性是一個基本的事實。這就是說,即使采用無偏差應(yīng)變計算模型,而接下來的進一步處理做與不做,如何做,最后所得到的結(jié)果仍然不一定相同。圖 2是采用式(4)計算川滇地區(qū)1999—2007年東、北向剪切應(yīng)變率并在不同的后處理情況下所得到的結(jié)果[12]。由此我們可以看到,計算結(jié)果的不同可能存在兩個方面的原因:其一是模型使用不恰當,其結(jié)果具有失真性;其二是模型使用恰當,但由于信息分離的不同而造成了差異,即應(yīng)變所包含的譜域范圍不同,它不屬于錯誤。圖 2(a)結(jié)果所包含的應(yīng)變譜域相對較寬,高頻應(yīng)變成分較為豐富,圖 2(b)結(jié)果所包含的應(yīng)變譜域相對較窄,主要側(cè)重于較低頻應(yīng)變信息的描述。

6 結(jié)語

在球面連續(xù)應(yīng)變分析中所存在的問題主要有兩方面:模型選用不恰當以及有的模型存在偏差。因此,使球面應(yīng)變分析產(chǎn)生了混亂。筆者舉薦球面無偏差計算模型式(4),該式不但沒有偏差,而且可包含各種頻譜的應(yīng)變信息 (依測站密度所能表征的),故可利用濾波方法分離出所需要的應(yīng)變信息[12];以三角形邊長相對變化為依托的耶格法在連續(xù)計算與分析中應(yīng)該拋棄;以式 (5)為依托的連續(xù)應(yīng)變計算是有偏差的,在低頻應(yīng)變分析中是不可忽略的,通過偏差修正則可得到無偏結(jié)果,可謹慎使用,Delaunay三角單元法 (斜面法)通過偏差修正后雖可進行的連續(xù)應(yīng)變分析,但在實際中難以滿足等邊條件,故形變信息比較混雜。此外,如果利用無偏應(yīng)變結(jié)果計算相對運動,則需按要求進行有偏處理。

圖 2 川滇地區(qū)東、北向間剪切應(yīng)變率Fig 2 Results of shear strain rate between E and N

致謝 感謝王敏研究員提供的運動場結(jié)果!

1 李延興,等.中國大陸活動塊體的運動與應(yīng)變狀態(tài)[J].中國科學(xué) (D輯),2003,46(增):65-81.(Li Yangxin,et al. Movement and strain conditions of active blocks in the Chinese mainland[J].Science in China(Series D),2003,46 (Supp):65-81)

2 李延興,等.中國大陸及周邊地區(qū)的水平應(yīng)變場[J].地球物理學(xué)報,2004,47(2):222-231.(Li Yangxing,et al. Horizontal strain field in the Chinese mainland and its surrounding areas[J].Acta Geophysics Sinica,2004,47(2): 222-231)

3 楊國華,江在森,王敏.印尼地震對我國川滇地區(qū)地殼水平活動的影響[J].大地測量與地球動力學(xué),2006,(1):9 -14.(Yang Guohua,Jiang Zaisen and WangMin.Effect of Indonesia earthquake on horizontal crustal movement in Sichuan-Yunnan region[J].Journal of Geodesy and Geodynamics,2006,(1):9-14)

4 江在森,等.GPS初步結(jié)果揭示的中國大陸水平應(yīng)變場與構(gòu)造變形[J].地球物理學(xué)報,2003,46(3):352-358. (Jiang Zaisen,et al.Horizontal strain field and tectonic deformation of China mainland revealed by preliminary GPS result[J].Acta Geophysics Sinica,2003,46(3):352-358)

5 張希,江在森.用最小二乘配置獲得地形變應(yīng)變場動態(tài)圖像的幾個問題研究[J].地殼形變與地震,1999,(3): 32-39.(Zhang Xi and Jiang Zaisen.Study on some questions of dynamic pictures of crustal deformation and strain flelds obtained by the least square collocation[J].Crustal Deformation and Earthquake,1999,(3):32-39)

6 Savage J C,Gan W and Svae J L.Strain accumulation and rotation in the Eastern California Shear Zone[J].J.Geophys.Res.,2001,106(B10):21 995-22 007

7 楊國華,等.由 GPS觀測結(jié)果推導(dǎo)中國大陸現(xiàn)今水平應(yīng)變場[J],地震學(xué)報,2002,24(4):337-347.(Yang Guohua,et al.Current horizontal strain field in Chinese mainland derived from GPS data[J].Acta Seismologica Sinica, 2002,24(4):337-347)

8 武漢測繪學(xué)院大地測量系地震測量教研組.大地形變測量學(xué)(下冊)[M].北京:地震出版社,1980.(Dynamic Geodesy Teaching and Research Group,GeodesyDepartment for Wuhan College of Geodesy and Cartography.Dynamic geodesy[M].Beijing:Seis mological press,1980)

9 翟偉,等.渭河盆地現(xiàn)今地殼水平運動及應(yīng)變特征[J].大地測量與地球動力學(xué),2009,(4):34-37.(Qu Wei,et al.Research on horizontal crustalmotion and strain features ofWeihe basin with GPS data[J].Journal of Geodesy and Geodynamics,2009,(4):34-37)

10 石耀霖,朱守彪.用 GPS位移資料計算應(yīng)變方法的討論[J].大地測量與地球動力學(xué),2006,(1):1-8.(Shi Yaolin and Zhu Shoubiao.Discussion onmethod of calculating strain with GPS data[J].Journal of Geodesy and Geodynamics,2006,(1):1-8)

11 劉序儼,黃聲明,梁全強.旋轉(zhuǎn)橢球面上的應(yīng)變與轉(zhuǎn)動張量表達 [J].地震學(xué)報,2007,29(3):240-249.(Liu Xuyan,Huang Shengming and Liang Quanqiang.Expression of strain and rotation tensor in geodetic coordinates [J].Acta Seismologica Sinica,2007,29(3):240-249)

12 楊博,張風(fēng)霜,韓月萍.球面水平應(yīng)變場無偏差計算的實現(xiàn)與濾波 [J].(待刊)(Yang Bo,Zhang Fengshuang and Han Yueping.No-deviation calculation and filtering of spherical horizontal strain field[J].(to be published))

SEVERAL PROBLEM S OF STRA IN CALCULATI ON AND ANALYSIS AND CORRECTI ON OF RELATED DEVIATI ON

Yang Guohua1),Wu Yanqiang2),Yang Bo1),Zhang Fengshuang1),
Han Yueping1),Zhou Haitao1)and Dong Yunhong1)
1)First CrustM onitoring and Application Center,CEA,Tianjin 300180
2)Institute of Earthquake Science,CEA,B eijing 100036

There existsmany kindsof calculation modelsof plane and spherical strain fields,but the results are different.It is analyzed the representative models,and points outwhich model has unbias aswell aswhich one has deviation,which models deviation can be corrected aswell aswhich one can only be used to the strain computing of unifor m medium but can not be extended and which model can be used the calculation and analysis of continuous strain field.Besides,it’s also provides the correction relationship for computing spherical difference movement (displacement)by strain results,and demonstrates themeaning of the non-differential ter m in spherical strain model.

GNSS;horizontal strain calculation;model applicability;deviation correction;strain filtering

P551

A

1671-5942(2010)04-0059-05

2009-12-12

中國地震局地震預(yù)測研究所基本科研業(yè)務(wù)專項(02092401);“十一五”國家科技支撐課題(2006BAC01B03-01-01);中國地震局“強化華北地區(qū)強震監(jiān)視跟蹤項目”

楊國華,男,1957年生,研究員,從事地殼形變及其動力學(xué)研究.E-mail:ygh5710@136.com