基于DEM的流域河網信息自動提取算法

劉永和張萬昌（ 河南理工大學資環學院，焦作 454000； 中國科學院遙感與數字地球研究所，北京 00094）

基于DEM的流域河網信息自動提取算法

劉永和1張萬昌2
（1 河南理工大學資環學院，焦作 454000；2 中國科學院遙感與數字地球研究所，北京 100094）

分布式陸面水文過程的模擬，除了經典水文模型所必需的流域氣象、水文信息外，還需要研究流域的詳細地形、水系信息，方可實現流域內產流和匯流的時空演算。以往通常需要借助商業軟件在研究區DEM上提取這些水系信息，不僅耗時、操作不便，還使得模型結構松散。在自主研發的一種分布式水文模型的基礎上，開發了一套用于提取水系信息的模塊，并與該水文模型以同一種語言緊密耦合，程序完全采用面向對象的方式自主編寫，大部分數據交換均在內存中直接進行，而無需占用磁盤空間，運行速度快，易于今后不斷完善和擴展。詳細介紹了模塊中填洼與平坦區域處理、流向與累積流向矩陣生成、Strahler河道等級的確定、子流域生成、匯流次序和流程長度信息的生成等多種流域河網信息提取的具體算法及程序實現。本模型系統完全采用自主方式開發，克服了以往使用商業軟件提取河網信息時的限制，使得分布式水文模式的流域模擬及分析功能更強大齊全和今后進一步完善和擴展。

分布式水文模型，數字高程模型，水系提取，耦合模塊

0　引言

模擬河流和土壤中水流的分布式水文模型被廣泛應用于水利工程及環境影響評價[1-5]。分布式模型可用于估計與水流有關屬性的空間變化，如化學物質遷移、土壤侵蝕、水分循環的陸面過程以及洪水狀況。由于流域水文狀況受天氣和氣候條件影響，人們期望在實現數值天氣預報和氣候變化預估的同時，能夠實現流域的水文預報或未來水文條件預估。由于水流嚴格地受制于流域的地形特征，只有預先提取出流域的水系信息，才能實現分布式模擬計算地表徑流和壤中流的匯流。基于數字高程模型（DEM）提取流域水系是目前一種較為高效且可靠的手段。目前通過現代雷達測量手段已能獲得高精度的地面高程信息，如美國地質調查局發布的SRTM DEM的精度已達到90m的水平分辨率，利用該數據能夠較為精準地提取水系河網信息。因此，從DEM中由計算機自動提取匯水網絡成為水文模擬研究中廣泛應用的手段[6]。

通常由DEM提取的水系信息包括像元間的流向、河段以及對應的子流域[7-8]。有關的研究已有很多文獻發表[3,7,9-13]。基于矩形規則格網的DEM是流域提取時應用最多的數據[7]，但在一些文獻中也會基于不規則三角網（TIN）DEM[14-15]。本文主要指基于規則格網結構的DEM。采用D8算法提取水系已經接近成熟，具體細節見文獻[6, 16]。目前在一些成熟的商業軟件中已集成了有關基于D8方法的水系提取工具，如ArcGIS中的空間分析工具、ENVI軟件的插件River Tools等。然而基于商業軟件的水系提取工具其源代碼不公開，人們無法知道其算法是否足夠高效且健壯，更重要的是商業軟件僅能作為一種數據處理工具，用戶無法對其處理過程和結果數據進行自由操控，限制了它在自主開發的分布式水文模型中的應用。 分布式水文模型通常表現為一種執行水文計算任務的計算機程序，由于缺少良好的人機交互界面而難以使用。目前一些重要的分布式水文模型（如SWAT）總是被集成于ArcGIS、GRASS等GIS軟件工具中，以充分借助GIS軟件所具有的DEM空間分析和可視化功能，增強其易用性。然而這種軟件集成方法對于分布式水文模型的研究者而言難度較大，影響了分布式水文模型的應用。因此，開發一套專門的水系信息提取工具，以解決當分布式水文模型應用中存在的問題和不足具有一定的實際意義。

基于上述原因，我們針對分布式水文模型自主開發了一套用于提供水系信息的耦合模塊。由于所建的分布式水文模型及水系提取模塊采用基于同一種語言（C#）的緊密耦合方式，大部分數據交換均在內存中直接進行，而無需占用磁盤空間，且運行速度快；程序完全采用面向對象的方式自主編寫，便于建立批量處理執行任務，易于不斷完善和擴展。

