李東臣 任俊杰 張志文 劉 亮

1)應(yīng)急管理部國家自然災(zāi)害防治研究院，北京 100085 2)復(fù)合鏈生自然災(zāi)害動力學(xué)應(yīng)急管理部重點(diǎn)實(shí)驗(yàn)室，北京 100085 3)中國地震局地殼動力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100085 4)應(yīng)急管理部國家減災(zāi)中心，北京 100124 5)應(yīng)急管理部衛(wèi)星減災(zāi)應(yīng)用中心，北京 100124

0 引言

大地震現(xiàn)場調(diào)查表明，震級>6.5級的地震通常會在地表形成幾千米至數(shù)十千米的地表破裂帶(鄧起東等，2008； 鐵瑞等，2016； 徐錫偉等，2018)。地表破裂的幾何結(jié)構(gòu)包含發(fā)震構(gòu)造的運(yùn)動學(xué)特征，不僅能為分析活動斷裂的同震位移空間分布規(guī)律和活動斷裂變形帶的寬度提供重要的定量資料(Klingeretal.，2005； Xuetal.，2009； 陳立春等，2010)，而且對于深入理解發(fā)震斷裂淺層地殼的破裂傳播規(guī)律與變形運(yùn)動學(xué)、 動力學(xué)機(jī)制等方面具有重要的科學(xué)意義(Brooksetal.，2017； Renetal.，2019； Rodriguez Padillaetal.，2021)。

目前獲取地震地表破裂主要采用2種方法： 1)野外地質(zhì)調(diào)查法。通過震后現(xiàn)場地表破裂實(shí)地調(diào)查獲得同震地表位移和破裂帶寬度等信息(徐錫偉等，2008； 李海兵等，2014； 張波等，2015)。2)影像目視解譯法?；谶b感影像，根據(jù)自身對地震地表破裂特征的理解，手動繪制地表破裂帶特征(徐岳仁等，2015； 劉超等，2021； 袁兆德等，2021)。這2種方法雖然可以獲得同震地表破裂的大致幾何結(jié)構(gòu)，但均具有一定的局限性。野外地質(zhì)調(diào)查法的可靠性高，但內(nèi)陸型地震多發(fā)生在地形地貌復(fù)雜的區(qū)域，很多地點(diǎn)難以到達(dá)，從而導(dǎo)致數(shù)據(jù)不完整，無法詳細(xì)繪制破裂帶的細(xì)節(jié)特征，同時野外調(diào)查往往費(fèi)時費(fèi)力(楊瑩輝等，2017)。影像目視解譯法雖然能解譯野外地質(zhì)調(diào)查無法到達(dá)區(qū)域的地表破裂，但結(jié)果精度容易受解譯者經(jīng)驗(yàn)的影響，整個解譯過程仍是時間和勞動密集型工作，且提取結(jié)果多為線型地表破裂，難以準(zhǔn)確獲取地表破裂寬度等精細(xì)特征(Renetal.，2019； Rodriguez Padillaetal.，2021； 張志文等，2021)。因此，地震地表破裂尤其是微破裂面精細(xì)結(jié)構(gòu)的自動提取是活動構(gòu)造研究中一項(xiàng)亟待解決的問題。

通過衛(wèi)星平臺獲取的遙感影像分辨率較低，且易受天氣因素影響，提取的地表破裂精細(xì)度不夠。而無人機(jī)平臺使用成本低，可低空飛行，不受云霧的影響，且可獲取厘米級分辨率的影像(郭慶華等，2021)，為大地震同震地表破裂精細(xì)結(jié)構(gòu)的提取提供了條件。本研究基于無人機(jī)航拍獲得的2021年青?，敹?.4級地震地表破裂的高分辨率影像，依據(jù)顏色空間色彩分割理論和面向?qū)ο蟮乃枷?，提出了一種面向?qū)ο蟮摹按址指?精提取”方法，可實(shí)現(xiàn)地震地表破裂帶特征的半自動提取。該研究為快速精確刻畫地震地表破裂提供了一種解決方案，有助于深入理解大地震發(fā)震斷裂地表變形的幾何學(xué)特征和運(yùn)動學(xué)機(jī)制。

1 面向?qū)ο蟮摹按址指?精提取”方法

航拍地表破裂得到的高分辨率無人機(jī)影像地面場景復(fù)雜，如影像上存在豐富的地物信息及大量噪聲，利用基于像元的機(jī)器學(xué)習(xí)方法提取地表破裂時，存在未充分利用影像信息、 分類效果差、 速度慢等局限(黃雪青，2008； 趙慧等，2012)。且目前尚未有公開的地表破裂數(shù)據(jù)集，無法在影像上應(yīng)用深度學(xué)習(xí)方法提取地表破裂。因此，本研究基于地表破裂在影像上的顏色和形狀等先驗(yàn)特征，提出了一種面向?qū)ο蟮摹按址指?精提取”方法。

