曾憲陽， 劉靜,*， 王偉， 姚文倩， 吳靜， 劉小利，韓龍飛， 王文鑫， 邢宇堃， 杜瑞林， 楊緒前

1 中國地震局地質(zhì)研究所地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 天津大學(xué)表層地球系統(tǒng)科學(xué)研究院，地球系統(tǒng)科學(xué)學(xué)院， 天津 300072 3 杭州五代通訊與大數(shù)據(jù)研究院， 杭州 310012 4 中國地震局地震研究所地震大地測量重點(diǎn)實(shí)驗(yàn)室， 武漢 430071

0 引言

地表破裂帶是地震破裂最直觀的現(xiàn)象.大地震的地表破裂帶精細(xì)填圖是詳細(xì)記錄地表破裂帶的幾何形態(tài)、同震位移及其空間分布規(guī)律的重要工作(比如，Ambraseys and Tchalenko, 1970; Tchalenko and Ambraseys, 1970；Sharp et al., 1982; King and Nábělek, 1985; Spotila and Sieh, 1995; Treiman et al., 2002; Haeussler et al., 2004; Liu-Zeng et al., 2009; Quigley et al., 2012; Choi et al., 2018; 韓龍飛等, 2022)，為探索地震發(fā)生機(jī)理、斷層破裂過程及運(yùn)動(dòng)方式提供重要基礎(chǔ)(Sieh et al., 1993; Haeussler et al., 2004; Oskin et al., 2012；Choi et al., 2018；Olson et al., 2019; DuRoss et al., 2020; Scholz, 2002, 2019).比如，有助于構(gòu)建和檢驗(yàn)大地震破裂擴(kuò)展模型,斷裂滑動(dòng)模型以及動(dòng)態(tài)破裂模擬和強(qiáng)地面運(yùn)動(dòng)建模等(Zachariasen and Sieh, 1995; Day and Ely,2002; Olsen et al., 2008; Xu et al., 2016; Lo Yi-Ching et al., 2018; Klinger et al., 2018; Hough et al., 2019; Pollitz et al., 2019; Goldberg et al., 2019; Lozos and Harris, 2020).此外，通過地表破裂觀測到的同震位錯(cuò)也是構(gòu)建地震位移與震級和破裂長度的經(jīng)驗(yàn)關(guān)系(鄧起東等，1992；Wells and Coppersmith,1994; Wesnousky, 2008)及研究古地震復(fù)發(fā)習(xí)性的重要參數(shù)(Schwartz and Coppersmith, 1984; Sieh, 1996; Liu-Zeng et al., 2006)，對預(yù)測未來地震的破裂尺度提供重要的理論基礎(chǔ)，為震后分析地震災(zāi)情、評估斷層的地震危險(xiǎn)性等工作提供重要依據(jù)(Mignan et al., 2015).

然而，對地震地表破裂的精細(xì)填圖并不容易.在以往對青藏高原內(nèi)部及周邊大地震的震后地表破裂調(diào)查中，研究人員一般通過野外實(shí)地考察的方式近距離對地表破裂進(jìn)行詳細(xì)觀測，雖然能夠在局部范圍內(nèi)獲取精度較高的地形地貌數(shù)據(jù)，但是會(huì)遺漏斷層外更細(xì)微的變形，同時(shí)需要大量的人力物力，也受到地區(qū)的限制，很難快速對整個(gè)地表斷裂進(jìn)行精細(xì)填圖.近年來，借助米級分辨率衛(wèi)星光學(xué)影像或無人機(jī)影像對單次地震的地表破裂或多次實(shí)踐累計(jì)地貌的特征進(jìn)行較詳細(xì)描述(Klinger et al., 2005; Xu et al., 2006; 劉靜等，2013；Li et al., 2016; Bi et al., 2017; Ren and Zhang, 2019).由于地表破裂往往比較細(xì)碎，米級分辨率衛(wèi)星影像不能對微小破裂成像，且對垂直位錯(cuò)很難進(jìn)行分辨.所以獲取的數(shù)據(jù)在精度和準(zhǔn)確度上仍然不足再現(xiàn)地表破裂的復(fù)雜度和豐富程度.得益于無人機(jī)在地球科學(xué)領(lǐng)域的廣泛普及和應(yīng)用，利用無人機(jī)攝影測量技術(shù)對強(qiáng)震造成的同震地表破裂進(jìn)行迅速和全面的影像采集是大勢所趨(Pierece et al., 2020; Ponti et al., 2020; Liu-Zeng et al., 2022； 王文鑫等, 2022)，此項(xiàng)技術(shù)極大推動(dòng)了地表破裂帶填圖的完整性和精細(xì)程度，尤其是針對一些無法涉足的區(qū)域.地震發(fā)生后第一時(shí)間對同震地表破裂進(jìn)行高精度無人機(jī)影像采集對于確定同震地表破裂的精確位置、范圍、最大同震位移等都起到了前所未有的推動(dòng)作用(Ponti et al., 2020; Liu-Zeng et al., 2022).

