基于無人機(jī)傾斜航空攝影三維點(diǎn)云測量同震傾滑變形研究
——以2021 年瑪多MS7.4 地震地表破裂為例1

李忠武 陳桂華

（中國地震局地質(zhì)研究所, 地震動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100029）

引言

地震地表破裂的定量化特征是刻畫整個(gè)破裂帶同震變形行為的基礎(chǔ)數(shù)據(jù)（Klinger 等，2006；Xu 等，2006；徐錫偉等，2010；Brooks 等，2017）。傳統(tǒng)的野外簡易測量僅可得到精度有限的長度和方位參數(shù)。利用全站儀等進(jìn)行測量可得到較高的測量精度，但效率較低，且受野外工作條件限制（Koehler 等，2021）。衛(wèi)星大地測量方法可獲取大范圍同震形變場（石峰等，2010；Jin 等，2021；He 等，2022），但強(qiáng)變形帶大變形可能導(dǎo)致失相干，無法獲取斷層強(qiáng)變形帶變形場和變形細(xì)節(jié)。具有差分定位功能的無人機(jī)航空攝影是近年來興起的高精度、高效率野外測繪技術(shù)，已大量應(yīng)用于大震現(xiàn)場調(diào)查和測量中（Kozac?等，2021；李智敏等，2021；潘家偉等，2021；張志文等，2021）。

無人機(jī)航空攝影測量一般獲取正射影像（DOM）和數(shù)字高程模型（DEM），用于地表形態(tài)獲取、地貌及變形測量與分析（魏占玉等，2014；曾洵等，2019；劉超等，2021；Little 等，2021；Padilla 等，2022）。同時(shí)，無人機(jī)航空攝影還可獲取高密度點(diǎn)云數(shù)據(jù)。在活動構(gòu)造三維建模與分析中，與DEM 相比，點(diǎn)云數(shù)據(jù)更具優(yōu)勢（高偉等，2017；Howell 等，2020；Wei 等，2021），這是因?yàn)辄c(diǎn)云數(shù)據(jù)更好地保留了原始數(shù)據(jù)，而非規(guī)則插值結(jié)果，且更好地保留了地表以下垂向變化信息。傳統(tǒng)的無人機(jī)航空攝影通過垂直地面下視以獲取地表DEM 為主，在垂向信息獲取方面存在一定不足。無人機(jī)傾斜攝影具有更好的三維建模優(yōu)勢（荊帥軍等，2019），為滿足相關(guān)測繪需求提供了新的解決方案。

對于具有傾滑位移的同震地表破裂，通常產(chǎn)生拉張裂縫和擠壓隆起等具有脆性的垂向不連續(xù)變形，第一時(shí)間采集野外數(shù)據(jù)并在測量拉張量或縮短量時(shí)充分考慮垂向變形信息至關(guān)重要。因此，引入無人機(jī)傾斜攝影測量技術(shù)，通過航空攝影點(diǎn)云數(shù)據(jù)對復(fù)雜、具有拉張性質(zhì)和擠壓性質(zhì)的地震地表破裂進(jìn)行精細(xì)三維建模，從而更準(zhǔn)確地測量同震變形量。

1 2021 年瑪多MS7.4 地震概況

據(jù)中國地震臺網(wǎng)中心測定，2021 年5 月22 日2 時(shí)4 分，青?，敹喟l(fā)生7.4 級地震，儀器震中位置（34.59°N，98.34°E），震源深度17 km。震源機(jī)制研究顯示，該地震發(fā)震斷層對應(yīng)的節(jié)面走向281°，傾角88°，滑動角1°，近乎左旋走滑破裂（張喆等，2021）。InSAR 觀測和野外調(diào)查顯示，地震發(fā)生于昆侖山口-江錯(cuò)斷裂的鄂陵湖-江錯(cuò)-昌麻河段（江錯(cuò)斷裂），產(chǎn)生約160 km 長的同震地表破裂（Chen 等，2021；Jin 等，2021；Ren 等，2021；李智敏等，2021；潘家偉等，2021；He 等，2022）。