圖 1 面向?qū)ο蟮摹按址指?精提取”方法Fig. 1 Object-oriented “rough segmentation-fine extraction” method.

該方法主要包括3個相互銜接的部分： 數(shù)據(jù)預(yù)處理、 地表破裂粗提取和地表破裂精提取，如圖 1 所示。該方法的具體流程為： 第1步，為了提高方法的提取精度和速度，要對原始的高分辨率無人機(jī)影像進(jìn)行預(yù)處理。本文首先設(shè)計(jì)了自定義不規(guī)則多邊形柵格裁剪方法，以快速去除大部分背景噪聲，獲得地表破裂感興趣區(qū)ROI(Region of Interest)。最后將ROI從RGB顏色空間轉(zhuǎn)換為HSV顏色空間。第2步，進(jìn)行地表破裂粗提取。首先，對ROI的HSV顏色空間進(jìn)行色彩分割，即通過設(shè)置不同的H、S、V值分割ROI圖像得到地表破裂候選區(qū)。由于人工選擇閾值分割圖像效率低，本研究設(shè)計(jì)了動態(tài)調(diào)整參數(shù)程序，實(shí)時顯示提取結(jié)果，以最大程度提取地表破裂并去除背景噪聲(即粗提取規(guī)則)，最終獲得了地表破裂候選區(qū)。第3步，采用數(shù)學(xué)形態(tài)學(xué)中的膨脹方法處理地表破裂候選區(qū)圖像。第4步，為了提高提取精度，基于面向?qū)ο蟮乃枷耄瑢Φ乇砥屏押蜻x區(qū)進(jìn)行子區(qū)域?qū)ο筇崛?。?步，進(jìn)行地表破裂精提取。首先，計(jì)算子區(qū)域?qū)ο蟮奶卣鲄?shù)，然后通過動態(tài)調(diào)整特征參數(shù)篩選子區(qū)域?qū)ο?，并?shí)時觀察精提取結(jié)果以獲得最好的精提取結(jié)果。最后，將精提取結(jié)果輸出為shp格式的矢量面數(shù)據(jù)。

2 半自動提取過程

基于以上方法，本研究采用Python語言設(shè)計(jì)了地震地表破裂半自動提取方法。下面將詳細(xì)介紹地震地表破裂半自動提取方法的基本原理和提取流程。

2.1 數(shù)據(jù)與分析

2021年5月22日，青?，敹喟l(fā)生了MS7.4 地震。震后快速調(diào)查表明，地震造成了長約160km的同震地表破裂帶，形成的左旋同震位錯量平均為1～2m(李智敏等，2021； 潘家偉等，2021)。本次地震的同震地表破裂帶由復(fù)雜的破裂類型組成，包括拉張裂縫、 剪切裂縫、 擠壓鼓包等(張志文等，2021)，這些現(xiàn)象被認(rèn)為與發(fā)震斷裂上的左旋走滑運(yùn)動相關(guān)的里德爾剪切作用有關(guān)(Renetal.，2021)。本次地震發(fā)生后，我們第一時間到地震現(xiàn)場進(jìn)行調(diào)查，利用大疆Phantom 4 RTK無人機(jī)對本次地震地表破裂的典型段落進(jìn)行了航拍，獲得了分辨率約2cm的高精度影像。本研究選取2021年瑪多地震地表破裂帶東段昌馬河附近的一段典型地表破裂開展實(shí)驗(yàn)研究(圖 2)。

圖 2 2021年青海瑪多 MS7.4 地震的地表破裂與余震分布Fig. 2 Aftershock distribution and coseismic rupture zone of the 2021 MS7.4 Maduo earthquake，Qinghai.斷層數(shù)據(jù)來源于文獻(xiàn)(鄧起東，2007)； 瑪多地震地表破裂來源于文獻(xiàn)(Ren et al.，2021)； 地震余震目錄來源于文獻(xiàn)(徐志國等，2021)； 底圖遙感衛(wèi)星影像來源于Google Earth； 研究區(qū)位于圖中藍(lán)色方框內(nèi)

2.2 自定義不規(guī)則多邊形柵格裁剪

柵格裁剪方式包括規(guī)則多邊形裁剪和不規(guī)則多邊形裁剪，為了盡可能多地去除背景噪聲，提高地表破裂的有效提取效率，首先要裁剪出影像上有效的地表破裂范圍?，F(xiàn)有的遙感圖像處理軟件雖然已經(jīng)實(shí)現(xiàn)了柵格裁剪功能，但需基于不規(guī)則多邊形矢量文件，不具有靈活性，無法實(shí)時得到裁剪結(jié)果。考慮到時間成本，本研究設(shè)計(jì)了自定義不規(guī)則多邊形柵格裁剪算法，基于輸入的影像可實(shí)時實(shí)現(xiàn)裁剪功能，具體算法如下：

