一種基于dem的自然水系自動提取方法

一種基于dem的自然水系自動提取方法

0確定高程數(shù)據(jù)點(diǎn)由于dem數(shù)據(jù)能夠反映特定分辨率的局部地形特征，因此可以根據(jù)地形的局部特征使用一定的算法來自動提取特定地理空間中的自然水系。算法方法如下：首先，根據(jù)水流沿斜面最陡方向的原理，確定dem中每個(gè)高數(shù)據(jù)點(diǎn)的水流方向。然后，根據(jù)高距離數(shù)據(jù)點(diǎn)的水流方向數(shù)據(jù)，計(jì)算每個(gè)高距離數(shù)據(jù)點(diǎn)的上游供水區(qū)，然后根據(jù)上游供水區(qū)的高距離數(shù)據(jù)，用閾值法確定水系的高距離數(shù)據(jù)點(diǎn)。最后，根據(jù)水流方向數(shù)據(jù)，從水系源中搜索整個(gè)水系。在國外的許多文獻(xiàn)中，討論了自然水系的自動提取過程，如yuan、smith、mats和lamper等。上述方法的適用條件是,在DEM中不能存在小平原、洼地等地形,且所有的復(fù)雜地形都必須由斜坡構(gòu)成,否則,所提取的水系通道是斷斷續(xù)續(xù)的.然而,在實(shí)際的DEM中,由于小平原和洼地的存在是普遍現(xiàn)象.為了使提取的自然水系,在流經(jīng)小平原和洼地部位時(shí),有一個(gè)明確的水系通道,因此在提取自然水系之前,需要對DEM中的小平原和洼地部位的高程數(shù)據(jù)進(jìn)行改造,以使洼地和小平原成為斜坡的延伸部分,經(jīng)過這樣的處理之后,DEM數(shù)據(jù)中的所有地形都由斜坡構(gòu)成.這樣才能保證從DEM數(shù)據(jù)中提取的自然水系是連續(xù)的.這種將DEM中的洼地和小平原改造成斜坡的處理過程稱為DEM數(shù)據(jù)的預(yù)處理.有關(guān)自然水系自動提取方面的研究資料在國內(nèi)尚不多見,為此,筆者總結(jié)了國外有關(guān)的研究資料,歸納出了一種從DEM中自動提取自然水系的實(shí)用方法,并用VB5.0編制了相應(yīng)的軟件系統(tǒng).該自然水系自動提取方法的工程流程如圖1所示.1斜坡預(yù)處理算法DEM數(shù)據(jù)預(yù)處理是為了將DEM中的洼地和小平原改造成斜坡,即,使DEM數(shù)據(jù)反映的地形特征均由斜坡構(gòu)成.預(yù)處理算法包括:①檢測DEM中是否有滲坑的算法;②洼地標(biāo)定算法;③洼地決口算法;④洼地填充算法;⑤平地起伏算法.1.1檢測dem中的滲坑根據(jù)水流總是沿斜坡最陡方向流動的原理,在DEM中,水流可以由相鄰柵格點(diǎn)流入某一柵格點(diǎn),再由該柵格點(diǎn)流向另一個(gè)相鄰的柵格點(diǎn).如果水流流入某柵格點(diǎn)后,不能流向其他相鄰柵格點(diǎn),那么該點(diǎn)即為DEM中的滲坑.可以使用如圖2所示的3×3模板來檢測DEM中的滲坑,即將模板沿掃描線方向移動,如果模板中心柵格值(指高程值)滿足條件:p0≤pk(k=1,2,…,8)則,p0即為滲坑.將以上操作針對DEM中的所有掃描線逐條進(jìn)行,即可以檢測出DEM中的所有滲坑.檢測DEM中是否有滲坑算法的VB偽碼如下:[用3×3模板檢測n行m列的DEM數(shù)據(jù)中的滲坑]1.2選取帶標(biāo)記的柵格點(diǎn)由于洼地(或小平原)總是圍繞滲坑分布,因此檢測出一個(gè)滲坑后,即可以通過以下算法確定洼地(或小平原)的范圍,其算法步驟如下:步驟1標(biāo)記滲坑(i0,j0),并以(i0,j0)為中心建立一個(gè)5×5窗口;步驟2用如下算法標(biāo)記窗口內(nèi)的滲坑供給區(qū):掃描窗口內(nèi)的所有柵格點(diǎn),將滿足下列條件的柵格點(diǎn)也標(biāo)上滲坑標(biāo)記:①與帶標(biāo)記的柵格點(diǎn)相鄰;②高程值不低于相鄰的帶標(biāo)記的柵格點(diǎn).該過程反復(fù)進(jìn)行,直到窗口內(nèi)不再有需要標(biāo)記的柵格點(diǎn)為止;步驟3用以下算法尋找洼地出水口:掃描窗口內(nèi)的所有柵格點(diǎn),將滿足下列條件的柵格點(diǎn)作為潛在出水口記錄下來:①柵格點(diǎn)帶有標(biāo)記;②與不帶標(biāo)記的柵格點(diǎn)相鄰;③高程值高于與其相鄰的不帶標(biāo)記的柵格點(diǎn).如果窗口內(nèi)不存在滿足條件①、②和③的柵格點(diǎn)或滿足條件①、②和③的柵格點(diǎn)其高程值高于滿足條件①和②,而不滿足條件③的柵格點(diǎn),則擴(kuò)大窗口,并重復(fù)步驟2和步驟3.步驟4當(dāng)窗口的大小可以覆蓋整個(gè)滲坑供給區(qū)時(shí),即可以先找到洼地的最低出水口;然后,檢查窗口內(nèi)所有帶標(biāo)記的柵格點(diǎn),如果存在高程值低于最低出水口高程值的柵格點(diǎn),則由帶有標(biāo)記,且高程值小于或等于最低出水口高程值的柵格點(diǎn)來構(gòu)成DEM中的封閉洼地;如果帶有標(biāo)記的柵格點(diǎn)中不存在高程值低于最低出水口高程值的柵格點(diǎn),則由高程值等于最低出水口的所有帶標(biāo)記的柵格點(diǎn)來構(gòu)成DEM中的小平原.1.3地面控制系統(tǒng)洼地決口算法主要用來處理封閉洼地.該算法通過模擬水流漫過水系通道中的狹窄障礙物,并最終在障礙物中侵蝕出一個(gè)缺口的過程來處理DEM中的洼地.這里,障礙物的厚度稱為決口長度,且一般決口長度小于或等于兩個(gè)DEM柵格.洼地決口算法的基本步驟如下:步驟1用以下算法尋找封閉洼地內(nèi)的潛在決口點(diǎn).掃描滲坑供給區(qū)內(nèi)的所有柵格點(diǎn),滿足下列條件的柵格點(diǎn)為潛在決口點(diǎn):①柵格點(diǎn)帶有標(biāo)記;②高程與最低出水口高程相同;③與不帶標(biāo)記的柵格點(diǎn)相鄰;④高程值高于與其相鄰的不帶標(biāo)記的柵格點(diǎn)高程;⑤高程值高于與其相鄰的洼地內(nèi)側(cè)的柵格點(diǎn)高程.如果找不到滿足上述條件的柵格點(diǎn),則轉(zhuǎn)移到洼地填充算法;如果找到滿足上述條件的柵格點(diǎn),則進(jìn)行步驟2.步驟2修改決口點(diǎn)的高程值,將決口點(diǎn)的高程值改為與其相鄰的不帶標(biāo)記的柵格點(diǎn)的高程值,即把決口點(diǎn)的高程值降低,使決口點(diǎn)成為封閉洼地的新的最低出水口.1.4出水口高程掃描滲坑供給區(qū)內(nèi)的所有柵格點(diǎn),將滿足下列條件的柵格點(diǎn)的高程值改為最低出水口的高程值:①帶有洼地標(biāo)記;②高程值低于最低出水口高程.通過上述的高程值修改步驟,滲坑供給區(qū)內(nèi)所有柵格點(diǎn)的高程值均大于或等于最低出水口高程.將滲坑供給區(qū)內(nèi)所有高程值等于最低出水口高程的柵格點(diǎn)標(biāo)上小平原標(biāo)記.1.5小原生的柵格點(diǎn)標(biāo)記平地起伏算法是利用在小平原上附加高程值的方法人為地將小平原改造成斜坡.人為的地形起伏是在原高程值基礎(chǔ)上附加小的高程增量而產(chǎn)生的.平地起伏算法的基本步驟如下:步驟1確定小平原邊界上不需要附加高程增量的柵格點(diǎn),并將其小平原標(biāo)記去掉:掃描滲坑供給區(qū)內(nèi)的所有柵格點(diǎn),將滿足下列條件的柵格點(diǎn)的小平原標(biāo)記去掉:①帶有小平原標(biāo)記;②與其相鄰的柵格點(diǎn)中,存在不帶小平原標(biāo)記的柵格點(diǎn);③相鄰柵格中,不帶小平原標(biāo)記的柵格點(diǎn)的高程值低于帶小平原標(biāo)記的柵格點(diǎn)高程.步驟2重新掃描滲坑供給區(qū)內(nèi)的所有柵格點(diǎn),將帶有小平原標(biāo)記的柵格點(diǎn)的高程值增加一個(gè)小的高程增量.返回到步驟1.重復(fù)步驟1和步驟2,直到滲坑供給區(qū)內(nèi)不存在帶有小平原標(biāo)記的柵格點(diǎn)為止.2水系提取自然水系自動提取的基本過程分以下3步:①計(jì)算DEM中每一個(gè)柵格點(diǎn)的水流方向;②計(jì)算每一個(gè)柵格點(diǎn)上游給水區(qū);③追索水系.2.1柵格點(diǎn)水流方向計(jì)算DEM中每一個(gè)柵格點(diǎn)的水流方向用與該柵格點(diǎn)相鄰的8個(gè)柵格點(diǎn)(即中心點(diǎn)四周8個(gè)方向上相鄰的柵格點(diǎn))中的某一個(gè)整數(shù)來表示.例如,當(dāng)水流方向數(shù)據(jù)為1時(shí),表示水流將由p0流向p1,此時(shí)則使用如圖2所示的3×3模板來計(jì)算柵格點(diǎn)水流方向,其算法的VB偽碼如下:[用3×3模板計(jì)算n行m列的DEM中柵格點(diǎn)的水流方向](注:程序中的“#”表示其前面的數(shù)為實(shí)數(shù),“*”為乘號)2.2柵格點(diǎn)上游給水區(qū)的確定根據(jù)DEM中柵格點(diǎn)的水流方向數(shù)據(jù),可以計(jì)算出DEM中每一個(gè)柵格點(diǎn)的上游給水區(qū).這里將一個(gè)柵格點(diǎn)的上游給水區(qū)定義為水流直接或間接地流向該柵格點(diǎn)的柵格點(diǎn)數(shù)目.計(jì)算DEM中每一個(gè)柵格點(diǎn)上游給水區(qū)算法的VB偽碼如下:[用3×3模板計(jì)算n行m列的DEM中柵格點(diǎn)的上游給水區(qū)]2.3自然水系的追求自然水系的追索過程分以下3個(gè)步驟:①用閾值法確定水系柵格點(diǎn);②篩選水系柵格點(diǎn);③進(jìn)行水系的追索與分級.2.3.1適當(dāng)選擇上游水由于DEM中某一個(gè)柵格點(diǎn)若能夠形成水系,則必須存在一定規(guī)模的上游給水區(qū),因此可以根據(jù)研究區(qū)氣候特征適當(dāng)選擇上游給水區(qū)高程閾值,并將柵格點(diǎn)上游給水區(qū)高程數(shù)值大于或等于給水區(qū)高程閾值的柵格點(diǎn)標(biāo)上水系柵格點(diǎn)的標(biāo)記.2.3.2追索水系通道篩選水系柵格點(diǎn)算法的基本步驟如下:步驟1找出DEM中的所有水系起始柵格點(diǎn).該水系起始柵格點(diǎn)滿足下列條件:①帶有水系柵格點(diǎn)標(biāo)記;②相鄰的水系柵格點(diǎn)水流方向均不指向它.步驟2從水系起始柵格點(diǎn)開始,追索水系通道,即根據(jù)水系柵格點(diǎn)的水流方向數(shù)據(jù),向下游逐步追索水系,并計(jì)算水系的累計(jì)長度;步驟3確定是否已追索到了水系交匯點(diǎn).每追索出一個(gè)水系柵格點(diǎn),均需檢查所有相鄰柵格點(diǎn),如果相鄰柵格點(diǎn)中有兩個(gè)或更多的水系柵格點(diǎn)的水流方向指向該水系柵格點(diǎn),則說明當(dāng)前柵格點(diǎn)為水系交匯點(diǎn),并將算法轉(zhuǎn)移到步驟4;步驟4根據(jù)追索出的水系通道的累積長度,決定該段水系是否需要刪除.如果累計(jì)長度大于預(yù)先給定的長度閾值,則不做任何操作,否則,將已經(jīng)追索出的屬于該水系的柵格點(diǎn)從水系柵格點(diǎn)中剔除.以上操作將針對每一個(gè)水系起始點(diǎn)逐個(gè)進(jìn)行,當(dāng)所有水系起始點(diǎn)均檢測完后,即結(jié)束水系柵格點(diǎn)的篩選過程.2.3.3建立健全水系分布及追索過程可以按Strahler的水系系統(tǒng)分級方案對整個(gè)水系系統(tǒng)進(jìn)行追索和分級.追索水系算法的基本步驟如下:步驟1追索Ⅰ級水系從水系起始點(diǎn)開始,根據(jù)水系柵格點(diǎn)的水流方向數(shù)據(jù),向下游逐步追索Ⅰ級水系,每追索出一個(gè)水系柵格點(diǎn),就將Strahler水系分級數(shù)據(jù)“1”賦給該柵格點(diǎn),當(dāng)追索到水系交匯點(diǎn)時(shí),交匯點(diǎn)的水系分級數(shù)據(jù)增加“1”,同時(shí),結(jié)束該段水系的追索.以上追索過程需針對每一個(gè)水系起始點(diǎn)逐個(gè)進(jìn)行.步驟2追索Ⅱ級及更高級水系從兩個(gè)Ⅰ級水系的交匯點(diǎn)開始,根據(jù)水系柵格點(diǎn)的水流方向數(shù)據(jù),向下游追索Ⅱ級水系,并將Strahler水系分級數(shù)據(jù)“2”賦給每一個(gè)追索出來的水系柵格點(diǎn).當(dāng)追索到水系交匯點(diǎn)時(shí),檢查與之交匯的水系是否為低級水系,如果是,就繼續(xù)向下游追索;否則,就結(jié)束該段水系的追索,并將當(dāng)前柵格點(diǎn)的水系分級數(shù)據(jù)增加“1”.以上追索過程需針對每一個(gè)Ⅱ級水系的起始點(diǎn)逐個(gè)進(jìn)行.重復(fù)與步驟2相似的追索過程,就可以追索出所有高級別的水系,當(dāng)所有的水系柵格點(diǎn)均被追索出來后,則結(jié)束水系的追索與分級過程.3dem數(shù)據(jù)的預(yù)處理筆者用VB5.0開發(fā)了上述算法的軟件系統(tǒng),并設(shè)計(jì)了軟件與MapInfo平臺的數(shù)據(jù)接口.該程序能自動將提取的自然水系按段編號、分級,并能計(jì)算各段長度,且每段的空間信息(河流上各個(gè)結(jié)點(diǎn)的坐標(biāo))和屬性信息(編號、級別、長度)均被寫入MapInfoInterchange文件中.為了測試上述算法的有效性,筆者選用小范圍內(nèi)的真實(shí)DEM數(shù)據(jù)進(jìn)行了試算.該DEM數(shù)據(jù)為54行×81列的點(diǎn)陣式數(shù)字高程數(shù)據(jù),采自我國1∶20萬敦化縣幅(編號:11-52-(3))地形圖東南角.DEM數(shù)據(jù)覆蓋面積為32.4km×21.6km.在地形圖上,數(shù)據(jù)網(wǎng)格大小為2.0mm×2.0mm,對應(yīng)的實(shí)際范圍是0.4km×0.4km的小網(wǎng)格.在DEM中,缺失的高程數(shù)據(jù)用1.70141×1038(用此數(shù)值充填,可使DEM數(shù)據(jù)格式與Surfer的輸入數(shù)據(jù)格式一致)充填.根據(jù)預(yù)處理后的DEM數(shù)據(jù),用本文所述算法來自動提取自然水系.圖3(a)為地形圖灰度圖象(*.srf),圖象中的黑線表示的是實(shí)際水系(以地形圖灰度圖象為底圖,在MapInfo中用“手工”方法繪制);圖3(b)為自動提取的自然水系圖(已做圓滑處理).在自然水系自動提取算法中,確定水系柵格點(diǎn)的上游給水區(qū)閾值為dq=30,水系