1　耦合系統的功能與耦合方案

1.1系統的功能

本文所述的流域河網信息自動提取系統及其與分布式水文模型的耦合系統具有以下三類功能：

1）類似于GIS軟件中各種位圖和簡單矢量圖顯示與漫游等可視化的功能。考慮到水文模型的主要任務是進行流域水文模擬，本系統既具有可顯示單幅DEM或其他地理資料矩陣的功能，也可任意打開并顯示用ARCGIS里所指定的ASCII二維數據格式或其他用ASCII碼矩陣方式存儲的二維數據。當光標在地圖上移動時，系統可自動計算并顯示當前光標處像元的值、行列號、經緯度、UTM坐標值。本系統可以作為查看任意兼容格式數據內容的基本工具，如查看分布式水文模擬的降水、蒸發、徑流的空間分布狀況。為使顯示的DEM底圖或其他屬性圖像與矢量地圖相聯系，該系統還提供了設置和顯示點類型和線類型矢量圖的功能。矢量圖顯示的數據源為ARCGIS軟件導出的文本矢量格式。

2）自動提取大量流域河網信息的功能。本系統中的流域河網提取模塊具有計算流向矩陣、累積流向矩陣、水系與河道等級矩陣、流域邊界、子流域劃分、匯流等級矩陣及匯流次序表、流程長度矩陣和匯流時間矩陣等功能，為分布式水文模擬提供必要的信息。

3）分布式水文模型的模擬功能。本系統還具有按逐日和逐時進行徑流和洪水模擬的功能，可按馬斯京根法實現河道洪水演算以及按等流時線滯時方法進行匯流。

1.2系統的耦合方案

水文模型在運行過程中所需要的氣象、土壤、植被、地形等空間分布數據均采用與所用DEM完全一致的分辨率。為方便使用該系統，在運行過程中系統將以程序配置文件中所記錄的河流出口斷面所在的像元位置（行號和列號）作為運行起始點，一次性完成所有河道信息提取過程。水文模擬過程中所需的所有河道信息以公用數據結構變量的形式保存在內存中，供水文模型中的產匯流模塊使用。

2　流域河網信息提取的算法

文中對流域河道信息的提取采用基于D8方法的思想。在每個DEM柵格上，待提取的信息包括流向、累積流向、水系與河道等級、流域邊界、子流域信息、匯流次序、匯流時間、流程長度矩陣。在計算流向前，需要對DEM進行填洼（depression filling）和處理平坦區域，以便能夠使DEM每個像元通過流向與主河道連接起來。

2.1填洼和平坦區域的處理

DEM中的洼地是指局部位置的像元高程低于外圍所有像元，使計算出的流向呈內流現象，即流向鏈條不能與外圍的主河道相連接。DEM上的洼地可能是地表上真實的洼地，如小的湖泊，也可能是由于測量誤差以及分辨率較低所造成的假性洼地。在水文研究中，無論流域中真實存在的局部洼地，還是假洼地，都被視作流域中排水網絡的組成部分，即必須通過流向將這些區域與水系主河道相連接。填洼曾經是DEM水系提取中最難解決的問題，有關處理該問題的算法較少，用傳統的一些算法（如Janson等[11]在1988年提出的方法）對一幅1000行×1000列的DEM進行填洼需花費數小時的時間，且因易出現錯誤而使運行不穩定，但Planchon等2002年提出的快速算法已解決了該問題[12]，本文使用了這種新算法，一般僅需數秒的時間即能完成，且運行穩定。

平坦區域中，由于所有像元的高程都相等，不存在任何流向，阻礙了水系的生成。目前常用做法是用不同程度的低于分辨率微小數值抬高平坦區域中像元的高程，使抬高的像元與其鄰接像元存在高程差，便于計算流向，并通過流向與外側水系相連接。

2.2水系生成與河道等級的確定

生成水系時需要完成流向計算、流向累積矩陣計算和確定河道等級。

