基于海底DEM的洋中脊火山錐自動(dòng)識(shí)別方法研究

黨 牛，余 星，韓喜球，陳安清

(1.成都理工大學(xué) 沉積地質(zhì)研究院,四川 成都 610059;2.自然資源部 海底科學(xué)重點(diǎn)實(shí)驗(yàn)室,浙江 杭州 310012;3.自然資源部 第二海洋研究所，浙江 杭州 310012)

0 引言

火山錐是中心式巖漿噴發(fā)的產(chǎn)物，廣泛分布于陸內(nèi)、島弧和洋盆等各個(gè)火山作用區(qū)域，錐體尺度一般為數(shù)十米至數(shù)十千米[1-4]。對(duì)于板塊擴(kuò)張邊界洋中脊，線性裂隙式噴發(fā)是最主要的巖漿產(chǎn)出形式。不過(guò)，從洋中脊的海底地形圖上常?？梢杂^察到錐狀火山的痕跡，直徑50～1 000 m。這些火山錐的形態(tài)特征、空間分布和化學(xué)成分與區(qū)域的構(gòu)造應(yīng)力、熔融深度、巖漿供應(yīng)和地殼厚度等因素密切相關(guān)[5-10]。通過(guò)認(rèn)識(shí)這些火山錐的空間分布信息，可以深入了解洋中脊巖漿作用的過(guò)程與特征。

火山錐識(shí)別通常采用人工解譯和機(jī)器自動(dòng)識(shí)別兩種方法[11-13]?；鹕藉F的形態(tài)在視覺(jué)上比較明顯，易于人工識(shí)別。人工提取方法的精度和可信度高，但耗時(shí)費(fèi)力、主觀性強(qiáng)、可重復(fù)性差[14-17]。機(jī)器自動(dòng)提取則效率高、成本低、可重復(fù)。隨著海洋調(diào)查的深入，數(shù)據(jù)規(guī)模不斷擴(kuò)大，人工解譯耗時(shí)費(fèi)力的局限性更加突顯，對(duì)海底多源數(shù)據(jù)的高效分析亟需自動(dòng)提取方法的有效應(yīng)用。

自動(dòng)提取方法主要包括閾值分類(lèi)、監(jiān)督分類(lèi)和非監(jiān)督分類(lèi)三種[18-20]。閾值分類(lèi)基于所提取要素的屬性臨界值進(jìn)行要素判別，主要適用于空間特征明顯的線性構(gòu)造識(shí)別[21-22]。監(jiān)督分類(lèi)與非監(jiān)督分類(lèi)通過(guò)建立判別函數(shù)來(lái)進(jìn)行圖像識(shí)別[23-26]。監(jiān)督分類(lèi)需要通過(guò)樣本訓(xùn)練來(lái)獲取先驗(yàn)知識(shí)，精度最高可達(dá)90%以上[27]，但受制于訓(xùn)練場(chǎng)地的選擇和先驗(yàn)知識(shí)的建立，多適用于多光譜遙感影像等多波段數(shù)據(jù)。非監(jiān)督分類(lèi)則不需要先驗(yàn)知識(shí)，只依據(jù)圖像集聚特征進(jìn)行分類(lèi)，精度一般低于監(jiān)督分類(lèi)，因而應(yīng)用較少。但非監(jiān)督分類(lèi)效率高，對(duì)于具有明顯形態(tài)特征的地質(zhì)要素，識(shí)別精度可達(dá)60%～90%[25,28]。本文嘗試基于卡爾斯伯格脊DEM(Digital Elevation Model,數(shù)字高程模型)數(shù)據(jù)，使用非監(jiān)督分類(lèi)方法進(jìn)行洋中脊火山錐的自動(dòng)識(shí)別研究，具有一定的開(kāi)創(chuàng)性。

1 地質(zhì)背景及數(shù)據(jù)來(lái)源

