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

（1. 山東建筑大學(xué) 測繪地理信息學(xué)院，山東 濟(jì)南 250101； 2. 可持續(xù)發(fā)展大數(shù)據(jù)國際研究中心，北京 100094；3. 中國科學(xué)院 空天信息創(chuàng)新研究院 數(shù)字地球重點(diǎn)實(shí)驗(yàn)室，北京 100094）

0 引言

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

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

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

1 方法介紹

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

1.1 泰森多邊形法

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

在眾多應(yīng)用領(lǐng)域中，泰森多邊形法的主要功能之一便是利用泰森多邊形的性質(zhì)來提取地理實(shí)體的中心線，并將該中心線作為地理實(shí)體的抽象以近似表示該實(shí)體［23］。在構(gòu)建泰森多邊形的過程中，算法首先會(huì)提取地理實(shí)體邊界線上的點(diǎn)作為離散點(diǎn)，根據(jù)離散點(diǎn)構(gòu)建泰森多邊形并對(duì)生成的泰森多邊形的邊進(jìn)行篩選，只保留位于地理實(shí)體內(nèi)部的部分作為該實(shí)體的抽象表示。

1.2 最短距離公式

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

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

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

表1 各參數(shù)值和單位Table 1 Parameter values and units

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

1.3 Dijkstra算法

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

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

2 應(yīng)用研究

2.1 研究區(qū)域與數(shù)據(jù)源

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

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

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

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

2.2 冰川中流線提取

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

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

泰森多邊形法生成的冰川中流線與冰川矢量邊界線上離散點(diǎn)采樣密度有關(guān)，即隨著離散點(diǎn)采樣密度的增加，泰森多邊形提取的中流線也會(huì)更加準(zhǔn)確（圖2），但是較高的采樣密度會(huì)降低模型的運(yùn)行效率。因此作為提取冰川中流線的第一步，模型首先會(huì)根據(jù)冰川外邊界線判斷冰川面積大小，并依據(jù)冰川面積采用經(jīng)驗(yàn)參數(shù)對(duì)模型進(jìn)行初始化。即對(duì)于大面積冰川（冰川面積大于5 km2），模型在構(gòu)建泰森多邊形時(shí)便采用100 m 的采樣密度，而對(duì)于小面積冰川（冰川面積小于5 km2）則采用10 m 的采樣密度以在保證準(zhǔn)確性的前提下提高模型的運(yùn)行效率。然而，由于部分冰川內(nèi)存在裸巖區(qū)即內(nèi)邊界線，因此為了保證模型生成的冰川中流線不穿過裸巖區(qū)，模型會(huì)根據(jù)冰川矢量邊界線數(shù)據(jù)對(duì)該冰川區(qū)域內(nèi)是否存在裸巖區(qū)進(jìn)行判斷。當(dāng)區(qū)域內(nèi)無裸巖區(qū)時(shí)，模型只需要提取外邊界線并對(duì)其進(jìn)行加密，而當(dāng)區(qū)域內(nèi)存在裸巖區(qū)時(shí)，模型則需要同時(shí)提取內(nèi)、外邊界線并對(duì)其進(jìn)行加密以保證結(jié)果滿足冰川中流線不穿過裸巖區(qū)的繪制原則。冰川矢量邊界線數(shù)據(jù)處理結(jié)果如圖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

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

圖4 局部最高點(diǎn)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）

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

3 應(yīng)用評(píng)價(jià)

對(duì)冰川中流線提取精度進(jìn)行評(píng)估最準(zhǔn)確的方法是將提取結(jié)果與實(shí)際冰川中流線進(jìn)行比較，然而，由于山地冰川地處偏遠(yuǎn)山區(qū)，很多冰川人力無法到達(dá)，因此很難實(shí)測冰川的真實(shí)中流線。但是為了對(duì)本文所提方法進(jìn)行精度評(píng)估，我們選擇采用與其他冰川中流線提取方法提取結(jié)果進(jìn)行對(duì)比的策略，以評(píng)估本文方法的可靠性和有效性。

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

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

為了進(jìn)一步比較兩方法的提取結(jié)果，本文引入了冰川長度比Ra/k用以評(píng)價(jià)本文使用的泰森多邊形法提取結(jié)果與使用Kienholz 方法提取結(jié)果的優(yōu)劣。具體計(jì)算公式如下：

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

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

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

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

圖7 三方法提取結(jié)果對(duì)比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 兩方法冰川中流線提取結(jié)果Fig. 8 The glacier centerline extraction results of two methods

4 討論與結(jié)論

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

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

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