本發(fā)明屬于水文技術(shù)領(lǐng)域，具體涉及一種基于數(shù)字高程模型計(jì)算地貌單位線初始概率的方法。

背景技術(shù)：

現(xiàn)代研究認(rèn)為，流域水文響應(yīng)是地貌擴(kuò)散和水動(dòng)力擴(kuò)散對(duì)降落在流域上具有一定時(shí)空分布的凈雨共同作用的結(jié)果。這種對(duì)于流域匯流機(jī)理的揭示，使得水文學(xué)家和自然地理學(xué)家試圖從理論上建立流域地形地貌特征與流域水文響應(yīng)之間的聯(lián)系成為可能。其中的一種方法就是以Horton-Strahler河流分級(jí)理論為基礎(chǔ)，借用統(tǒng)計(jì)物理學(xué)中大量“粒子”運(yùn)動(dòng)宏觀表現(xiàn)的方法來建立地貌單位線。

在建立基于Horton-Strahler河流分級(jí)的地貌瞬時(shí)單位線中重要的一個(gè)步驟就是計(jì)算水滴隨機(jī)降落到流域上各級(jí)坡面的概率(即初始概率)，以及水滴從低級(jí)別河流到高級(jí)別河流的轉(zhuǎn)移的概率(即轉(zhuǎn)移概率)。在1979年Rodiguze-Iturbe等人指出初始概率均與水系結(jié)構(gòu)有關(guān)，并在Horton地貌參數(shù)的基礎(chǔ)上運(yùn)用Smart定律對(duì)其進(jìn)行估算，從而得到水滴從降落在流域面后選擇某一條匯流路徑運(yùn)動(dòng)至流域出口的概率(即路徑概率)，進(jìn)而計(jì)算得到流域匯流時(shí)間的概率分布密度，即地貌瞬時(shí)單位線。但是Smart定律本身就是一種基于數(shù)據(jù)統(tǒng)計(jì)而的得到的規(guī)律，運(yùn)用其估算初始概率時(shí)往往并不能得到很好地精度，從而限制了地貌瞬時(shí)單位線的應(yīng)用，不利于地貌瞬時(shí)單位線的發(fā)展。

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

為了解決現(xiàn)有技術(shù)中存在的不足，本發(fā)明提供了一種基于數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，具有數(shù)據(jù)來源穩(wěn)定可靠、計(jì)算效率高、結(jié)果客觀合理等優(yōu)點(diǎn)，有利于地貌單位線初始概率的直接求取。

為解決上述問題，本發(fā)明具體采用以下技術(shù)方案：

一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，其特征在于，包括以下步驟：

步驟1，提取流域數(shù)字高程模型(DEM)數(shù)據(jù)；

步驟2，利用流域DEM數(shù)據(jù)提取流域分級(jí)河道柵格；

步驟3，基于分級(jí)河道柵格提取河流交匯點(diǎn)；

步驟4，依據(jù)河流交匯點(diǎn)提取各子流域出口點(diǎn)；

步驟5，運(yùn)用子流域出口點(diǎn)提取各子流域面積；

步驟6，以各子流域面積為基礎(chǔ)分析計(jì)算得到地貌單位線初始概率。

前述的一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，其特征在于，所述步驟1中利用地形信息系統(tǒng)軟件ArcGIS構(gòu)建地理處理工作流，獲取流域數(shù)字高程模型(DEM)數(shù)據(jù)，包括以下步驟：

1.1，填洼，得到填洼之后的柵格Fill；

1.2，計(jì)算流向，得到流向柵格Dir；

1.3，計(jì)算匯流流量，得到累計(jì)流量柵格Acc；

1.4，確定流域出口站點(diǎn)StationPiont；

1.5，提取目標(biāo)流域AimWatershed。

前述的一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，其特征在于，所述步驟2中利用流域DEM數(shù)據(jù)提取流域分級(jí)河道柵格的步驟包括：

2.1，由流域DEM數(shù)據(jù)經(jīng)過填洼，計(jì)算流向，計(jì)算匯流流量得到累計(jì)流量柵格Acc；

2.2，設(shè)置匯流閾值threshold，依據(jù)匯流閾值threshold對(duì)累計(jì)流量柵格Acc進(jìn)行重分類，低于匯流閾值threshold的柵格賦空值NODATA，高于匯流閾值threshold的柵格賦值為1，提取出流域河道柵格RiverRaster；

