薛樹強，黨亞民，秘金鐘，劉紀平，董 春，吳 波，王世進

1.中國測繪科學(xué)研究院，北京100830；2.長安大學(xué)地質(zhì)工程與測繪學(xué)院，陜西 西安710054

1 引 言

地理國情監(jiān)測是新時期空間信息科學(xué)的重要使命［1］，是地理信息更深層次的應(yīng)用。地理國情包括3種信息：感知信息、統(tǒng)計信息和分析信息。許多國家和組織都開展了地理國情監(jiān)測方面的項目或工程，以便更好地服務(wù)于本國或地區(qū)的資源、環(huán)境、能源、社會等領(lǐng)域［2－4］。

數(shù)字高程模型（digital elevation model，DEM）是表示地面高程的一種實體地面模型［5］，在測繪、水文、氣象、地貌、地質(zhì)等科學(xué)領(lǐng)域應(yīng)用廣泛。地表高斯投影面積、橢球面積和地表面積是3類常用統(tǒng)計量［6－7］。文獻［8］討論了高斯投影中央子午線變化和投影面高程變化引起的面積變形，以及橢球面面積與高斯面上面積的差異。地球旋轉(zhuǎn)橢球和地球三軸橢球是兩種常用的地球形狀簡化模型。已有文獻討論多邊形區(qū)域的旋轉(zhuǎn)橢球或三軸橢球面積計算的嚴密公式［3，9－10］。在許多應(yīng)用中，地表面積比投影面積或橢球面積更具有參考價值［11－12］。構(gòu)造地表三角網(wǎng)并建立地表面積柵格數(shù)據(jù)是地表面積統(tǒng)計的常用方法［11］。

對于相同的投影面積，地表面積會因地形因子不同而異。利用DEM提取地形因子可分為數(shù)值分析方法和局部曲面擬合方法［14－16，30］。坡度和坡向的常用數(shù)值分析方法主要包括簡單差分、二階差分、三階差分、Frame算法等［28］。應(yīng)用表明，三階不帶權(quán)差分算法具有較高精度［28］；局部曲面擬合方法則通過最小二乘曲面擬合再進行求導(dǎo)獲取地形的坡度和坡向等信息［30］?；跀?shù)值積分的復(fù)化公式法則需要高分辨率的微觀地形因子數(shù)據(jù)［13］。文獻［17］基于二維小波的多分辨分析理論，提出并證明了不同分辨率下地面坡度的面積與復(fù)雜地貌形態(tài)復(fù)雜度及尺度之間的變化關(guān)系。對于具有特定分辨率的DEM數(shù)據(jù)，提高地形因子計算精度需根據(jù)地形起伏情況設(shè)計合適的窗口，最大限度消除隨機觀測誤差的影響。利用準確的地形因子及其變化率信息構(gòu)造地表三角網(wǎng)可提高地表面積的統(tǒng)計精度。

本文研究地表面積的理論面積及其近似方法，討論影響地表面積統(tǒng)計精度的主要因素，提出了一種顧及任意地形因子的地表面積統(tǒng)計方法。該方法首先利用泰勒級數(shù)逼近原理對地形因子進行最小二乘估計，然后利用這些地形因子對DEM進行加密，最后利用加密后的DEM構(gòu)建地表三角網(wǎng)統(tǒng)計地表面積。試驗表明，在局部地形因子顯著的山區(qū)或丘陵地區(qū)，不同分辨率下地表面積統(tǒng)計結(jié)果存在較大的差異，而利用本文提出的顧及地形局部地形因子的地表面積統(tǒng)計方法可明顯提高低分辨率DEM地表面積統(tǒng)計精度。

2 地表面積統(tǒng)計原理與精度分析

2.1 地表的理論面積及其近似

數(shù)字高程模型是一種采用二維高斯投影坐標附加正（常）高實現(xiàn)的一種三維空間地理信息，是反映離散地表曲面模型。如圖1所示，對于地表封閉區(qū)域Ω，設(shè)投影區(qū)域ΩG的邊界線為A1、A2、…、Ap。由于（似）大地水準面和參考橢球面均為曲面參考基準，且高斯投影變換存在面積變形，當利用DEM數(shù)據(jù)計算地表面積時，可將DEM據(jù)中恢復(fù)為地形表面的三維空間模型［23－24］。

圖1 投影DEM格網(wǎng)劃分Fig.1 Projection DEM grid division

