本發(fā)明涉及一種內(nèi)流流域提取方法，特別是一種基于priority-flood的內(nèi)流流域提取方法。

背景技術(shù)：

內(nèi)流區(qū)指的是它的集水區(qū)與外部海洋環(huán)境沒有水文聯(lián)系的流域，內(nèi)流區(qū)域占地球表面的20％，分布范圍廣泛，是地球地貌、地表生態(tài)環(huán)境的重要組成部分。針對內(nèi)流流域自身的地形、水文特點(diǎn)，進(jìn)行流域單元的劃分和提取，是研究內(nèi)流區(qū)及內(nèi)流流域的水文過程、地貌發(fā)育、生物演化等特征重要基礎(chǔ)。

基于數(shù)字高程模型(dem)的水文分析方法是目前流域劃分提取和研究分析的常用方法。其中流向模型算法是首要步驟且是最重要的基礎(chǔ)部分，目前有單流向算法和多流向算法。d8(deterministiceight-neighbors)作為典型的單流向算法，將水流方向可能性抽象為平面上的八個(gè)方向，限制了對水流實(shí)際上無窮多種運(yùn)動(dòng)方向的描述精度。后續(xù)也提出來諸如rho8、rho4、demon和dinf等諸多改進(jìn)方法，但單流向算法將所有水量匯入統(tǒng)一下游單元的設(shè)計(jì)必然導(dǎo)致它不適合模擬在地勢平坦的坡面上水流散漫流動(dòng)的情況，在平面區(qū)域會(huì)出現(xiàn)大量不合理的平行河網(wǎng)。為解決該問題，人們提出來多流向算法，可以使中心格網(wǎng)單元有多個(gè)流出單元，并按一定的權(quán)重分配水量。但合理參數(shù)值的選取成為難以確定的問題。不管基于何種流向算法，基于dem的流域劃分提取方法通常包括三個(gè)步驟：計(jì)算流向—提取河網(wǎng)—?jiǎng)澐至饔颉?/p>

在現(xiàn)有的流向算法中，計(jì)算流向的必要前提步驟是填洼，因?yàn)橥莸氐拇嬖跁?huì)令水流無法流出，導(dǎo)致接下來的河網(wǎng)及流域提取失敗。但是針對于內(nèi)流區(qū)，此種方法顯然是錯(cuò)誤的。內(nèi)流區(qū)本身就是作為洼地存在，這也是內(nèi)流區(qū)與外界沒有水文聯(lián)系的原因。按常規(guī)的方法填充洼地之后，喪失了整個(gè)內(nèi)流區(qū)的地形地貌信息，整個(gè)內(nèi)流區(qū)作為一個(gè)平地存在，后續(xù)提取匯流網(wǎng)絡(luò)時(shí)會(huì)產(chǎn)生大量的平行河網(wǎng)，所劃分和提取內(nèi)流流域的結(jié)果也是錯(cuò)誤的。這也就導(dǎo)致了現(xiàn)有的方法無法用于內(nèi)流區(qū)流域的劃分與提取。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明所要解決的技術(shù)問題是提供一種基于priority-flood的內(nèi)流流域提取方法，它能夠自動(dòng)化、高效的提取內(nèi)流流域。

為解決上述技術(shù)問題，本發(fā)明所采用的技術(shù)方案是：

一種基于priority-flood的內(nèi)流流域提取方法，其特征在于包含以下步驟：

步驟一：定義存儲(chǔ)數(shù)據(jù)及判斷柵格訪問順序的三個(gè)隊(duì)列及流向柵格：包括平地隊(duì)列、洼地隊(duì)列、坡地隊(duì)列，以及用來存儲(chǔ)流向結(jié)果的柵格數(shù)據(jù)；