2.3，在2.2中提取的流域河道柵格RiverRaster的基礎(chǔ)上進(jìn)行河流分級(jí)，將流域河道柵格RiverRaster依據(jù)Strahler分級(jí)法進(jìn)行河流分級(jí)，得到分級(jí)河道柵格。

前述的一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，其特征在于，所述步驟3中基于分級(jí)河道柵格提取河流交匯點(diǎn)的步驟包括：

3.1，將分級(jí)河道柵格轉(zhuǎn)換成矢量格式的線型河道LineRiver；

3.2，將分級(jí)河道柵格轉(zhuǎn)換成矢量格式的點(diǎn)型河道PointRiver；

3.3，對(duì)于點(diǎn)型河道PointRiver，按線型河道LineRiver中的鏈與其空間位置相交的原則，逐一提取出點(diǎn)型河道PointRiver中與線型河道LineRiver中的水系的各條鏈相交的點(diǎn)集RallyPoint：

RallyPoint＝LineRiver*PointRiver (式1)

其中：*為相交運(yùn)算；

3.4，將點(diǎn)集RallyPoint之間作相交運(yùn)算，提取出點(diǎn)集RallyPoint之間的交點(diǎn)，即為河網(wǎng)中各條河流鏈之間的交匯點(diǎn)IntersectPoint _ij：

IntersectPoint_ij＝RallyPoint_i*RallyPoint_j(i，j＝1,2……n，且i≠j) (式2)

其中：*為相交運(yùn)算，n為點(diǎn)集RallyPoint的總個(gè)數(shù)，i、j為RallyPoint的序列號(hào)，且i≠j。

前述的一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，其特征在于，所述步驟4中依據(jù)河流交匯點(diǎn)提取各子流域出口點(diǎn)的步驟包括：

4.1，提取所用柵格的分辨率a，依據(jù)分辨率a確定搜索距離d：

4.2，對(duì)步驟3中提取出來的各條河流鏈之間的交匯點(diǎn)進(jìn)行鄰域分析，提取出距每一個(gè)交匯點(diǎn)在搜索距離d之內(nèi)的柵格單元RasterUnit；

4.3，對(duì)4.2中提取出的每一個(gè)的柵格單元RasterUnit進(jìn)行分析，判斷其是否為子流域的出口點(diǎn)，找出柵格單元中的分支河流鏈的注入點(diǎn)即各子流域出口點(diǎn)InjectPoint：

(InjectPoint＝Con(RasterUnit＞threshold)and(RasterUnit＜max(RasterUnit))) (式4)

其中：Con為ArcGIS中條件判斷工具函數(shù)，threshold為匯流閾值。

前述的一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，其特征在于，所述步驟5中運(yùn)用子流域出口點(diǎn)提取各子流域面積的步驟包括：

5.1，在步驟2的2.2中提取出的流域河道柵格RiverRaster基礎(chǔ)上，將流域河道柵格RiverRaster轉(zhuǎn)換為矢量格式Rivershp：

Rivershp＝Rastertoshp(RiverRaster) (式5)

其中：Rastertoshp為ArcGIS中柵格轉(zhuǎn)換為矢量的工具函數(shù)；

5.2，將4.3中提取出來的各子流域出口點(diǎn)InjectPoint與5.1中提取出來的矢量格式Rivershp做相交計(jì)算，提取出各子流域出口點(diǎn)InjectPoint所在的矢量格式Rivershp的小分塊即子流域出口分塊InjectChunk，并讀取子流域出口分塊InjectChunk的值，然后乘以每一個(gè)柵格單元RasterUnit的面積a²，得到相對(duì)應(yīng)的各子流域的面積AChunk：

InjectChunk＝InjectPoint*Rivershp (式6)

其中：*為相交運(yùn)算；

AChunk＝Value(InjectChunk)×a² (式7)

其中：Value(InjectChunk)為子流域出口分塊InjectChunk小塊中的值。

前述的一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，其特征在于，所述步驟6中以各子流域面積為基礎(chǔ)分析計(jì)算得到地貌單位線初始概率的步驟包括：

6.1，將子流域出口分塊InjectChunk與經(jīng)過其中的河流做連接匯總分析，將河流的級(jí)別屬性賦予InjectChunk；

