山地冰川中流線自動提取方法研究

張一新， 周建民， 桑文剛， 李 震， 黃 磊， 魯安新

（1. 山東建筑大學 測繪地理信息學院，山東 濟南 250101； 2. 可持續發展大數據國際研究中心，北京 100094；3. 中國科學院 空天信息創新研究院 數字地球重點實驗室，北京 100094）

0 引言

冰川是氣候條件與地形條件相結合的產物，是氣候變化的指示器，對氣候變化具有重要的反饋作用［1-2］。根據第二次冰川編目，中國是世界上中低緯度地區山地冰川最豐富的國家，共有冰川48 571條，面積約5.18×104km2［3］。冰川在自身重力的作用下會由高海拔區域向低海拔方向運動進而形成冰川冰流［4］，冰川冰的最大流量軌跡稱為冰川主流線［5］。當前，冰川主流線的提取方法是利用水文分析通過計算匯水線來實現，然而該方法計算的匯水線在地形平坦地區無法追索到冰川頂端邊界［6］。因此，學者們提出了一種利用冰川中流線替代冰川主流線的提取方法［7-8］，冰川中流線是冰川主流線的中心線［5］，是反映冰川幾何形態的重要參數之一，在測量冰川長度隨時間變化［9］、分析冰川運動速度變化特征［10-11］、估算冰川體積［12］及構建一維冰川模型等應用中都有著重要的意義［13］。

近年來國內外學者針對冰川中流線的提取開展了大量的研究工作［5］，一些自動與半自動的提取方法也被陸續提出來。Le Bris等［7］率先將冰川軸線引入到冰川中流線提取研究中，通過連接不同高程帶等高線中點的方式獲得冰川中心線實現冰川中流線的提取，但該方法只能提取單一中心線，且對于面積較大、形狀復雜的冰川其提取結果并不一定為主冰流的中心線。Machguth等［8］提出了一種結合冰川表面寬度和坡度的冰川中流線提取方法，其成功率達到了95%～98%，但該方法需要較多的人工參與，進行模型參數的調整，對于數字高程模型數據也有較高的質量要求，且無法提取冰川支流的中心線。Kienholz等［14］將成本距離思想引入到冰川中流線提取過程中，利用歐式距離和高程對冰川中流線進行約束從而提高結果的準確性，該方法彌補了多分支冰川只生成單一中心線的缺點，但算法過程復雜，參數較多，且部分冰川中流線仍需要人工調整。姚曉軍等［15］在GIS軟件的支持下提出了針對單一盆地與單一出口、復式盆地與單一出口、冰帽三種類型冰川中流線的自動提取方案，其優化方案使提取結果更加合理，但該方法僅實現了對單一盆地與單一出口類型冰川中流線的自動提取，對于復式盆地與單一出口和冰帽類型的冰川中流線的提取仍需要專家知識的輔助，未實現完全自動化處理。楊佰義等［6］綜合利用冰川中心線法和冰川主流線法獲取了冰川長度線，在一定程度上提高了結果的準確性，但對于大面積冰川，由于受到地形起伏和數字高程模型數據質量的影響，使得獲取的冰川長度線無法從海拔最高點流向最低點，需要人工進行修正。Ji 等［16］在ArcGIS 環境下利用Le Bris 和Pual 提出的方法，結合數字高程模型數據和冰川矢量邊界線數據實現了冰川中流線的提取，然而該研究在高程最值選取等過程仍需要人工參與。Zhang 等［5］基于歐幾里得分配和冰川表面地形特征，在ArcGIS環境下實現了對冰川中流線的自動提取，但該模型結構復雜，對于部分多出口坡面冰川和未分割的單一出口冰川的提取效果不理想，致使所得的冰川長度不準確。Hansen 等［17］通過計算數字高程模型數據的局部線性回歸梯度下降來提取冰川中流線，利用組合對應于不同窗口大小的冰川中流線來提高結果的準確性，但是，此方法計算過程較為復雜且對于冰川中流線的提取仍需要人工提供初始點，而未完全實現自動化。

鑒于此，本文提出了基于泰森多邊形法的山地冰川中流線自動提取方法，該方法以冰川編目數據中的冰川矢量邊界線數據與數字高程模型數據作為輸入，通過泰森多邊形法生成冰川中心線，然后對中心線進行篩選從而獲得冰川中流線。本方法在無需人工參與的情況下，不僅實現了對冰川中流線的自動化快速提取，還能準確地對冰川各支流的中流線進行提取，解決了基于冰川軸線方法只能提取單一中心線的問題，提高了冰川中流線提取效率。