對于在地固地心坐標系中給出的地表曲面（X，Y，Z），其尺度來自于地心坐標系基準的定義。將地表曲面（X，Y，Z）表示為函數(shù)式Z＝f（X，Y），并設(shè)其為分段光滑曲面，地表區(qū)域Ω在XOY平面上的投影為ΩO，則地表區(qū)域Ω的面積由以下曲面積分給出［13，22］

式中，fX、fY為曲面函數(shù)Z＝f（X，Y）的一階偏導(dǎo)數(shù)。如圖1所示，對于給定的高斯投影邊界區(qū)域ΩG，則分別率為g的DEM格網(wǎng)對ΩG的劃分為

式中，UG，i和ΦG，i分別為完整單元格和破碎單元格。

基于高斯投影劃分式（2），地表區(qū)域Ω的面積可由表示為

式中

表示高斯投影區(qū)域UG，i和ΦG，i映射為地表區(qū)域（由坐標變換公式和基準定義給出）。若子區(qū)域內(nèi)地表為平面或起伏不大，則式（4）中子區(qū)域的面積可近似為

式中，T（UO，i）和T（ΦO，j）分別為子區(qū)域Ui和Φj在XOY平面上的投影面積，fX、fY在相應(yīng)子區(qū)域內(nèi)的一階偏導(dǎo)數(shù)?？梢姡檬剑?）計算地形表面積需要地形表面的嚴格曲面方程，而利用式（4）近似地表面積則需要高分辨率DEM數(shù)據(jù)。

2.2 局部地形因子對地表面積計算的影響

以高斯投影子區(qū)域UG，i邊界上的任意DEM格網(wǎng)點（xG，0，yG，0）按高斯投影逆變換映射至橢球面的o點，并將該映射點o為原點在橢球面上建立地平空間直角坐標系o－xyH，則在o－xyH系下，地表曲面方程為

式中

式中，x、y為高斯投影坐標。顧及大地高和正（常）高的轉(zhuǎn)換關(guān)系，可得

式中，h（x，y）為正常高或正高；ζ（x，y）為高程異?；虼蟮厮疁拭娌罹啵籯i為高斯投影長度變形比。在空間直角坐標系o－xyH下，類似于式（3），地表區(qū)域Ui的地表面積為

式中，積分區(qū)域UE，i為地表區(qū)域Ui在地平面上的投影區(qū)域。

式 中，α∈ ［0，1］

對于邊長為kig的高斯投影正方形區(qū)域UG，i，結(jié)合式（3）和式（13），其表面積為

式（17）第1項為顧及地形函數(shù)一階偏導(dǎo)數(shù)的地表區(qū)域Ui的近似面積，第2項為地形函數(shù)H＝f（x，y）二階偏導(dǎo)數(shù)對地表面積計算的影響，即利用線性地形因子（如坡度或地表函數(shù)一階偏導(dǎo)數(shù)）計算地表面積的截斷誤差。由第2項可知，當且僅當DEM分辨率足夠高，即g足夠小，或者足夠小，即地表函數(shù)二階偏導(dǎo)數(shù)足夠小時，利用式（4）計算地表面積才是嚴格的。為此，本文將地表曲面在地平坐標系下的各偏導(dǎo)數(shù)統(tǒng)稱為局部地形因子，以便討論和揭示地形起伏對地表面積統(tǒng)計的影響。

3 局部地形因子的最小二乘估計

對于給定的DEM數(shù)據(jù)，提高地表面積統(tǒng)計精度需提取DEM數(shù)據(jù)中所蘊含的局部地形因子。由式（11）可得以下微分關(guān)系

式中，Hx、Hy，Hxx、Hyy、Hxy分別為函數(shù)式（9）的一階、二階。在小區(qū)域內(nèi)，ζx和ζy的幾何意義為大地水準面相對于橢球面的夾角正切值。

由式（18）可知，曲面（x，y，H）的偏導(dǎo)數(shù)信息可由DEM 曲面（xG，yG，h）的偏導(dǎo)數(shù)間接計算。局部地形因子空間函數(shù)近似服從Tobler定理，即距離較近的點比距離較遠的點其特征值具有更大的相似性［21，25］。使用窗口大小為N×N實測數(shù)據(jù)構(gòu)建觀測方程，利用最小二乘原理計算窗口中心DEM格網(wǎng)點處的各階偏導(dǎo)數(shù)。下面基于DEM 數(shù)據(jù)討論地形函數(shù)h＝h（x）（x：＝［xGyG］∈R2為高斯投影坐標）在已知點處的各階偏導(dǎo)數(shù)計算方法。在已知DEM點xi處，利用窗口內(nèi)的其他已知點可構(gòu)造以下觀測方程［26］