6.2，將同一級(jí)的子流域出口分塊InjectChunk的值相加得到areaA_w，并用高一級(jí)的河流的匯流面積areaA_w+1減去低一級(jí)的河流的匯流面積areaA_w計(jì)算得到每一級(jí)河流的匯水面積A_w：

A_w＝areaA_w+1-areaA_w (式8)

6.3，基于每一級(jí)河流的匯水面積計(jì)算其初始概率θ_w：

θ_w＝A_w/A (式9)

其中：A為流域總的面積。

本發(fā)明的有益效果：本發(fā)明提供的一種利用數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，包括以下主要步驟：提取流域數(shù)字高程模型(DEM)數(shù)據(jù)；利用流域DEM數(shù)據(jù)提取流域分級(jí)河道柵格；基于分級(jí)河道柵格提取河流交匯點(diǎn)；依據(jù)河流交匯點(diǎn)提取各子流域出口點(diǎn)；運(yùn)用子流域出口點(diǎn)提取各子流域面積；以各子流域面積為基礎(chǔ)分析計(jì)算得到地貌單位線初始概率。本發(fā)明具有數(shù)據(jù)來源穩(wěn)定可靠、計(jì)算效率高、結(jié)果客觀合理等優(yōu)點(diǎn)，有利于地貌單位線初始概率的直接求取。

附圖說明

圖1是本發(fā)明的計(jì)算流程示意圖。

圖2為本發(fā)明提取目標(biāo)流域柵格的地理信息系統(tǒng)工作流示意圖。

圖3為本發(fā)明提取出的流域DEM柵格示意圖。

圖4為本發(fā)明提取出的流域河道柵格示意圖。

圖5為本發(fā)明中由河道柵格轉(zhuǎn)換成矢量格式的線型河道LineRiver示意圖。

圖6為本發(fā)明中由河道柵格轉(zhuǎn)換成矢量格式的點(diǎn)型河道PointRiver示意圖。

圖7為本發(fā)明提取河流鏈之間的交匯點(diǎn)IntersectPoint的示意圖。

圖8為本發(fā)明提取出的距每一個(gè)交匯點(diǎn)在搜索距離d之內(nèi)的柵格單元RasterUnit示意圖。

圖9為本發(fā)明提取出的子流域出口分塊InjectChunk示意圖。

具體實(shí)施方式

下面結(jié)合附圖和具體實(shí)施例對(duì)本發(fā)明作進(jìn)一步描述。

如圖1至圖9所示，本發(fā)明提供的一種基于數(shù)字高程模型計(jì)算地貌單位線初始概率的方法，包括以下步驟：

S1、利用地形信息系統(tǒng)軟件ArcGIS構(gòu)建地理處理工作流如圖2，提取目標(biāo)流域的數(shù)字高程模型(DEM)柵格數(shù)據(jù)如圖3，其中包括：

1)填洼，得到填洼之后的柵格Fill；

2)計(jì)算流向，得到流向柵格Dir；

3)計(jì)算匯流流量，得到累計(jì)流量柵格Acc；

4)確定流域出口站點(diǎn)StationPiont；

5)提取目標(biāo)流域AimWatershed；

S2、利用目標(biāo)流域DEM數(shù)據(jù)提取流域分級(jí)河道柵格，包括：

1)由流域DEM數(shù)據(jù)經(jīng)過填洼，計(jì)算流向，計(jì)算匯流流量得到累計(jì)流量柵格Acc；

2)設(shè)置匯流閾值threshold，依據(jù)匯流閾值threshold對(duì)累計(jì)流量柵格Acc進(jìn)行重分類，低于匯流閾值threshold的柵格賦空值NODATA，高于匯流閾值threshold的柵格賦值為1，提取出流域河道柵格RiverRaster；

3)在2.2中提取的流域河道柵格RiverRaster的基礎(chǔ)上進(jìn)行河流分級(jí)，將流域河道柵格RiverRaster依據(jù)Strahler分級(jí)法進(jìn)行河流分級(jí)，得到分級(jí)河道柵格。

S3、基于分級(jí)河道柵格提取河流交匯點(diǎn)，包括：

1)將分級(jí)河道柵格轉(zhuǎn)換成矢量格式的線型河道LineRiver；

2)將分級(jí)河道柵格轉(zhuǎn)換成矢量格式的點(diǎn)型河道PointRiver；

3)對(duì)于點(diǎn)型河道PointRiver，按線型河道LineRiver中的鏈與其空間位置相交的原則，逐一提取出點(diǎn)型河道PointRiver中與線型河道LineRiver中的水系的各條鏈相交的點(diǎn)集RallyPoint：