1 方法介紹

本文所用到的核心算法是泰森多邊形法、最短距離公式和Dijkstra 算法，利用泰森多邊形法生成冰川內部的中心線作為相應冰川中流線的預選線段并使用最短距離公式對具有多起始點的中心線進行篩選，刪除多余起始點，然后使用Dijkstra 算法沿冰川內部中心線進行最短路徑選擇，完成冰川中流線的提取。

1.1 泰森多邊形法

泰森多邊形又稱為Voronio圖或Dirichlet圖，最初是由俄國數學家Georgy Fedoseevich Voronio 提出的一種關于空間鄰近的算法，并在1900年由荷蘭科學家A. H. Thiesseny 引入到降雨量的空間分析中，并將該方法正式稱為泰森多邊形法。該方法實質上是對空間平面的劃分，即泰森多邊形內的點至相應離散點的距離最近且位于泰森多邊形上的點至其兩邊離散點的距離相等［18-19］。泰森多邊形的性質使該方法在圖像幾何、計算機圖形處理、計算機視覺等領域有著廣泛的應用，同時隨著對泰森多邊形法研究的逐漸深入，也使其成為解決空間分析問題的重要方法［20-22］。

在眾多應用領域中，泰森多邊形法的主要功能之一便是利用泰森多邊形的性質來提取地理實體的中心線，并將該中心線作為地理實體的抽象以近似表示該實體［23］。在構建泰森多邊形的過程中，算法首先會提取地理實體邊界線上的點作為離散點，根據離散點構建泰森多邊形并對生成的泰森多邊形的邊進行篩選，只保留位于地理實體內部的部分作為該實體的抽象表示。

1.2 最短距離公式

由于本文需要對冰川所有中流線進行提取，因此需要確定每條冰川中流線的起始點。本文對冰川中流線起始點的確定主要分為兩步：（1）沿冰川矢量邊界線結合數字高程模型數據提取邊界線上各折點的高程，并獲取該邊界線上高程最高點與最低點。（2）將邊界線上的折點與邊界線上鄰近折點（每個方向取20個鄰近點）的高程進行比較，若該折點高程高于所有鄰近折點，同時大于相應冰川邊界線上各折點高程分布的三分之一，則該點即為當前冰川中流線的起始點。由于冰川起始點通常位于高海拔地區，采用每個冰川三分之一高度作為閾值，降低了提取的冰川中流線起始于地勢低洼地區的可能性，同時鄰近點的選取也保證了所提取的起始點即為局部高程最值點。

然而，部分冰川中流線存在多個起始點的問題，為了解決該問題，本文引入了最短距離公式對上述過程提取到的冰川中流線起始點進行篩選［14］，具體計算公式如下：

式中：r為兩起始點間最短距離（m）；S為冰川面積（m2）；q1、q2分別為方程的兩個參數，其值如表1所示。

表1 各參數值和單位Table 1 Parameter values and units

如果范圍r中有多個起始點，方法只保留高程最高的起始點并刪除其余起始點，同時方法限定r的取值范圍最小為500 m。

1.3 Dijkstra算法

本文方法提取的冰川中流線結果是利用Dijkstra 算法搜索起始點與終點間距離最短的連線即最短路徑實現的。Dijkstra 算法是目前廣泛應用的求解最短路徑算法之一，它是由荷蘭科學家E.D.Dijkstra 提出的一種典型的單源最短路徑算法。算法采用貪心算法的策略，通過計算起始點至終點不同路徑下的權重，篩選出權重值最小的路線作為兩點間的最短路徑［24］。

Dijkstra算法的基本思想如下：對任意一個帶權無向圖，取S 為頂點集合，P 為各點間權重的集合，以兩點間連線的距離作為權重，若兩點間無直線連接則將權重設為無窮大。在算法開始時將頂點集合分為兩部分，第一部分為S1代表已求出最短路徑的頂點，此時該部分僅包含起始點，在之后每求得一個最短路徑便將對應頂點放入該數組中，直到所有頂點遍歷完畢。第二部分為S2代表待確定最短路徑頂點集，算法將遍歷起始點的所有鄰近頂點并從集合P 中查詢兩點間權重值，篩選權重值最小的連線作為當前頂點的最短路徑，并將該頂點存儲至S1中，當S2中頂點遍歷完成后，集合S1中的所有頂點的連線即為起始點至目標點間的最短路徑。