傳統(tǒng)處理高分辨率無人機(jī)影像通常采取人工繪制地表破裂的方法，非常耗時(shí)，同時(shí)高度依賴個(gè)人經(jīng)驗(yàn)，往往會(huì)造成漏畫、過度解譯等情況(Rodriguez Padilla et al., 2022; Mattéo et al., 2021).隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，機(jī)器學(xué)習(xí)(Machine learning)為處理這類遙感影像數(shù)據(jù)提供了更多的可能 (Bergen et al., 2019).近些年來，機(jī)器學(xué)習(xí)已經(jīng)被逐步廣泛地應(yīng)用在地球科學(xué)領(lǐng)域的多個(gè)方面(Tingdahl and De Rooij, 2005; Meier et al., 2007; Adeli and Panakkat, 2009; Lary et al., 2016; 周永章等，2018；龔健雅和季順平, 2018)，如地震震級預(yù)測(Adeli and Panakkat, 2009)、大洋海底地形中海山的自動(dòng)識別和填圖(Valentine et al., 2013)、地形特征識別(Li et al., 2017, 2020; Du et al., 2019)、河道形態(tài) (Beechie and Imaki, 2014)、滑坡陡坎結(jié)構(gòu)(Dunning et al., 2009)、滑坡易發(fā)性研究 (Yilmaz, 2010；Brenning, 2005)、冰川流向識別(Smith et al., 2016)、冰緣地貌(Wainwright et al., 2015)、海岸地貌和海洋沉積物特性(Li et al., 2011; Martin et al., 2015)、地表變形和表面過程的復(fù)雜性和相互作用(Hu et al., 2021)等.目前機(jī)器學(xué)習(xí)在震后地表破裂自動(dòng)識別上的研究相對較少，僅有個(gè)別學(xué)者對陸地?cái)鄬雍土严兜葮?gòu)造要素的識別進(jìn)行了嘗試(Mattéo et al.，2021).因此，我們基于2021年5月22日瑪多MW7.4地震發(fā)生后無人機(jī)航拍采集的地表破裂的高精度影像數(shù)據(jù)，采用機(jī)器學(xué)習(xí)的方法——基于卷積神經(jīng)網(wǎng)絡(luò)(CNN)(LeCun et al., 1999)的有監(jiān)督深度學(xué)習(xí)算法，構(gòu)建了地表破裂識別模型；通過遷移學(xué)習(xí)訓(xùn)練模型，使用模型進(jìn)行運(yùn)算，完成了對研究區(qū)地表破裂的快速識別；通過對比自動(dòng)識別的結(jié)果與人工填圖的結(jié)果，分析了當(dāng)前模型的優(yōu)點(diǎn)和不足，為后續(xù)利用該方法進(jìn)行地表破裂的快速自動(dòng)識別提供了新的優(yōu)化策略.

1 研究背景

2021年5月22日，在青海省果洛藏族自治州瑪多縣發(fā)生了MW7.4地震，造成大量房屋倒塌、部分橋梁等基礎(chǔ)設(shè)施受損.瑪多地震發(fā)生在此前未被識別出來的江錯(cuò)斷裂上，位于東昆侖斷裂以南70 km處.震中位于(34.59°N，98.34°E)，地震震源機(jī)制解為一條近EW走向、具有顯著拉張分量的左旋走滑錯(cuò)動(dòng)事件，發(fā)震斷層走向N102°E，傾向S，傾角81°，滑動(dòng)角-11°(中國地震臺(tái)網(wǎng)中心，2021)(圖1).震后我們對瑪多地震進(jìn)行野外調(diào)查，同時(shí)于5月24日—6月15日對地表破裂全段進(jìn)行了高精度無人機(jī)影像(UAV)數(shù)據(jù)采集，發(fā)現(xiàn)瑪多地震震后形成了長達(dá)158 km的地表破裂帶 (王文鑫等，2022；姚文倩等，2022).地表破裂帶西起鄂陵湖南(97.6135°E，34.7459°N)，向E終止于昌麻河鄉(xiāng)(99.2698°E，34.5102°N)，依據(jù)其走向變化自西往東依次可劃分為鄂陵湖段、野馬灘北段、野馬灘南段、黃河鄉(xiāng)段、東哦段、昌瑪河段和江錯(cuò)分支段等8段(姚文倩等，2022；劉小利等，2022；Liu-Zeng et al., 2022)，主要分布在海拔4200～4600 m的山區(qū)、草原、沼澤、沖積扇和部分沙丘地區(qū)，主要沿N105°E走向展布，由一系列張裂縫、張剪裂縫、剪切裂縫以及擠壓鼓包、擠壓脊和裂陷等雁行狀組合而成(李智敏等，2021；潘家偉等，2021；邵延秀等，2022)，近地表植被主要以低矮的高寒沼澤灌叢—嵩草為主(李寧云，2018)，因此植被對遙感圖像上裂縫的識別影響較少.在影像上，這些地表破裂大都表現(xiàn)為深灰色細(xì)條帶，其寬度和長度分別為幾十厘米至 1～10 m和20～500 m(邵延秀等，2022；姚文倩等，2022；劉小利等, 2022).

圖1 瑪多地震構(gòu)造背景圖Fig.1 Seismogenic tectonic background map of the Madoi earthquake

此次瑪多MW7.4地震是繼汶川地震之后我國發(fā)生的震級最高的一次地震事件，同時(shí)也是獲得地表破裂帶厘米級無人機(jī)影像資料最為詳細(xì)的地震.瑪多地震地表破裂的高精度填圖和特征研究為震后評估地震斷層危害性等工作提供重要依據(jù)，對于地震破裂過程的精細(xì)研究具有重要的科學(xué)意義.

2 方法

2.1 無人機(jī)圖像采集與正射影像(DOM)的生成

我們于5月28日—6月8日用縱橫大鵬CW-15垂直起降固定翼無人機(jī)沿著地表破裂帶進(jìn)行了UAV圖像采集(圖2)，無人機(jī)可以垂直起降，續(xù)航時(shí)間3 h，續(xù)航里程180 km，最高起降海拔>4500 m，抗風(fēng)能力強(qiáng)(王文鑫等，2022).無人機(jī)在CW Commander平臺(tái)上操作，配備JOUAV CA103全畫幅正射相機(jī)(4200萬像素)，操作中使用CW Commander航空成像應(yīng)用軟件進(jìn)行飛行規(guī)劃和監(jiān)測，以及采集航空照片，其中航空照片的航向、旁向重疊率均為80%，地面比例尺大多設(shè)置為1∶300(在地形復(fù)雜的破裂東端設(shè)置為1∶600).同時(shí)在研究區(qū)內(nèi)設(shè)置臨時(shí)基站并使用RTK (Real-time Kinematic)測量其地理坐標(biāo)，定位精度為厘米級.數(shù)據(jù)采集后，使用Photoscan軟件對UAV數(shù)據(jù)進(jìn)行整理，生成數(shù)字正射影像(DOM).影像的質(zhì)量取決于航拍區(qū)域地形、植被、作業(yè)時(shí)的風(fēng)向、飛行高度、相鄰照片的重疊度、特定區(qū)域照片數(shù)量等外界條件.尤其是在高原地區(qū)，地勢起伏大，地形復(fù)雜，天氣多變，不同飛行區(qū)域的無人機(jī)飛行高度不同，以上諸多因素導(dǎo)致影像分辨率有所差異(Liu-Zeng et al., 2022).最終生成的DOM(圖1)空間分辨率為3～7 cm/pix，包含紅、綠、藍(lán)三個(gè)光譜帶，其中研究區(qū)東端山地地形較為復(fù)雜，起伏度大，分辨率在5～7 cm/pix.之后使用ArcGIS10.2提取相應(yīng)的山影圖、坡度圖等，通過對各要素圖層進(jìn)行拉伸、濾波、變換等處理，以便更好的識別地表破裂，提升數(shù)據(jù)解譯能力.