(1)掩膜圖像生成： 在窗口顯示的影像上手動繪制任意多邊形，根據(jù)該多邊形生成掩膜圖像(黑白二值圖)，其中掩膜圖像中多邊形區(qū)域內(nèi)的像素值為255，區(qū)域外值為0。

(2)待掩膜圖像生成： 計(jì)算出手畫多邊形的四至范圍(四至范圍即手畫多邊形上下左右的像素坐標(biāo))，根據(jù)四至范圍裁剪原圖像得到待掩膜圖像。

(3)掩膜過程實(shí)現(xiàn)： 基于掩膜圖像和待掩膜圖像通過卷積運(yùn)算生成最終的掩膜結(jié)果ROI，卷積運(yùn)算公式為

(1)

其中，mask表示步驟(1)中的掩膜圖像，image表示步驟(2)中的待掩膜圖像。

通過對影像進(jìn)行自定義不規(guī)則多邊形柵格裁剪，得到感興趣區(qū)ROI。整幅高分辨率無人機(jī)影像數(shù)據(jù)及裁剪結(jié)果如圖 3 所示，獲取的ROI用于后續(xù)的地表破裂粗提取。

圖 3 自定義不規(guī)則多邊形柵格裁剪結(jié)果Fig. 3 Raster clipping result of custom irregular polygon.a 高分辨率無人機(jī)影像數(shù)據(jù)； b 裁剪結(jié)果。綠框?yàn)樽远x不規(guī)則多邊形

2.3 顏色空間轉(zhuǎn)換

顏色空間不同，反映的圖像特性和人眼感知特性也不同。為更好地基于ROI進(jìn)行彩色圖像分割，需要對ROI進(jìn)行顏色空間轉(zhuǎn)換。ROI的顏色空間是RGB，顏色空間分布不均勻，人眼視覺感知和像素距離存在較大差異，直接基于RGB顏色空間進(jìn)行彩色圖像分割的效果不佳。而HSV顏色空間分布均勻，符合人眼視覺感知，常被用于顏色識別、 顏色跟蹤、 彩色對象提取。因此，本研究將ROI的RGB顏色空間轉(zhuǎn)化為HSV顏色空間。

HSV又名六角錐體模型(Hexcone Model)，由3種感知變量組成： 色調(diào)(Hue)、 飽和度(Saturation)、 明度(Value)。其中，色調(diào)(H)表示顏色的差別，飽和度(S)表示顏色的飽和程度，明度(V)表示顏色的明亮程度。Smith(1978)根據(jù)顏色的直觀特性提出了HSV顏色空間，并給出了RGB與HSV之間的轉(zhuǎn)換關(guān)系，如式(2)、 (3)和(4)所示。其中，R、G、B分別是三通道上歸一化后的像素值，即R，G，B∈[0，1]； max＿Value和min＿Value分別是三通道上(即R、 G、 B)像素值的最大值和最小值；m是最大值和最小值的差值，即m=max＿Value-min＿Value； G＿minus＿B是G和B通道對應(yīng)位置上像素值的差值，即G＿minus＿B=G-B。

由于地表破裂由偏黑色的暗像元組成，背景噪聲由偏亮的像元組成?；谝陨戏治隹芍琀SV顏色空間適用于影像粗分割以提取地表破裂候選區(qū)(王博等，2017)。

(2)

(3)

V=max＿Value

(4)

2.4 地表破裂粗提取

基于HSV顏色空間分割影像，首先要根據(jù)H、S、V三者在圓錐模型中的相互變化關(guān)系量化地震地表破裂的色彩模板。為有效提取地表破裂候選區(qū)，需要將顏色空間量化為地表破裂區(qū)與背景噪聲區(qū)(胡焯源等，2017)。由于三分量H、S、V之間彼此獨(dú)立，根據(jù)六角椎體模型的顏色分布特性(王琪，2011)，可以確定H、S、V三分量之間關(guān)系，即當(dāng)V非常高或非常低時，H值的變化對于圖像視覺特性幾乎無意義； 當(dāng)S非常低時，H不穩(wěn)定，同一物體對象可能對應(yīng)不同的H值； 當(dāng)V非常高或非常低時，S值的變化對于圖像視覺特性幾乎無意義。因此，為了最大程度提取地表破裂并減少背景噪聲，一般不改變S(飽和度)，先調(diào)整V(亮度)，均衡地表破裂和背景噪聲，再適當(dāng)調(diào)整H值，減少背景噪聲，最終得到地表破裂候選區(qū)。