式中，j＝1，2，…，n，j≠i，Δxj：＝ ［Δxj，Δyj］＝xj－xi，xj為已知DEM點xi處的鄰近已知點；rk（Δxj）為泰勒級數(shù)展開殘余項，n＝N×（N－1）為窗口內(nèi)可構(gòu)造的虛擬觀測方程數(shù)目，k為泰勒級數(shù)展開的階數(shù)。將方程（19）簡記為

式中

類似的，可繼續(xù)寫出h（x）在xi處的高階偏導(dǎo)數(shù)。將方程（20）中的殘余項r視為隨機誤差，則該問題的最小二乘解為

采用局部地形因子信號的最小二乘估計，可平滑DEM量測誤差，提高局部地形因子信號估值的精度和可靠性。當N＝3，方程（20）存在8個方程，即n＝8。若令k＝1，即僅估計式（28）給出的兩個地形參數(shù)，此時

則由式（31）可得

式（33）即為坡度坡向計算數(shù)學(xué)模型之三階不帶權(quán)差分公式［28］。事實上，若對虛擬觀測量Δyi賦予適當權(quán)重，則可導(dǎo)出其他配權(quán)坡度坡向計算數(shù)學(xué)模型，如使用反距離配權(quán)模型。值得注意的是，文獻［28］在坡度坡向算法精度分析中指出，在實際生產(chǎn)中，因DEM誤差占主導(dǎo)因素，三階不帶權(quán)差分算法精度最高，推薦使用三階不帶權(quán)差分坡度計算模型。文獻［26］指出：若平差模型式（20）僅估計空間函數(shù)的一階偏導(dǎo)數(shù)，將空間函數(shù)非線性信號的高階偏導(dǎo)數(shù)信息歸為平差隨機模型是有條件的。

4 DEM加密與地表面積統(tǒng)計

當DEM分辨率不夠高或地表非線性起伏顯著時，可利用地表函數(shù)的高階偏導(dǎo)數(shù)信息提高地表面積統(tǒng)計精度。如前所述，利用式（4）計算地表面積則需要地表離散數(shù)據(jù)，且DEM數(shù)據(jù)的分辨率越高，式（4）近似地表理論面積的精度越高。下面給出利用式（31）所得地表函數(shù)偏導(dǎo)數(shù)估值對DEM進行加密方法，在通過提高DEM分辨率和地面三角網(wǎng)的分辨率，進而提高地表面積統(tǒng)計精度。

地表區(qū)域邊界一般由高斯投影坐標給出，且區(qū)域邊界的頂點一般為非整數(shù)格網(wǎng)點，需要推估這些多邊形邊界區(qū)域的高程。當利用低分辨率DEM數(shù)據(jù)統(tǒng)計地表面積時，可利用局部地形因子對DEM進行加密，以提高地表面積統(tǒng)計的精度。設(shè)距離待估點x處最近的整數(shù)格網(wǎng)點為xi，顧及整數(shù)格網(wǎng)點xi處的各階偏導(dǎo)數(shù)，則待估點x的高程推估模型可表示為

式中，各階偏導(dǎo)數(shù)由式（31）給出。

當整數(shù)格網(wǎng)點xi處的各階偏導(dǎo)數(shù)為0時，則式（34）給出的高程推估退化為鄰近點插值法。式（34）是一種顧及局部地形因子的鄰近點高程插值法，因其使用一個已知點進行高程推估，可將該模型稱為單點推估模型。基于空間函數(shù)高階偏導(dǎo)數(shù)的多點推估方法可參考文獻［26］。

如圖1所示，根據(jù)加密DEM數(shù)據(jù)的分辨率，依次計算統(tǒng)計封閉多邊形區(qū)域邊界與整數(shù)格網(wǎng)的交點，并利用式（34）計算這些交點的高程值。對封閉多邊形邊界點和區(qū)域內(nèi)單元格網(wǎng)按DEM等間距格網(wǎng)線逐行逐列構(gòu)建地表規(guī)則三角網(wǎng)（保證每個三角形的頂點坐標互差均小于等于1倍的整數(shù)格網(wǎng)間距）。地表三角網(wǎng)構(gòu)造方法并不唯一，地表不規(guī)則三角網(wǎng)構(gòu)建方法可參考文獻［27］。將地表三角網(wǎng)點坐標轉(zhuǎn)化為大地空間直角坐標后，累加地表三角網(wǎng)中三角形的面積，可得地表面積的三角網(wǎng)近似為

