李春曉

1 中國科學院云南天文臺，昆明市官渡區(qū)，650216

MORVEL板塊運動模型估計了全球25 個板塊（占全球面積比為97.2%）的相對運動參數(shù)及其不確定度［1］?；贛ORVEL 模型，Donald將剩余的31 個小板塊（由Bird［2］定義，占全球面積比為2.8%）統(tǒng)一到MORVEL 模型中，從而構(gòu)成NNR－MORVEL56［3］絕對板塊運動模型。NNR－MORVEL56描述了全球56個板塊相對于無旋轉(zhuǎn)（no－net－rotation，NNR）參考架的絕對運動參數(shù)，每個板塊的幾何參數(shù)、運動參數(shù)及其不確定度可參考文獻［1，3］。

ITRF2008［4］是對4種空間大測量技術(shù)（VLBI，SLR，GPS，DORIS）的觀測數(shù)據(jù)進行重新解算獲得的，考慮到它的解較以往國際地球參考架有了很大提高，本文以ITRF2008 參考架下225個GPS測站的速度場作為計算板塊運動參數(shù)（歐拉矢量）的原始數(shù)據(jù)。

1 數(shù)學模型

若將地球近似為球體，那么總可以用歐拉定理將板塊在地球表面的運動描述為整個板塊圍繞一個穿過地心的旋轉(zhuǎn)軸作轉(zhuǎn)動。定義CM（center of mass）為由激光測距解算的地球質(zhì)心（包括地球、海洋、大氣的質(zhì)量）位置，CF（center of figure）為固態(tài)地球的形心位置［4－5］，板塊的運動可用以下方程表述：

其中，ve、vn分別為測站水平速度的東向和北向分量，λ和φ分別為測站的大地測量經(jīng)緯度，ωx、ωy、ωz分別為ω的3 個分量。用以下公式將求得的歐拉矢量ω用大地測量坐標（λ，φ）和模量ω表示為：

其中，φ′為地心緯度，f為地球扁率，f＝1／298.257 22。若要求解式（2）中的3 個未知量，則至少需要兩個測站的速度場。通常情況下一個板塊上會含有多個測站，從而式（2）為超定方程。利用加權(quán)最小二乘法求解超定方程可以得到板塊的歐拉矢量，而歐拉矢量的估計誤差可以通過誤差傳播定律來求得。假設(shè)已知測站在地心坐標系中的三維速度及其不確定度，那么可以利用以下方程得到測站的水平速度及其不確定度：

其中，cov（ven）為測站地平坐標系下水平速度的協(xié)方差矩陣，cov（vxyz）為ITRF2008參考架下測站三維速度的協(xié)方差矩陣，R為?ven／?vxyz，即從地心坐標系到測站地平坐標系的旋轉(zhuǎn)矩陣，RT為轉(zhuǎn)置矩陣。本文假定各個測站之間的水平速度場互不影響［7］，則由加權(quán)最小二乘法得到的擬合參數(shù)為：

其中，Ai為式（2）中的設(shè)計矩陣，表示矩陣的轉(zhuǎn)置，vi為測站i的水平速度，n代表板塊上總的測站數(shù)目，Wi為加權(quán)矩陣：

其中，rei為測站i東向速度殘差，rni為測站i北向速度殘差，n為板塊上測站總數(shù)。擬合參數(shù)的協(xié)方差矩陣可由式（9）來計算：

歐拉矢量在球面坐標下的協(xié)方差矩陣及其誤差橢圓也可以用誤差傳播定律來推導。首先，已知歐拉矢量的直角坐標（ωx，ωy，ωz）與球面坐標（ω，λ，φ′）之間有如下關(guān)系：

然后求出雅可比矩陣及其逆矩陣，也即

最后求出球面坐標下歐拉矢量的協(xié)方差矩陣：

對上面3×3協(xié)方差矩陣cov（ωωλφ′）提取出歐拉矢量經(jīng)度和緯度之間的2×2 協(xié)方差矩陣cov（ωλφ′），然后利用誤差傳播定律推導出歐拉矢量的經(jīng)度和大地測量緯度之間的協(xié)方差矩陣為：