在確定基于HSV顏色空間設(shè)置色彩模板提取地表破裂候選區(qū)的算法后，通過分析地表破裂的顏色特性，初步確定H、S、V三變量的默認(rèn)閾值，S的最小值Sat＿min=0，S的最大值Sat＿max=255，H的最小值Hue＿min=0，V的最小值Val＿min=0，則H的最大值Hue＿max與V的最大值Val＿max將直接決定地表破裂候選區(qū)中地表破裂的精細(xì)程度。 若閾值設(shè)置過高(Hmax=58，Vmax=105)，分割結(jié)果會出現(xiàn)大量噪聲(圖4a)； 若閾值設(shè)置過低(Hmax=58，Vmax=39)，噪聲能夠基本消除，但部分有效破裂也受到了抑制(圖4b)。 因此需要確定合適的Hue＿max與Val＿max，以最大限度保留地表破裂，同時有效抑制噪聲。

圖 4 HSV顏色空間色彩分割結(jié)果Fig. 4 The results of HSV color space color segmentation.a 閾值過高時的分割結(jié)果； b 閾值過低時的分割結(jié)果。 白色區(qū)域即分割得到的地表破裂候選區(qū)

傳統(tǒng)的閾值選擇方法是通過人工不斷輸入有限且離散的閾值，得到多個提取結(jié)果，并從中選擇最優(yōu)結(jié)果所對應(yīng)的閾值作為最佳閾值。但這種方法步驟繁瑣，難以快速修改閾值，閾值覆蓋不夠全面，且無法實(shí)時可視化提取結(jié)果，操作效率較低。因此，本文將提取結(jié)果的輪廓疊加在輸入圖像上，并在顯示窗口上設(shè)置閾值動態(tài)調(diào)整滑動條，通過動態(tài)調(diào)整閾值(圖 5)，實(shí)時觀察提取結(jié)果，經(jīng)前后對比選出最優(yōu)提取結(jié)果，最終從彩色影像中快速分割出地表破裂候選區(qū)(圖 6)。

圖 5 粗分割動態(tài)調(diào)整界面Fig. 5 The dynamic adjustment interface of coarse extraction.

圖 6 地表破裂候選區(qū)Fig. 6 The ROI area of coseismic surface ruptures.a 原始影像； b 地表破裂的提取結(jié)果(紅色即為提取結(jié)果的輪廓)。 底圖為原始影像灰度圖

具體分割算法實(shí)現(xiàn)如下：

(1)顏色空間轉(zhuǎn)換： 根據(jù)式(2)—(4)可知，H∈[0，360]，S∈[0，1.0]，V∈[0，1.0]，但影像數(shù)據(jù)類型一般為8bit的無符號整型，為了匹配影像數(shù)據(jù)類型8bit的特性，調(diào)整H、S、V的范圍，使得H∈[0，180]，S∈[0，255]，V∈[0，255]。

(2)參數(shù)動態(tài)調(diào)整： 基于研究區(qū)影像，通過動態(tài)調(diào)整H和V，最終確定V∈[0，85]和H∈[0，32]，最大程度保留了地表破裂，減少了背景噪聲，至此實(shí)現(xiàn)了地表破裂粗提取(圖 5)。

2.5 子區(qū)域?qū)ο筇崛?/h3>

數(shù)學(xué)形態(tài)學(xué)(mathematical morphology)是一種非線性濾波方法，分為二值數(shù)學(xué)形態(tài)學(xué)和灰度數(shù)學(xué)形態(tài)學(xué)2種運(yùn)算。 本研究使用的是二值數(shù)學(xué)形態(tài)學(xué)運(yùn)算，將待處理的二值圖像當(dāng)作一個集合，并用結(jié)構(gòu)元素在該集合上進(jìn)行“探測”。結(jié)構(gòu)元素是一個可以在圖像上平移且小于圖像尺寸的集合。所用的基本運(yùn)算是使結(jié)構(gòu)元素在圖像范圍內(nèi)進(jìn)行平移，并施加交、 并等集合運(yùn)算，本質(zhì)是通過圖像集合與結(jié)構(gòu)元素間的相互作用來提取有意義的圖像信息(羅朝陽等，2020)，不同的結(jié)構(gòu)元素可以提取不同層面的圖像信息，以滿足圖像分析、 識別、 提取的目的。膨脹是二值數(shù)學(xué)形態(tài)學(xué)的基本運(yùn)算之一，對于二值圖像中的亮區(qū)域具有膨脹(即擴(kuò)大)和平滑作用，基本思想是： 當(dāng)目標(biāo)區(qū)域包含結(jié)構(gòu)元素的任一點(diǎn)時，結(jié)構(gòu)元素原點(diǎn)對應(yīng)的像素點(diǎn)屬于該區(qū)域； 反之，則不屬于該區(qū)域(陸可等，2021)。由于粗提取的地表破裂輪廓疊加在原始影像上顯得更小(即結(jié)果輪廓顯示在破裂邊界內(nèi)側(cè))，平滑性和連續(xù)性也較差，因此采用膨脹運(yùn)算膨脹和平滑粗提取結(jié)果。