式中，Ti為第i個地表三角形的面積。

5 算 例

5.1 算例1：DEM分辨率對地表面積統(tǒng)計精度的

為驗證DEM分辨率對地表面積統(tǒng)計精度的影響，在5m分辨率實測DEM數(shù)據(jù)的基礎(chǔ)上，提取了25m分辨率DEM數(shù)據(jù)，使用這兩種分辨率DEM構(gòu)造地表三角網(wǎng)計算地表面積。四川茂縣東興鄉(xiāng)（山區(qū)）和山東東營六戶鎮(zhèn)（平原地區(qū)）的地表面積計算結(jié)果表明，對于山區(qū)，地表面積和橢球面積之間的差異明顯，不同DEM分辨率計算地表面積結(jié)果差異也較大，例如，利用5m分辨率DEM計算四川茂縣東興鄉(xiāng)地表面積比利用25m分辨率DEM計算所得表面大9.93%。在平原地區(qū)，地表面積和橢球面積相差不大，不同DEM分辨率統(tǒng)計所得地表面積差異也較小，例如，不同DEM統(tǒng)計山東東營六戶鎮(zhèn)地表面積差異僅為0.236%。

5.2 算例2：實測DEM與內(nèi)插DEM統(tǒng)計地表面積精度對比

DEM分辨率和地形起伏程度共同決定了地表三角網(wǎng)的地表面積統(tǒng)計精度。為驗證本文的方法，采用式（31）提取區(qū)域地形因子，并采用式（34）將實測25m分辨率DEM內(nèi)插為5m分辨率DEM后，構(gòu)建地表三角網(wǎng)并統(tǒng)計地表面積，即按式（35）計算地表面積。本算例采用以下7種方案（表1）。

以實測5m分辨率DEM計算所得地表面積為參考，評價內(nèi)插5m分辨率DEM對地表面積統(tǒng)計的精度。結(jié)果表明，在地形起伏較大的茂縣，上述7種DEM加密方案中，當考慮地表函數(shù)的二階或三階偏導(dǎo)數(shù)時，地表面積統(tǒng)計精度最高，其25m內(nèi)插5m計算結(jié)果和實測5mDEM統(tǒng)計結(jié)果相差僅為0.5%，相比于9.93%，地表面積計算精度提高顯著。然而，若僅提取線性地形因子，采用三階差分算法具有最高的地表面積統(tǒng)計精度，而簡單差分算法精度最低（其差異高達39.48%），二階差分算法次之，這與文獻［28］討論坡度信息提取算法時的結(jié)論一致。而在局部地形起伏較小的平原，各種方法對地表面積統(tǒng)計結(jié)果均有改進，其中基于簡單差分算法提取地形因子對地表面積幾乎沒有改進，而基于三階差分算法提取地形因子對地表面積的改進最為顯著?？梢?，地形因子的提取精度及其分辨率是最終決定地表面積計算精度的兩個重要因素。對于特定分辨率的DEM，將本文給出的公式（31）作為局部地形因子的統(tǒng)一計算模型，對于山區(qū)或丘陵地帶，可采用相對大的窗口提取非線性地形因（如二階偏導(dǎo)數(shù)），而對于平原地區(qū)，可直接小窗口使用三階差分算法僅提取線性地形因子。

表1 顧及地形非線性起伏統(tǒng)計地表面積Tab.1 Surface area Computation with regard to the nonlinear terrain factors

6 結(jié) 論

地表理論面積計算由地表曲面積分給出，其計算需要附加地表曲面函數(shù)假設(shè)。在實踐中，可利用基于局部地形因子的最小二乘估計，對地表曲面進行剖分和面積近似。影響地形表面區(qū)域面積計算精度的因素由局部地形因子、DEM精度和分辨率。當僅利用坡度信息計算地表面積時，其截斷誤差與DEM分辨率的三次方成正比，而與地表空間函數(shù)的二次階偏導(dǎo)數(shù)成正比。在實際應(yīng)用中，可利用DEM格網(wǎng)點的對稱性導(dǎo)出地形因子的實用計算公式，避免復(fù)雜的矩陣求逆運算。利用DEM提取所得地形因子有望提高地表面積的統(tǒng)計精度。當僅考慮地形線性地形因子時，本文建議采用三階差分算法提取地形因子。在地形起伏顯著的地區(qū)，窗口不應(yīng)小于5×5，地形參數(shù)不應(yīng)少于5個，即應(yīng)考慮坡度變率對地表面積計算的影響。鑒于我國5m分辨率DEM覆蓋率不足50%，建議采用非線性插值對低分辨率DEM數(shù)據(jù)進行加密或根據(jù)需要重新采集生產(chǎn)高分辨率DEM數(shù)據(jù)，以提高地表面積統(tǒng)計精度。

