基于聚類分析的川滇地區(qū)塊體劃分與應變研究

張小青　許才軍　王帥

摘要 本文研究川滇地區(qū)塊體劃分與區(qū)域地殼應變.首先，對川滇地區(qū)GPS水平速度場采用K-medoids聚類法直接聚類;然后，基于改進的K-medoids方法，利用GPS站點地理位置與水平速度場對初步聚類結(jié)果聯(lián)合約束進行第2次聚類分析，并結(jié)合F統(tǒng)計檢驗法確定各塊體邊界;最后，運用整體旋轉(zhuǎn)線性應變模型計算各塊體的應變參數(shù)，分析水平運動、應變場空間分布與區(qū)域構(gòu)造變形的關(guān)系.研究結(jié)果表明：川滇地區(qū)的GPS聚類分析結(jié)果與地質(zhì)學提出的板塊劃分結(jié)果相吻合，為塊體劃分提供另一種方式;整體旋轉(zhuǎn)線性應變模型計算的應變結(jié)果顯示川滇菱形塊體表現(xiàn)為最大剪應變和面應變率梯度高值區(qū)，表明該區(qū)域積累了一定的彈性應變能，具有發(fā)生地震的危險性.

關(guān)鍵詞 GPS速度場;川滇地區(qū);聚類分析;應變

中圖分類號 P315

文獻標志碼 A

0 引言

川滇地區(qū)位于青藏高原東緣，為印度板塊向青藏高原東北向擠壓、碰撞作用的前沿地帶，區(qū)內(nèi)構(gòu)造活動顯著和強震頻發(fā)，是中國大陸構(gòu)造變形最復雜和最強烈的區(qū)域之一.近年來，隨著GPS技術(shù)的進一步發(fā)展以及GPS觀測網(wǎng)絡的密集布設，利用GPS觀測數(shù)據(jù)研究塊體運動和地殼變形，已經(jīng)成為地殼形變研究的熱點[1-5].利用GPS 觀測到的地殼形變速度場，結(jié)合地質(zhì)、地球物理學等觀測資料，通過構(gòu)建地表變形的場源數(shù)學物理模型，可進一步求解地殼應力-應變參數(shù)，進而深入掌握區(qū)域地殼形變的內(nèi)在機制特征[1-5].黨亞民等[1]根據(jù)兩期GNSS實測速度場，運用平移-旋轉(zhuǎn)-應變模型，研究了川滇地區(qū)地震前后塊體運動與應變特征.溫揚茂等[2]聯(lián)合GPS資料、重力資料反演了川滇地區(qū)的現(xiàn)今地殼形變，結(jié)果表明川滇地區(qū)地殼運動模式更趨向于連續(xù)變形.丁開華等[3]根據(jù)川滇地區(qū)1998—2004 年間的GPS觀測結(jié)果，基于地質(zhì)學給出的塊體劃分模型，反演了川滇地區(qū)各活動塊體的現(xiàn)今地殼運動和變形參數(shù)，認為川滇地區(qū)的地殼變形模式以連續(xù)和漸變的剪切為主.依據(jù)劉峽等[4]、黨亞民等[1]和廖思佩等[6]對川滇地區(qū)活動塊體劃分的研究結(jié)果，川滇地區(qū)活動塊體可分為巴顏喀拉塊體（Ⅰ）、川滇菱形塊體（Ⅱ）、滇東塊體（Ⅲ）、藏東塊體（Ⅳ）和滇西南塊體（Ⅴ）.其中巴顏喀拉塊體（Ⅰ）以龍日壩斷裂為界又分為阿壩（Ⅰ1）和龍門山（Ⅰ2）兩個次級塊體;川滇菱形塊體（Ⅱ）以麗江—小金河斷裂為界可進一步劃分為川西北次級塊體（Ⅱ1）和滇中（Ⅱ2）兩個次級塊體，各板塊分布如圖1所示.