傳統(tǒng)的圖像分類、 分割或?qū)ｎ}信息的提取主要基于像元的光譜特性，這種方式僅孤立考慮單個像元的灰度特征，并未充分利用地物形狀、 紋理、 空間上下文關(guān)系等對象特征。因此，僅在噪聲少的高質(zhì)量影像或中低分辨率影像上進(jìn)行分類、 分割或?qū)ｎ}信息提取的效果顯著。而高分辨率影像場景復(fù)雜，含有豐富的地物信息，且像元一般是非混合像元，只是某一地物的局部，使用傳統(tǒng)方法所得結(jié)果極易包含大量噪聲，如果用濾波抑制噪聲又會降低影像的分辨率，影像就會失去高分辨率的意義。而面向?qū)ο蟮乃枷敫先祟惸恳暯庾g的思維邏輯，能夠較好地克服“同物異譜”或“異物同譜”，可以極大提高高分辨率影像分類分割或?qū)ｎ}信息提取的精度(曾舉等，2011)。

地表破裂粗提取是基于像素光譜特征進(jìn)行分割的過程。通過觀察發(fā)現(xiàn)，粗提取結(jié)果是由眾多分散的子區(qū)域組成，為了將子區(qū)域作為對象進(jìn)行處理，必須提取子區(qū)域的外輪廓(李政國，2008)，外輪廓可根據(jù)輪廓跟蹤算法計(jì)算得到(Suzukietal.，1985)，該算法的核心思想是： 基于鏈碼原理，在二值圖像上搜索連通域的輪廓點(diǎn)及其方向，并用鏈碼表示，從中標(biāo)記出內(nèi)外輪廓。粗提取結(jié)果同樣也是二值圖像，因此，可通過輪廓跟蹤算法處理地表破裂候選區(qū)圖像，得到子區(qū)域的外輪廓，即子區(qū)域?qū)ο蟆?/p>

2.6 面向?qū)ο蟮牡乇砥屏丫崛?/h3>

以計(jì)算得到的子區(qū)域?qū)ο鬄榛締卧诿嫦驅(qū)ο蟮乃枷雽Φ乇砥屏堰M(jìn)行精提取。地表破裂的形態(tài)有扁長型、 近圓型等，而背景噪聲的形態(tài)也復(fù)雜多樣，通過周長、 面積、 矩形度、 圓形度、 最小外接矩形的長寬比、 形狀指數(shù)等形狀特征無法很好地將兩者區(qū)分。但微調(diào)形狀特征參數(shù)可去除部分背景噪聲，且大部分背景噪聲呈近圓形，因此在算法中選擇最小外接矩形的長寬比輔助地表破裂的精提取。

光譜特征為地物經(jīng)過吸收與發(fā)射后在特定傳感器上成像時的數(shù)值特征，一般同類地物光譜特征相同或相近，故可用光譜特征提取某種地物(鞠明明，2013)。無人機(jī)影像分辨率高，背景噪聲復(fù)雜多樣且相對于地表破裂更勻質(zhì)，通過多次實(shí)驗(yàn)發(fā)現(xiàn)光譜標(biāo)準(zhǔn)差可以很好地區(qū)分背景噪聲和地表破裂。又由于背景的光譜均值較大，在算法中可利用均值有效濾除背景噪聲。綜上所述，為高效精提取地表破裂，選用子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差作為主要特征，光譜均值、 最小外接矩形的長寬比作為輔助特征提取地表破裂。

面向?qū)ο蟮牡乇砥屏丫崛∷惴ǖ木唧w實(shí)現(xiàn)過程包括：

(1)對象特征計(jì)算： 基于提取的子區(qū)域?qū)ο蠛脱芯繀^(qū)高分影像，分別計(jì)算子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差、 光譜均值、 最小外接矩形的長寬比。

(2)動態(tài)參數(shù)調(diào)整： 將提取的地表破裂輪廓作為前景與研究區(qū)影像疊加顯示，通過動態(tài)調(diào)整子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差，實(shí)時觀察提取效果，直到最大程度提取地表破裂和剔除背景噪聲，再適當(dāng)調(diào)整輔助特征，盡可能多地去除背景噪聲。精提取動態(tài)參數(shù)調(diào)整界面如圖 7 所示。