誤差橢圓的半長軸、半短軸以及長軸方向可通過對角化cov（ωλφ）求得，分別對應矩陣的兩個特征值和較大特征值對應的特征向量。

2 測站篩選

測站的篩選主要遵循以下5個準則：1）測站的觀測時間段不少于3a，主要用來避免季節(jié)的周期性對測站速度的影響；2）測站的速度場精度優(yōu)于3mm／a，也即測站的東向和北向的速度不確定度都要小于3 mm／a；3）測站的第二不變應變率小于10－14／a，因為第二不變應變率的大小可以間接衡量板塊內(nèi)測站的形變大小，從而判定測站是否位于板塊的剛性區(qū)域，這樣避免了用測站是否位于板塊邊界地帶作為判定依據(jù)的不可靠性；4）剔除受冰期回彈［8－10］的影響而產(chǎn)生附加速度較大的測站［6］，本文選取附加水平速度小于0.5mm／a作為限值［5］，也可以用附加垂直速度小于2.5mm／a作為限值；5）保留速度殘差小于3mm／a的測站。

2.1 全球應變率模型

由于板塊并非嚴格的剛體，當板塊受到外界驅(qū)動（擠壓、分離、滑動）時，其整體會發(fā)生形變，而且形變在空間上的分布也各不相同。板塊上空間點的形變表現(xiàn)為應變。為了衡量全球應變的大小，相應的全球應變率模型相繼建立。模型建立的基本方法如下［11］：1）獲取全球測站的速度場，主要為GPS速度場，以及中亞地區(qū)第四紀斷層滑動率和全球震源較淺（小于40km）的地震矩張量等數(shù)據(jù)；2）根據(jù)Kostrov理論將斷層滑動率和地震矩張量轉(zhuǎn)換為應變率［12］，并利用雙三次Bessel插值［13］將旋轉(zhuǎn)矢量函數(shù)（實際上該旋轉(zhuǎn)矢量類似于歐拉矢量，不同的是歐拉矢量應用在整個板塊上，為常矢量，而旋轉(zhuǎn)矢量函數(shù)則為經(jīng)度和緯度的函數(shù)）展開為球面坐標的函數(shù)，其中16個控制點作為未知參數(shù)；3）通過最小二乘法對基于模型和觀測的速度場和應變率進行擬合，得到16個控制點，從而得到全球分布的旋轉(zhuǎn)矢量函數(shù)；4）基于歐拉定理，利用得到的旋轉(zhuǎn)矢量函數(shù)計算全球速度場，然后利用速度場計算速度場梯度張量場，根據(jù)速度場梯度張量場可以得到第二不變應變率。

本文使用的應變率模型為GSRMv1.2，該模型采用了全球4 281個GPS測站的5 170個水平速度場，第二不變應變率的分辨率為0.2°×0.2°，覆蓋范圍為68°S～80°N。圖1為基于全球應變率模型GSRMv1.2 的第二不變應變率等值線圖。從圖中可以看出，板塊邊界地區(qū)的應變率比較大，同時有些板塊的內(nèi)部也存在較大形變，如青藏高原地區(qū)的應變圖表現(xiàn)出歐亞板塊內(nèi)部存在一定的形變，說明板塊并非嚴格的剛體。

圖1 基于全球應變率模型GSRMv1.2的第二不變應變率等值線圖Fig.1 Contour line for second invariant strain rate based on GSRMv1.2

通過對全球應變率模型［14－15］GSRMv1.2 中分辨率為0.2°×0.2°的網(wǎng)格節(jié)點第二不變應變率進行雙三次二維平面插值得到測站的應變率，并根據(jù)測站篩選的第3 條準則選取了310 個GPS測站。

2.2 冰期回彈模型