在本系統中，我們采用一種新的鄰域位置編碼方法（圖1）。在這種編碼方法中，編碼值排列為三行。從左到右，第一行中的編碼設為0、1、2，第二行中的編碼設為3、4、5，第三行的編碼設為6、7、8。需要注意的是，在被設為4的3×3的編碼窗口中心的編碼不代表任何鄰域位置，而是代表了當前指定像元的位置。由這些方向編碼很容易算出當前指定像元的鄰域像元的行列號。假定當前像元的位置為（i，j），該像元處的水流流向編碼為m，則其下游鄰接像元的行列號為（i+m/3-1，j+m%3-1），其中符號%表示取余運算。如果一個流向編碼為n的鄰域像元的流向指向本像元，則其流向編碼應為8-n。可見，這種流向編碼方法非常有用，且很直接，這是因為當前感興趣像元的流出和流入像元的位置能夠高效計算出來。

本系統中采用單流向算法來計算流向，算法步驟是：（1）首先對DEM邊界上像素的流向設為指向DEM外部；（2）對除邊界外的所有內部像素，比較其與周圍8個鄰接像元的高程差構成的梯度，選擇高程低于當前像素且梯度最大的鄰接像元所在的方向作為流向。這里的梯度是以當前像元與鄰接像元的高程差除以像元距離D來計算。對上左下右四個緊鄰像元，像元距離D為1，而對角線上的四個像元，像元距離D為。

流向累積矩陣中每個像元上的累積流向值被定義為流過該像元的所有上游像元的總數，它代表了流過指定像元的上游匯水面積。常用的統計當前像元累積流向的方法有兩種，一種是遞歸找出當前像元的上游像元，統計這些找到的像元總數，即為累積流向值。另一種方法是逐個遍歷DEM中的每個像元，當遍歷到某個像元時，該像元的所有下游像元的累積流向值加1。當DEM中的像元全部遍歷一遍時，便得到了累積流向矩陣。

在得到流向矩陣和累積流向矩陣的基礎上，可進一步對河流進行基于Strahler分級方法[17]的水系分級，即把位于頂端的不再有分支的河流稱為第1級河流，由兩個以上的同級河流匯合組成更高一級河流，如果多條河流同時匯聚，但最高級別的支流只有一條，則匯合后的河流級別與其最高級別的支流的級別相同。文獻[18]對Strahler河流分級方法進行了較詳細的描述，采用從最上游的像素開始向下游搜索，較為可行，但可能有一定的搜索難度。而本文則提出了另一種不需要從最上游像素開始搜索的方法。具體步驟如下：

1）給定一個累積流向（匯流面積）閾值，認為大于該閾值的像元位于河流上，不大于該閾值的不屬于河流；

2）在DEM中找出所有累積流向值（匯流面積）大于指定閾值的所有像元，即找出所有位于河網上的像元。為便于快速檢索，將這些河網像元放入一個以行號和列號為關鍵字的字典集合中；

3）對字典集合中的所有河網像元逐個遍歷，對找到的每個像元執行第（4）～（8）步操作；

4）若當前像元（用C1表示）已被分配子流域標志，則返回第（3）步，即跳過當前像元，對遍歷到的下一個像元進行處理；

5）從當前像元的8個鄰接像元中找出屬于C1上游的所有像元，若沒有找到這樣的像元，則將C1的河流級別設為1（屬初級河流）并返回第（3）步執行下一趟循環。若找到了這樣的像元，執行后續步驟；

6）對C1的上游像元，統計最大河流級別值m，以及擁有最大河流級別值的像元個數n，并判斷是否C1的全部上游像元已被標記過河流級別，如果已被標記則執行后續步驟，否則返回第（3）步執行下一趟循環；

7）如果在前面已找出的C1的上游像元個數為1時，則C1的河流級別設為與其上游像元相同的值；如果C1的上游像元多于1個時，且有兩個或更多像元級別具有相同的上游最大級別m時，C1的河流級別設為m+1；如果C1的上游像元多于1個，但只有一個上游像元具有最大級別m，則C1的河流級別設為m；

8）返回第（3）步；

9）如果第（2）～（8）步中沒有任何像元被標記過，則退出；否則從字典集合中刪除已標記過的像素記錄，并返回第（2）步繼續執行；

10）結束后退出。

2.3子流域劃分

子流域一般是指流域面積超過一定面積的分支流域。文獻[19]中提出了一種子流域劃分方法，但算法不夠具體。文獻[20]中提出了具體的算法，但略顯復雜。本文采用思路與文獻[20]基本一樣，但采用下面的算法來完成：