(3)矢量結(jié)果輸出： 將提取的地震地表破裂輸出為shp格式的矢量面數(shù)據(jù)。

圖 7 精提取動態(tài)調(diào)整界面Fig. 7 The dynamic adjustment interface of fine extraction.

3 結(jié)果分析

本研究以高分辨率無人機(jī)影像為數(shù)據(jù)源，實(shí)現(xiàn)了地震地表破裂帶的半自動提取，在實(shí)驗(yàn)區(qū)提取的地震地表破裂結(jié)果如圖 8 所示，其中共包含1125個矢量面。通過對比全局提取效果發(fā)現(xiàn)(圖8a，b)，提取結(jié)果雖仍含有少量背景噪聲，但研究區(qū)內(nèi)地表破裂面的提取效果非常好。為了更直觀地展示提取結(jié)果的精細(xì)結(jié)構(gòu)，將提取結(jié)果輪廓與原始影像疊加顯示(圖8c—k)，結(jié)果表明： 1)受到破壞的地表破裂提取效果不佳(圖8c)，這是因?yàn)樵撎幍乇砥屏咽艿秸鸷蠛恿鳑_刷的影響，與未受破壞的地表破裂相比，其顏色偏亮且邊界不明顯，地表破裂特征并不顯著。另外，該結(jié)果是使用全局閾值從原始影像中提取的，因此提取效果存在一定的局限性。2)本文提出的方法能夠滿足半自動提取地表破裂面精細(xì)結(jié)構(gòu)的要求(圖8d—k)。提取結(jié)果中雖仍含有少量背景噪聲，但其輪廓與原始影像上的地表破裂邊界基本吻合，提取結(jié)果中甚至包含某些次級地表破裂(圖8j 中的黃框)。同時需要指出的是，提取結(jié)果的某些細(xì)節(jié)存在問題(圖8f 中的黃框)，這是由于在本研究所使用的分辨率約2cm的影像中，這些次級破裂很小，很多為主地表破裂內(nèi)部的細(xì)節(jié)。3)僅以原始影像為數(shù)據(jù)源，本文方法甚至可提取單像素寬且連續(xù)性差的地表破裂面(圖8g，h)。

綜上所述，利用本文方法從高分辨率無人機(jī)影像中提取的地表破裂面與影像中的地表破裂邊界基本一致，表明該方法具有有效性和可靠性。

圖 8 地表破裂的提取結(jié)果Fig. 8 Extraction results of surface rupture.a 地表破裂所在范圍的研究區(qū)原始影像； b 地表破裂提取結(jié)果，底圖為原始影像的灰度圖； c—k 地表破裂輪廓線的細(xì)部結(jié)構(gòu)

圖 9 影像分塊后的提取結(jié)果Fig. 9 Extraction results of image segmentation.a 河流處地表破裂的原提取效果； b 分塊后河流處地表破裂的提取效果； c 圖8f處的提取效果； d 分塊后圖8f處的提取效果，其中綠色矢量面為原提取結(jié)果，紅色輪廓為分塊后的提取結(jié)果

4 討論

4.1 分塊提取地表破裂

由于受河流沖刷破壞的地表破裂提取效果不佳，為驗(yàn)證本研究方法的可靠性，對影像進(jìn)行分塊以使用局部閾值提取地表破裂(圖9a 中的影像底圖即為分塊后的影像)。先用局部閾值對影像進(jìn)行粗分割，然后進(jìn)行精提取，提取效果如圖9b 所示。提取結(jié)果的輪廓與影像中的地表破裂邊界基本吻合，且輪廓線較圓滑。雖然該地表破裂受到了河流的沖刷，使其與原生地表破裂呈現(xiàn)出不同的特征，無法通過全局閾值對其進(jìn)行提取，但將影像分塊后，采用局部閾值可提取該地表破裂，使得提取效果有了極大提升。

另外，部分特殊的噪聲也可能影響提取的精度。圖9c 紫框中提取結(jié)果的輪廓內(nèi)含有次級地表破裂，將影像進(jìn)行分塊后，再利用本研究的方法提取地表破裂，并與原提取結(jié)果進(jìn)行對比(圖9d)。結(jié)果表明，在保證大部分地表破裂能被提取的情況下，排除了所有背景噪聲，并提取出次級地表破裂。

綜上所述，在本研究方法的基礎(chǔ)上，設(shè)置分塊提取地表破裂的功能，使用局部閾值改進(jìn)全局閾值的不足之處，可以更為準(zhǔn)確地提取地震地表破裂。

4.2 提取的目標(biāo)需求