合理有效地劃分構(gòu)造活動塊體是構(gòu)造運動及地震活動性研究的重要基礎.許多研究者致力于板塊，尤其是板內(nèi)塊體及構(gòu)造亞板塊的劃分[3-10].石耀霖等[7]對基于GPS 觀測資料計算的歐拉極進行聚類分析，將中國大陸劃分為 9 個活動塊體，該方法在計算歐拉極時存在研究范圍大、計算精度低的不足，且研究區(qū)域內(nèi)GPS測站較少，分布不均勻，板塊劃分難于精細.Savage等[8]對加利福尼亞東南部莫哈韋地區(qū)的GPS水平速度進行聚類分析，將莫哈韋地區(qū)分成4個聚類地塊，分類結(jié)果與Meade[9]的劃分結(jié)果類似.

本文首先利用K-medoids聚類算法對川滇地區(qū)的GPS水平速度場進行聚類分析，獲得初始聚類中心;然后，聯(lián)合GPS站點地理位置與速度場數(shù)據(jù)對初始聚類結(jié)果進行約束，重新判斷各點與初始聚類中心的距離，進行第2次聚類，并運用F統(tǒng)計檢驗方法進一步篩選GPS測站，得到最終的聚類分析結(jié)果;最后，基于塊體劃分結(jié)果，利用整體旋轉(zhuǎn)線性應變模型計算川滇地區(qū)地殼運動與變形參數(shù).

1 聚類算法與數(shù)據(jù)

1.1 K-medoids聚類算法

K中心點算法（K-medoids）的基本思想是：對含有n個數(shù)據(jù)的數(shù)據(jù)集，欲分k個簇，首先任意選擇k個不同的數(shù)據(jù)作為初始簇中心，再根據(jù)每個數(shù)據(jù)對象到簇中心的距離，將它們重新分配到距離最近的簇;然后通過反復地用非中心點代替初始簇中心，而選取的標準就是當該樣本點成為新的中心點后能提高類簇的聚類質(zhì)量，使得類簇更緊湊.K-medoids相比k-means算法采用均值計算法獲得中心點，能改善k-means的“噪聲”敏感特性[11].

1.2 川滇地區(qū)GPS水平運動速度場

本文利用Liang等[12]、Gan 等[13]公布的川滇地區(qū)1999—2013年間GPS觀測數(shù)據(jù).數(shù)據(jù)處理主要有以下4步：

1）利用GAMIT軟件處理GPS數(shù)據(jù)，獲得站坐標和衛(wèi)星軌道的單日松弛解;

2）將獲得的單日松弛解與SOPAC提供的全球解用GLOBK軟件進行融合;

3）用QOCA軟件根據(jù)卡爾曼濾波方法計算全球參考框架ITRF2000站坐標和速度;

4）將基于全球參考框架ITRF2000的測站速度轉(zhuǎn)換到穩(wěn)定歐亞參考框架的速度.

川滇地區(qū)的水平運動速度場由北向南、西向東逐漸減弱，呈順時針渦旋式運動[1-3]（圖1）.

2 聚類結(jié)果

2.1 K-medoids直接聚類

首先，利用K-medoids聚類分析方法對川滇地區(qū)的水平運動速度場進行初步聚類分析.在聚類分析過程中，設定K取不同的值，以顧及不同聚類數(shù)對聚類分析結(jié)果的影響.圖2a—2d分別是聚類數(shù)K取3、4、5和6時的聚類分析結(jié)果，其中同一集群的GPS站點用相同的符號表示.

圖2a為聚類數(shù)K取3時的分析結(jié)果，可以看出川滇菱形塊體與滇東塊體能夠很好地被劃分出來;但巴顏喀拉塊體中的龍門山塊體與滇東塊體被劃分到了一起，其中滇東塊體與滇西南塊體也被分在一起.值得說明是，龍門山塊體中的一個點被錯誤地劃分到川滇菱形塊體中.

圖2b為聚類數(shù)分成4類的分析結(jié)果，可以看出與3類聚類分析結(jié)果相似，川滇菱形塊體與滇東塊體能夠被很好地識別并劃分出來.除此之外，滇東塊體與滇西南塊體也被成功地劃分出來，巴顏喀拉塊體中的阿壩塊體與藏東塊體被分在一起，并且龍門山斷裂帶未被檢測到，玉樹-甘孜斷裂帶也未被檢測到.其中滇東塊體的3個點被錯誤地劃分到滇西南塊體中（圖2b粉色三角形表示的點），同時巴顏喀拉塊體和滇西塊體各有1個點和2個點被錯誤地分到滇東塊體中（藍色方塊點所示）.