卡爾斯伯格脊是位于印度洋西北角的一小段洋中脊，全長(zhǎng)約1 400 km，全擴(kuò)張速率為24～26 mm/a，屬于典型的慢速擴(kuò)張洋中脊[29-33]。北部邊界為歐文斷裂帶，與亞丁灣的希巴洋脊分隔；南部邊界為酈道元斷裂帶，與中印度洋中脊分隔。卡爾斯伯格脊從北到南共分為5個(gè)一級(jí)洋脊段，本研究選取第4段61°E—65°E洋脊段為自動(dòng)識(shí)別方法的試驗(yàn)對(duì)象。該段洋脊為典型的對(duì)稱(chēng)擴(kuò)張洋脊段，軸部水深較淺，兩側(cè)地形起伏較小，且分布有大量的火山錐(圖1)。海底地形數(shù)據(jù)來(lái)源于中國(guó)大洋24航次的船載多波束測(cè)深資料，利用Teledyne Caris HIPS和SIPS軟件將原始多波束數(shù)據(jù)生成連續(xù)平滑的曲面，使用克里金插值方法生成分辨率約為55 m的DEM數(shù)據(jù)[34-35]。

圖1 研究區(qū)位置及典型火山錐地形特征Fig.1 Location of the study area and topographical features of typical volcanic cones(圖a為研究區(qū)位置，圖b為研究區(qū)地形特征，圖c、d、e為典型火山錐三維地形特征。紅色線指示擴(kuò)張中心，黑色線指示斷層，藍(lán)色方框?yàn)閳D4a范圍。)(Fig.a shows the location of the study area.Fig.b shows the topographical features of the study area.Fig.c,d,and e show the 3D terrains of typical volcanic cones.The red line indicates the spreading center.The black lines indicate the faults.The blue square shows the location of Fig.4a.)

2 自動(dòng)識(shí)別方法

計(jì)算機(jī)直接基于原始DEM圖像開(kāi)展自動(dòng)識(shí)別難度較大，需要對(duì)原始圖像進(jìn)行特征變換，使火山錐具有機(jī)器可讀的更加明顯的空間特征。本研究使用ISO(Iterative Self Organizing,迭代自組織)聚類(lèi)非監(jiān)督分類(lèi)方法對(duì)特征變換圖像進(jìn)行自動(dòng)聚類(lèi)，選擇包含火山錐的聚類(lèi)進(jìn)行面積與形狀篩選，完成火山錐自動(dòng)提取。具體思路及處理流程如圖2 所示。

圖2 基于DEM提取火山錐的流程圖Fig.2 The flow chart of volcanic cone recognition based on DEM data

2.1 圖像特征變換

研究區(qū)火山錐多呈圓形或近圓形，半徑400～1 000 m，高度約120～250 m。火山錐底部往上坡度顯著增大，而在坡頂坡度明顯減小，火山錐的邊界對(duì)應(yīng)坡腳坡度開(kāi)始增大的位置。地表粗糙度反映地表的起伏，是地表實(shí)際面積與投影面積的比值。從火山錐底部至頂部地表粗糙度先增大后減小，與坡度正相關(guān)。火山錐的邊界也對(duì)應(yīng)坡腳地表粗糙度開(kāi)始增加的位置(圖3)。

圖3 研究區(qū)坡度與地表粗糙度計(jì)算結(jié)果Fig.3 The calculated slope and surface roughness for the study area

火山錐像元值在原始圖像中過(guò)于分散，不利于機(jī)器自動(dòng)提取。通過(guò)對(duì)原始DEM數(shù)據(jù)進(jìn)行數(shù)學(xué)運(yùn)算，使火山錐空間信息產(chǎn)生集聚效果，像元值集中于一定區(qū)間，使計(jì)算機(jī)聚類(lèi)時(shí)與其他地物類(lèi)別有效區(qū)分。特征變換公式如下

(1)

(2)

式中：xz為特征變換后像元值，x為原始像元值，xrough為地表粗糙度圖像像元值，xmean為滑動(dòng)窗口方法求取的平均像元值；n為滑動(dòng)窗口大小，多為奇數(shù)。

由于火山錐一般為近圓形目標(biāo)，適合采用圓形滑動(dòng)窗口。取平滑窗口為15個(gè)像元對(duì)水深圖像進(jìn)行平滑[22,36-37]，將原始圖像與平滑圖像相減生成正負(fù)地形[38]，使區(qū)域凸起地形與低洼地形分離，從而所有的火山錐均具有正值或?qū)?yīng)凸起區(qū)域。將正負(fù)地形與地表粗糙度相乘，縮小火山邊緣與坡頂之間像元值的差距，增強(qiáng)火山錐邊緣特征，并保持火山錐邊緣的完整性(圖4)。然后移除不包含火山錐的無(wú)效區(qū)域，完成圖像特征變換，火山錐對(duì)應(yīng)像素高值區(qū)域，并且特征變換后火山錐的聚類(lèi)特征明顯增強(qiáng)，形態(tài)保持完好(圖5)。