圖2 瑪多地震地表破裂帶無人機(jī)圖像采集區(qū)Fig.2 UAV image acquisition area of the surface rupture zone of the Madoi earthquake

2.2 訓(xùn)練數(shù)據(jù)集準(zhǔn)備

深度學(xué)習(xí)中的監(jiān)督學(xué)習(xí)算法需要大量的標(biāo)簽數(shù)據(jù)來訓(xùn)練網(wǎng)絡(luò)，這些數(shù)據(jù)泛化性質(zhì)量的好壞決定著模型魯棒性的強(qiáng)弱.因此我們使用ArcGIS在DOM上對地表裂縫采用勾畫折線的形式進(jìn)行了大量的手工標(biāo)注，獲得鄂陵湖段西端A區(qū)和野馬灘北段B區(qū)(圖2)約40 km2的無人機(jī)圖像數(shù)據(jù)與之對應(yīng)的約800條裂縫標(biāo)注.為了減少遺漏、避免重復(fù)、方便復(fù)查，對整個(gè)測區(qū)采取統(tǒng)一的解譯標(biāo)志，將影像分為等大小格網(wǎng)，依次對每個(gè)格網(wǎng)采取1∶1000、1∶500、1∶250等比例尺顯示，保障最小裂縫識別長度<0.1 m，按照移動(dòng)窗口法拉網(wǎng)式、高密度人工解譯最終獲得沿總體走向長度約160 km、累積破裂長度超310 km的精細(xì)地表破裂數(shù)據(jù).解譯期間，各解譯人員采取統(tǒng)一的解譯標(biāo)志，并及時(shí)溝通.并于9月在野外對標(biāo)注數(shù)據(jù)進(jìn)行現(xiàn)場確認(rèn)以及調(diào)整，減少錯(cuò)標(biāo)誤標(biāo)，減輕了這些因素帶來的不利影響.數(shù)據(jù)現(xiàn)場確認(rèn)后，將獲取到的瑪多震后遙感影像與人工標(biāo)注信息清洗制作為PASCAL VOC2007格式的數(shù)據(jù)集.

PASCAL VOC2007數(shù)據(jù)集在檢測、分割等任務(wù)上催生了很多優(yōu)秀的計(jì)算機(jī)視覺模型.瑪多震后遙感影像的地表裂縫識別可以劃分為二分類的分割任務(wù)，一類為背景，另一類為地表裂縫.在計(jì)算機(jī)編程語言中，可以利用地理空間數(shù)據(jù)抽象庫(GDAL)讀取矢量數(shù)據(jù)中的標(biāo)注信息，根據(jù)標(biāo)注信息生成標(biāo)簽圖片(圖3).由于人工標(biāo)注是在ArcGIS軟件中采用勾畫折線段的形式展示，而非多邊形.因此利用標(biāo)注數(shù)據(jù)直接生成的標(biāo)簽是由連續(xù)的單點(diǎn)像素組成，較為狹窄，無法完全覆蓋地表裂縫.我們選擇利用形態(tài)學(xué)中的膨脹算法處理，對標(biāo)注的線膨脹40個(gè)像素元，把點(diǎn)序列形成的線轉(zhuǎn)成裂縫區(qū)域，這樣的處理方式可以使當(dāng)前不能識別的數(shù)據(jù)變成可以利用的數(shù)據(jù)集，同時(shí)可以使標(biāo)簽更好的擬合裂縫的延展方向與空間展布.但是利用膨脹算法處理可能會(huì)導(dǎo)致圖像裂縫中一些局部紋理細(xì)節(jié)丟失嚴(yán)重，從而影響某些像素點(diǎn)準(zhǔn)確定位.最后切割為(512×512)大小的標(biāo)簽圖片，形成2135份數(shù)據(jù)集，并把數(shù)據(jù)集按7∶2∶1的比例劃分為訓(xùn)練集、測試集、驗(yàn)證集.

圖3 地表裂縫數(shù)據(jù)集(a) 原圖； (b) 標(biāo)簽圖，為8位彩圖.Fig.3 Surface fracture data(a) Original drawing; (b) Label map, 8-bit color map.

2.3 使用卷積神經(jīng)網(wǎng)絡(luò)識別裂縫區(qū)域

2.3.1 卷積神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)