步驟二：將所有邊界柵格壓入坡地隊(duì)列中：首先選取dem數(shù)據(jù)的邊界柵格，然后根據(jù)制定的邊界柵格流向規(guī)則計(jì)算每一個(gè)柵格的流向，并將流向柵格對應(yīng)位置賦以相應(yīng)的流向編碼，然后將邊界柵格標(biāo)記為已處理，最后將所有柵格壓入坡地隊(duì)列中；

步驟三：循環(huán)處理三個(gè)隊(duì)列中的元素，直到三個(gè)隊(duì)列都為空：根據(jù)元素選取順序決策規(guī)則，從三個(gè)隊(duì)列中選擇彈出一個(gè)元素進(jìn)行處理：循環(huán)處理彈出柵格所有未處理的鄰域柵格；

步驟四：根據(jù)高程判斷待處理柵格所屬類型：通過比較待處理柵格高程與溢出高程的大小，確定鄰域柵格所屬的隊(duì)列；

步驟五：根據(jù)待處理柵格類型選擇不同的流向計(jì)算方法：通過比較待處理柵格高程與中心柵格高程的大小，根據(jù)不同的情況選擇相應(yīng)的流向計(jì)算方法，然后將流向柵格對應(yīng)位置賦以相應(yīng)的流向編碼；

步驟六：判斷待處理柵格點(diǎn)是否為洼地出水口點(diǎn)：確立洼地出水口的特征，判斷待處理柵格是否為洼地出水口點(diǎn)，若是，標(biāo)記此柵格點(diǎn)；更新洼地面積：判斷待處理柵格是否為洼地，若是，更新其所屬洼地面積大小，然后將其壓入所屬的隊(duì)列中，并將待處理柵格標(biāo)記為已處理；

步驟七：篩選洼地，歸并流域：選擇合適的洼地面積閾值，篩選出符合要求的洼地及其洼地出水口點(diǎn)，然后遞歸處理將所有流向出水口點(diǎn)的柵格標(biāo)記為同一流域。

進(jìn)一步地，所述步驟一中三個(gè)柵格要素隊(duì)列包括平地隊(duì)列為普通的先進(jìn)先出隊(duì)列、洼地隊(duì)列也為普通的先進(jìn)先出隊(duì)列和坡地隊(duì)列為最小優(yōu)先隊(duì)列。

進(jìn)一步地，所述步驟二中，邊界柵格流向規(guī)則為：按照d8編碼順序遍歷該邊界柵格的鄰域柵格，找到第一個(gè)在邊界外的或者值為nodata的柵格，邊界柵格即流向該柵格。

進(jìn)一步地，所述步驟三中，元素選取決策順序規(guī)則為：首先判斷洼地隊(duì)列是否為空，若不為空，彈出洼地隊(duì)列隊(duì)首元素作為中心柵格；否則，繼續(xù)判斷平地隊(duì)列是否為空，若平地隊(duì)列不為空，則彈出平地隊(duì)列隊(duì)首元素作為中心柵格；否則，則繼續(xù)判斷坡地隊(duì)列是否為空，若坡地隊(duì)列不為空，彈出坡地隊(duì)列隊(duì)首元素作為中心柵格，并將溢出高程設(shè)為該中心柵格高程值。

進(jìn)一步地，所述步驟四中，根據(jù)高程判斷待處理柵格所屬類型的方法為：若待處理柵格高程小于溢出高程，則待處理柵格屬于洼地；若待處理柵格高程等于溢出高程，則待處理柵格屬于平地；若待處理柵格高程大于溢出高程，則待處理柵格屬于坡地。

進(jìn)一步地，所述步驟五中，選擇不同的流向計(jì)算方法為：若待處理柵格高程值等于中心柵格高程值，則待處理柵格流向中心柵格；若待處理柵格高程值大于中心柵格，則找出待處理柵格鄰域中所有已處理柵格，在這些柵格中找到高程值比待處理柵格高程值小且坡降最大的柵格，流向該柵格；若待處理柵格小于中心柵格，則找出待處理柵格鄰域中所有已處理柵格，在這些柵格中找到高程值比待處理柵格高程值大且坡降最小的柵格，流向該柵格。