圖5 研究區(qū)特征變換后的圖像(a)以及特征變換前后的典型火山錐剖面圖(b和c)Fig.5 Image after feature transformation(a),profiles of typical volcanic cones before and after feature transformation(b,c)

2.2 ISO聚類(lèi)非監(jiān)督分類(lèi)

完成圖像特征變換后，使用ISO聚類(lèi)非監(jiān)督分類(lèi)方法對(duì)圖像進(jìn)行自動(dòng)分類(lèi)[39-41]。ISO聚類(lèi)非監(jiān)督分類(lèi)方法包括ISODATA算法(迭代自組織數(shù)據(jù)分析算法)與最大似然法分類(lèi)兩部分。ISODATA算法基于多元分析屬性空間中的像元值來(lái)計(jì)算聚類(lèi)中心，通過(guò)遷移平均值法，多次迭代計(jì)算各個(gè)聚類(lèi)的平均值，得到新的聚類(lèi)中心，直到完成迭代次數(shù)或兩次迭代差異小于2%[42]，獲得最終聚類(lèi)結(jié)果。具體算法如下[43]。

(1)任意選擇初始聚類(lèi)中心，使用公式(3)在聚類(lèi)中心間分配像素：

x∈iif|ω(x)-ωi|<|ω(x)-ωj|for allj≠i

(3)

式中：ωi和ωj分別是聚類(lèi)i和j的聚類(lèi)中心，ω(x)是像素x的特征向量。

(2)通過(guò)遷移平均值法計(jì)算新的聚類(lèi)中心：

(4)

式中：K為聚類(lèi)個(gè)數(shù)，Qi是類(lèi)別i包含的總像素個(gè)數(shù)。

(3)像素被分配到鄰近的聚類(lèi)，并計(jì)算新聚類(lèi)的平均值和協(xié)方差。

(4)兩次迭代差異小于2%或完成迭代次數(shù)結(jié)束運(yùn)算，否則重復(fù)第三步與第四步。

通過(guò)ISO聚類(lèi)可得到像元與聚類(lèi)之間的所屬關(guān)系、聚類(lèi)的平均值、方差以及協(xié)方差等多元統(tǒng)計(jì)信息。依據(jù)ISO聚類(lèi)統(tǒng)計(jì)特征文件，使用最大似然法分類(lèi)將圖像分為若干個(gè)聚類(lèi)，創(chuàng)建分類(lèi)圖像。保留自動(dòng)分類(lèi)圖像中包含火山錐的聚類(lèi)作為有效區(qū)域，進(jìn)行下一步提取分析(圖6)。

圖6 ISO聚類(lèi)非監(jiān)督分類(lèi)圖像中包含火山錐的有效聚類(lèi)Fig.6 The effective clusters containing volcanic cones within the ISO clustering unsupervised classification image

2.3 景觀形狀指數(shù)篩選

將ISO聚類(lèi)非監(jiān)督分類(lèi)的有效聚類(lèi)轉(zhuǎn)化為shape格式的圖斑，這些圖斑不僅包含了需要提取的火山錐，也包括同為高值聚類(lèi)的斷層、山脊等干擾地貌類(lèi)型。根據(jù)火山錐的獨(dú)特形貌特征，對(duì)圖斑進(jìn)行面積與形狀篩選，可以得到自動(dòng)識(shí)別的火山錐。通過(guò)計(jì)算各圖斑的周長(zhǎng)、面積等形態(tài)參數(shù)，利用景觀形狀指數(shù)公式進(jìn)行形狀篩選[44-45]：

(5)

式中：E為斑塊周長(zhǎng)；A為斑塊面積；LSI為景觀形狀指數(shù)(圓形)，景觀生態(tài)學(xué)中常用于計(jì)算區(qū)域內(nèi)某斑塊形狀與相同面積的圓的偏離程度，LSI值為1時(shí)，代表標(biāo)準(zhǔn)的圓形，偏離程度越大LSI值越大。

由于區(qū)域內(nèi)火山錐半徑多為0.4～1.0 km，將斑塊面積為0.35～4.00 km2的區(qū)域視為有效區(qū)域(圖7a)。滿(mǎn)足面積條件后，依據(jù)火山錐形狀特征，選取景觀形狀指數(shù)小于1.25的斑塊為最終的火山錐區(qū)域(圖7b)。