訓(xùn)練神經(jīng)網(wǎng)絡(luò)往往需要大量的訓(xùn)練集樣本，而我們的訓(xùn)練集僅僅只有約1495份，在訓(xùn)練集較少的情況下易產(chǎn)生過擬合的問題.因此我們使用卷積神經(jīng)網(wǎng)絡(luò)VGG-Unet作為識別裂縫區(qū)域的深度學(xué)習(xí)網(wǎng)絡(luò)，U-net采用了跳躍連接(skip-connection)結(jié)構(gòu)對編碼器和解碼器中的各級分辨率特征進(jìn)行連接，這種結(jié)構(gòu)能使模型更好的利用圖像的低級紋理特征，該方法曾在數(shù)據(jù)集樣本較少的醫(yī)療影像結(jié)構(gòu)分割任務(wù)中有很好的表現(xiàn).本文使用VGG-16圖像分類網(wǎng)絡(luò)的圖像特征提取模塊作為U-net的主干特征提取器.VGG-16是一個(gè)經(jīng)典的圖像分類網(wǎng)絡(luò)，在使用中舍棄VGG中的分類模塊，在保留了高效的特征提取能力的同時(shí)，避免了VGG龐大的分類模塊對整體性能的影響.同時(shí)考慮到裂縫提取是一個(gè)紋理特征提取任務(wù)，我們使用U-net對提取到的特征進(jìn)行解碼分類，U-net的優(yōu)勢在于能很好兼顧紋理特征和語義特征，對于紋理類的目標(biāo)識別有很好的效果(圖4).

圖4 基于VGG-16模型的U-net Encoder模塊Fig.4 U-net Encoder module based on VGG-16 model

首先，我們將切分后的RGB圖像輸入模型(圖5)，圖像經(jīng)由編碼器模塊(Encoder Module)的卷積層(Conv)濾波和池化層(Pooling)下采樣得到不同尺寸的特征圖，這些特征圖被保存在Features Map中，小尺寸的特征圖會(huì)在Decoder Module中經(jīng)過上采樣層轉(zhuǎn)換成與較大尺寸的特征圖相同的大小，然后與較大尺寸的特征圖進(jìn)行特征融合.這種做法可以讓網(wǎng)絡(luò)在小尺寸的特征圖中通過高級語義特征，以較大的感受野搜索裂縫區(qū)域，同時(shí)融合大尺寸的特征圖中的低級紋理特征，補(bǔ)全裂縫區(qū)域的細(xì)節(jié).最終，融合后的特征通過一個(gè)Conv 2d層進(jìn)行特征壓縮，將特征壓縮為兩個(gè)通道，用于分別表示每個(gè)像素點(diǎn)分別為裂縫或非裂縫的權(quán)重.

圖5 基于上采樣和skip-connection的U-net解碼器模塊Fig.5 U-net decoder module based on upsampling and skip-connection

2.3.2 評價(jià)標(biāo)準(zhǔn)

(1)

2.3.3 訓(xùn)練

預(yù)訓(xùn)練已經(jīng)在許多實(shí)驗(yàn)中被證明是一種有效提高訓(xùn)練結(jié)果穩(wěn)定性和加快網(wǎng)絡(luò)參數(shù)擬合的方法.對于Encoder部分，我們使用經(jīng)由ImageNet數(shù)據(jù)集訓(xùn)練的VGG網(wǎng)絡(luò)，在舍棄全連接層參數(shù)后作為遷移學(xué)習(xí)預(yù)訓(xùn)練網(wǎng)絡(luò)，而對于Decoder部分，我們選擇隨機(jī)初始化參數(shù)重新訓(xùn)練.為了防止隨機(jī)參數(shù)Decoder產(chǎn)生的較大LCE值使網(wǎng)絡(luò)反向傳播時(shí)對預(yù)訓(xùn)練網(wǎng)絡(luò)參數(shù)造成破壞，在訓(xùn)練的前50個(gè)Epoch，我們鎖定Encoder模塊的參數(shù)，訓(xùn)練Decoder使其獲得基本的識別能力以抑制LCE值.在該過程中，由于Adam優(yōu)化器具有計(jì)算高效、占內(nèi)存少、收斂快等優(yōu)點(diǎn)，我們使用Adam優(yōu)化器、批數(shù)據(jù)量8條、學(xué)習(xí)率0.0001和學(xué)習(xí)率下降率0.92作為訓(xùn)練參數(shù).在第50個(gè)Epoch后，我們解鎖Encoder模塊的參數(shù)，對網(wǎng)絡(luò)整體參數(shù)進(jìn)行微調(diào).此時(shí)由于Encoder需要參與網(wǎng)絡(luò)參數(shù)更新，導(dǎo)致網(wǎng)絡(luò)參數(shù)數(shù)據(jù)增加，同時(shí)因顯存容量有限，我們將批數(shù)據(jù)量改為4條，其余訓(xùn)練參數(shù)不變.

2.4 后處理

在通過深度學(xué)習(xí)網(wǎng)絡(luò)識別到裂縫區(qū)域后，為了將預(yù)測結(jié)果放進(jìn)ArcGIS中進(jìn)一步觀察裂縫分布與走向，可在識別到的裂縫區(qū)域內(nèi)對裂縫進(jìn)行近一步精細(xì)刻畫.注意到裂縫所在位置像素灰度值快速變化的特征，我們使用Canny算法在裂縫區(qū)域中提取裂縫實(shí)體，該方法主要分為以下四個(gè)步驟：(1) 首先利用高斯濾波，將圖像轉(zhuǎn)化為灰度圖并使用高斯算子對灰度圖像進(jìn)行濾波，以消除高斯噪聲，防止高斯噪聲在后續(xù)過程中被誤識別為裂縫. (2) 使用一階離散微分算子在濾波后的灰度圖像x和y方向分別計(jì)算近似梯度，用該值表示圖像在該點(diǎn)x和y方向上灰度值變化的劇烈程度.為了獲得清晰、明確的裂縫，我們需要找到圖像中灰度變化最快的方向，即近似梯度最大的方向.為此，需要通過x方向的近似梯度Gx，y方向上的近似梯度Gy，來計(jì)算梯度幅值G和梯度方向θ，計(jì)算公式如下：

(2)

(3)