對地表破裂提取的目標(biāo)需求不同，其最終提取效果也不同。1)當(dāng)提取需求為最大程度提取地表破裂并盡可能減少背景噪聲時，提取方法需要兼顧較大破裂和較小破裂(單像素寬且連續(xù)性差的地表破裂)，提取結(jié)果會含有一些背景噪聲(圖10a)，共包括1125個矢量面。2)當(dāng)提取需求是可忽略部分較小破裂但需盡可能減少背景噪聲時，提取結(jié)果中的背景噪聲明顯減少，其中共包括692個矢量面(圖10b)。對比不同需求下提取結(jié)果的精細(xì)結(jié)構(gòu)可以發(fā)現(xiàn)(圖10c—f)，當(dāng)需求改變后，提取結(jié)果中的較小破裂減少了一部分，同時背景噪聲也明顯減少。因此，提取目標(biāo)需求不同，地表破裂提取結(jié)果也不同，這也表明本研究提出的方法具有靈活性。

圖 10 需求改變后的提取結(jié)果Fig. 10 Extraction results after demand changes.a 本文的提取結(jié)果； b 需求改變后的提取結(jié)果； c、 e 本文提取結(jié)果輪廓的精細(xì)結(jié)構(gòu)； d、 f 需求改變后提取結(jié)果輪廓的精細(xì)結(jié)構(gòu)

4.3 數(shù)學(xué)形態(tài)學(xué)——膨脹

為了提高地表破裂提取的準(zhǔn)確性，本研究設(shè)計(jì)的算法使用了二值數(shù)學(xué)形態(tài)學(xué)中的膨脹運(yùn)算。為了驗(yàn)證膨脹運(yùn)算對于地表破裂提取的重要性，采用未使用膨脹運(yùn)算和使用膨脹運(yùn)算的算法分別對高分辨率無人機(jī)影像進(jìn)行地震地表破裂半自動提取。圖11a 為采用含有膨脹運(yùn)算算法的提取結(jié)果，可見要素更具有連續(xù)性，數(shù)目較少。而未采用膨脹運(yùn)算的算法提取結(jié)果中的要素則更離散，數(shù)目相對較多(圖11b)。為進(jìn)一步對比膨脹運(yùn)算的效果，我們比較2種算法提取結(jié)果的細(xì)部結(jié)構(gòu)(圖11c—h)，發(fā)現(xiàn)未采用膨脹運(yùn)算算法的提取結(jié)果輪廓在影像中地表破裂邊界的內(nèi)側(cè)且提取結(jié)果更破碎，對較小破裂提取效果不佳(圖11f)。采用膨脹運(yùn)算后的提取結(jié)果輪廓和影像中的地表破裂邊界基本吻合，輪廓鋸齒少(即更平滑)，可很好地提取小破裂。因此，在算法中加入膨脹運(yùn)算可有效提高地表破裂提取的準(zhǔn)確性。

圖 11 采用膨脹運(yùn)算與否的提取結(jié)果Fig. 11 Extraction results with or without expansion operation.a 采用膨脹運(yùn)算后的提取結(jié)果； b 未采用膨脹運(yùn)算的提取結(jié)果； c、 e、 g 采用膨脹運(yùn)算后提取結(jié)果輪廓的精細(xì)結(jié)構(gòu)； d、 f、 h 未采用膨脹運(yùn)算提取結(jié)果輪廓的精細(xì)結(jié)構(gòu)

4.4 特征參數(shù)的選擇

為了深入探究用于地表破裂精提取的特征參數(shù)主次關(guān)系，下面將選擇不同的特征參數(shù)組合精提取地震地表破裂，并將提取結(jié)果與本研究的提取結(jié)果進(jìn)行對比。

經(jīng)對比發(fā)現(xiàn)，當(dāng)參數(shù)閾值相同時，僅使用子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差為約束參數(shù)提取地表破裂，其結(jié)果和本文提取結(jié)果中的地表破裂總體形態(tài)一致，但含有更多背景噪聲(圖12a，b)。其中共包括1536個矢量面，而本文的提取結(jié)果中包含1125個矢量要素。由此可知，在此基礎(chǔ)上又選擇子區(qū)域?qū)ο蟮墓庾V均值和最小外接矩形的長寬比作為輔助約束參數(shù)，可進(jìn)一步減少背景噪聲。分析不同特征組合對于地表破裂精細(xì)結(jié)構(gòu)提取的影響(圖12c—f)，發(fā)現(xiàn)2種情況下的提取結(jié)果包含的地表破裂形態(tài)、 數(shù)目均一致。不同的是僅使用子區(qū)域的光譜標(biāo)準(zhǔn)差作為約束參數(shù)時，提取結(jié)果包含更多的背景噪聲； 而另一種情況可進(jìn)一步減少結(jié)果中的背景噪聲，從而降低了錯分率。因此，本研究算法中以子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差為精提取主特征、 以光譜均值和最小外接矩形的長寬比為輔助特征是合理的。