分5類的結(jié)果如圖2c所示，與分4類聚類分析結(jié)果相比，巴顏喀拉塊體中的阿壩塊體與滇西塊體已經(jīng)被成功分開，并且能夠檢測到玉樹-甘孜斷裂帶，但龍門山塊體與滇東塊體被分在了一起.雖然能夠識別出滇西塊體與川滇菱形塊體，但滇東塊體中有一些點被錯誤地劃分到滇西南塊體中.

相比前幾類分析結(jié)果，分6類的聚類分析能夠很好地識別出川滇地區(qū)的各個塊體，如圖2d，川滇菱形塊體可分為川西北和滇中兩個次級塊體.

2.2 K-medoids聚類結(jié)果分析

因K-medoids算法聚類中心選擇的隨意性，即使兩次聚類分析取相同的聚類數(shù)，聚類分析的結(jié)果也存在較大的差異，如圖3a—3b所示.所以，需要一定的指標來評定聚類分析的結(jié)果的可靠性.目前聚類分析效果的評定指標主要有輪廓系數(shù)、戴維森堡丁指數(shù)和偏差估計等[14].本文采用輪廓系數(shù)對聚類分析結(jié)果進行評定.將待分類數(shù)據(jù)分成k個簇，對簇中的各點分別計算輪廓系數(shù)S（i）：

其中，a（i）為樣本i到同簇其他樣本的平均距離，b（i）為樣本i到其他簇的所有樣本的平均距離.根據(jù)式（1），簇內(nèi)距離越小，簇間距離越大，也即Sk值越大，聚類效果越好.圖3中，聚類數(shù)為5的圖3a和3b中，圖3a聚類輪廓系數(shù)小于3b的輪廓系數(shù)值，也就是說圖3b聚類效果更好;聚類數(shù)是6的圖3c和3d中，圖3c輪廓系數(shù)小于圖3d，顯然圖3d聚類效果更好.

2.3 改進的K-medoids聚類分析

由于對GPS速度場直接聚類時未考慮站點的地理空間分布連貫性，分類結(jié)果中仍有一些點被錯分.為獲得更好的聚類結(jié)果，本文對K-medoids聚類算法進行改進.首先利用K-medoids進行初步聚類，獲得初始聚類中心;然后聯(lián)合各點地理位置與速度場數(shù)據(jù)對初始聚類中心進行約束，重新判斷各點與原始聚類中心的距離，進行第2次聚類.利用改進后的K-medoids算法獲得的聚類結(jié)果如圖4所示，可以看出聚類效果得到了明顯的改善，不再有錯分的點.對比圖1與圖4b，可看出改進的聚類算法獲得的結(jié)果與基于地質(zhì)學提出的塊體劃分結(jié)果基本吻合.

2.4 GPS 速度解的篩選

值得注意的是，雖然改進的K-medoids算法能夠明顯改善聚類的結(jié)果，但仍有部分測站的速度值與周邊測站明顯不符，出現(xiàn)這種速度奇異點的原因很復雜，可能是由觀測期間測站的不穩(wěn)定性等原因造成的.在板塊構(gòu)造劃分和地殼形變研究中，常使用F統(tǒng)計檢驗法來判定點的相符性[15]，從而研究塊體間的獨立性.采用同時考慮塊體剛性運動和內(nèi)部應變的整體旋轉(zhuǎn)模型來計算各塊體的應變參數(shù)[16]，進而根據(jù)塊體模型和應變參數(shù)預測、估計各測站速度，擬合速度值見圖1中紅色箭頭，黑色箭頭為GPS觀測速度.針對塊體邊界點采用F統(tǒng)計檢驗，若該測站與其他測站相互獨立，則將此站點從該塊體去除，劃到地理相鄰的塊體中.F統(tǒng)計量：

式中n為所選用的測站數(shù)，χ2n為所有測站速度共同運動的歐拉極擬合后的速度殘差，

χ2n-1

為去除某一測站擬合后的速度殘差.

3 塊體運動特征分析

3.1 塊體應變分析