1）生成一個空隊列Q，將代表河流出口處的像元設為子流域最低編號（如0或1），并將該像元的記錄加入隊列；

2）若隊列不空，從隊列中推出一個像元C1，執行后續步驟，否則退出執行；

3）從像元C1的8個鄰接方向中逐個查找C1的上游像元；

4）若找到C1的一個上游像元為C2，將C2加入隊列；

5）判斷C2的累積流向（流域面積）值大于給定閾值，且C2與C1之間的累積流向值之差也大于給定閾值的條件是否成立，若成立則將C2標記為新子流域（其標記值為C1的標記值加1），若不成立，則C2的標記與C1相同；

6）返回第（2）步繼續執行。

上面算法第（5）步中的條件十分關鍵。C2作為C1的上游像元，其上游流域面積必須大于閾值才能成為一個子流域。另外，還需要保證C1存在除C2外的較大分支才能使C2成為新子流域的像元，這是因為，如果除C2外的其他分支的匯流面積過小，小分支無法獨立成為子流域，而與C2、C1屬同一子流域。步驟（5）中的條件“C2與C1之間的累積流向值之差也大于給定閾值”正是為了保證C1存在除C2外的較大分支。

2.4匯流等級矩陣、匯流次序表、流程長度與匯流時間矩陣的生成

在分布式水文模型的洪水模擬中，需要根據流域中每個像元間的流向關系從上游向下游采用馬斯京根方法進行河道匯流。僅僅根據流向矩陣無法得出正確的匯流結果，仍需要根據匯流路徑的遠近對所有柵格劃分匯流順序，從而得到整個流域或子流域的匯流等級矩陣和匯流次序表。對于同樣屬最低級別的河流像元，模擬時距離流域出口像元遠的像元應該先計算，距離流域出口近的像元應后計算。匯流等級矩陣中的元素值為理論上根據到流域出口的距離來確定的匯流次序，等級值小的像元要先計算匯流，相同匯流等級的多個像元應同時計算匯流。但從實際計算的角度來看，匯流等級相同的多個像元仍是分別來計算匯流的，但計算順序可是任意的。為方便，一般要導出匯流次序表，用以代替匯流等級矩陣。匯流次序表以匯流的實際計算次序存放流域中所有的像元的行列號，而不是采用矩陣形式存放。

流程長度矩陣是與流域DEM完全對應的二維矩陣（數組），其每個元素存放流域中對應像元通過流向拓撲關系確定的通向流域出口的空間距離。匯流時間矩陣與流程長度矩陣類似，只是其每個元素存放的是對應像元的“水流”流往流域出口所花費的時間。匯流時間矩陣可以通過將流程長度矩陣直接除以平均流速粗略得到，也可以通過將各個像元的水流流速與相鄰上下游像元間的梯度建立關系來精確導出。

匯流次序表可以通過匯流等級矩陣或流向矩陣導出。從流向矩陣導出的方法是從河流出口逆向上溯的方法搜索。由于流程長度矩陣的計算以及通過梯度精確計算匯流時間矩陣時也需要采用同樣的逆向上溯方法來計算，可將三個計算任務用同一個算法來完成。算法仍借助隊列來實現，具體步驟如下：

1）將流域出口像元入隊列Q，并往匯流次序表中添加該像元，對應該像元的流程長度矩陣設為1（對匯流時間矩陣可采用類似的方法賦值）；

2）當隊列不空時，從中推出一個像元C，否則跳到第（6）步執行；

3）從8個鄰接像元中查找屬于C的上游流入像元；

4）對找到的所有上游流入像元的流程長度值進行賦值，對斜對角方向的流入像元，流程長度值為C的流程長度加1.414，對上下左右四方向的流入像元，流程長度值為C的流程長度加1。用類似的方法可對匯流時間矩陣進行賦值。把找到的上游流入像元按任意順序添加入匯流次序表的末尾；

5）返回第（2）步執行；

6）將匯流次序表反轉；退出。

3　其他地面信息的生成

除與水系有關的信息外，分布式水文模型仍需要借助地面坡度和坡向信息來計算匯流速度和水分蒸散發量。計算坡度和坡向的方法可參見文獻[21]。

