一種用于數字地形分析建模知識案例化自動處理方法
【專利摘要】本發明涉及一種用于數字地形分析建模知識案例化自動處理方法，其中對數字地形分析(DTA)建模知識的案例化是將數字地形分析所應用的研究區定量表達為一個包含柵格分辨率、面積、地形起伏度、高程?坡度累積頻率二維表、面積高程積分值等5種定量指標的五元組，以柵格DEM為輸入數據，通過調用開源的柵格空間數據轉換庫GDAL的讀柵格數據功能以及地理信息系統軟件ArcGIS的部分地形因子計算模塊，實現對上述5個定量指標進行自動處理，所得到的數字地形分析建模知識案例可用于對某一研究區內已有數字地形分析建模知識在其他研究區中可用性的定量評價。本發明減少人為干預，能夠為DTA建模知識的案例化及可適用范圍推測提供技術支持。
【專利說明】
一種用于數字地形分析建模知識案例化自動處理方法
技術領域
[0001] 本發明涉及一種用于數字地形分析建模知識案例化的定量指標自動處理方法，屬 于基于案例推理的數字地形分析應用技術領域。
【背景技術】
[0002] 數字地形分析(DTA，Digital terrain analysis)是針對地表形態基于數字高程 模型（DEM，Digital Elevation Model)進行地形屬性計算及特征提取的數字信息處理技 術，DTA應用建模在自然地理學領域應用廣泛。然而，DTA應用建模過程較為復雜，涉及到大 量DTA專業知識。在實際應用中，算法的選擇和參數的設置等細節往往對整個模型的運行結 果有很大影響，需要根據區域和數據特征對模型進行配置，由于區域和數據特征本身較為 復雜，配置模型時需要使用大量的專業知識。由于這類知識難以形式化表達，現有的DTA工 具對此類知識缺乏利用，從而導致普通用戶在應用DTA時建模困難。
[0003] 基于案例的知識表達和推理方法適用于對非系統化的DTA專業知識進行形式化表 達和推理應用，能夠降低用戶在應用DTA時建模難度。基于案例的知識表達和推理方法的基 本思想是通過參考已有的相似案例的解決方案應用于新的問題。其中案例的知識表達是將 案例的屬性進行提取并描述的過程，它是進行案例推理的前提。一個完整的案例通常由3個 部分組成:案例問題、解決方案、案例輸出。其中案例問題是指案例發生時與案例有關因素 的狀態，該部分內容直接參與案例推理;案例的解決方案是作為案例推理的結果加以應用， 不參與案例推理;案例輸出是指案例發生后與案例相關因素的狀態，該部分內容是可選的， 不參與案例推理。
[0004] DTA建模知識的案例表達主要是對案例問題進行描述，即確定影響DTA算法選擇和 參數設置的主要因素（例如區域和數據特征等），進而確定適合于案例推理應用的定量指 標，以有效體現DTA應用建模相關的因素的狀態。在描述DTA建模知識的案例問題時，主要考 慮的因素有案例的數據特征和研究區特征，目前尚未有一套對相應的定量描述指標的自動 處理方法。
【發明內容】
[0005] 本發明的技術解決問題:克服現有技術的不足，提供一種用于數字地形分析建模 知識案例化自動處理方法，針對DTA建模知識的案例化表達，在案例指標確定的基礎上，通 過調用開源的柵格空間數據轉換庫GDAL的讀柵格數據功能以及地理信息系統軟件ArcGIS (其他地理信息系統軟件如GRASS GIS、SuperMap等也可以）的部分地形因子計算模塊，實現 案例定量描述指標的自動處理。
[0006] 本發明的技術解決方案為:一種用于數字地形分析建模知識案例化的定量指標自 動處理方法，原理是:通過調用GDAL庫的讀柵格功能以及ArcGIS的部分地形因子計算模塊， 自動獲取數字地形分析所應用研究區的柵格分辨率、面積、地形起伏度、高程-坡度累積頻 率表、面積-高程積分值，得到定量指標五元組，步驟如下：
[0007] (1)采用柵格DEM(Digital Elevation Model)為輸入數據，利用開源的柵格空間 數據轉換庫GDAL讀取柵格DEM數據，利用GetGeoTransform函數獲取柵格分辨率GridSize;
[0008] (2)利用GetNoDataValue函數獲取柵格DEM的Nodata值，遍歷整個柵格DEM數據，統 計非Nodata的柵格總數，用柵格總數乘以分辨的平方得到面積Area;
[0009] (3)遍歷整個柵格DEM數據，統計非Nodata柵格DEM數據的最大值和最小值，用最大 值減最小值得到地形起伏度數據Relief;
[0010] (4)采用柵格DEM(Digital Elevation Model)為輸入數據，利用Python腳本調用 地理信息系統軟件ArcGIS(其他地理信息系統軟件如GRASS GIS、SuperMap等也可以）計算 坡度，得到坡度數據;遍歷柵格DEM數據和坡度數據，統計各個高程等級里各坡度等級的頻 率，生成一張10X7的二維表;將各坡度等級在高程等級上進行累加運算，生成高程-坡度累 積頻率二維表。
[0011] (5)遍歷柵格DEM數據，統計每個高程值所對應的柵格數（由于面積與柵格數是成 比例的，為了計算方便用柵格數代替面積）；將高程值數組標準化(公式1)得到標準化后的 高程值數組;將柵格數數組標準化(公式2)得到標準化后的柵格數數組;遍歷標準化后的高 程值數組和柵格數數組，依次相乘并累加得到面積高程積分值Hypsolnt。
[0012] y = (H-Hmin)/Relief (公式 1)
[001 3] 式中，y表示面積高程積分曲線中高程所占比例;H表示DEM中的實際高程值;Hmir^ 示EffiM的最低高程值;Re 1 i ef?表示EffiM的高程差。
[0014] x = count/sumCount (公式 2)
[0015]式中，x表示面積高程積分曲線中面積所占比例；count表示某一高程值所對應的 柵格數;sumCount表示DEM的總柵格數。
[0016] (6)將上述5個步驟封裝成獨立函數，將得到的定量指標五元組輸出記錄，可實現 對數字地形分析建模知識案例化的自動化處理。
[0017] 所述步驟（1)中如何利用GetGeoTransform函數獲取柵格分辨率GridSize的過程 如下：
[0018]首先，注冊所有已知的驅動(這里也可以選擇只注冊指定驅動），然后將DEM文件數 據讀取到內存，調用GetGeoTransform函數獲取仿射變換系數矩陣，該系數矩陣的第二個元 素代表東西方向一個像素對應的距離，即水平分辨率，若不加特殊說明，一般所說的柵格分 辨率就是水平分辨率。
[0019] 所述步驟(4)中將各坡度等級在高程等級上進行累加運算，生成高程-坡度累積頻 率二維表的過程如下：在每個高程等級上，根據坡度從小到大進行累加，若累加前的高程-坡度頻率分布表上高程為i等級坡度為j等級的頻率為 aij，累加之后高程為i等級坡度為j等 級的累加頻率為bij:
[0020] bij= 2n=i.. jain
[0021] 式中，i = l，2,3,4,5,6,7,8,9，10;j = l，2,3,4,5,6,7;n 為整數。
[0022] 所述步驟(4)中的地理信息系統軟件包括ArcGIS、GRASS GIS、SuperMap。
[0023] 所述步驟(4)中生成一張10 X 7的二維表。
[0024]本發明與現有技術相比的優點在于:本發明對數字地形分析(DTA)建模知識的案 例化是將數字地形分析所應用的研究區定量表達為一個包含柵格分辨率、面積、地形起伏 度、高程-坡度累積頻率二維表、面積高程積分值等5種定量指標的五元組，以柵格數字高程 模型(DEM)為輸入數據，通過調用開源的柵格空間數據轉換庫GDAL的讀柵格數據功能以及 地理信息系統軟件ArcGIS的部分地形因子計算模塊，實現對上述5個定量指標進行自動處 理，所得到的數字地形分析建模知識案例可用于對某一研究區內已有數字地形分析建模知 識在其他研究區中可用性的定量評價。本發明的優勢是僅需要基本的DEM輸入數據，能夠自 動計算案例指標，完成相應研究區中數字地形分析建模知識的自動案例化，減少人為干預， 本發明能夠為DTA建模知識的案例化及可適用范圍推測提供技術支持。
【附圖說明】
[0025]圖1為本發明的實現流程圖；
[0026]圖2為以福建省朱溪河小流域柵格DEM為輸入數據、用Java語言編程輸出的結果截 圖。
【具體實施方式】
[0027] 在描述DTA建模知識的案例問題時，為定量刻畫案例的數據特征和研究區特征，采 用的定量指標分別為:DEM柵格分辨率(考慮到目前實際應用中主要采用柵格DEM，本說明目 前針對柵格DTA應用）、面積、地形起伏度、坡度分布、發育特征。其中DEM柵格分辨率用來刻 畫數據特征;面積、地形起伏度、坡度分布、發育特征這4個指標用來刻畫研究區特征。柵格 分辨率、面積、地形起伏度這3個指標采用直觀、常用的單一值描述;坡度和發育特征分別采 用高程-坡度累積頻率二維表和面積高程積分曲線進行描述。其中坡度的分級根據坡度對 土壤侵蝕的影響分7級（0~3°、3~8°、8~15°、15~25°、25~35°、35~45°、45~90°)，高程 按等間隔分為10個等級，生成一張10X7的二維表。由于面積高程積分曲線在案例推理過程 中是以面積高程積分值來計算的，所以其形式化是以面積高程積分值的形式來體現的。最 終得到一個五元組〈柵格分辨率、面積、地形起伏度、高程-坡度累積頻率二維表、面積高程 積分值〉，實現DTA建模知識的案例表達，通過該五元組表達的案例可用于對某一研究區內 已有數字地形分析建模知識在其他研究區中可用性的定量評價。其中的關鍵技術就是是對 柵格分辨率、面積、地形起伏度、高程-坡度累積頻率二維表、面積高程積分值這5個定量指 標的自動計算。
[0028] 下面對本發明采用的相關術語和含義等在此進行一下簡要說明。
[0029] GetGeoTransform函數是⑶AL的庫函數，返回的是一個指向六個元素的double類 型指針，其中六個元素的意義分別為：左上角橫坐標、東西方向一個像素對應的距離、旋轉 (0表示上面為正北）、左上角縱坐標、旋轉(0表示上面為正北）、南北方向一個像素對應的距 離。GetNoDataValue函數也是GDAL的庫函數，用于獲取數據的Nodata值。Nodata值指的是柵 格數據中"無意義值"，柵格數據通常為矩形，但其所表達的數據區域可能為一個不規則形 狀，且其中部分柵格不具有實際數值，這些情況下需要對一些有效數據區域之外的柵格填 充以"無意義值"，在實際應用中，無意義值通常不參與計算。
[0030] 下面結合如圖1所示的流程圖，以Java編程語言為例，說明本發明的具體實施方 法：
[0031] (1)配置好Java以及GDAL庫的運行環境。
[0032] (2)用Java語言編寫讀數據函數readRaster(DEMfile)(括號中表示該函數要輸入 的參數，下同）：調用GDAL庫中的gdal.open(DEMfile)函數將柵格數據讀取到全局變量 Dataset中。具體做法是:首先，注冊所有已知的驅動(這里也可以選擇只注冊指定驅動），然 后調用GDAL的庫函數open()函數將DEM數據讀取到內存。
[0033] (3)編寫獲取柵格分辨率函數getResolution(Dataset):調用 Dataset.GetGeoTransform()函數獲取柵格數據分辨率。具體做法是：調用 GetGeoTransformO函數獲取仿射變換系數矩陣，該系數矩陣的第二個元素代表水平分辨 率，返回該系數矩陣的第二個元素值，單位是米。
[0034] (4)編寫計算面積函數calculateArea(Dataset):利用GDAL庫函數 GetNoDataValue ()函數獲取柵格DEM的Nodata值，遍歷整個柵格DEM數據，統計非Nodata的 柵格總數，用柵格總數乘以分辨的平方得到面積Area:
[0035] Area = SumCount XGridSize2 X 10-6
[0036] 式中，SumCount表示DEM中非Nodata值的柵格數量，GridSize表示分辨率(單位是 米），面積單位是平方千米。
[0037] (5)編寫計算地形起伏度函數calculateRelief (Dataset):利用GDAL庫函數 GetNoDataValue ()函數獲取柵格DEM的Nodata值，遍歷整個柵格DEM數據，統計非Nodata柵 格DEM數據的最大值和最小值，用最大值減最小值得到地形起伏度數據Re 1 i ef:
[0038] Re1i e f = Vmax-Vm i n
[0039] 式中，Vmax表示DEM高程值的最大值，Vmin表示DEM高程值的最小值，單位均是米。
[0040] (6)用Python腳本編寫坡度計算程序calculateSlope ? py，調用arcpy里面的函數 計算坡度（需要安裝ArcGIS)。具體做法為：先調用arcpy .CheckOutExtension(〃Spatial〃） 函數獲取擴展模塊Spatial的許可，再調用arcpy ? Slope (DEMfile，"DEGREE"，1)函數計算坡 度，最后調用arcpy .save (SlopeFile)函數保存坡度數據。其中arcpy .Slope函數的第一個 參數是DEM原始數據文件，第二個參數"DEGREE"表示坡度的單位是角度(根據具體需要可以 改成弧度等），第三個參數默認值是1，表示DEM數據中x、y單位和z單位采用相同的測量單 位，DEM中x、y單位也就是分辨率的單位，z單位也就是高程單位，都是以米為單位，所以這里 采用默認值即可。
[0041 ] (7)編寫計算高程-坡度累積頻率二維表函數calculateDemSlp(Dataset):調用 calculateSlope. py腳本，計算坡度數據，調用讀數據函數readRaster讀取坡度數據，利用 GDAL庫函數GetNoDataValueO函數獲取柵格DEM和坡度數據的Nodata值，遍歷柵格DEM數據 和坡度數據，統計非Nodata數據中各個高程等級里各坡度等級的頻率aij，將各坡度等級在 高程等級上進行累加運算b 1J;
[0042] bij= 2n=i-jain
[0043] 式中，i = l，2,3,4,5,6,7,8,9，10;j = l，2,3,4,5,6,7;n 為整數。
[0044] (8)編寫計算面積-高程積分值函數calculateHypsoInt(Dataset):利用GDAL庫函 數GetNoDataValue ()函數獲取柵格DEM的Nodata值，遍歷柵格DEM數據，統計每個非Nodata 高程值所對應的柵格數，分別將高程值數組和柵格數數組標準化(根據公式1和公式2分別 對高程值數組和柵格數數組進行標準化），遍歷標準化后的高程值數組y和柵格數數組X，依 次相乘并累加；
[0045] HypsoInt= 2n=o-(len-i)xnyn
[0046] 式中，x，y代表標準化后的高程值數組和柵格數數組，Xn，yn代表數組中第n+1個元 素的值，其中n = 0，1，2，3，…，1 en-1，1 en代表數組的長度。
[0047] (9)在主函數中依次運行readRaster、getResolution、calculateArea、 calculateRelief、calculateDemSlp、calculateHypsoInt函數，設置輸入和輸出，運行得到 結果。
[0048]如圖2所示，以福建省朱溪河小流域柵格DEM為輸入數據，經過以上程序計算的到 的指標值分別為：分辨率5m，面積5.22km2，地形起伏度556m，面積-高程積分值0.21，高程-坡度累積頻率二維表：
[0050] 本發明說明書中未作詳細描述的內容屬于本領域專業技術人員公知的現有技術。
[0051] 所有上述僅是本發明的優選實施方式，應當指出，對于本技術領域的普通技術人 員來說，在不脫離本發明原理的前提下，還可以做出若干改進和潤飾，這些改進和潤飾也應 視作本發明的保護范圍。
【主權項】
1. 一種用于數字地形分析建模知識案例化自動處理方法，其特征在于:實現步驟如下： (1) 采用柵格DEM(Digital Elevation Model，數字高程模型）為輸入數據，采用開源的 柵格空間數據轉換庫GDAL讀取柵格DEM數據，利用GetGeoTransform函數獲取柵格分辨率 GridSize； (2) 利用GetNoDataValue函數獲取柵格DEM的Nodata值，遍歷整個柵格DEM數據，統計非 Nodata的柵格總數，用柵格總數乘以格分辨率GridSize得到面積Area; (3) 遍歷整個柵格DEM數據，統計非Nodata柵格中DEM數據的最大值和最小值，用最大值 減去最小值得到地形起伏度數據Relief; (4) 采用柵格DEM(Digital Elevation Model)為輸入數據，利用Python腳本調用地理 信息系統軟件計算坡度，得到坡度數據;遍歷柵格DEM數據和坡度數據，統計各個高程等級 里各坡度等級的頻率，生成一張二維表;將各坡度等級在高程等級上進行累加運算，生成高 程-坡度累積頻率二維表； (5) 遍歷柵格DEM數據，通過面積Area統計每個高程值所對應的柵格數，面積與柵格數 是成比例的，為了計算方便用柵格數代替面積;將高程值數組標準化公式1得到標準化后的 高程值數組;將柵格數數組標準化公式2得到標準化后的柵格數數組;遍歷標準化后的高程 值數組和柵格數數組，依次相乘并累加得到面積高程積分值HypsoInt， y= (H-Hmin)/Relief (1) 式中，y表示面積高程積分曲線中高程所占比例;H表示DEM中的實際高程值;Hmin表示 DEM的最低高程值;Relief表示DEM的高程差； X = count/SumCount (2) 式中，X表示面積高程積分曲線中面積所占比例；count表示某一高程值所對應的柵格 數;SumCount表示DEM的總柵格數； (6) 將上述5個步驟封裝成獨立函數，得到的定量指標五元組輸出記錄，實現對數字地 形分析建模知識案例化的自動化處理。2. 根據權利要求1所述的用于數字地形分析建模知識案例化自動處理方法，其特征在 于:所述步驟(1)中如何利用GetGeoTransform函數獲取柵格分辨率GridSize的過程如下： 首先，注冊所有已知的驅動，也可以選擇只注冊指定驅動，然后將DEM文件數據讀取到 內存，調用GetGeoTransfom函數獲取仿射變換系數矩陣，該系數矩陣的第二個元素代表東 西方向一個像素對應的距離，為水平分辨率，即柵格分辨率。3. 根據權利要求1所述的用于數字地形分析建模知識案例化自動處理方法，其特征在 于:所述步驟(4)中將各坡度等級在高程等級上進行累加運算，生成高程-坡度累積頻率二 維表的過程如下：在每個高程等級上，根據坡度從小到大進行累加，若累加前的高程-坡度 頻率分布表上高程為i等級坡度為j等級的頻率為 aij，累加之后高程為i等級坡度為j等級的 累加頻率為bij: bij- Σ η=1. . jclin 式中，1 = 1，2,3,4,5,6,7,8,9，10;]_ = 1，2,3,4,5,6,7;11為整數。4. 根據權利要求1所述的用于數字地形分析建模知識案例化自動處理方法，其特征在 于:所述步驟(4)中的地理信息系統軟件包括ArcGIS、GRASS GIS、SuperMap。5. 根據權利要求1所述的用于數字地形分析建模知識案例化自動處理方法，其特征在 于:所述步驟(4)中生成一張10 X 7的二維表。
【文檔編號】G06F17/30GK105893590SQ201610213395
【公開日】2016年8月24日
【申請日】2016年4月7日
【發明人】吳雪薇, 秦承志, 盧巖君, 朱阿興
【申請人】中國科學院地理科學與資源研究所