通過上述聚類算法與測站篩選法確定各塊體邊界之后，采用同時考慮塊體剛性運動和內(nèi)部應變的整體旋轉(zhuǎn)模型來計算各塊體的歐拉矢量（表1）與應變參數(shù)（圖5），藏東、巴顏喀拉、滇西南、滇中和滇東塊體的旋轉(zhuǎn)速率與肖卓輝等[17]的對應區(qū)域大致相同.

由于2008年汶川地震后川滇地區(qū)震后形變機制復雜，本文使用的1999—2013年GPS速度場未考慮2008以后川滇地區(qū)的形變，因此本文只對2008年汶川地震之前的應變狀態(tài)進行分析.圖5a為川滇地區(qū)最大剪應變結(jié)果，結(jié)果表明紅河斷裂帶、安寧河斷裂帶、小江斷裂帶和麗江—小金河斷裂帶處的最大剪應變比較顯著，尤其在鮮水河斷裂帶北段以及玉樹—甘孜斷裂帶最為顯著，與江在森[18]的研究結(jié)果相吻合，表明川滇地區(qū)受到較明顯的外力.

從圖5b可以看出，川滇地區(qū)內(nèi)部以面膨脹為主，而沿川滇塊體邊界以面壓縮為主，并且在紅河斷裂帶、安寧河斷裂帶、小江斷裂帶和麗江—小金河斷裂帶處的面應變比較顯著，與江在森[18]以及魏文薪等[19]的研究結(jié)果相吻合.

圖5c是面應變梯度結(jié)果，可以看出川滇地區(qū)西部的面應變梯度比東部大，同樣表現(xiàn)在紅河斷裂帶、安寧河斷裂帶、小江斷裂帶和麗江—小金河斷裂帶，面應變梯度最大可達4.8×10-10 km.龍門山斷裂帶面應變梯度為中值區(qū)，應變特征不明顯.滇西南地區(qū)的面應變梯度高值區(qū)轉(zhuǎn)換為梯度中值區(qū)，這是由于震后地殼的松弛和變形調(diào)整所致.面應變梯度高低值交匯處主要沿區(qū)內(nèi)主干走滑斷裂帶分布，1900年以來發(fā)生的中強以上地震多出現(xiàn)在面應變梯度高低值交匯處.

3.2 塊體間相對運動

在獲取塊體歐拉運動矢量后，本文以GPS點的外包絡線作為塊體邊界，計算了邊界點的相對運動速度，計算結(jié)果如表2所示.從表2可看出，巴顏喀拉塊體與藏東塊體、滇東塊體以擠壓為主，巴顏喀拉塊體與川西北塊體以擴張為主，川西北塊體與藏東塊體、滇中塊體以擠壓為主，滇中塊體與滇西南塊體以擠壓為主，滇中與滇東塊體以擴張為主，這與黨亞民等[1]的研究結(jié)果基本吻合.

4 結(jié)論

本文基于1999—2004年以及1999—2008年間川滇地區(qū)的GPS觀測資料，在K-medoids初步聚類的基礎上，聯(lián)合各點地理位置與速度場數(shù)據(jù)對初始

聚類中心進行約束，進行第2次聚類，能有效改善被錯分的站點，最終將川滇地區(qū)劃分為6個次級塊體，包括巴顏喀拉塊體、川滇菱形塊體、滇東塊體、藏東塊體和滇西南塊體.其中川滇菱形塊體以麗江—小金河斷裂為界可進一步劃分為川西北和滇中兩個次級塊體.塊體劃分結(jié)果與地質(zhì)上塊體劃分模型相吻合，表明聚類分析用于塊體劃分的合理性，可與地質(zhì)上塊體劃分方法形成互補.值得一提的是，利用GPS觀測資料進行聚類分析，塊體劃分精度依賴于GPS站點的密集度.由于各站點速度代表鄰近區(qū)域的平均速度，因此，GPS站點越密集，其所代表的鄰近區(qū)域越小，塊體劃分精度越高.