此次瑪多地震震中位于黃河鄉(xiāng)附近，向兩側(cè)破裂，西至鄂陵湖南側(cè)，東達(dá)昌麻河鄉(xiāng)東側(cè)河谷。野外調(diào)查顯示，地震沿江錯(cuò)斷裂鄂陵湖-江錯(cuò)-昌麻河段形成4 個(gè)相對連續(xù)的地表破裂段（圖1）。根據(jù)地震地表破裂幾何結(jié)構(gòu)、地質(zhì)地貌特征，此次瑪多地震地表破裂可通過震中附近黃河鄉(xiāng)與格波隆格恰階區(qū)分為2 段（圖1）。2 段呈左行左階排列，階區(qū)表現(xiàn)為以沼澤為特征的山間盆地，發(fā)育大量砂土液化帶和裂縫帶。黃河鄉(xiāng)以西地表破裂帶展布于三疊系基巖山與第四系堆積盆地邊界，而格波隆格恰以東地表破裂帶切割三疊系基巖山和山間谷地。地表破裂最西段走向近東西，分布在鄂陵湖南側(cè)山腳，長約25 km，以走滑為主，兼具正斷變形（張扭）。進(jìn)入野馬灘盆地后，沿盆地北部邊界向東延伸至江錯(cuò)附近，長約20 km，以走滑為主，兼具逆斷變形（壓扭）。在江錯(cuò)至黃河鄉(xiāng)段未觀察到有明顯斷錯(cuò)位移的地表破裂。格波隆格恰以東的拉木草至東湖段，可見連續(xù)的地表破裂切過山坡和溝谷，長約23 km，繼續(xù)向東延伸至沙漠中。朗瑪哦爾-昌麻河段同樣切割基巖山地、山間盆地和昌麻河谷，長約25 km，在其南側(cè)存在1 個(gè)分支，破裂特征與主破裂帶類似，但變形更小。同時(shí)，衛(wèi)星形變測量和余震分布顯示瑪多地震地表變形連續(xù)地集中在較窄的變形帶內(nèi)（Chen 等，2021；王未來等，2021；He 等，2022），可通過變形帶內(nèi)的變形細(xì)節(jié)測量進(jìn)一步分析斷層同震變形。

圖1 2021 年瑪多MS7.4 地震地表破裂展布Fig. 1 Surface rupture zone of the 2021 Maduo MS7.4 earthquake

整個(gè)地震地表破裂帶雖以左旋走滑為主，但在最西段具有一定程度的張扭性質(zhì)，而在野馬灘段具有一定程度的壓扭性質(zhì)，格波隆格恰以東無顯著的傾滑變形。分別在朗瑪加合日和野馬灘西選取2 個(gè)典型觀測點(diǎn)，基于無人機(jī)傾斜攝影測量技術(shù)獲取的高密度點(diǎn)云數(shù)據(jù)對拉張量和縮短量進(jìn)行測量。

2 數(shù)據(jù)采集與處理

2.1 無人機(jī)傾斜攝影數(shù)據(jù)采集與點(diǎn)云數(shù)據(jù)處理

采用大疆M300RTK 型無人機(jī)進(jìn)行破裂帶2 個(gè)典型觀測點(diǎn)航空攝影，搭載具有4 500 萬像素的禪思P1相機(jī)，配置35 mm 鏡頭。航空攝影過程中采用大疆D-RTK2 型移動站進(jìn)行實(shí)時(shí)差分定位（RTK），系統(tǒng)標(biāo)稱的照片定位精度為水平1 cm+1 ppm 和垂直1.5 cm+1 ppm。航空攝影的航向重疊率為72%，旁向重疊率為55%。采用智能擺動模式，從前視、后視、左視、右視及下視角度對目標(biāo)進(jìn)行攝影。室內(nèi)使用Agisoft Metashape 軟件進(jìn)行三維建模獲取點(diǎn)云、正射影像（DOM）、數(shù)字高程模型（DEM）數(shù)據(jù)。觀測點(diǎn)無人機(jī)航空攝影航高、地面分辨率、范圍、建模誤差和點(diǎn)云密度如表1 所示。

表1 觀測點(diǎn)無人機(jī)航空攝影參數(shù)與建模誤差Table 1 UAV aerial photography parameters and modeling errors for two observation sites

2.2 同震傾滑位移測量方法