為使流域位置上DEM的地形狀況更加直觀，可生成漫反射光照明暗圖。在計算機圖形學中，漫反射光照采用Lambert余弦定理來實現[22]。Lambert余弦定理是根據入射光向量與反射面的法向來計算的，因此先要計算出DEM中的每個像元所處位置上的法向。由幾何學知道，空間中的平面法向為通過平面內的任意兩個向量的向量積。顯然，只要已知平面中的三點，就能確定該平面的法向，因此三角形是最適于光照處理的面片。而規則格網DEM中每相鄰四個像元都能拼成兩個三角形（圖2）。需要指出的是，在拼三角形時需要保證三角形的頂點順序都為逆時針或都為順時針。設三角形的三個頂點分別為則三角形某兩邊上的矢量

當每個DEM像元的三角形法向被算出后，再給定一束來自無窮遠處的虛擬入射光，即可算出像元上的光強。

4　流域河網信息提取實例

本文截取了山東省境內的沂河流域范圍的DEM數據，分辨率為240m，在此基礎上生成水文模擬時所用的各種與流域和河網有關的信息。圖3顯示的是系統的主窗體，其主要工作區中顯示了DEM圖像以及該區域內的河流、縣界及水文測站的矢量信息。圖4顯示了DEM原圖及由其導出的Strahler河流等級、子流域分區、流程長度、匯流次序等信息的圖像。圖4b顯示的由DEM導出的河網與現有的矢量主河道基本一致，這表明流域河網信息提取模塊中的填洼及計算累積流向的算法是較為可靠的。圖4c顯示了以沂河流域內臨沂水文站為出口點的面積閾值超過2000像素的子流域劃分，由圖可見，劃分的子流域基本合理，但仍存在不同子流域面積差別較大的問題。圖4d顯示了由DEM導出的沂河流域內各像素至河流出口處的流程長度，圖4e顯示了由DEM導出的各像元的匯流次序，這兩幅圖表明對應的算法計算結果是正確的。圖4f顯示了由DEM直接生成的虛擬化光照明暗圖，該圖能夠反映研究區域地形起伏的狀況。

5　結語

本文介紹了在分布式水文模型中基于DEM提取流域河網信息的耦合模塊的算法設計以及耦合系統的功能和耦合方案，對水系生成和河道等級的確定、子流域生成、匯流次序及流程長度信息的生成以及光照明暗圖的生成等算法進行了詳細介紹。耦合模塊采用基于內部數據結構的緊密耦合方式與分布式水文模型耦合成一體化系統，對水文模型所需的各種流域河網信息采用自動批量式生成，節省了對磁盤空間的占用，提高了運行速度，同時也簡化了對水文模擬時流域信息的準備過程，能夠較大幅度的提高使用效率。耦合模塊中的各種算法結果基本合理或正確，運行速度快，十分穩定。

本水文模型系統的開發完全采用自主方式，完全避免了使用商業GIS軟件在基于DEM準備各種河網信息時的限制和不足，更有利于今后的進一步完善和擴展。

[1]Yang X J. Use of LIDAR elevation data to construct a highresolution digital terrain model for an estuarine marsh area. International Journal of Remote Sensing, 2005, 26(23): 5163-5166.

[2]Gumbo B, Munyamba N, Sithole G, et al. Coupling of digital elevation model and rainfall-runoff model in storm drainage network design. Physics and Chemistry of the Earth, 2002, 27(11-22): 755-764.

[3]Wang X H, Yin Z Y. A comparison of drainage networks derived from digital elevation models at two scales. Journal of Hydrology, 1998, 210(1-4): 221-241.

[4]Wolock D M, Mccabe G J. Differences in topographic characteristics computed from 100-and 1000-m resolution digital elevation model data. Hydrological Processes, 2000, 14(6): 987-1002.

[5]Dodov B A, Foufoula-Georgiou E. Floodplain morphometry extraction from a high-resolution digital elevation model: A simple algorithm for regional analysis studies. IEEE Geoscience and Remote Sensing Letters, 2006, 3(3): 410-413.

[6]Tribe A. Automated recognition of valley lines and drainage networks from grid digital elevation models - A review and a new method. Journal of Hydrology, 1992, 139(1-4): 263-293.

[7]Turcotte R, Fortin J P, Rousseau A N, et al. Determination of thedrainage structure of a watershed using a digital elevation model and a digital river and lake network.Journal of Hydrology, 2001, 240(3-4): 225-242.