2 應用研究

2.1 研究區域與數據源

青藏高原及周邊地區的平均海拔超過4 000 m，是除南極和北極以外冰川最為發育的地區，

冰川面積約為9.8×104km2，也是中國境內山地冰川最豐富的地區［25-26］。因此，本文選擇青藏高原東南部地區作為研究區域，該區域內廣泛分布有山地冰川1 014 條，冰川形狀各異，且地形條件極其復雜［27］，為本文方法進行山地冰川中流線的自動提取研究提供了優異的實驗區域。

本文選擇Randolph Glacier Inventory（RGI）v6.0 數據集作為研究區域的冰川數據源，該數據集是由美國國家冰雪數據中心負責管理和更新的全球范圍的冰川編目數據集。RGI v6.0 數據集將青藏高原及其周邊地區劃分為了喜馬拉雅地區、西藏南部和東部地區、橫斷山脈地區及昆侖山和祁連山地區等共計13個子區域［28］，并分別記錄了它們的面積、坡度和坡向等空間信息，為進行山地冰川中流線的自動提取研究提供了條件。

本文使用的數字高程模型是美國航天航空局（NASA）及地理信息情報局等聯合提供的SRTM1（Shuttle Radar Topography Mission1）數字高程模型，該模型由美國航天飛機“奮進號”歷時11天采集的60° N 到56° S 之間的所有高程數據組成［29-30］，覆蓋面積達1.19×108km2，占地球面積的80%。同時，由于其在地球大部分區域內具有統一的分辨率和精度而被廣泛應用于各個研究領域［31］。相較于ASTER DEM 等數字高程模型數據，SRTM1 數字高程模型有良好的垂直精度且該數據受坡度、土地利用類型及地貌等因素的影響小，保證了其即使在復雜的山地地區也能提供良好的高程數據［32-33］。

2.2 冰川中流線提取

本方法以冰川矢量邊界線數據和數字高程模型數據作為輸入，使用Python 作為基礎編程語言架構山地冰川中流線自動提取模型，同時結合GDAL、Scipy、Shapely 等第三方空間數據處理庫，實現對冰川矢量邊界線數據和數字高程模型數據的讀取、泰森多邊形的構建、模型結果輸出等功能。該模型的主要處理過程如圖1所示。

圖1 冰川中流線提取技術流程Fig. 1 The technical process of glacier centerline extraction

泰森多邊形法生成的冰川中流線與冰川矢量邊界線上離散點采樣密度有關，即隨著離散點采樣密度的增加，泰森多邊形提取的中流線也會更加準確（圖2），但是較高的采樣密度會降低模型的運行效率。因此作為提取冰川中流線的第一步，模型首先會根據冰川外邊界線判斷冰川面積大小，并依據冰川面積采用經驗參數對模型進行初始化。即對于大面積冰川（冰川面積大于5 km2），模型在構建泰森多邊形時便采用100 m 的采樣密度，而對于小面積冰川（冰川面積小于5 km2）則采用10 m 的采樣密度以在保證準確性的前提下提高模型的運行效率。然而，由于部分冰川內存在裸巖區即內邊界線，因此為了保證模型生成的冰川中流線不穿過裸巖區，模型會根據冰川矢量邊界線數據對該冰川區域內是否存在裸巖區進行判斷。當區域內無裸巖區時，模型只需要提取外邊界線并對其進行加密，而當區域內存在裸巖區時，模型則需要同時提取內、外邊界線并對其進行加密以保證結果滿足冰川中流線不穿過裸巖區的繪制原則。冰川矢量邊界線數據處理結果如圖3所示。

圖2 不同采樣密度下泰森多邊形法提取的冰川中流線Fig. 2 Glacier centerlines extracted by Thiessen polygon method with different sampling densities （200 m and 50 m）

圖3 冰川矢量邊界線的提取Fig. 3 The extraction of glacier outlines