在獲得幅值G和梯度方向θ后，將θ離散為8個(gè)方向，對應(yīng)每個(gè)像素周圍8個(gè)像素方向.并根據(jù)離散后θ指向的方向，找到該方向幅值G的最大值，最大值所在的像素點(diǎn)即判斷為候選裂縫.(3) 最后為了排除一些灰度變化不明顯的紋理噪聲，同時(shí)保證裂縫的連續(xù)和裂縫末端的完整，我們采用雙閾值的滯后閾值處理進(jìn)行較長邊緣裂縫的鏈接.通過對梯度幅值G設(shè)置高門限G1為120、低門限G2為80.將候選裂縫進(jìn)行檢測，其中檢測值大于G1的定為強(qiáng)裂縫；檢測值在G2與G1之間，且與確定為強(qiáng)裂縫的像素鄰接的定為弱裂縫；檢測值小于G1的定為非裂縫.最終將強(qiáng)裂縫、與強(qiáng)裂縫鄰接的弱裂縫作為識別結(jié)果.

為了能夠在ArcGIS軟件中觀察到裂縫識別結(jié)果的幾何位置信息，需要對裂縫識別結(jié)果的圖片數(shù)據(jù)進(jìn)行計(jì)算并轉(zhuǎn)化為文檔數(shù)據(jù)，最后結(jié)合投影文件信息，形成用點(diǎn)來描述裂縫幾何位置的ESRI Shapefile文件.

2.5 評價(jià)指標(biāo)

因?yàn)榱芽p數(shù)據(jù)是以矢量的方式保存，直接判斷裂縫識別的矢量結(jié)果缺乏客觀性，因此我們將矢量數(shù)據(jù)柵格化后膨脹的裂縫標(biāo)簽來驗(yàn)證預(yù)測裂縫的準(zhǔn)確度.我們將DOM上的裂縫標(biāo)注和裂縫預(yù)測結(jié)果分別進(jìn)行膨脹，兩者的膨脹結(jié)果分別視作真實(shí)的裂縫區(qū)域和預(yù)測的裂縫區(qū)域.我們根據(jù)預(yù)測結(jié)果和實(shí)際結(jié)果的一致性將識別結(jié)果以像素級分為以下四類.TP(真正例)：正確預(yù)測裂縫區(qū)域的像素?cái)?shù)；FP(假正例)：預(yù)測為裂縫區(qū)域但非真實(shí)裂縫的像素?cái)?shù)；TN(真反例)：正確預(yù)測非裂縫區(qū)域的像素?cái)?shù)；FN(假反例)：預(yù)測為非裂縫但實(shí)際是裂縫區(qū)域的像素?cái)?shù).在此基礎(chǔ)上，我們計(jì)算查準(zhǔn)率(Precision)和召回率(Recall)，前者為預(yù)測為裂縫區(qū)域中預(yù)測正確的比重，后者為真實(shí)裂縫區(qū)域中預(yù)測正確的比重.查準(zhǔn)率體現(xiàn)了模型捕獲裂縫區(qū)域的能力，召回率體現(xiàn)了模型捕獲裂縫區(qū)域的質(zhì)量.查準(zhǔn)率和召回率的計(jì)算方式如下：

(4)

(5)

為了減少非裂縫區(qū)域誤識別的同時(shí)盡可能獲取所有的裂縫區(qū)域，我們希望模型同時(shí)擁有較好的查準(zhǔn)率指標(biāo)和查全率指標(biāo)，因此我們使用Dice系數(shù)指標(biāo)作為最終全面評估識別效果的指標(biāo).Dice系數(shù)的計(jì)算方式如公式(6)所示：

(6)

在Dice系數(shù)計(jì)算公式中，由于準(zhǔn)確率和召回率取值為0～1，因此無論其中一個(gè)指標(biāo)有多好的表現(xiàn)，如果另一個(gè)指標(biāo)表現(xiàn)過差都會(huì)使Dice系數(shù)指標(biāo)變差.同時(shí)Dice系數(shù)指標(biāo)在分子乘上系數(shù)2，保證Dice系數(shù)指標(biāo)取值范圍在0～1之間.

這些評價(jià)指標(biāo)為我們的工作提供了重要的幫助，其中查準(zhǔn)率和召回率為我們調(diào)整訓(xùn)練方法和修改超參數(shù)提供了重要的方向指導(dǎo)，Dice評分為評估我們的方法結(jié)果的總體準(zhǔn)確性提供了一個(gè)定量的標(biāo)定.

3 結(jié)果

3.1 實(shí)驗(yàn)結(jié)果

本實(shí)驗(yàn)中需要對113 GB的數(shù)據(jù)進(jìn)行計(jì)算，因此要求實(shí)驗(yàn)設(shè)備具有強(qiáng)大的計(jì)算能力.在實(shí)驗(yàn)硬件搭配中，采用了1 TB固態(tài)硬盤儲(chǔ)存數(shù)據(jù)與RTX 3060 12 GB顯卡做計(jì)算力支持.我們采用隨機(jī)梯度下降的方式訓(xùn)練網(wǎng)絡(luò)，設(shè)置總Epoch為300.訓(xùn)練過程中平均每個(gè)Epoch耗時(shí)2分52秒，總共耗時(shí)12小時(shí)36分.為了衡量模型的性能，在每個(gè)Epoch分別計(jì)算出訓(xùn)練集小批次與驗(yàn)證集小批次的平均Dice系數(shù)指標(biāo)及損失值(圖6).

圖6 圖像分割評測指標(biāo)Dice、損失值變化曲線Fig.6 Image segmentation evaluation index Dice，LCE value curve

從圖6可以看出模型參數(shù)擬合良好，最優(yōu)性能模型對驗(yàn)證集進(jìn)行預(yù)測的Dice約為87.5%.選取在驗(yàn)證集上表現(xiàn)最優(yōu)性能的模型參數(shù)(第263 Epoch)，并在驗(yàn)證集上進(jìn)行識別將標(biāo)注區(qū)域與機(jī)器識別區(qū)域結(jié)果區(qū)域勾繪出來(圖7c)，緊接著對模型的識別結(jié)果進(jìn)行預(yù)測區(qū)域和地面實(shí)況區(qū)域之間相似性數(shù)據(jù)統(tǒng)計(jì)(圖7d)，其中TP的像素點(diǎn)個(gè)數(shù)為79430，F(xiàn)N的像素點(diǎn)個(gè)數(shù)為3292，F(xiàn)P的像素點(diǎn)個(gè)數(shù)為13051，由此計(jì)算出Precision=0.85、Recall=0.96、Dice=0.9.