進(jìn)一步地，所述步驟六中，判斷待處理柵格點(diǎn)是否為洼地出水口點(diǎn)和更新洼地面積的方法為：當(dāng)洼地隊(duì)列為空時(shí)，第一個(gè)比溢出高程小的點(diǎn)即為新的洼地出水口點(diǎn)；更新洼地面積的方法：根據(jù)第四步中判斷的待處理柵格所屬的類型，若待處理柵格點(diǎn)屬于洼地且不是洼地出水口點(diǎn)，則將當(dāng)前洼地面積大小加一，并將其壓入洼地隊(duì)列中；若屬于平地，則將其壓入平地隊(duì)列；若屬于坡地，則將其壓入坡地隊(duì)列中。

進(jìn)一步地，所述步驟七中，篩選洼地及歸并流域的方法為：通過合適的面積閾值篩選洼地并尋找出水口點(diǎn)，并將流向出水口點(diǎn)的柵格壓入隊(duì)列中，彈出隊(duì)列中元素，將流向該元素的柵格繼續(xù)壓入隊(duì)列中，循環(huán)處理，直到隊(duì)列為空。

本發(fā)明與現(xiàn)有技術(shù)相比，具有以下優(yōu)點(diǎn)和效果：

1、本發(fā)明采用先漫水示蹤法，通過模擬水流注滿洼地并溢出的過程，獲取洼地內(nèi)正確的匯流方向并同時(shí)獲得水溢出每個(gè)洼地時(shí)的“出水口點(diǎn)”以及該洼地的面積；然后根據(jù)洼地面積大小篩選符合要求的洼地；根據(jù)流向數(shù)據(jù)，將流向出水口點(diǎn)的所有柵格歸并為一個(gè)區(qū)域，即完成內(nèi)流流域的劃分與提取。

2、本發(fā)明適用于基于dem進(jìn)行內(nèi)流區(qū)的流域提取，無需填充洼地，且只需對所有柵格遍歷一次，能夠自動(dòng)化、高效的提取內(nèi)流流域，為下一步內(nèi)流流域的地形地貌和水文分析奠定重要的基礎(chǔ)。

3、改變了傳統(tǒng)利用dem計(jì)算流向時(shí)必須進(jìn)行洼地填充的方法，使得dem水文分析擴(kuò)展到內(nèi)流區(qū)域，打破了傳統(tǒng)方法只能用于外流區(qū)域的局限；為進(jìn)一步分析內(nèi)流區(qū)流域的水文、地貌特征等奠定了基礎(chǔ)。

附圖說明

圖1是本發(fā)明的基于priority-flood的內(nèi)流流域提取方法的流程框圖。

圖2是本發(fā)明的流向算法流程圖。

圖3是本發(fā)明的洼地篩選及流向歸并流程圖。

圖4是本發(fā)明的單流向d8流向編碼示意圖。

圖5是本發(fā)明的優(yōu)先漫水算法(priority-flood)的示意圖。

圖6是本發(fā)明的柵格類型示意圖(深色柵格表示邊界，外圍淺色柵格表示nodata)。

圖7是本發(fā)明的柵格流向計(jì)算示意圖。

圖8是本發(fā)明的整體流向示意圖。

圖9是本發(fā)明的洼地及內(nèi)流流域示意圖。

圖10是本發(fā)明的實(shí)施例的青藏高原內(nèi)流區(qū)內(nèi)流流域劃分圖。

具體實(shí)施方式

下面結(jié)合附圖并通過實(shí)施例對本發(fā)明作進(jìn)一步的詳細(xì)說明，以下實(shí)施例是對本發(fā)明的解釋而本發(fā)明并不局限于以下實(shí)施例。

如圖所示，本發(fā)明的基于priority-flood的內(nèi)流流域提取方法，包含以下步驟：