當冰川矢量邊界線加密完成后，模型便引入數字高程模型數據沿加密后冰川外邊界線上各點獲取高程值，再依據1.2 節中的步驟選取高程同時滿足大于其附近點高程并大于相應冰川邊界線上各點高程分布的三分之一的點作為冰川中流線的起始點，結果如圖4所示。同時，以加密后外邊界線上的點作為離散點繪制泰森多邊形，結果如圖5（a）所示，由于位于冰川外部的泰森多邊形的邊無法應用于冰川中流線的提取，所以我們只保留位于冰川內部的邊作為冰川中流線的預選線段，結果如圖5（b）所示。然后利用最短距離公式，根據冰川面積計算冰川中流線起始點間的最短距離，若最短距離內存在多個起始點，則需根據數字高程模型數據提取各起始點高程并進行比較，保留高程最高的起始點。

圖4 局部最高點Fig. 4 The local highest point

圖5 利用泰森多邊形法提取冰川中流線Fig. 5 Extraction of glacier centerlines by Thiessen polygon method： voronoi diagram based on scattered points （a），voronoi polygons clipped by glacier vector boundary lines （b）， glacier centerlines extracted using Dijkstra （c）， and glacier length line （d）

在完成冰川中流線多起始點處理后，模型便使用Dijkstra 算法以各中流線起始點與其鄰居頂點之間連線的距離作為該直線的權重計算起始點到終點之間的最短路徑，完成所有中流線的選取工作，結果如圖5（c）所示。最后，選擇最長冰川中流線作為冰川長度線以近似計算該冰川長度［6］，選取結果圖5（d）所示。

3 應用評價

對冰川中流線提取精度進行評估最準確的方法是將提取結果與實際冰川中流線進行比較，然而，由于山地冰川地處偏遠山區，很多冰川人力無法到達，因此很難實測冰川的真實中流線。但是為了對本文所提方法進行精度評估，我們選擇采用與其他冰川中流線提取方法提取結果進行對比的策略，以評估本文方法的可靠性和有效性。

在現有的研究中，Kienholz 等［14］提出的方法同樣采用了冰川矢量邊界線數據和數字高程模型數據作為輸入，為兩方法間同數據源的對比分析提供了條件。本文通過計算平均長度和長度比對兩方法的提取結果進行了初步的比較。由表2 可以看出，利用Kienholz 等［14］的方法計算得到的冰川平均長度為3.156 km，使用本文方法計算得到的冰川平均長度為3.13 km，兩種方法的平均吻合長度為3.11 km 約占總長度的99.4%，兩種方法所提取到的冰川長度平均值幾乎相同。為了便于表述，本文將Kienholz等［14］提出的方法定義為Kienholz方法。

表2 本研究計算的長度與Kienholz等［14］的比較Table 2 Comparison between the length calculated in this study and Kienholz et al［14］

為了進一步比較兩方法的提取結果，本文引入了冰川長度比Ra/k用以評價本文使用的泰森多邊形法提取結果與使用Kienholz 方法提取結果的優劣。具體計算公式如下：

式中：La為本文方法所提取的冰川長度；Lk為Kienholz 方法得到的冰川長度。由表2 可知兩種方法的長度比平均值為0.98，表明兩種方法之間存在較小的負向偏差。

本文根據冰川面積將研究區域劃分為0.1～0.5 km2、0.5～2.5 km2、2.5～10 km2及大于10 km2四組并分別計算兩方法的長度比，以比較兩方法對不同面積冰川的提取結果，計算結果如圖6 所示。結果表明，共有90.2%的數據長度比分布在0.9～1.1之間，4.4%的數據長度比大于1.1，5.4%的數據長度比低于0.9。離散較大的數據多集中于小面積冰川中，當冰川面積大于0.5 km2時兩方法的長度比值向1收斂。造成兩種方法長度比數據離散的原因與冰川長度的不確定性有關，在通常情況下小面積冰川的不確定性約為20%，大面積冰川的不確定性約為15%，且隨著數字高程模型數據質量的提高而逐漸降低［8］。

圖6 不同面積冰川長度比Fig. 6 Length ratio of glaciers of different areas： histogram of calculated results for both methods （a），box plot of calculated results for both methods （b）