在聚類分析的基礎上，采用F統(tǒng)計檢驗法能有效去除部分與周邊測站的速度值明顯不符的測站，估算得到的各塊體歐拉矢量與其他研究者得到的結(jié)果相吻合.給出的應變特征表明：紅河斷裂帶、安寧河斷裂帶、小江斷裂帶和麗江—小金河斷裂帶處的最大剪應變比較顯著，尤其在鮮水河斷裂帶處最為顯著，川滇地區(qū)西部是最大剪應變作用最廣泛、最強烈的區(qū)域.川滇地區(qū)以面膨脹為主，整個川滇菱形塊體為應變能和面應變率梯度高值區(qū)，為不同活動方式的構(gòu)造斷裂共同作用的結(jié)果.

參考文獻

References

[1] 黨亞民，楊強，梁詩明，等.川滇區(qū)域活動塊體運動與應變特征地震影響分析[J].測繪學報，2018，47（5）：559-566

DANG Yamin，YANG Qiang，LIANG Shiming，et al.Block movement and strain characteristics effected by earthquake in Sichuan-Yunnan region[J].Acta Geodaetica et Cartographica Sinica，2018，47（5）：559-566

[2] 溫揚茂，許才軍.聯(lián)合GPS與重力資料反演分析川滇地區(qū)現(xiàn)今地殼形變[J].武漢大學學報（信息科學版），2009，34（5）：568-572

WEN Yangmao，XU Caijun.Contemporary crustal deformation in Sichuan-Yunnan region from GPS and gravity data[J].Geomatics and Information Science of Wuhan University，2009，34（5）：568-572

[3] 丁開華，許才軍，鄒蓉，等.利用GPS分析川滇地區(qū)活動地塊運動與應變模型[J].武漢大學學報（信息科學版），2013，38（7）：822-827

DING Kaihua，XU Caijun，ZOU Rong，et al.Crustal movement and strain model of active blocks analyzed by GPS in Sichuan-Yunnan region[J].Geomatics and Information Science of Wuhan University，2013，38（7）：822-827

[4] 劉峽，馬瑾，杜雪松，等.川滇主要斷裂帶近期運動變化及與地震活動關(guān)聯(lián)性[J].中國科學（地球科學），2016，59（6）：1267-1282

LIU Xia，MA Jin，DU Xuesong，et al.Recent movement changes of main fault zones in the Sichuan-Yunnan region and their relevance to seismic activity[J].Science China Earth Sciences，2016，59（6）：1267-1282

[5] 金雙根，朱文耀.太平洋板塊運動和形變及其邊緣現(xiàn)今相對運動[J].大地測量與地球動力學，2002，22（2）：57-60

JIN Shuanggen，ZHU Wenyao.The motion and deformation of pacific plate and contemporary relative motion at its boundaries[J].Crustal Deformation and Earthquake，2002，22（2）：57-60

[6] 廖思佩，侯強，杜永超.基于GPS形變資料的川滇地區(qū)應力場數(shù)值模擬研究[J].大地測量與地球動力學，2016，36（7）：645-649

LIAO Sipei，HOU Qiang，DU Yongchao.Numerical simulation of the stress field in Sichuan-Yunnan region with GPS deformation data[J].Journal of Geodesy and Geodynamics，2016，36（7）：645-649

[7] 石耀霖，朱守彪.利用GPS觀測資料劃分現(xiàn)今地殼活動塊體的方法[J].大地測量與地球動力學，2004，24（2）：1-5

SHI Yaolin，ZHU Shoubiao.Method for division of present active crustal blocks by gps survey data[J].Crustal Deformation and Earthquake，2004，24（2）：1-5

[8] Savage J C，Simpson R W.Clustering of GPS velocities in the mojave block，southeastern california[J].Journal of Geophysical Research（Solid Earth），2013，118（4）：1747-1759

[9] Meade B J.Block models of crustal motion in southern California constrained by GPS measurements[J].Journal of Geophysical Research，2005，110（B3）：B03403

[10] Aranganayagi S，Thangavel K.Clustering categorical data using silhouette coefficient as a relocating measure[C]∥International Conference on Conference on Computational Intelligence & Multimedia Applications.IEEE，2007

[11] 朱純，吳建華，潘毅.關(guān)于優(yōu)化K-medoids聚類算法搜索策略研究[J].計算機仿真，2016，33（10）：244-248，277

