基于多面函數法的青藏高原應變特征分析*

吳嘯龍 楊志強 武繼峰 胡 洋

(長安大學地質工程與測繪學院，西安 710054)

基于多面函數法的青藏高原應變特征分析*

吳嘯龍 楊志強 武繼峰 胡 洋

(長安大學地質工程與測繪學院，西安 710054)

將曲線擬合法解算大區域低頻應變分布與地球參考橢球面上的應變張量矩陣相結合，提出基于多面函數模型解算在給定參考橢球面上的均勻視應變場方法。并根據GPS實測速度場資料解算出青藏高原的均勻視應變場，解算結果與青藏高原各主要構造活動特征以及地震震源機制解所得P、T軸優勢方向分布與變化規律吻合較好。認為基于旋轉橢球面解算大范圍均勻視應變場方法能更加精確地獲取大范圍內地殼形變的總體變化趨勢。關鍵詞 多面函數;參考橢球;應變率;震源機制解;青藏高原

1 引言

WGS84大地坐標系是GPS觀測結果所采用的基本參考框架之一，加之參考橢球幾何模型更加接近地球的真實形態，因此，利用 GPS觀測結果在WGS84定義的參考橢球面上解算地殼應變率，不僅在理論上更加嚴密，也可以進一步減小模型解算過程中的系統誤差。尤其是在大區域地應變解算時，其解算結果將更加逼近真實情況［1］。

青藏高原是世界上海拔最高，地殼最厚、年代最新、面積最廣的高原，具有十分復雜的地質結構與構造特征。自20世紀90年代至今，我國先后在青藏高原及其鄰域布設的一系列GPS觀測站，已經累積了豐富的觀測資料。為了更加精確地獲取大范圍內地殼形變的總體變化趨勢，本文利用現有的GPS觀測資料，基于多面函數理論模型和WGS84所定義的參考橢球，計算了青藏高原的應變張量，獲得了青藏高原整體在參考橢球面上的地殼應變率場分布，并對應變率場分布特征進行了分析。

2 旋轉橢球面上的應變張量

應變張量矩陣是一個對稱矩陣，在二維情況下為［2－4］:

式(2)就是基于球面的應變張量解。由此可知，應變張量的球面解實際上是橢球面解的一種特殊情況。本文主要研究橢球面上的應變場解算方法，忽略式(1)中的大地高及高程變化，并令uh=0，h=0，可得:

3 基于多面函數理論模型求解應變參數

3.1 多面函數模型［5］

多面函數基本數學模型為:

式中，n為所選建立模型的結點的個數，R(x，y;xj，yj)為核函數。經驗上在對地殼運動速率場的擬合中，選取下述函數作為核函數:

式中，(x，y)為推估點的坐標，(xj，yj)為用于建立多面函數模型所選取的測點資料，δ為光滑因子，β為冪數，通過取不同的值來構建不同的核函數，一般取1/2(正雙曲型)，3/2(三次曲面型)，－1/2(負雙曲型)。

為了考察光滑因子的不同取值與核函數的選取對模型構建的有效性。本文采取交叉驗證法對建立的模型進行評價。并分別嘗試利用正雙曲面、倒雙曲面和三次曲面進行建模。選擇50、100、145個數據作為結點，其余觀測點作為檢核點，對高原地區GPS實測速率場進行擬合分析。每次建模時，令光滑因子δ在區間［0，2］內隨機取值。并將東向、北向速率各自的擬合精度求平均數，來評價當前模型的擬合精度。通過反復試驗，當δ=1.14時，得到的擬合結果最好(表1);進一步選取雙三次曲面(β=)作為核函數，多面函數建模效果最好。

表1 多面函數擬合精度Tab.1 Fitting accuracy of multi-surface function

3.2 應變參數解算方法

對GPS速度場的東向分量Ve與北向分量Vn分別建立多面函數模型為:

對式中兩項分別求偏微分得:

4 應變場解算結果及分析

分析數據來源于CMONOC工程，數據處理主要分為四步［6］:1)基于 GAMIT軟件處理載波相位以獲取測站坐標和衛星軌道的單日松馳解;2)利用GLOBK將IGS產出的全球IGS跟蹤站單日松馳解作為約束條件，解算得到包含區域點和IGS站的單日松馳解;3)基于ITRF2000框架下的速度解，采用QOCA軟件對全部獲得的單日松馳解進行平差，估算研究區域內測站坐標和速度;4)基于ITRF2000求解歐亞板塊歐拉矢量，并以此為基礎進一步解算研究區域各測站相對于穩定歐亞參考框架下的相對運動速率場。

首先選定青藏高原及其鄰域75°E ～105°E，23°N～39°N的范圍內建立1°×1°的格網節點;然后利用多面函數模型建立研究區域均勻分布的速率場。如圖1(藍色箭頭代表擬合多面函數模型擬合結果，黑色箭頭代表實測速度場)所示，圖中的實測速度場與擬合速度場二者符合甚好，說明利用多面函數法對青藏高原GPS觀測速度場進行擬合是一種十分有效的方法。但是利用多面函數法擬合研究區域速率場類似于低通濾波，過濾掉了局部變化或者偏差較大的數據，進而使得到的速率場的整體趨勢性特征更加明顯。從圖1中看出，青藏高原分布的速度場整體上從喜馬拉雅碰撞帶的主邊界中心區域向北東方向偏轉，并呈扇形輻射打開。整體上，高原南部的速率場遠大于北部的，同時青藏高原的東部邊緣區域發生順時針偏轉。這顯示出了青藏高原內部形變的兩個基本特征:1)高原內部最為顯著的特點是印度板塊在喜馬拉雅沖撞帶與歐亞大陸發生碰撞，NNE向主壓應力促成了青藏高原整體北東向運移，在青藏高原的北緣受到一系列穩定塊體的阻擋，進而形成高原現今南北向的壓縮;2)高原南部在側向擠壓和重力拖拽的作用下，塑性的下地殼物質流拖動了上地殼脆性層發生東向逃逸。

基于上述建立的函數模型，按應變參數解算方法求解區域均勻視應變。得到研究區域面膨脹值分布、最大剪應變分布以及最大、最小主應變率場(圖2～4)。

圖2顯示青藏高原南北邊界都處于壓縮狀態，在青藏高原南緣的主邊界沖斷帶附近，地殼被劇烈壓縮。這一近NS向的推擠作用通過高原內部各塊體的相互作用向北傳遞，最終在高原北部邊界受到塔里木及阿拉善等穩定塊體的阻擋形成了南強北弱的地殼壓縮。另外在青藏高原中部存在大范圍的引張區。這與該區域正斷型地震多發以及一系列的南北走向的拉張地塹分布十分吻合。許志勤等［7］研究認為在這一區域下地殼存在大規模的地幔上升流可能是造成高原中部相對較大區域的膨脹作用的深部原因。青藏高原內部存在另一個規模宏大、以川滇塊體的邊界走滑斷裂為界的面膨脹區。結合該地區大量的火山分布，推測其上地殼擴展的原因可能是下地殼軟流圈的上涌以及高原物質的大規模東流共同作用結果。

圖2 面膨脹Fig.2 Distribution of surface dilatation rate

圖3顯示出高原內最大剪應變作用最廣泛、強烈的區域同樣位于青藏高原南部的喜馬拉雅主沖斷帶上。鮮水河斷裂與阿薩姆地區的剪應變率次之。最大剪應變分布與高原內部各深大走滑斷裂系分布均可以較好吻合。最大剪應變分布能夠清晰地顯示出阿爾金斷裂帶、祁連山斷裂帶以及鮮水河斷裂帶的基本輪廓，且最大剪應變幾乎全部傾向于于青藏高原內部。說明青藏高原相對于其周圍穩定地體，具有明顯的側向擠出趨勢。如此也驗證了基于橢球面解算大區域地殼均勻視應變場分布的方法不僅在理論上更加嚴密，也可以得到十分精細的結果。