圖7 面積特征提取結(jié)果(a)與形狀特征提取結(jié)果(b)Fig.7 Features extracted from size screening(a)and shape screening(b)

3 自動(dòng)識(shí)別結(jié)果與評(píng)價(jià)

通過(guò)ISO聚類(lèi)非監(jiān)督分類(lèi)方法和景觀形態(tài)指數(shù)篩選，共自動(dòng)識(shí)別火山錐25處。而利用研究區(qū)DEM數(shù)據(jù)，基于山地陰影、坡度、地表曲率、地表粗糙度等衍生參數(shù)進(jìn)行人工目視解譯，共提取28處火山錐(圖8a)，兩者存在一定差異。

圖8 自動(dòng)提取結(jié)果與人工提取結(jié)果對(duì)比Fig.8 The comparison of automatic extraction and manual extraction results

評(píng)價(jià)機(jī)器自動(dòng)識(shí)別的效果可以參考真陽(yáng)性、假陽(yáng)性、偽陰性三個(gè)指標(biāo)[36]。真陽(yáng)性是指提取結(jié)果與實(shí)地相符，即正確識(shí)別；假陽(yáng)性是指提取結(jié)果與實(shí)地不符，即錯(cuò)誤識(shí)別；偽陰性是指實(shí)際存在而未被識(shí)別，即遺漏識(shí)別。最后通過(guò)正確率和重疊率評(píng)價(jià)自動(dòng)識(shí)別效果，正確率和重疊率計(jì)算公式如下

(6)

(7)

對(duì)比自動(dòng)識(shí)別與人工提取結(jié)果，發(fā)現(xiàn)機(jī)器識(shí)別結(jié)果呈真陽(yáng)性20處，假陽(yáng)性5處，偽陰性8處。機(jī)器識(shí)別正確率達(dá)到0.8，重疊率約0.7，精度較好。

Kappa系數(shù)是衡量分類(lèi)效果和一致性檢驗(yàn)的常用指標(biāo)，系數(shù)越高代表一致性越好，數(shù)據(jù)范圍一般為0～1[23]。對(duì)比自動(dòng)提取與人工提取的火山錐位置、形態(tài)等參數(shù)，發(fā)現(xiàn)兩者的Kappa值為0.996，一致性很好。自動(dòng)提取與人工提取的火山錐圖斑的面積重疊率為0.816，提取的幾何形態(tài)也較吻合。

4 結(jié)論與展望

使用ISO聚類(lèi)非監(jiān)督分類(lèi)方法對(duì)特征變換圖像進(jìn)行聚類(lèi)分析，結(jié)合形態(tài)特征進(jìn)行進(jìn)一步篩選，完成火山錐自動(dòng)提取，可以達(dá)到較為滿(mǎn)意的識(shí)別效果。提取的火山錐三維形態(tài)信息可用于定量統(tǒng)計(jì)分析，有助于深入了解洋中脊構(gòu)造-巖漿活動(dòng)的細(xì)節(jié)過(guò)程。非監(jiān)督分類(lèi)適用于大面積地形數(shù)據(jù)的快速處理和特征信息提取，在自動(dòng)提取獨(dú)立分布且形態(tài)完整的火山錐時(shí)效果較好。由于洋中脊地形復(fù)雜，火山錐常常與其他線狀火山脊交疊，影響了自動(dòng)提取的效果。本次研究中遺漏識(shí)別8處，錯(cuò)誤識(shí)別5處，影響了整體的識(shí)別正確率。

與人工解譯相比，雖然目前機(jī)器自動(dòng)識(shí)別的精確度仍有待進(jìn)一步提高，但其高效、易重復(fù)的特點(diǎn)，將使其在大批量地形數(shù)據(jù)的快速自動(dòng)識(shí)別方面發(fā)揮重要作用。未來(lái)通過(guò)大量樣本的先驗(yàn)訓(xùn)練進(jìn)行監(jiān)督分類(lèi)，將極大地提高識(shí)別的準(zhǔn)確度。將非監(jiān)督分類(lèi)和監(jiān)督分類(lèi)相結(jié)合，有望實(shí)現(xiàn)不同效率、不同精度要求的機(jī)器自動(dòng)識(shí)別。

致謝感謝大洋24航次“李四光號(hào)”科考船的全體工作人員對(duì)地形數(shù)據(jù)采集提供的幫助，感謝審稿人的寶貴建議。