步驟一：各種變量的定義及初始化。首先，讀取原始dem數(shù)據(jù)，定義一個(gè)與原始dem大小相同的柵格存儲(chǔ)流向數(shù)據(jù)。其次，定義存儲(chǔ)柵格單元元數(shù)據(jù)的結(jié)構(gòu)體，包括柵格行號、柵格列號及柵格高程值。該結(jié)構(gòu)體即為之后三個(gè)隊(duì)列中所要存儲(chǔ)的數(shù)據(jù)。接著定義平地隊(duì)列，用來存儲(chǔ)標(biāo)記為平地類型的柵格；定義坡地隊(duì)列，用來存儲(chǔ)標(biāo)記為坡地類型的柵格；定義洼地隊(duì)列，用來存儲(chǔ)標(biāo)記為洼地類型的柵格。其中，平地隊(duì)列和洼地隊(duì)列為普通的先進(jìn)先出隊(duì)列，坡地隊(duì)列為最小優(yōu)先隊(duì)列(即屬性值最小的元素?fù)碛凶罡叩膬?yōu)先級，會(huì)最優(yōu)先被彈出處理)。接著定義洼地出水口標(biāo)識為false，定義溢出高程為0。

步驟二：將所有邊界柵格壓入坡地隊(duì)列中。首先，循環(huán)遍歷所有柵格，對于所有的邊界柵格(即柵格有效數(shù)據(jù)范圍的邊界，如果一個(gè)柵格的鄰域柵格中存在空值或無數(shù)據(jù)，即判斷該柵格為邊界柵格)，按照邊界柵格流向算法，計(jì)算該柵格的流向，修改對應(yīng)位置的流向柵格的值。然后將該柵格壓入坡地隊(duì)列中。以上流向均采用d8流向編碼。計(jì)算邊界柵格流向的算法為：按照d8編碼順序遍歷一個(gè)邊界柵格的鄰域，柵格流向第一個(gè)處于數(shù)據(jù)范圍外或值為nodata的柵格。d8算法計(jì)算坡降的計(jì)算公式如下：

式中，slpc為待求中心柵格的最大等高線加權(quán)高差鄰域單元的坡度，hc為中心柵格單元的高程值，hi為第i號鄰域單元的高程值，li為第i號鄰域單元的等高線加權(quán)值。li的定義如下：

步驟三：處理三個(gè)隊(duì)列元素。首先，判斷三個(gè)隊(duì)列是否都為空，若是，則流向計(jì)算結(jié)束，跳轉(zhuǎn)到第七步；否則繼續(xù)完成本步計(jì)算和接下來第四、五、六步。依次按照洼地柵格、平地柵格和坡地柵格的順序，找到第一個(gè)不為空的隊(duì)列，并將彈出隊(duì)列中元素，將其稱之為中心柵格。如果彈出的是坡地隊(duì)列，還要設(shè)置溢出高程為中心柵格高程。如果洼地隊(duì)列為空，需將洼地出水口標(biāo)識設(shè)置為false，以便開始標(biāo)記下一塊洼地出水口。

步驟四：判斷中心柵格鄰域柵格所屬類型。在第三步確定中心高程后，遍歷處理每一個(gè)鄰域柵格(稱為待處理柵格)，若待處理柵格高程小于溢出高程，則待處理柵格屬于洼地，在完成第五、六步計(jì)算后，將其壓入洼地隊(duì)列；若待處理柵格高程等于溢出高程，則待處理柵格屬于平地，在完成第五、六步計(jì)算后，將其壓入平地隊(duì)列；若待處理柵格高程大于溢出高程，則待處理柵格屬于坡地，在完成第五、六步處理后，將其壓入坡地隊(duì)列。

步驟五：計(jì)算待處理柵格的流向。若待處理柵格高程等于中心柵格高程，則待處理柵格流向中心柵格即可；若待處理柵格高程小于中心柵格高程，則先找到待處理柵格鄰域中已處理的柵格，再從這些柵格中找到高程值大于待處理柵格且坡降最小的柵格，流向該柵格即可；若待處理柵格高程大于中心柵格高程，則先找到待處理柵格鄰域中已處理的柵格，再從這些柵格中找到高程值小于待處理柵格且坡降最大的柵格，流向該柵格即可。