圖 12 選取不同特征參數(shù)時的提取結(jié)果Fig. 12 Extraction results with different characteristic parameters.a 本文的提取結(jié)果； b 僅使用子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差為約束參數(shù)的提取結(jié)果； c、 e 本文提取結(jié)果輪廓的精細(xì)結(jié)構(gòu)； d、 f 僅使用子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差為約束參數(shù)的提取結(jié)果輪廓的精細(xì)結(jié)構(gòu)

4.5 方法分析

本研究所述方法存在以下優(yōu)點(diǎn)： 1)未使用濾波處理影像以去除噪聲，不會破壞影像真值，從而可充分利用影像高分辨率的特性，有效地提取出地表破裂的精細(xì)結(jié)構(gòu)。2)通過設(shè)計(jì)自定義不規(guī)則多邊形柵格裁剪功能以去除部分背景噪聲，提高了本文方法的效率和提取效果。3)基于面向?qū)ο蟮乃枷?，利用地表破裂對象的光譜特征和形狀特征，提高了地表破裂的提取精度。4)設(shè)置滑動條動態(tài)調(diào)整參數(shù)代替人工選擇閾值分割影像以實(shí)現(xiàn)實(shí)時粗分割和精提取過程，可快速提取地表破裂。5)與以往影像目視解譯的線性地表破裂不同，本文方法能夠半自動提取地表破裂矢量面，從而可進(jìn)一步準(zhǔn)確計(jì)算寬度等定量參數(shù)，為理解發(fā)震斷裂運(yùn)動學(xué)機(jī)制奠定基礎(chǔ)。6)目前深度學(xué)習(xí)方法需要大量的訓(xùn)練標(biāo)注數(shù)據(jù)，雖然一些自動化標(biāo)注工具能夠?qū)崿F(xiàn)自動化標(biāo)注，但是仍然會受到使用場景的限制，標(biāo)注之后的樣本需要再次進(jìn)行人工處理，成本很高。而本研究的方法可作為地表破裂樣本數(shù)據(jù)集的生成工具，通過人工適當(dāng)干預(yù)，自適應(yīng)得到標(biāo)注結(jié)果，從而能夠?qū)崟r、 準(zhǔn)確地滿足不同的標(biāo)注需求。

但該方法也存在2個不足之處： 1)未能完全自動化，需要人為調(diào)整滑動條以確定相關(guān)參數(shù)，提取結(jié)果仍然會受到提取者主觀經(jīng)驗(yàn)的影響； 2)對于一些特征不明顯或受后期破壞的地表破裂的整體提取效果不佳，需要利用分塊提取以提高精度。

5 結(jié)論

針對大地震地表破裂面自動化提取的難題，本研究提出了一種面向?qū)ο蟮摹按址指?精提取”方法，基于2021年青海瑪多MS7.4 地震地表破裂的分辨率約2cm的無人機(jī)影像，執(zhí)行以下步驟提取地表破裂： 首先，通過自定義不規(guī)則柵格裁剪方法裁剪原始圖像得到感興趣區(qū)ROI。其次，將ROI顏色空間轉(zhuǎn)換為HSV，再利用亮度值和色調(diào)將ROI的HSV顏色空間分割為地表破裂候選區(qū)，然后采用膨脹運(yùn)算處理地表破裂候選區(qū)； 再利用輪廓跟蹤方法將地表破裂候選區(qū)對象轉(zhuǎn)化為一系列子區(qū)域?qū)ο蟆Ｗ詈?，利用子區(qū)域?qū)ο蟮墓庾V標(biāo)準(zhǔn)差、 光譜均值、 最小外接矩形的長寬比為特征參數(shù)對地表破裂進(jìn)行精提取，并輸出地表破裂矢量面結(jié)果。實(shí)驗(yàn)結(jié)果表明，除受影像分辨率的限制、 地表破裂受河流沖刷破壞導(dǎo)致的幾處提取效果不佳的區(qū)域外，提取結(jié)果具有較高的可靠性，可用于大地震地表破裂定量參數(shù)的提取和地表破裂數(shù)據(jù)集(用于深度學(xué)習(xí))的制作。而為了進(jìn)一步驗(yàn)證方法的有效性并提高提取的準(zhǔn)確率，又結(jié)合不同的提取目標(biāo)需求，采用分塊提取，最終構(gòu)建了一套完整的基于高分辨率無人機(jī)影像的地震地表破裂精細(xì)結(jié)構(gòu)半自動提取技術(shù)框架。