綜合來看，面膨脹值的大范圍引張地區與最大剪應變分布的峰值區在川滇地塊十分吻合，這可能預示著青藏高原的東向擠出在高原東部受到揚子地臺的阻擋之后，形成的以“東構造節”為中心的滇藏漩渦。“漩渦構造”的旋轉主體正是青藏高原地區剪應變最為顯著區域之一的川滇菱形塊體。不難看出，川滇菱形塊體邊界及內部的一系列斷裂均顯示出明顯的剪切走滑特征。而對應嘉黎斷裂、理塘斷裂的共軛剪切和鮮水河-小江斷裂系的左旋走滑，可能是由于旋轉體內外圈旋轉速度差異造成的。地震波層析資料表明，這一區域分布的大型活動走滑斷裂多為巖石圈-地幔剪切帶。這似乎也暗示著青藏高原“東向逃逸”實際上是以川滇菱形地塊為通道，以大規模軟塑性物質流拖動上地殼脆性巖石層的形式發生的。

由圖4可看出，青藏高原均勻分布的最大、最小主應變場具有類似于GPS速度場的空間變化趨勢:主應變方向自西向東存在順時針旋轉，且以南北向的擠壓和東西向的拉張為其主要特征。最大主壓應變速率的峰值出現在喜馬拉雅主沖撞帶附近，其值約(5～9)×10－8/a;其次，高值分布區域主要出現在阿薩姆地區至川滇一帶。對最大、最小主應變軸優勢方向的特征分布與空間變化進行歸納統計，結果如圖5所示。青藏高原內部主應變率方位分布統計特征及變化趨勢與徐紀人等［8］根據區域內發生的中、強型地震震源機制解得出的地殼巖石圈應力場P、T軸分布方位具有較好的一致性(表2)。

整體來看，來自印度次大陸的NNE向的推擠作用產生的構造應力為現今青藏高原地區地殼巖石圈形變的主要影響。高原內部主壓應變方位集中于近南北向分布:自南向北，青藏高原內主壓應變方位由其南端的密集分布，向北過渡逐漸向近NS方向的均勻打開。圖6所示為青藏高原內部自西向東的主張應變方位分布。不難看出，青藏高原內部普遍存在近東西向的拉張。張應變優勢方位自高原內部過渡到東部邊緣時，發生的明顯的順時針偏轉。

表2 震源機制解P、T軸主方向與主應變率軸方位對比Tab.2 Contrast between the main direction of P，T-axis of the focal mechanism solutions and the principal strain rate axis

5 結論

本文利用多面函數法建立了青藏高原及其鄰近區域地殼水平形變函數模型，得到了1°×1°的格網節點上建立了擬合速率場。通過對比實驗，認為當光滑因子δ取1.14，核函數為三次曲面時，建模效果最優。根據此多面函數模型進一步導出了基于旋轉橢球面上的應變張量公式，并解算得到了青藏高原及其臨區的均勻視應變場。

對青藏高原應變場分析發現，青藏高原內部存在南北向的擠壓與東西向的拉伸。最大主壓應變的峰值主要集中于喜馬拉雅沖斷帶附近。高原內部東西向的拉張與其間分布的一系列南北向正斷層與地塹較好對應。面膨脹值分布特征清晰地顯示出高原的南北擠壓邊界、側向擠出通道與順時針旋轉等高原形變的主要特征。最大剪應變與高原南北邊界及內部一系列大型走滑斷層空間分布非常吻合。計算得到的高原總體以及各細部地殼主應變率方位均與地震震源機制解得到的P、T軸方向可以較好吻合。主壓應變方向在高原南部相對集中且近似與NS分布，向北逐漸扇形打開且趨于東向偏轉。主張應變優勢方位自西向東發生明顯的順時針轉動，說明本文方法具有較好的精確度與可靠性。