當(dāng)同震地表變形具有傾向位移時(shí)，通常會形成不同形式的斷層陡坎（于貴華等，2010）。斷層傾滑變形測量方法包括：①直接找到可靠的標(biāo)識點(diǎn)獲取三維運(yùn)動矢量數(shù)據(jù)，如上、下盤相交的線性標(biāo)識（陳桂華等，2009）；②測量傾滑位移的垂直分量，根據(jù)斷層傾角計(jì)算傾滑變形；③通過詳細(xì)測量變形細(xì)節(jié)，解析同震變形帶的垂直位移和水平伸縮位移?，敹嗟卣鸬乇砥屏褞Т蟛糠謪^(qū)域無較好的變形標(biāo)識，難以獲得穩(wěn)定可靠的斷層面傾角等信息。因此，通過無人機(jī)傾斜攝影測量獲取的點(diǎn)云數(shù)據(jù)，分別測量傾滑位移的垂直分量和水平分量，進(jìn)而計(jì)算斷層傾滑位移和斷層面傾角。

2.2.1 傾滑位移垂直分量

斷層傾滑的結(jié)果是其中一盤相對另一盤抬升，在地貌上形成斷層陡坎等（圖2）。通過測量陡坎兩側(cè)同一地貌面高度差h得到斷層的垂直位移量，即傾滑位移的垂直分量（圖2（a）、（b））。在獲取的點(diǎn)云數(shù)據(jù)中，垂直斷層提取地表點(diǎn)生成地形剖面。地形剖面在斷層變形帶兩側(cè)延伸一定距離，以確定未變形的斷層兩盤原始地形。通過線性方程擬合斷層陡坎兩側(cè)地貌面，分別測量變形帶兩端地貌面的高差值，以其平均值作為斷層垂直位移量，以2 個(gè)高差值與其平均值的差值作為斷層垂直位移量的誤差。

2.2.2 傾滑位移水平分量

張扭性同震地表破裂帶通常形成斜列的裂縫（圖2（b）），壓扭性同震地表破裂帶出現(xiàn)鼓包、擠壓脊等微地貌（圖2（a）），通過變形前后的地形對比確定傾滑位移的水平分量（水平拉張量或水平縮短量）。

對于壓扭性同震地表破裂帶，選擇未逆沖疊覆的位置測量水平縮短量（圖2（a）），利用地震前后兩點(diǎn)間距離變化表示水平縮短量。圖2（a）中，地表由于剪切作用出現(xiàn)裂縫，切割形成相對獨(dú)立的小條塊，條塊由于斜向的擠壓作用出現(xiàn)縮短，原為直線的PQ 線段在剖面上整體表現(xiàn)為上凸的弧形彎曲，距離縮短。沿?cái)鄬佣缚财拭鏈y量P、Q 點(diǎn)變形前的距離，再與變形后P、Q 點(diǎn)在地貌面上的距離lPQ作差，得到水平縮短量。需注意的是，當(dāng)剖面中出現(xiàn)裂縫時(shí)，剖面的縮短量等于鼓包部分縮短量減去裂縫部分的拉張量。

圖2 斷層傾滑變形量測量方法示意Fig. 2 The measurement method of dip-slip

對于傾斜地貌面上發(fā)育的壓扭性破裂，測量P、Q 點(diǎn)變形后距離時(shí)，需測量P、Q 點(diǎn)在原始未變形坡面上的距離，而非直線距離。因此，當(dāng)基于最小二乘法擬合陡坎兩側(cè)地貌面時(shí)，應(yīng)將P、Q 點(diǎn)距離分別投影到陡坎下地貌面擬合線P1Q1，得到lPQ1，并將P、Q 點(diǎn)距離分別投影到陡坎上地貌面擬合線P2Q2，得到lPQ2，最后分別計(jì)算2 條擬合線得到的縮短量及誤差。

張扭性同震地表破裂帶傾向變形測量如圖2（b）所示，以未發(fā)生坍塌的地表變形前后剖面長度差作為傾滑位移水平拉張量。具體測量過程為：首先基于破裂空間結(jié)構(gòu)確定斷層帶整體走向、垂直斷層走向、垂直方向，利用點(diǎn)云數(shù)據(jù)生成地形剖面L2；然后通過剖面判斷裂縫是否存在坍塌，如果存在坍塌，則平移剖面位置生成新的剖面；最后在剖面L2 識別裂縫端點(diǎn)M、N，測量點(diǎn)M、N 在陡坎上、下地貌面擬合線上的投影距離lMN，即可得到傾向水平拉張量。當(dāng)剖面中出現(xiàn)擠壓鼓包時(shí)，需減去由鼓包調(diào)節(jié)的水平縮短量。