圖7 模型識別結(jié)果(a) 數(shù)字正射影像； (b) 人工識別地表破裂，其中紅線為人工標(biāo)注的地表破裂； (c) 標(biāo)注區(qū)域與機(jī)器識別區(qū)域?qū)Ρ?，其中紅線為標(biāo)簽標(biāo)注區(qū)域， 藍(lán)線為機(jī)器識別區(qū)域； (d) 機(jī)器識別結(jié)果，其中綠色區(qū)域?yàn)檎嬲?TP)、紅色區(qū)域?yàn)榧俜蠢?FN)、藍(lán)色區(qū)域?yàn)榧僬?FP)．Fig.7 Model recognition results(a) DOM image； (b) Artificially identified surface rupture map，where the red line is the artificially marked surface rupture; (c) Comparison of marked area and machine recognition area, where the red line is the labeling area, and the blue line is the machine recognition area; (d) Machine recognition result map,where the green area is TP, the red area is FN, the blue area is FP.

按照上述識別率統(tǒng)計(jì)分析方法，隨機(jī)在不同區(qū)域的驗(yàn)證集中抽取10張圖進(jìn)行識別與數(shù)據(jù)統(tǒng)計(jì)，結(jié)果如表1所示，各驗(yàn)證集的Dice系數(shù)在0.7～0.91之間，這也說明了我們的結(jié)果魯棒性較好.與DOM上的裂縫進(jìn)行對比(圖7)表明，大多數(shù)裂縫能夠很好地被識別.

表1 驗(yàn)證集樣例識別數(shù)據(jù)統(tǒng)計(jì)Table 1 Verification set sample identification data statistics

3.2 測試結(jié)果

選取B-2(19839×16384)的測試集區(qū)域(即未在訓(xùn)練集與驗(yàn)證集中出現(xiàn)的區(qū)域)進(jìn)行地震裂縫識別，為了觀察模型效果，對B-2區(qū)域進(jìn)行了人工標(biāo)注，并將識別結(jié)果與人工標(biāo)注進(jìn)行對比計(jì)算，得出Precision、Recall、Dice計(jì)算結(jié)果如圖8.

圖8 B-2區(qū)域識別結(jié)果數(shù)據(jù)Fig.8 B-2 test area identification map

根據(jù)圖8中的數(shù)據(jù)，該方法對B-2區(qū)域的識別結(jié)果為Precision=0.417、Recall=0.663、Dice=0.512.為了更好地觀察識別效果與地震裂縫的識別位置，將識別結(jié)果按不同顏色進(jìn)行展現(xiàn)，可以更直觀地觀察混肴矩陣中TP、TN、FP、FN各個(gè)樣例的分布(圖9).由于誤差不可避免，在這種分辨率高的圖片進(jìn)行人工標(biāo)注工作時(shí)，會(huì)產(chǎn)生漏標(biāo)、錯(cuò)標(biāo)等錯(cuò)誤因素，從而影響數(shù)據(jù)統(tǒng)計(jì)結(jié)果，使Precision、Recall、Dice的計(jì)算值變小.從圖9(a、b)中可以看出藍(lán)色FP區(qū)域?yàn)槿斯ぢ?biāo)區(qū)域，而模型預(yù)測結(jié)果為裂縫，說明模型具有一定的泛化性能.圖9(c、d)中紅色FN區(qū)域的裂縫紋理雜亂呈斑塊狀，模型錯(cuò)誤預(yù)測為背景.因此，性能優(yōu)良的模型依賴于良好數(shù)據(jù)集的建立，需要盡可能多、地質(zhì)紋理種類多、標(biāo)簽正確的數(shù)據(jù)集.盡管我們期望模型在其他數(shù)據(jù)中的表現(xiàn)如同訓(xùn)練集上一樣良好，但是顯然，模型只能盡可能地趨近于人類期望結(jié)果，我們依舊可以根據(jù)檢測結(jié)果，觀察裂縫分布.

圖9 B-2測試區(qū)域標(biāo)簽圖與識別結(jié)果對比圖(a,b,c,d) 分別為4個(gè)不同的區(qū)域，紅線為標(biāo)簽標(biāo)注區(qū)域，綠色區(qū)域?yàn)檎嬲?TP),紅色區(qū)域?yàn)榧俜蠢?FN),藍(lán)色區(qū)域?yàn)榧僬?FP) ．Fig.9 B-2 test area label map and identification map comparison chart(a,b,c,d) are four different areas , the red line is the labeling area, the green area is TP, the red area is FN, and the blue area is FP.

4 討論

4.1 假反例分析

召回率(Recall)的大小受真正例(TP)、假反例(FN)的影響，其主要因素為真實(shí)標(biāo)簽的質(zhì)量，此外高分辨率無人機(jī)影像每張約1 GB的儲(chǔ)存容量，圖片較大使得標(biāo)注工作量也十分巨大，這需要我們在進(jìn)行人工標(biāo)注地表破裂時(shí)更加細(xì)致，并對標(biāo)注進(jìn)行清洗.

我們對A-4區(qū)選取的數(shù)據(jù)集識別(圖10)，從圖10a中可看出FN(假反例)區(qū)域主要覆蓋在一些紋理細(xì)小、長度短的裂縫上，圖10b中FN(假反例)區(qū)域并沒有覆蓋在真實(shí)裂縫上，這與人工標(biāo)注時(shí)將非裂縫標(biāo)記為裂縫所導(dǎo)致的誤差有關(guān)，圖10c中FN(假反例)區(qū)域主要覆蓋在一些紋理線路雜亂、呈斑塊狀的裂縫上.因此，裂縫紋理特征不明顯或聚集成斑塊狀、真實(shí)標(biāo)簽標(biāo)注不細(xì)致等因素，模型將裂縫預(yù)測為背景的概率會(huì)變大，從而導(dǎo)致FN數(shù)量增多，召回率偏低.