RallyPoint＝LineRiver*PointRiver(式1)

其中：*為相交運(yùn)算；

4)將點(diǎn)集RallyPoint之間作相交運(yùn)算，提取出點(diǎn)集RallyPoint之間的交點(diǎn)，即為河網(wǎng)中各條河流鏈之間的交匯點(diǎn)IntersectPoint _ij：

IntersectPoint_ij＝RallyPoint_i*RallyPoint_j(i，j＝1,2……n，且i≠j) (式2)

其中：*為相交運(yùn)算，n為點(diǎn)集RallyPoint的總個(gè)數(shù)，i、j為RallyPoint的序列號(hào)，且i≠j。

S4、依據(jù)河流交匯點(diǎn)提取各子流域出口點(diǎn)，包括：

1)提取所用柵格的分辨率a，依據(jù)分辨率a確定搜索距離d：

2)對(duì)步驟3中提取出來的各條河流鏈之間的交匯點(diǎn)進(jìn)行鄰域分析，提取出距每一個(gè)交匯點(diǎn)在搜索距離d之內(nèi)的柵格單元RasterUnit；

3)對(duì)4.2中提取出的每一個(gè)的柵格單元RasterUnit進(jìn)行分析，判斷其是否為子流域的出口點(diǎn)，找出柵格單元中的分支河流鏈的注入點(diǎn)即各子流域出口點(diǎn)InjectPoint：

(InjectPoint＝Con(RasterUnit＞threshold)and(RasterUnit＜max(RasterUnit))) (式4)

其中：Con為ArcGIS中條件判斷工具函數(shù)，threshold為匯流閾值。

S5、運(yùn)用各子流域出口點(diǎn)提取相對(duì)應(yīng)的各子流域面積，包括：

1)在步驟2提取出的流域河道柵格RiverRaster基礎(chǔ)上，將流域河道柵格RiverRaster轉(zhuǎn)換為矢量格式Rivershp：

Rivershp＝Rastertoshp(RiverRaster) (式5)

其中：Rastertoshp為ArcGIS中柵格轉(zhuǎn)換為矢量的工具函數(shù)；

2)將4.3中提取出來的各子流域出口點(diǎn)InjectPoint與5.1中提取出來的矢量格式Rivershp做相交計(jì)算，提取出各子流域出口點(diǎn)InjectPoint所在的矢量格式Rivershp的小分塊即子流域出口分塊InjectChunk，并讀取子流域出口分塊InjectChunk的值，然后乘以每一個(gè)柵格單元RasterUnit的面積a²，得到相對(duì)應(yīng)的各子流域的面積AChunk：

InjectChunk＝InjectPoint*Rivershp (式6)

其中：*為相交運(yùn)算；

AChunk＝Value(InjectChunk)×a² (式7)

其中：Value(InjectChunk)為子流域出口分塊InjectChunk、中的值。

S6、以各子流域面積為基礎(chǔ)分析計(jì)算得到地貌單位線初始概率，包括：

1)將子流域出口分塊InjectChunk與經(jīng)過其中的河流做連接匯總分析，將河流的級(jí)別屬性賦予InjectChunk；

2)將同一級(jí)的子流域出口分塊InjectChunk的值相加得到areaA_w，并用高一級(jí)的河流的匯流面積areaA_w+1減去低一級(jí)的河流的匯流面積areaA_w計(jì)算得到每一級(jí)河流的匯水面積A_w：

A_w＝areaA_w+1-areaA_w (式8)

3)基于每一級(jí)河流的匯水面積計(jì)算其初始概率θ_w：

θ_w＝A_w/A (式9)

其中：A為目標(biāo)流域總的面積。

以上顯示和描述了本發(fā)明的基本原理、主要特征及優(yōu)點(diǎn)。本行業(yè)的技術(shù)人員應(yīng)該了解，本發(fā)明不受上述實(shí)施例的限制，上述實(shí)施例和說明書中描述的只是說明本發(fā)明的原理，在不脫離本發(fā)明精神和范圍的前提下，本發(fā)明還會(huì)有各種變化和改進(jìn)，這些變化和改進(jìn)都落入要求保護(hù)的本發(fā)明范圍內(nèi)。本發(fā)明要求保護(hù)范圍由所附的權(quán)利要求書及其等效物界定。