2.2.3 傾滑位移

通過上述方法獲取斷層傾滑位移垂直分量h與水平伸縮量s后，可通過式（1）得到傾滑位移量D：

3 野馬灘西觀測點(diǎn)壓扭性破裂同震傾滑量測量

野馬灘西觀測點(diǎn)位于長條形野馬灘盆地北側(cè)邊界（圖1），本次地震地表破裂切割山前洪積扇，形成了北高南低的斷層陡坎，發(fā)育大量擠壓鼓包和斜列的地震裂縫（圖3（a）、（b）），變形集中于斷層上升盤，具有擠壓環(huán)境下逆斷層變形特征，是壓扭性破裂帶?；邳c(diǎn)云數(shù)據(jù)，通過垂直破裂帶整體走向提取剖面，獲取斷層傾滑位移。

同震地表破裂的精細(xì)解譯結(jié)果如圖3（c）所示，確定破裂帶的整體走向113°，垂直破裂帶整體走向提取了YMT1、YMT2 剖面（圖3（d）），另外提取平行和垂直裂縫方向的YMT3、YMT4 剖面，以分析剖面方向?qū)鄬觾A滑位移的影響。為充分獲取斷層兩盤變形信息，在統(tǒng)計(jì)研究區(qū)2 個(gè)觀測點(diǎn)附近破裂帶最大寬度為18 m 的基礎(chǔ)上，向兩側(cè)延伸至50 m 提取剖面，然后利用最小二乘法擬合陡坎兩側(cè)線性地形線，確定震前地貌面形態(tài)。

測線剖面如圖4 所示，YMT1 得到的斷層垂直位移為（0.69±0.06）m，YMT2 得到的斷層垂直位移為（1.01±0.16）m，YMT3 得到的斷層垂直位移為（0.92±0.03）m，YMT4 得到的斷層垂直位移為（0.73±0.09）m。由等高線（圖3（d））可知，該觀測點(diǎn)破裂帶走向與洪積扇坡向近于垂直，斷層發(fā)生水平位移時(shí)不會產(chǎn)生顯著的垂直視位移。

圖4 野馬灘西觀測點(diǎn)點(diǎn)云剖面圖Fig. 4 The point cloud profiles of Yematan

根據(jù)前述傾滑位移水平縮短量測量方法，選擇剖面中無逆沖疊覆及塌陷的位置測量野馬灘西觀測點(diǎn)破裂帶水平縮短量。其中垂直斷層剖面YMT1 縮短量為（0.17±0.01）m，垂直斷層剖面YMT2 縮短量為（0.41±0.01）m，而垂直裂縫的剖面YMT3 縮短量為（-0.27±0.01）m，平行裂縫的剖面YMT4 縮短量為（0.54±0.01）m。垂直斷層剖面縮短量代表的是走滑位移在剖面方向的分量，而垂直、平行裂縫的剖面縮短量代表的是水平縮短和走滑位移在剖面方向上的矢量和。

得到斷層傾滑位移水平分量和垂直分量后，基于同震滑移的矢量關(guān)系（圖2（c）），計(jì)算得到垂直斷層剖面YMT1、YMT2 傾滑位移分別為0.71、1.09 m，斷層傾角分別為76°、68°。

4 朗瑪加合日觀測點(diǎn)張扭性破裂同震傾滑量測量

朗瑪加合日觀測點(diǎn)展布于近東西向谷地的南側(cè)（圖1），同震地表破裂如圖5 所示，破裂帶整體走向東西，裂縫和擠壓鼓包相間發(fā)育，但鼓包主要出現(xiàn)在裂縫首尾交接處，且變形整體集中于斷層下降盤，具有拉張環(huán)境下正斷層變形特征，是張扭性破裂帶。跨破裂帶提取4 條地形剖面，剖面LJR1、LJR2、LJR4 與破裂帶垂直，剖面LJR3 與局部裂縫平行。沿剖面位置從點(diǎn)云數(shù)據(jù)中生成地形剖面，如圖6 所示。

圖6 朗瑪加合日觀測點(diǎn)點(diǎn)云剖面圖Fig. 6 The point cloud profiles of Langma Gaheri