圖10 A-4測試區(qū)域標(biāo)簽圖與識別結(jié)果對比圖(a,b,c) 分別為3個(gè)不同的區(qū)域，紅線為標(biāo)簽標(biāo)注區(qū)域，綠色區(qū)域?yàn)檎嬲?TP),紅色區(qū)域?yàn)榧俜蠢?FN),藍(lán)色區(qū)域?yàn)榧僬?FP) ．Fig.10 A-4 test area label map and identification map comparison chart(a,b,c) are three different areas, the red line is the labeling area, the green area is TP, the red area is FN, and the blue area is FP.

4.2 假正例分析

地表破裂的精確率(Precision)與假反例數(shù)量主要受到訓(xùn)練集樣本影響，為了使模型具有更好的準(zhǔn)確性與泛化性，往往需要準(zhǔn)備盡量多、且包含不同紋理地形的訓(xùn)練數(shù)據(jù)集.相反，當(dāng)訓(xùn)練數(shù)據(jù)集樣例較少或樣本不均衡時(shí)，模型不能學(xué)習(xí)到最佳的參數(shù)并獲得更好的泛化能力，會(huì)使得假正例預(yù)測結(jié)果增多.

從圖11a中可看到FP(假正例)區(qū)域覆蓋在標(biāo)簽的周圍，模型預(yù)測區(qū)域比標(biāo)簽區(qū)域廣，該種情況可視為模型預(yù)測正確.而圖11b中可看出FP(假正例)區(qū)域覆蓋在真正裂縫上，這是由于人工漏標(biāo)導(dǎo)致的誤差，可通過細(xì)致勾繪標(biāo)簽減少誤差.圖11c中可看出FP(假正例)區(qū)域覆蓋在類似于對比圖中的裂縫紋理上且未在真實(shí)地裂縫周圍，這種情況視為模型錯(cuò)誤預(yù)測，可在訓(xùn)練數(shù)據(jù)集中放入一批類似于裂縫紋理的圖片，標(biāo)簽設(shè)為全背景進(jìn)行訓(xùn)練.顯然，上述3種情況都是影響精確率數(shù)值變低與假正例樣例增多的主要因素.

圖11 B-3測試區(qū)域標(biāo)簽圖與識別結(jié)果對比圖(a,b,c) 分別為3個(gè)不同的區(qū)域，紅線為標(biāo)簽標(biāo)注區(qū)域，綠色區(qū)域?yàn)檎嬲?TP),紅色區(qū)域?yàn)榧俜蠢?FN),藍(lán)色區(qū)域?yàn)榧僬?FP) ．Fig.11 B-3 test area label map and identification map comparison chart(a,b,c) are three different areas, the red line is the labeling area, the green area is TP, the red area is FN, and the blue area is FP.

4.3 深度學(xué)習(xí)在地表破裂識別中的優(yōu)勢與局限

同震地表破裂帶的精細(xì)填圖可以為理解斷裂帶復(fù)雜幾何結(jié)構(gòu)、動(dòng)態(tài)破裂過程和古地震復(fù)發(fā)習(xí)性等提供重要的基礎(chǔ)信息(Oskin et al., 2012；Klinger et al., 2018; Choi et al., 2018; 劉小利等，2022；姚文倩等，2022).而隨著海量高精度地形數(shù)據(jù)的便捷化日益獲取，如何采取更高效可靠的(半)自動(dòng)化定量處理和分析也逐漸引發(fā)了眾多科研人員的關(guān)注和投入(Palamara et al., 2007；Walker et al., 2016；姚文倩等，2019；韓龍飛等，2019；Mattéo et al.，2021).相較于傳統(tǒng)的人工識別和專家識別系統(tǒng)而言，本次研究基于深度學(xué)習(xí)的算法對瑪多地震同震破裂實(shí)現(xiàn)了自動(dòng)識別和提取，拓展了利用深度學(xué)習(xí)在同震地表破裂研究領(lǐng)域上的廣泛應(yīng)用.總體而言，深度學(xué)習(xí)在構(gòu)造地貌自動(dòng)識別研究中具備以下優(yōu)勢：(1)相較于傳統(tǒng)的人工識別方法，基于深度學(xué)習(xí)的自動(dòng)識別技術(shù)可以提升數(shù)十倍效率，且可以降低標(biāo)注過程中漏檢誤檢的概率.(2)深度學(xué)習(xí)的自動(dòng)識別結(jié)果是基于嚴(yán)格統(tǒng)一的判別標(biāo)準(zhǔn)，相對于人工識別更多依賴于判讀者的主觀認(rèn)知來說，后者難以標(biāo)準(zhǔn)一致，因而深度學(xué)習(xí)的結(jié)果具有較為穩(wěn)定的識別標(biāo)準(zhǔn)和和更為優(yōu)異的質(zhì)量.(3)傳統(tǒng)手動(dòng)識別中要求研究人員具有一定的專業(yè)背景知識，相對而言，深度學(xué)習(xí)網(wǎng)絡(luò)可以具有較高的自主學(xué)習(xí)能力，因此不需要過多強(qiáng)調(diào)研究人員的圖像知識、地質(zhì)知識等來設(shè)計(jì)優(yōu)秀的特征模式.

