利用ITRF2008分析全球板塊運(yùn)動(dòng)

李春曉

1 中國(guó)科學(xué)院云南天文臺(tái)，昆明市官渡區(qū)，650216

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

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

1 數(shù)學(xué)模型

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

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

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

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

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

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

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

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

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

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

誤差橢圓的半長(zhǎng)軸、半短軸以及長(zhǎng)軸方向可通過(guò)對(duì)角化cov（ωλφ）求得，分別對(duì)應(yīng)矩陣的兩個(gè)特征值和較大特征值對(duì)應(yīng)的特征向量。

2 測(cè)站篩選

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

2.1 全球應(yīng)變率模型

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

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

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

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

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）預(yù)測(cè)的地殼垂直運(yùn)動(dòng)速率如圖2（a）所示，圖中顯示冰期回彈對(duì)歐亞板塊的北歐地區(qū)、北美板塊和南極洲板塊的部分區(qū)域影響最顯著，垂直運(yùn)動(dòng)速率可達(dá)到1cm／a甚至更大，而且地殼上升的最大運(yùn)動(dòng)速率比地殼下降的運(yùn)動(dòng)速率要大。地殼水平運(yùn)動(dòng)速率東向分量和北向分量分別如圖2（b）和2（c）所示，圖中顯示北美地區(qū)相對(duì)于其他地區(qū)有較大的地殼水平運(yùn)動(dòng)速率。對(duì)比不同地區(qū)的垂直和水平地殼運(yùn)動(dòng)速率看出，垂直運(yùn)動(dòng)速率一般比水平運(yùn)動(dòng)速率大4～5 倍，但是也存在例外，因?yàn)椴糠值貐^(qū)具有很小的垂直運(yùn)動(dòng)速率，但卻有很大的水平運(yùn)動(dòng)速率，這可能與不同地區(qū)地殼對(duì)荷載的反應(yīng)能力不同有關(guān)。

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

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

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

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

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

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

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

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

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

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

表2 ITRF2008 參考架下全球13個(gè)板塊對(duì)的相對(duì)運(yùn)動(dòng)參數(shù)Tab.2 Kinematic parameters for 13pairs of plates in ITRF2008

4 結(jié) 語(yǔ)

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

［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