[8]Paz A R, Collischonn W. River reach length and slope estimates for large-scale hydrological models based on a relatively hill highresolution digital elevation model. Journal of Hydrology, 2007, 343(3-4): 127-139.

[9]Colombo R, Vogt R V, Soille P, et al. Deriving river networks and catchments at the European scale from medium resolution digital elevation data. CATENA, 2007, 70(3): 296-305.

[10]Ahamed T, Rao K G, Murthy J. Automatic extraction of tank outlets in a sub-watershed using digital elevation models. Agricultural Water Management, 2002, 57(1): 1-10.

[11]Jenson S K, Domingue J O. Extracting topographic structure from digital elevation data for geographic information-system analysis. Photogrammetric Engineering and Remote Sensing, 1988, 54(11): 1593-1600.

[12]Planchon O, Darboux F. A fast, simple and versatile algorithm to fill the depressions of digital elevation models. CATENA, 2002, 46(2-3): 159-176.

[13]Jana R, Reshmidevi T V, Arun P S, et al. An enhanced technique in construction of the discrete drainage network from low-resolution spatial database. Computers & Geosciences, 2007, 33(6): 717-727.

[14]Tucker G E, Lancaster S T, Gasparini N M, et al. An objectoriented framework for distributed hydrologic and geomorphic modeling using triangulated irregular networks. Computers & Geosciences, 2001, 27(8): 959-973.

[15]劉學軍, 王永君, 任政, 等. 基于不規則三角網的河網提取算法.水利學報, 2008, 39(1): 27-34.

[16]Mcmaster K J. Effects of digital elevation model resolution on derived stream network positions. Water Resources Research, 2002, 38(4): 13-1-13-8.

[17]Strahler A N. Quantitative analysis of watershed geomorphology. Transactions AGU, 1957, 38(6): 913-920.

[18]張渭軍, 王文科, 孔金玲, 等. 基于數字高程模型的水系自動生成. 大地測量與地球動力學, 2006, 26(4): 41-44.

[19]沈中原, 李占斌, 李鵬, 等. 基于DEM的流域數字河網提取算法研究. 水資源與水工程學報, 2009, 20(1): 20-28.

[20]葉愛中, 夏軍, 王綱勝, 等. 基于數字高程模型的河網提取及子流域生成. 水利學報, 2005, 36(5): 531-537.

[21]湯國安, 劉學軍, 閭國年. 數字高程模型及地學分析的原理與方法. 北京: 科學出版社, 2005.

[22]唐澤圣, 周嘉玉, 李新友. 計算機圖形學. 北京: 清華大學出版社, 1994.

A DEM Based Drainage Networks Extraction Module Seamlessly Integrated in a Distributed Hydrological Model

Liu Yonghe1, Zhang Wanchang2
(1 Institute of Resources and Environment, Henan Polytechnic University, Jiaozuo 454000 2 Institute of Remote Sensing and Digital Earth (RADI), Chinese Academy of Sciences, Beijing 100094)

For simulation of the hydrology in a watershed by distributed hydrological models, besides the meteorological and hydrological data needed, the detailed information such as topography and watershed networks is also necessary. Traditionally, the drainage networks were derived from digital elevation models (DEM) using commercial software, which is time consuming and inconvenient for operation in hydrological modeling. Based on a distributed hydrological model (DHM) developed by us, an automatic watershed information extraction module was developed which is able to be integrated into the DHM seamlessly using the same computer language of C#. Therefore, most of the data transferring can be finished in memory with no occupation on hard disks, so it is highly efficient when running. The algorithms for the watershed information extraction were introduced, such as the algorithms of removing of depression and flatten areas, generation of flow direction and accumulating flow direction, obtaining Strahler river order, dividing of sub watershed and calculation the sequence of flow order and water-drainage length. This system overcomes the limitations of using commercial software for extracting watershed information, and also favors the modeling system for convenience of modifying and extending the model easily in the future.

distributed hydrological models, digital elevation models, drainage networks extraction, coupled module

10.3969/j.issn.2095-1973.2015.02.009

2013年12月12日；

2014年3月11日

劉永和（1976—），Email: sucksis@163.com

資助信息：國家自然科學基金項目（41105074和41275108）；中科院數字地球重點實驗室開放基金項目（2011LDE010）；河南理工大學博士基金（B2011-038）