冰期回彈模型包含冰模型和地球結(jié)構(gòu)模型，冰模型主要包含冰在全球范圍內(nèi)的分布模型和冰的演化模型，地球結(jié)構(gòu)模型主要包含地球的分層結(jié)構(gòu)以及分層結(jié)構(gòu)的物理參數(shù)，如密度、粘彈性等。目前最新的冰模型為ICE－5G，廣泛采用的地球結(jié)構(gòu)模型為VM2。本文所采用的冰模型為ICE－3G，地球模型為VM2。冰期回彈模型ICE－5G（VM2）預測的地殼垂直運動速率如圖2（a）所示，圖中顯示冰期回彈對歐亞板塊的北歐地區(qū)、北美板塊和南極洲板塊的部分區(qū)域影響最顯著，垂直運動速率可達到1cm／a甚至更大，而且地殼上升的最大運動速率比地殼下降的運動速率要大。地殼水平運動速率東向分量和北向分量分別如圖2（b）和2（c）所示，圖中顯示北美地區(qū)相對于其他地區(qū)有較大的地殼水平運動速率。對比不同地區(qū)的垂直和水平地殼運動速率看出，垂直運動速率一般比水平運動速率大4～5 倍，但是也存在例外，因為部分地區(qū)具有很小的垂直運動速率，但卻有很大的水平運動速率，這可能與不同地區(qū)地殼對荷載的反應能力不同有關(guān)。

基于冰期回彈模型［16］ICE3G－VM2，用TABOO［17－18］軟件對§2.1中選取的310個測站的垂直和水平位移速率進行估算，根據(jù)測站篩選的第4條準則從中選取242個測站。

圖2 基于冰期回彈模型ICE－5G（VM2）的地殼垂直和水平運動速率的等值線圖Fig.2 Contour line for rate of vertical and horizontal motion of solid earth based on ICE－5G（VM2）

3 計算結(jié)果與分析

利用上述4個準則進行篩選后，GPS測站剩余242個（圖1中紅色圓點所示）。通過將這些測站分派到全球56個板塊中發(fā)現(xiàn)，只有16個板塊含有GPS測站，其中揚子板塊、加勒比板塊和安納托利亞板塊各只有1 個GPS 測站，分別為WUHN、CRO1和ANKR。由于解算板塊的運動參數(shù)至少需要2個測站，所以本文通過最小二乘法剔除速度殘差較大（residual＞3σ）的測站并不斷迭代，直到測站的數(shù)目和殘差的均方根收斂為止，最后位于所有13個板塊上的GPS測站數(shù)目為225個。

基于NNR－MORVEL56的板塊邊界劃分情況，估計全球13個板塊相對于ITRF2008的運動參數(shù)，建立了ITRF2008板塊運動模型，這13個板塊分別為阿穆尼亞板塊、南極洲板塊、阿拉伯板塊、澳大利亞板塊、歐亞板塊、印度板塊、北美板塊、納茲卡板塊、努比亞板塊、太平洋板塊、南美板塊、索馬里板塊和巽他古陸板塊。表1給出了本文和Altamimi［6］計算的ITRF2008 參考架下全球13個板塊的運動參數(shù)、中誤差和水平速度殘差的均方根，表2給出了ITRF2008板塊運動模型和MORVEL板塊運動模型的相對運動參數(shù)。

通過對比本文和Altamimi計算的結(jié)果看出，兩者得到的歐拉矢量基本一致，但有幾個板塊在歐拉矢量的位置上達到6°的偏差：阿穆尼亞板塊在沿緯度方向上存在6.95°的偏差；北美板塊沿經(jīng)度方向存在6.01°的偏差；巽他古陸板塊在沿經(jīng)度方向上存在6.44°的偏差。對比兩者歐拉矢量的中誤差可以看出，本文擬合的中誤差比Altamimi擬合的中誤差要小。尤其是對于巽他古陸板塊，兩者都選擇了GETI和NTUS 兩個GPS測站，但是擬合的結(jié)果卻存在很大的差異。本文計算的歐拉矢量轉(zhuǎn)速為0.346 4°／Ma±0.002 0°／Ma，歐拉極的經(jīng)度和緯度為－89.207°±1.922°、50.880°±4.217°，中誤差較??；Altamimi計算的歐拉極轉(zhuǎn)速為0.388°／Ma±0.308°／Ma，歐拉極的經(jīng)度和緯度為－87.309°±22.186°、44.436°±44.892°，中誤差較大。對比東向速度和北向速度的均方根差可以看出，兩者在數(shù)量級上保持一致，沒有大的偏差。