剖面LJR1 測得的垂直位移為（0.34±0.06）m，剖面LJR2 測得的垂直位移為（0.54±0.03）m，剖面LJR3 測得的垂直位移為（0.35±0.02）m，剖面LJR4 測得的垂直位移為（0.44±0.02）m，其中方向斜交的剖面LJR1 和LJR3 垂直位移近乎相等。

剖面LJR1 測得的水平拉張量為（1.99±0.01）m，剖面LJR2 測得的水平拉張量為（2.08±0.01）m，剖面LJR3 測得的水平拉張量為（2.58±0.01）m。其中剖面LJR1、LJR3 處的裂縫是瑪多MS7.4 級地震地表破裂帶中最寬的裂縫，在地震過程中有牛、羊掉入其中。垂直裂縫走向的剖面LJR3 較垂直破裂帶整體走向的剖面LJR1 有更大的拉張量。由點(diǎn)云地形剖面可知，剖面LJR4 裂縫發(fā)生了垮塌，因此不能用于水平拉張量的測量。

傾滑位移水平拉張量明顯大于垂直分量。由圖5（d）可知，破裂帶整體走向東西，而地貌面傾向北西，兩者斜交，破裂帶的左旋走滑產(chǎn)生北盤上升的垂直視位移，使傾滑位移的垂直分量測量值等于真值減去走滑產(chǎn)生的垂直視位移。

圖5 朗瑪加合日觀測點(diǎn)觀測結(jié)果Fig. 5 Observation of Langma Gaheri

5 討論

本文通過跨破裂帶剖面的精細(xì)化測量得到野馬灘觀測點(diǎn)傾滑位移為0.71～1.09 m，而該點(diǎn)附近走滑位移為1.8 m，與已有震源機(jī)制研究揭示的近乎純左旋走滑存在差異（張喆等，2021）。

描述斷層形成機(jī)制的安德森斷層模式指出，斷層的傾滑、走滑運(yùn)動和應(yīng)力方向與斷層走向關(guān)系密切，當(dāng)斷層走向與主壓應(yīng)力方向低角度斜交時(shí)，以走滑運(yùn)動為主；而當(dāng)斷層走向與主壓應(yīng)力方向交角增大時(shí)，則可能出現(xiàn)傾滑運(yùn)動（Anderson，1951）。已有學(xué)者通過地表破裂調(diào)查、InSAR 及震源機(jī)制研究等獲取了瑪多Ms7.4 級地震地表破裂帶整體走向?yàn)?76°～285°（華俊等，2021；潘家偉等，2021；張喆等，2021），利用震源機(jī)制反演得到最大主壓應(yīng)力方向?yàn)?40°（張建勇等，2022），符合安德森斷層模式。而本次測量得到的野馬灘觀測點(diǎn)附近破裂帶整體走向?yàn)?00°，因此，初步推測由于破裂帶走向發(fā)生變化，使野馬灘觀測點(diǎn)產(chǎn)生傾滑運(yùn)動分量。

6 結(jié)論

2021 年5 月22 日瑪多MS7.4 地震形成了長達(dá)160 km 的地震地表破裂帶，發(fā)震斷層以左旋運(yùn)動為主，但在西段朗瑪加合日段和野馬灘段處分別具有張扭性質(zhì)及壓扭性質(zhì)。通過無人機(jī)傾斜攝影測量可以獲得高精度三維點(diǎn)云數(shù)據(jù)并用于破裂帶變形定量分析。在2 個(gè)觀測點(diǎn)獲得了同震地表破裂的垂直位移、傾滑位移水平伸縮量。測量結(jié)果顯示，野馬灘觀測點(diǎn)斷層垂直位移為0.69～1.01 m，傾向水平縮短量為0.17～0.41 m，傾滑位移為0.71～1.09 m，基于安德森斷層模式初步推測為破裂帶走向發(fā)生變化，使該點(diǎn)產(chǎn)生傾滑運(yùn)動；朗瑪加合日觀測點(diǎn)斷層垂直位移為0.34～0.54 m，傾向水平拉張量為1.99～2.08 m，由于作為標(biāo)識的地貌面傾向與破裂帶整體走向斜交，垂直位移可能因走滑位移產(chǎn)生的垂直視位移而減小。

致謝 感謝審稿專家提出的寶貴修改意見，感謝編輯部的高效處理，感謝地震科考過程中李濤、蘇鵬、郭鵬、孫浩越、哈廣浩、袁兆德同志的幫助。