基于海底DEM的洋中脊火山錐自動識別方法研究

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

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

0 引言

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

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

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

1 地質背景及數據來源

卡爾斯伯格脊是位于印度洋西北角的一小段洋中脊，全長約1 400 km，全擴張速率為24～26 mm/a，屬于典型的慢速擴張洋中脊[29-33]。北部邊界為歐文斷裂帶，與亞丁灣的希巴洋脊分隔；南部邊界為酈道元斷裂帶，與中印度洋中脊分隔?？査共窦箯谋钡侥瞎卜譃?個一級洋脊段，本研究選取第4段61°E—65°E洋脊段為自動識別方法的試驗對象。該段洋脊為典型的對稱擴張洋脊段，軸部水深較淺，兩側地形起伏較小，且分布有大量的火山錐(圖1)。海底地形數據來源于中國大洋24航次的船載多波束測深資料，利用Teledyne Caris HIPS和SIPS軟件將原始多波束數據生成連續平滑的曲面，使用克里金插值方法生成分辨率約為55 m的DEM數據[34-35]。

圖1 研究區位置及典型火山錐地形特征Fig.1 Location of the study area and topographical features of typical volcanic cones(圖a為研究區位置，圖b為研究區地形特征，圖c、d、e為典型火山錐三維地形特征。紅色線指示擴張中心，黑色線指示斷層，藍色方框為圖4a范圍。)(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 自動識別方法

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

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

2.1 圖像特征變換

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

圖3 研究區坡度與地表粗糙度計算結果Fig.3 The calculated slope and surface roughness for the study area

火山錐像元值在原始圖像中過于分散，不利于機器自動提取。通過對原始DEM數據進行數學運算，使火山錐空間信息產生集聚效果，像元值集中于一定區間，使計算機聚類時與其他地物類別有效區分。特征變換公式如下

(1)

(2)

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

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

圖5 研究區特征變換后的圖像(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聚類非監督分類

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

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

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

(3)

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

(2)通過遷移平均值法計算新的聚類中心：

(4)

式中：K為聚類個數，Qi是類別i包含的總像素個數。

(3)像素被分配到鄰近的聚類，并計算新聚類的平均值和協方差。

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

通過ISO聚類可得到像元與聚類之間的所屬關系、聚類的平均值、方差以及協方差等多元統計信息。依據ISO聚類統計特征文件，使用最大似然法分類將圖像分為若干個聚類，創建分類圖像。保留自動分類圖像中包含火山錐的聚類作為有效區域，進行下一步提取分析(圖6)。

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

2.3 景觀形狀指數篩選

將ISO聚類非監督分類的有效聚類轉化為shape格式的圖斑，這些圖斑不僅包含了需要提取的火山錐，也包括同為高值聚類的斷層、山脊等干擾地貌類型。根據火山錐的獨特形貌特征，對圖斑進行面積與形狀篩選，可以得到自動識別的火山錐。通過計算各圖斑的周長、面積等形態參數，利用景觀形狀指數公式進行形狀篩選[44-45]：

(5)

式中：E為斑塊周長；A為斑塊面積；LSI為景觀形狀指數(圓形)，景觀生態學中常用于計算區域內某斑塊形狀與相同面積的圓的偏離程度，LSI值為1時，代表標準的圓形，偏離程度越大LSI值越大。

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

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

3 自動識別結果與評價

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

圖8 自動提取結果與人工提取結果對比Fig.8 The comparison of automatic extraction and manual extraction results

評價機器自動識別的效果可以參考真陽性、假陽性、偽陰性三個指標[36]。真陽性是指提取結果與實地相符，即正確識別；假陽性是指提取結果與實地不符，即錯誤識別；偽陰性是指實際存在而未被識別，即遺漏識別。最后通過正確率和重疊率評價自動識別效果，正確率和重疊率計算公式如下

(6)

(7)

對比自動識別與人工提取結果，發現機器識別結果呈真陽性20處，假陽性5處，偽陰性8處。機器識別正確率達到0.8，重疊率約0.7，精度較好。

Kappa系數是衡量分類效果和一致性檢驗的常用指標，系數越高代表一致性越好，數據范圍一般為0～1[23]。對比自動提取與人工提取的火山錐位置、形態等參數，發現兩者的Kappa值為0.996，一致性很好。自動提取與人工提取的火山錐圖斑的面積重疊率為0.816，提取的幾何形態也較吻合。

4 結論與展望

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

與人工解譯相比，雖然目前機器自動識別的精確度仍有待進一步提高，但其高效、易重復的特點，將使其在大批量地形數據的快速自動識別方面發揮重要作用。未來通過大量樣本的先驗訓練進行監督分類，將極大地提高識別的準確度。將非監督分類和監督分類相結合，有望實現不同效率、不同精度要求的機器自動識別。

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