導致兩者之間個別板塊存在較大差異的主要原因是模型建立的出發(fā)點不一樣。Altamimi建立模型（參見式（1））的過程中考慮了原點偏差率，而且擬合方法為同時對所有板塊進行擬合，以得到所有13個板塊的歐拉矢量和一個原點偏差率，也即一次性擬合45個未知參數(shù)，那么各個板塊之間的歐拉矢量必然存在聯(lián)系，各個板塊歐拉矢量之間的協(xié)方差系數(shù)不會近似為0，同時這些板塊的歐拉矢量與原點偏差率之間也會存在一定的相關(guān)性。本文在建立模型的過程中，不考慮整體原點的偏差率，對13個板塊的歐拉矢量分別進行擬合，每次擬合3個參數(shù)，消除了各個板塊之間的相關(guān)性。由于板塊之間的相關(guān)性，導致巽他古陸板塊歐拉矢量的中誤差很大，因為該板塊僅存在兩個測站，相對其他板塊，這兩個測站所占的權(quán)重較小，因而該板塊的歐拉矢量受其他板塊的影響較大。

測站的篩選準則也會對擬合結(jié)果造成影響。Altamimi采用以遠離板塊邊界至少100km 作為判斷測站是否位于剛性板塊的判據(jù)之一，但是由于板塊邊界具有彌散性和選擇100km 作為邊界距離值具有較大的主觀性，并且在精確計算測站到板塊邊界的過程中，需要計算點到球面多邊形（將板塊近似為球面多邊形）的最短距離，運算量較大，故本文采用全球應變模型作為測站是否位于剛性板塊的判據(jù)。造成偏差的另外一個原因是在最小二乘擬合的過程中，本文采用了3σ準則，即測站的東向和北向速度殘差都要小于3σ，而Altamimi則采用的是3σ－3mm／a準則，即測站水平速度的殘差既要小于3 倍水平速度的不確定度，也要小于3 mm／a。由于測站篩選的準則不同，每個板塊上所篩選出的測站也會不同，最終擬合的歐拉矢量也會有偏差。

比較本文計算結(jié)果與MORVEL板塊運動模型發(fā)現(xiàn)，努比亞－南極洲板塊對的歐拉極在沿經(jīng)度方向上的偏差比較大，達到15.8°；索馬里－南極洲板塊對的歐拉極在沿經(jīng)度和緯度方向的偏差分別達到15.3°和13.5°；歐亞－北美板塊對的歐拉極在沿經(jīng)度和緯度方向的偏差分別達到10.3°和14.4°；努比亞－北美板塊對的歐拉極在沿緯度方向上的偏差非常大，高達30.5°；印度－索馬里板塊對的歐拉極在沿緯度方向上的偏差為9.9°，在轉(zhuǎn)速上的偏差較大，為0.118°／Ma。這些偏差的出現(xiàn)，可能說明努比亞、索馬里、歐亞和印度板塊在1～3 Ma 期間的運動確實發(fā)生了變化。由于ITRF2008板塊運動模型與MORVEL 模型建立的時間尺度不一樣，前者采用了近十幾a的大地測量資料，為短期的結(jié)果，而后者則主要采用近幾Ma的地質(zhì)資料，為長期平均的結(jié)果，兩者之間的偏差項可能意味著驅(qū)動板塊運動的合力矩在零點發(fā)生了波動。另外，測站在空間上的幾何分布不均會對計算結(jié)果造成一定的偏差，并且本文沒有按照傳統(tǒng)的測站篩選準則對測站進行篩選，而是以全球應變率模型為基準進行測站篩選，這也可能是造成偏差的原因，尤其是北美板塊。

表1 ITRF2008 參考架下全球13個板塊的運動參數(shù)Tab.1 Kinematic parameters for 13plates in ITRF2008

表2 ITRF2008 參考架下全球13個板塊對的相對運動參數(shù)Tab.2 Kinematic parameters for 13pairs of plates in ITRF2008

4 結(jié) 語