然而需要注意的是，目前深度學(xué)習(xí)在地學(xué)領(lǐng)域的運(yùn)用仍存在一些問題，主要體現(xiàn)在：(1)深度學(xué)習(xí)網(wǎng)絡(luò)的推廣性尚且有限.深度學(xué)習(xí)網(wǎng)絡(luò)的識別能力是基于對已有數(shù)據(jù)的學(xué)習(xí)，當(dāng)構(gòu)造地貌類型發(fā)生其他變化時(shí)，需要重新收集相關(guān)的高精度地形數(shù)據(jù)對網(wǎng)絡(luò)進(jìn)行二次訓(xùn)練，需要追加一定的成本.(2)盡管深度學(xué)習(xí)網(wǎng)絡(luò)總體上可以實(shí)現(xiàn)較高準(zhǔn)確度的識別，但在較小的粒度上仍有不同程度的誤識別和漏識別現(xiàn)象.因此，目前深度學(xué)習(xí)更適用于宏觀構(gòu)造地貌的精細(xì)刻畫.為了克服以上問題，未來我們需要在以下幾個(gè)方面進(jìn)行改進(jìn)和創(chuàng)新：(1)盡可能采用最先進(jìn)的算法，嘗試突破已經(jīng)訓(xùn)練好的模型的應(yīng)用局限性；(2)繼續(xù)改進(jìn)算法，利用大型圖像處理器更加快速地處理覆蓋較大區(qū)域的數(shù)據(jù)集，減少運(yùn)算時(shí)間.(3)結(jié)合遞歸神經(jīng)網(wǎng)絡(luò)和CNN，利用多波段信息，達(dá)到更準(zhǔn)確識別的效果.

4.4 展望

深度學(xué)習(xí)方法也可拓展到其他地貌和地表過程等定量研究中.近年來，移動(dòng)攝影測量技術(shù)(SfM)以其方便快捷、低成本獲得數(shù)公里到數(shù)十公里范圍的厘米級三維地形數(shù)據(jù)和光學(xué)影像(魏占玉等, 2015; 畢海蕓等, 2017)，激光雷達(dá)掃描(LiDAR)能快速、高精度地獲得更大范圍的三維地形數(shù)據(jù)，這些都極大地推進(jìn)了活動(dòng)構(gòu)造和構(gòu)造地貌的量化研究(Hudnut et al., 2002; Gold et al., 2009; Cowgill et al., 2012; Oskin et al., 2012; 劉靜等, 2013; 任治坤等, 2014).隨著大規(guī)模LiDAR和SfM技術(shù)的應(yīng)用，海量密集地貌點(diǎn)云和影像數(shù)據(jù)獲得越來越便捷，數(shù)據(jù)呈現(xiàn)指數(shù)級增長.對于海量高精度地形描述數(shù)據(jù)的有效挖掘，一方面在三維數(shù)據(jù)可視化方面嘗試應(yīng)用虛擬現(xiàn)實(shí)技術(shù)研究目標(biāo)區(qū)域的構(gòu)造和地貌特征(Cowgill et al., 2012; Fowler et al., 2012; 王鵬等，2020)，另一方面基于高分辨率地形數(shù)據(jù)對地表地貌和地質(zhì)巖性的識別也逐漸從定性描述向定量分析轉(zhuǎn)變.

以往研究人員通常以感性差異的裸眼經(jīng)驗(yàn)進(jìn)行人工地貌解譯，或以簡單的計(jì)算方法和指標(biāo)(如地形起伏度和粗糙度)反映不同期次的地貌單元發(fā)育特征 (韓龍飛等，2019; Frankel and Dolan, 2007; Regmi et al., 2014)，進(jìn)行不同地貌面特征的參數(shù)表征和自動(dòng)化提取后，開展半人工填圖.機(jī)器學(xué)習(xí)為基于高精度海量數(shù)據(jù)的大范圍地貌定量研究提供新的機(jī)會(huì).深度學(xué)習(xí)網(wǎng)絡(luò)更容易更快速地將遙感影像和DEM、坡度等其他多類型的數(shù)據(jù)集進(jìn)行有效地挖掘，提高地貌識別性能，以幫助更好地了解地表過程(Du et al., 2019).

總之，隨著算法和模型不斷的改進(jìn)和創(chuàng)新，機(jī)器學(xué)習(xí)會(huì)在地貌和地表過程的研究中有著更廣泛的運(yùn)用.機(jī)器學(xué)習(xí)技術(shù)的應(yīng)用近年來已經(jīng)涉及地學(xué)的各個(gè)領(lǐng)域，為地球深淺部過程研究提供了直接的手段，開辟了新視野.未來在算法、精度、應(yīng)用場景和批處理便利性等方面還有更大的發(fā)展空間.

5 結(jié)論

我們在MW7.4 瑪多地震發(fā)生后第一時(shí)間沿整個(gè)地表破裂帶進(jìn)行了高分辨率無人機(jī)影像數(shù)據(jù)采集，獲取的數(shù)據(jù)集能夠很好的識別地表破裂，為探索地震發(fā)生機(jī)理、斷層破裂過程及運(yùn)動(dòng)方式提供基礎(chǔ)信息.本文提出了一種基于機(jī)器學(xué)習(xí)算法，從高分辨率無人機(jī)影像數(shù)據(jù)中自動(dòng)檢測同震地表破裂的能力，采用卷積神經(jīng)網(wǎng)絡(luò)(CNN)的有監(jiān)督深度學(xué)習(xí)算法構(gòu)建了地表破裂識別模型，并應(yīng)用于鄂陵湖段和野馬灘北段約40 km2的無人機(jī)航片數(shù)據(jù)集，快速、高效地繪制了該段的地表破裂.盡管存在一些誤報(bào)和漏識別，但與人工解譯結(jié)果對比表明，深度學(xué)習(xí)算法與自動(dòng)化處理系統(tǒng)相結(jié)合能夠有效地繪制地震所產(chǎn)生的地表破裂，為震后分析地震災(zāi)情、評估地震斷層危害性等工作提供重要依據(jù)，也為未來研究大地震地表破裂提供新思路.本次采用的卷積神經(jīng)網(wǎng)絡(luò)模型為后續(xù)利用深度學(xué)習(xí)進(jìn)行地表破裂的快速自動(dòng)識別提供了一種可能的優(yōu)化策略.同時(shí)展示了機(jī)器學(xué)習(xí)在地震地質(zhì)、地表過程和地貌等定量研究中的巨大優(yōu)勢和廣闊前景.

致謝感謝評審專家以及編委對于本文提出的問題和建議；感謝青海省地震局、成都縱橫自動(dòng)化技術(shù)有限公司在數(shù)據(jù)采集過程中給予的支持和幫助.在此一并表示感謝！