通過兩方法冰川長度對比結果可知，兩種方法所提取到的冰川平均吻合長度達到了總長度的99.4%，表明兩方法提取的冰川平均長度幾乎相同，而對于不同面積的冰川，本文分別提取了其冰川長度并與Kienholz 方法提取的相應長度進行比較，計算的長度比平均值達到了0.98，與完全一致（Ra/k=1）相差較小。整體而言，兩種方法所提取到的冰川長度具有很高的一致性，表明本文方法提取的冰川長度具有可靠性和有效性。

對于不同類型冰川的中流線提取，本文又引入了Zhang 等［5］的方法與本文方法和Kienholz 方法的提取結果進行對比，結果如圖7所示，同時為了便于表述本文定義Zhang等［5］的研究方法為Zhang方法。從圖中可以看出對于不同類型的冰川，各方法對冰川中流線的提取結果有一定的差異。對于結構較為簡單的單式山谷冰川，本文方法的提取結果與Kienholz 方法和Zhang 方法的提取結果差異較小，對于結構較為復雜的復式山谷冰川和坡面冰川，本文方法所提取到的冰川中流線數量多于Kienholz方法但少于Zhang方法，而對于冰帽類型冰川，本文方法所提取到的冰川中流線數量均多于Kienholz方法和Zhang 方法的提取結果。圖8 顯示了本文方法與Kienholz方法對研究區域內冰川中流線的提取結果對比，由圖8 可以發現，本文方法提取出了Kienholz方法所沒有提取的冰川支流的中流線，如方框區所示，表明較之Kienholz 方法本文方法提取結果覆蓋更全。

圖7 三方法提取結果對比Fig. 7 Comparison of extraction results of three methods， the extraction results of single-type valley glaciers （a），compound-type valley glaciers （b）， slope-type glaciers （c）， and ice caps （d）

圖8 兩方法冰川中流線提取結果Fig. 8 The glacier centerline extraction results of two methods

4 討論與結論

本文提出了一種基于泰森多邊形法的山地冰川中流線自動提取方法，該方法結合冰川編目數據中的冰川矢量邊界線數據和冰川區數字高程模型數據，通過計算最低點和局部最高點，利用泰森多邊形法提取研究區域內所有冰川中心線，并使用Dijkstra 算法得到冰川中流線。解決了基于冰川軸線方法進行多支流冰川中流線提取過程中只能提取單一中流線的問題，在無人工參與的情況下，實現了山地冰川中流線的自動化快速提取。通過對青藏高原研究區內1 014 條山地冰川進行的中流線提取實驗，本方法共提取了冰川中流線2 114 條，平均長度為3.13 km。通過與現有方法的提取結果進行對比分析，本文方法通過自動化的提取模型獲得的冰川長度與現有方法具有很高的一致性，且本文方法提取的冰川中流線結果覆蓋更全。

然而本文方法在實際應用過程中依然存在著部分缺陷。比如，本方法對于冰川矢量邊界線上離散點采樣密度的選取依然采用經驗參數，這導致了在形狀復雜冰川的中流線提取過程中部分局部最高點無法有效識別。同時，與現行的通用自動化冰川中流線提取算法相同，本文算法在對多出口類型冰川中流線提取過程中同樣以冰川矢量邊界線上最低點作為冰川出口，這可能會導致所提取到的冰川中流線長度出現偏差等問題。本文提出的冰川中流線提取算法，需要以冰川矢量邊界線數據作為輸入。但由于該輸入數據自身存在部分誤差，導致冰川中流線提取結果受到影響，像RGI v6.0 冰川編目數據中部分冰川矢量邊界線呈現鋸齒狀形態，這會使得基于冰川矢量邊界線提取的冰川局部最高點和最低點出現誤差，同時鋸齒狀的冰川矢量邊界線也會使本文方法提取的冰川中流線出現明顯的偏折。

針對上述的缺陷，在后續研究工作中可以采取以下方法提高提取結果的精度：（1）根據冰川類型和復雜程度進一步完善冰川矢量邊界線上離散點采樣密度及局部最高點和最低點的選取方法，以提高方法的適用性。（2）在提取冰川中流線的過程中加入輔助線將多出口類型冰川分割為單一起始點單一出口類型的冰川和多起始點單一出口類型的冰川，在降低其復雜程度的同時實現對多出口類型冰川中流線的提取。（3）在模型中加入預處理算法對存在鋸齒形態的冰川矢量邊界線進行平滑處理，以避免因鋸齒狀冰川矢量邊界線導致的冰川中流線提取錯誤的問題。