步驟六：判斷出水口點(diǎn)并更新洼地面積。若待處理柵格點(diǎn)的高程小于中心柵格點(diǎn)，且洼地標(biāo)識為false，則將該待處理柵格點(diǎn)標(biāo)記為洼地出水口點(diǎn)，修改洼地標(biāo)識為true，同時(shí)新建一條洼地記錄，包括洼地出水口點(diǎn)的坐標(biāo)位置、高程信息和洼地面積大小。若待處理柵格高程小于中心柵格高程但不是洼地出水口點(diǎn)，則該待處理柵格點(diǎn)屬于洼地，將上一條洼地記錄的面積加一。

步驟七：篩選洼地及歸并流域。根據(jù)實(shí)地情況確定合適洼地面積閾值，選擇符合要求的洼地，生成一個(gè)洼地列表。從列表中取出一塊洼地，獲得出水口點(diǎn)柵格，將其加入一個(gè)臨時(shí)隊(duì)列中，并賦以一個(gè)標(biāo)記值(代表某一流域)。從臨時(shí)隊(duì)列中取出一個(gè)柵格元素，判斷該柵格所有鄰域柵格，如果鄰域柵格流向該柵格且鄰域柵格不屬于洼地列表中洼地的出水口點(diǎn)，則將該鄰域柵格標(biāo)記為同一流域，并加入臨時(shí)隊(duì)列中。循環(huán)處理臨時(shí)隊(duì)列，直到臨時(shí)隊(duì)列為空。

下面以青藏高原的內(nèi)流區(qū)(面積)為例對本發(fā)明進(jìn)行說明。

1)獲取數(shù)據(jù)。本例中使用分辨率90米的srtm數(shù)據(jù)。根據(jù)青藏高原內(nèi)流區(qū)范圍截取相應(yīng)的dem數(shù)據(jù)。

2)隊(duì)列及變量初始化。讀入dem數(shù)據(jù)，新建與dem大小相等的流向數(shù)據(jù)，初始化相關(guān)的隊(duì)列、變量等。

3)數(shù)據(jù)邊界壓入隊(duì)列。識別數(shù)據(jù)邊界，將邊界存儲(chǔ)到優(yōu)先隊(duì)列中。

4)隊(duì)列數(shù)據(jù)處理。根據(jù)選擇順序選取不同的隊(duì)列元素，處理元素的鄰域柵格，計(jì)算各柵格的流向，判斷洼地出水口點(diǎn)，計(jì)算洼地面積，并將其壓入所屬的隊(duì)列中，循環(huán)處理直到三個(gè)隊(duì)列都為空。最后，獲得柵格流向，識別除青藏高原內(nèi)流區(qū)內(nèi)的所有洼地，以及每一塊洼地的面積大小和相應(yīng)的出水口點(diǎn)位置。

5)確定面積閾值并篩選。根據(jù)青藏高原樣區(qū)的實(shí)際情況及研究目的，選擇合適的洼地面積閾值，并根據(jù)閾值的大小篩選識別出來的洼地。

6)歸并流域。根據(jù)所識別出來的洼地出水口的位置，將所有匯入該點(diǎn)的流向柵格歸并為同一流域，合并小流域。最后得到青藏高原內(nèi)流流域劃分提取結(jié)果。

本說明書中所描述的以上內(nèi)容僅僅是對本發(fā)明所作的舉例說明。本發(fā)明所屬技術(shù)領(lǐng)域的技術(shù)人員可以對所描述的具體實(shí)施例做各種修改或補(bǔ)充或采用類似的方式替代，只要不偏離本發(fā)明說明書的內(nèi)容或者超越本權(quán)利要求書所定義的范圍，均應(yīng)屬于本發(fā)明的保護(hù)范圍。