ZHU Chun，WU Jianhua，PAN Yi.Research about optimizing search strategy of K-medoids clustering algorithm[J].Computer Simulation，2016，33（10）：244-248，277

[12] Liang S M，Gan W J，Shen C Z，et al.Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements[J].Journal of Geophysical Research（Solid Earth），2013，118（10）：5722-5732

[13] Gan W J，Zhang P Z，Shen Z K，et al.Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements[J].Journal of Geophysical Research，2007，112（B8）：B08416

[14] 賈虎，鄧力琿.基于流線聚類人工智能方法的水驅(qū)油藏流場識別[J].石油勘探與開發(fā)，2018，45（2）：312-319

JIA Hu，DENG Lihui.Oil reservoir water flooding flowing area identification based on the method of streamline clustering artificial intelligence[J].Petroleum Exploration and Development，2018，45（2）：312-319

[15] 郭東美，許厚澤.基于統(tǒng)計檢驗法研究Amurian板塊相對獨立性[J].地球物理學報，2013，56（4）：1246-1252

GUO Dongmei，HSU Houtse.Research on the existence of Amurian plate based on statistical tests[J].Chinese Journal of Geophysics，2013，56（4）：1246-1252

[16] 李延興，張靜華，何建坤，等.由空間大地測量得到的太平洋板塊現(xiàn)今構(gòu)造運動與板內(nèi)形變應變場[J].地球物理學報，2007，50（2）：437-447

LI Yanxing，ZHANG Jinghua，HE Jiankun，et al.Current-day tectonic motion and intraplate deformation-strain field obtained from space geodesy in the Pacific Plate[J].Chinese Journal of Geophysics，2007，50（2）：437-447

[17] 肖卓輝，許才軍，江國焰，等.汶川地震前十年間龍門山區(qū)域顧及斷層閉鎖的地殼應變場[J].地球物理學報，2017，60（3）：953-961

XIAO Zhuohui，XU Caijun，JIANG Guoyan，et al.Crustal strain in the Longmenshan region considering fault locking during ten years before the 2008 Wenchuan earthquake[J].Chinese Journal of Geophysics，2017，60（3）：953-961

[18] 江在森.利用動態(tài)大地測量資料研究中國大陸構(gòu)造形變與強震關(guān)系[D].武漢：武漢大學，2012

JIANG Zaisen.Research on the relationship between the tectonic deformation and the strong earthquake in the Chinese mainland using dynamic geodetic data[D].Wuhan：Wuhan University，2012

[19] 魏文薪，江在森，劉曉霞，等.川滇地區(qū)應變率場分布及其變化特征研究[J].地震，2015，35（4）：11-20

WEI Wenxin，JIANG Zaisen，LIU Xiaoxia，et al.Distribution and variation characteristics of strain rate field in the Sichuan：Yunnan region[J].Earthquake，2015，35（4）：11-20

Block partitioning and strain study in Sichuan-Yunnan

region based on cluster analysis

ZHANG Xiaoqing1 XU Caijun2 WANG Shuai2

1 Department of Hydraulic Engineering，F(xiàn)ujian College of Water Conservancy and Electric Power，Yongan 366000

2 School of Geodesy and Geomatics，Wuhan University，Wuhan 430079

AbstractThis paper focuses on block partitioning and strain calculation in the Sichuan-Yunnan region.First，the K-medoids clustering method is used to cluster the GPS horizontal velocity field in Sichuan-Yunnan.Next，the initial classification results are further clustered based on constraints from the geographical location and velocity field of GPS stations by using the improved K-medoids method.Subsequently，the F statistical test method is used to determine the boundary of each block.Finally，the strain parameters of each block are calculated using an integral rotation linear strain model and the relationship between horizontal motion，spatial distribution of strain field，and regional tectonic deformation is analyzed.The results of the velocity field clustering analysis of GPS stations in the Sichuan-Yunnan region are consistent with the results of the alternative block division methodology based on geological theory of plate division.The strain results，which are calculated using the integral rotation linear strain model，show that the Sichuan-Yunnan rhombic block has the maximum shear strain and the highestsurface strain rate gradient.In addition，this region has accumulated some elastic strain energy，an indicator of earthquake risk.

Key wordsGPS velocity field;the Sichuan-Yunnan region;cluster analysis;strain