本文采用基于全球應變模型來判斷測站是否位于形變區(qū)內(nèi)作為測站篩選的準則之一，把這種方法應用到ITRF2008板塊運動模型中，并與Altamimi建立的板塊運動模型和MORVEL 板塊運動模型相比較，結(jié)果表明本文方法是可行的。地球表面的運動及變化本質(zhì)上是由地球內(nèi)部作用力驅(qū)動的，如果已知板塊的運動參數(shù)，則從宏觀上為研究地球內(nèi)部的驅(qū)動力提供了方便。同時，以板塊作為參考基準，可以為研究區(qū)域性地殼應變場和應力場提供借鑒。

［1］Demets C，Gordon R G，Argus D F.Geologically Current Plate Motions［J］.Geophysical Journal International，2010，181（1）：1－80

［2］Bird P.An Updated Digital Model of Plate Boundaries［J］.Geochemistry，Geophysics，Geosystems，2003，4（3）：101－112

［3］Argus D F，Gordon R G，Demets C.Geologically Current Motion of 56Plates Relative to the No－Net－Rotation Reference Frame［J］.Geochemistry Geophysics Geosystems，2011，12（11）：75－87

［4］Altamimi Z，Collilieux X，Métivier L.ITRF2008：An Improved Solution of the International Terrestrial Reference Frame［J］.Journal of Geodesy，2011，85（8）：457－473

［5］Argus D F，Gordon R G，Heflin M B，et al.The Angular Velocities of the Plates and the Velocity of Earth’s Centre from Space Geodesy［J］.Geophysical Journal International，2010，180（3）：913－960

［6］Altamimi Z，Métivier L，Collilieux X.ITRF2008 Plate Motion Model［J］.Journal of Geophysical Research：Solid Earth，2012，117（B7）：47－56

［7］Goudarzi M A，Cocard M，Santerre R.EPC：Matlab Software to Estimate Euler Pole Parameters［J］.GPS Solutions，2014，18（1）：153－162

［8］Peltier W R.Global Glacial Isostasy and the Surface of the Ice－Age Earth：The ICE－5G（VM2）Model and Grace［J］.Earth and Planetary Sciences，2004，32：111－149

［9］Argus D F，Peltier W R.Constraining Models of Postglacial Rebound Using Space Geodesy：A Detailed Assessment of Model ICE－5G（VM2）and Its Relatives［J］.Geophysical Journal International，2010，181（2）：697－723

［10］Klemann V，Martinec Z，Ivins E R.Glacial Isostasy and Plate Motion［J］.Journal of Geodynamics，2008，46（3－5）：95－103

［11］Haines A J，Holt W E.A Procedure for Obtaining the Complete Horizontal Motions within Zones of Distributed Deformation from the Inversion of Strain Rate Data［J］.Journal of Geophysical Research：Solid Earth，1993，98（B7）：12 057－12 082

［12］Holt W E，Chamot－Rooke N，Pichon X L，et al.Velocity Field in Asia Inferred from Quaternary Fault Slip Rates and Global Positioning System Observations［J］.Journal of Geophysical Research，2000，105（B8）：19 185－19 209

［13］David Salomon.Curves and Surfaces for Computer Graphics［M］.New York：Springer，2007

［14］Kreemer C，Haines J，Holt W E，et al.On the Determination of a Global Strain Rate Model［J］.Earth Planets ＆Space，2000，52（10）：765－770

［15］Kreemer C，Holt W E，Haines A J.An Integrated Global Model of Present－Day Plate Motions and Plate Boundary Deformation［J］.Geophysical Journal International，2003，154（1）：8－34

［16］Spada G，Galassi G.New Estimates of Secular Sea Level Risefrom Tide Gauge Data and GIA Modelling［J］.Geophysical Journal International，2012，191（3）：1 067－1 094

［17］Giorgio Spada.The Theory behind TABOO［EB／OL］.http：／／samizdat.mines.edu／taboo／teoria.pdf，2003

［18］Giorgio Spada.TABOO User Guide［EB／OL］.http：／／samizdat.mines.edu／taboo／user guide.pdf，2003