［1］ LI Deren，SUI Haigang，SHAN Jie.Discussion on Key Technologies of Geographic National Conditions Monitoring［J］.Geomatics and Information Science of Wuhan University，2012（5）：505－512.（李德仁，眭海剛，單杰.論地理國情監(jiān)測的技術(shù)支撐［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2012（5）：505－512.）

［2］ CHEN Junyong.Study Notes on National Condition Monitoring［J］.Acta Geodaetica et Cartographica Sinica，2012，41（5）：633－635.（陳俊勇.地理國情監(jiān)測的學(xué)習札記［J］.測繪學(xué)報，2012，41（5）：633－635.）

［3］ CAO Zhenyu，TAN Mingjian，SHEN Yanmin.Study on Polygon Area Calculation Methods in Geodetic Coordinate System［J］.Bulletin of Surveying and Mapping，2012（S1）：570－572.（曹振宇，譚明建，沈艷敏.大地坐標系中多邊形面積計算方法研究［J］.測繪通報，2012（S1）：570－572.）

［4］ QIAO Chaofei.Review of Overseas National Condition Monitoring［J］.Bulletin of Surveying and Mapping，2011（11）：81－83.（喬朝飛.國外地理國情監(jiān)測概況與啟示［J］.測繪通報，2011（11）：81－83.）

［5］ American Society for Photogrammetry and Remote Sensing.Digital Elevation Model Technologies and Applications：The DEM Users Manual［M］.Bethesda ASPRS Publications，2007.

［6］ SJ?BERG L E.Determination of Areas on the Plane，Sphere and Ellipsoid［J］.Survey Review，2006.38（301）：583－593.

［7］ GILLISSEN I.Area Computation of a Polygon on an Ellipsoid［J］.Survey Review，1993，32（248）：92－98.

［8］ WANG Jiexian，YU Zhenwu.Area Calculation Error Caused by Gauss Projection［J］.Bulletin of Surveying and Mapping，2003（4）：5－6.（王解先，俞振武.高斯投影引起的面積計算誤差［J］.測繪通報，2003（4）：5－6.）

［9］ PEDZICH P，KU Z＇MA M.Application of Methods for Area Calculation of Geodesic Polygons on Polish Administrative Units［J］.Geodesy and Cartography，2012，61（2）：105－115.

［10］ JI Bing，BIAN Shaofeng.Calculation of Surface Area of Triaxial Ellipsoid［J］.Journal of Naval University of Engineering，2008（4）：21－24.（紀兵，邊少鋒.三軸橢球表面積的計算［J］.海軍工程大學(xué)學(xué)報，2008（4）：21－24.）

［11］ ZHANG Y.Surface Area Processing in GIS for Different Mountain Regions［J］.Forestry Studies in China，2011，13（4）：311－314.

［12］ JENNESS J S.Calculating Landscape Surface Area from Digital Elevation Models［J］.Wildlife Society Bulletin，2004，32（3）：829－839.

［13］ JIANG Fan，LU Xiaohua，WANG Zhonglan.The Analysis of DEM Surface Area Based on Composed Calculus Algorithms［J］.Journal of Institute of Surveying and Mapping，2005（4）：263－265.（江帆，呂曉華，王仲蘭.基于復(fù)化公式的DEM表面積算法分析［J］.測繪學(xué)院學(xué)報，2005（4）：263－265.）

［14］ JENSON S，DOMINGUE J.Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis［J］.Photogrammetric Engineering and Remote Sensing，1988，54（11）：1593－1600.

［15］ HOECHSTETTER S.Effects of Topography and Surface Roughness in Analyses of Landscape Structure—a Proposal to Modify the Existing Set of Landscape Metrics［J］.Landscape Online，2008，3：1－14.