1 賴錫安，等.中國大陸現今地殼運動［M］.北京:地震出版社，2004.(Lai Xian，et al.Present-day crustal movement in China constrained［M］.Beijing:Seismological Press，2004)

2 劉序儼，黃聲明，林巖釗.地形旋轉張量探討［J］.大地測量與地球動力學，2010，(5):57－63.(Liu Xuyan，Huang Shengming and Lin Yanzhao.Research on rotation tensor of crustal deformation［J］.Journal of Geodesy and Geodynamics，2010，(5):57 －63)

3 劉序儼，黃聲明，梁全強.旋轉橢球面上的應變與轉動張量表達［J］.地震學報，2007，29(3):240 －249.(Liu Xuyan，Huang Shengming and Liang Quangqiang.Expression of Strain and rotation tensor in geodetic coordinates［J］.Acta Seismologica Sinica，2007，29(3):240 －249)

4 陳健，晁定波.橢球大地測量學［M］.北京:測繪出版社，1989.(Chen Jian and Chao Dingbo.Ellipsoidal geodesy［M］.Beijing:Surveying and Mapping Press，1989)

5 石耀林，朱守彪.用GPS位移資料計算應變方法的討論［J］.大地測量與地球動力學，2006，(1):1 －8.(Shi Yaoling and Zhu Shoubiao.Discussion on method of calculating strain with GPS displacement data［J］.Journal of Geodesy and Geodynamics，2006，(1):1 －8)

6 王琪，張培震，馬宗晉.中國大陸現今構造變形GPS觀測數據與速度場［J］.地學前緣，2002，9(2):415－429.(Wang Qi，Zhang Peizhen and Ma Zongjin.Chinese contemporary tectonic deformation of the GPS data and velocity field［J］.Earth Science Frontiers，2002，9(2):415 －429)

7 許志勤，等.青藏高原與大陸動力學——地體拼合、碰撞造山及高原隆升的深部驅動力［J］.中國地質，2006，33(2):221 －238.(Xu Zhiqin，et al.The Qinghai-Tibet Plateau and continental dynamics:a review on terrain tectonics，collisional orogenesis，and progress and mechanisms for the rise of the plateau［J］.Geodesy in China，2006，33(2):221－238)

8 徐紀人，趙志新，石川有三.中國大陸地殼應力場與構造運動區域特征研究［J］.地球物理學報，2008，51(3):770－781.(Xu Jiren，Zhao Zhixin and Ishikawa Yuzo.Regional characteristics of crustal stress field and tectonic motions in and around Chinese mainland［J］.Chinese J Geophys.，2008，51(3):770 －781)

ANALYSIS ON QINGHAI-TIBET PLATEAU STRAIN FIELD BASED ON MULTI-SURFACE FUCTION

Wu Xiaolong，Yang Zhiqiang，Wu Jifeng and Hu Yang
(Academy of Geological Engineering and Surveying，Chang’an University，Xi’an710054)

Combined the method of curve fitting when calculating a lower-frequency strain field in a large area with the strain tensor matrix on the earth reference ellipsoid，we put forward a method of calculating strain parameters on reference ellipsoid by using multi-surface function.According to the GPS data，we get the Qinghai-Tibet plateau strain field.The result is in line with both characters of the major tectonic activity within this area and the main direction of P，T-axis from focal mechanism solutions.

multi-surface function;reference ellipsoid;strain rate;focal mechanism solution;Qinghai-Tibet plateau

P315.72+5

A

1671-5942(2013)04-0017-05

2013-01-12

吳嘯龍，男，1988年生，在讀博士生，研究方向為地殼形變監測與動態大地測量學.E－mail:xlong_wu@foxmail.com