［16］ BEASOM S L，WIGGERS E P，GIARDINO J R.A Technique for Assessing Land Surface Ruggedness［J］.The Journal of Wildlife Management，1983，47（4）：1163－1166.

［17］ QIN Yanhong，LU Lin，WU Jitao.A Study on the Multi－resolution Representation of the Areas Based on Two－dimension Wavelet Transformation ［J］. Mathematical Theory and Applications，2005（2）：43－48.（秦艷紅，盧林，吳紀桃.地形坡度面積的二維多分辨分析研究［J］.數(shù)學(xué)理論與應(yīng)用，2005（2）：43－48.）

［18］ AGUILAR F J.Effects of Terrain Morphology，Sampling Density，and Interpolation Methods on Grid DEM Accuracy［J］.Photogrammetric Engineering and Remote Sensing，2005，71（7）：805.

［19］ WISE S M.The Effect of GIS Interpolation Errors on the Use of Digital Elevation Models in Geomorphology［J］.Landform Monitoring，Modelling and Analysis，1998，99：139－164.

［20］ CHILDS C.Interpolating Surfaces in ArcGIS Spatial Analyst［R］.［S.l.］：ESRI Education Services，2004.

［21］ LI Xin，CHENG Guodong，LU Ling.Comparison of Spatial Interpolation Methods［J］.Advance in Earth Sciences，2000，15（3）：260－265.（李新，程國棟，盧玲.空間內(nèi)插方法比較［J］.地球科學(xué)進展，2000，15（3）：260－265.）

［22］ THOMAS G B，F(xiàn)INNEY R L.Thomas’Calculus［M］.［S.l.］：Addison Wesley Publishing Company，2002.

［23］ DANG Yamin，CHENG Yingyan，XUE Shuqiang.Geodetic Coordinate System and Its Applications［M］.Beijing：Surveying and Mapping Press，2010.（黨亞民，成英燕，薛樹強.大地坐標系統(tǒng)及其應(yīng)用［M］.北京：測繪出版社，2010.）

［24］ ZHU Huatong.Establishment of Geodetic Coordinate System［M］Beijing：Surveying and Mapping Press，1986.（朱華統(tǒng).大地坐標系的建立［M］.北京：測繪出版社，1986.）

［25］ ZHU Huiyi，LIU Shulin，JIA Shaofeng.Problems of the Spatial Interpolation of Physical Geographical Elements［J］.Geographical Research，2004，23（4）：425－432.（朱會義，劉述林，賈紹鳳.自然地理要素空間插值的幾個問題［J］.地理研究，2004，23（4）：425－432.）

［26］ XUE Shuqiang，YANG Yuanxi.Generalized Inverse Distance Weighting Method for Spatial Interpolation［J］.Geomatics and Information Science of Wuhan University，2013，38（12）：1435－1439.（薛樹強，楊元喜.廣義反距離加權(quán)空間推估法［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2013，38（12）：1435－1439.）

［27］ ZHU Qing，CHEN Chujiang.Quick Generation of TIN and Its Dynamic Updating［J］.Journal of Wuhan Technical University of Surveying and Mapping，1998（3）：18－21.（朱慶，陳楚江.不規(guī)則三角網(wǎng)的快速建立及其動態(tài)更新［J］.武漢測繪科技大學(xué)學(xué)報，1998（3）：18－21.）

［28］ LIU Xuejun，GONG Jianya，ZHOU Qiming，et al.A Study of Accuracy and Algorithms for Calculating Slope and Aspect Based on Grid Digital Elevation Model（DEM）［J］.Acta Geodaetica et Cartographica Sinica，2004，33（3）：258－263.（劉學(xué)軍，龔健雅，周啟鳴，等.基于DEM坡度坡向算法精度的分析研究［J］.測繪學(xué)報，2004，33（3）：258－263）.

［29］ LIU Xuejin，REN Zhifeng，WANG Yanfang，et al.Slope Model at Arbitrary Direction Derived from Grid－based DEM［J］.Areal Research and Development，2009（4）：139－141.（劉學(xué)軍，任志峰，王彥芳，等.基于DEM的任意方向坡度計算方法［J］.地域研究與開發(fā)，2009（4）：139－141.）

［30］ HE Dan.Window Analysis for Extraction of Terrain Factors from Different Resolution DEM［D］.Xi’an：Northwest University，2012.（賀丹.不同分辨率DEM提取地形因子的適宜分析窗口研究［D］.西安：西北大學(xué)，2012.）