DEM提取坡向變率中的誤差特征與消除方法

謝軼群，湯國安，2＊，江 嶺，2

（1.南京師范大學(xué)地理科學(xué)學(xué)院，江蘇 南京 210023；2.南京師范大學(xué)虛擬地理環(huán)境教育部重點實驗室，江蘇 南京 210023）

0 引言

坡向變率（Slope of Aspect，SOA）［1］是反映地形在水平方向變化特征的關(guān)鍵指標(biāo)，也是地形定量指標(biāo)體系的重要組成部分，在基于DEM的數(shù)字地形分析中具有特殊意義［2－10］。坡向變率的提取一般采用對DEM坡向矩陣提取坡度的方法獲得。但由于坡向值的取值區(qū)間［0，360）無法以漸變的方式表達由360到0的過渡，尤其會夸大在正北方向上的變化率，對該方法的提取結(jié)果若不加修正則會出現(xiàn)大量誤差（后文中簡稱SOA誤差），不利于后續(xù)分析及應(yīng)用。因此，研究SOA誤差的特征，并在此基礎(chǔ)上提出完善的消除方法有重要意義。

湯國安等提出了利用正反DEM結(jié)合消除誤差的方法［11］，并得到廣泛應(yīng)用［5－10］。但仍存在兩個問題：1）誤差特征的概括尚不全面，事實上，誤差不只存在于正北坡位置，只要提取SOA時坡向矩陣中兩坡向值之差的絕對值大于180，則會出現(xiàn)誤差；2）由于在提取坡度時需要采用相應(yīng)算法，而原有的誤差消除方法并沒有考慮到坡度提取算法對結(jié)果造成的影響，因此，無法取得理想的效果。本文在對SOA誤差成因研究的基礎(chǔ)上，通過分析，揭示了正反DEM方法中存在的缺陷，并提出了基于“六分法”的誤差消除方法。

1 誤差成因與解決途徑

1.1 SOA誤差的成因

對于地面任何一個微分單元，坡度為過該點的切面的傾角，而坡向為該切面法線在水平面的投影方向。坡向線為矢量線，所以，鄰坡坡向線的夾角應(yīng)為兩向量的夾角。由向量夾角的定義［12］可知，兩向量夾角的取值范圍為［0°，180°］，因此，相鄰兩面坡的坡向線的夾角不應(yīng)大于180°；而由于坡向值的取值范圍是［0，360），致使相鄰坡向值之差的取值范圍為［0，360），超出了正確的范圍。至于北坡誤差，其實為一種達到最大值的特例。

1.2 正反DEM消除誤差方法的缺陷

提取坡度時，三階反距離平方權(quán)差分算法的精度較高［13］，且在諸多GIS平臺軟件中得到了廣泛應(yīng)用，公式如下：

其中，z代表鄰域內(nèi)的柵格值，i、j為鄰域內(nèi)的行、列號，分別自上至下、自左至右逐漸增大。

當(dāng)用該算法提取SOA時，柵格值為坡向值，因此，結(jié)合式（2）、式（3）及本文對SOA誤差特征的概括，可以看出當(dāng)提取SOA時，若fx或fy分子式中對應(yīng)的三對坡向差值的絕對值存在大于180的情況，則會造成誤差。

正反DEM的誤差消除方法［11］，可以理解為在每個柵格處分別根據(jù)正反DEM提取SOA，并取最小值作為最終結(jié)果。圖1所示為正反DEM分別對應(yīng)的坡向，其值相差180，而由于坡向值不能存在負值和大于等于360的情況，因此，正反DEM坡向矩陣中同一位置坡向值的關(guān)系如下：

其中，aspectrev、aspectobv分別為正、反DEM對應(yīng)柵格的坡向值。

圖1 正反地形坡向示意Fig.1 Aspects of obverse and reverse terrains

正反DEM的誤差修正方法主要存在以下缺陷：（1）求算fx或fy時，忽略了正反DEM對應(yīng)的式（2）、式（3）分子中3組差值正負關(guān)系的不一致性。圖2所示為基于正反DEM在相同位置分別提取的局部坡向矩陣，由正、反DEM提取的SOA值分別為SOA1和SOA2。圖2a（左）中坡向矩陣對應(yīng)的X、Y方向上6對差值中，不存在絕對值大于180的錯誤情況，而圖2a（右）的反DEM坡向矩陣中則出現(xiàn)了這種錯誤。但是，由于X方向上3對差值間的正負關(guān)系發(fā)生了變化，使差值結(jié)合后反而小于正確情況，導(dǎo)致SOA1大于SOA2。根據(jù)其SOA求算公式［6］，SOA2被誤作正確結(jié)果。

（2）求算fx和fy時，忽略了正反DEM對應(yīng)的式（2）、式（3）分子中6對差值正誤關(guān)系的不一致性。由于采用三階差分算法提取SOA時，會同時使用鄰域窗口內(nèi)X、Y方向上的6對差值，因此，必須保證6對差值都在正確范圍內(nèi)（絕對值小于180），并且相互之間的正負關(guān)系正確。雖然借助反DEM，修正了提取SOA時正DEM坡向矩陣中存在的錯誤差值，但同時可能生成新的錯誤差值（如南坡變?yōu)楸逼拢?，因此，若提取SOA時，新生成的錯誤差值與修正的差值處于同一個鄰域窗口中，則無法求算正確的SOA。圖2b為基于正反DEM在相同位置分別提取的鄰域坡向矩陣，其中，左圖在前兩行對應(yīng)的差值存在錯誤；經(jīng)過反地形處理后，原有誤差得以消除，但在第3行與第2列生成了新誤差。因此，由兩矩陣提取得到的SOA均為錯，無法對誤差進行修正。

圖2 正反地形DEM鄰域坡向矩陣Fig.2 Neighborhood aspects in obverse and reverse DEMs

1.3 SOA誤差消除方法——六分法

1.3.1 三階反距離平方權(quán)差分算法的拆分 通過正反地形DEM結(jié)合消除SOA誤差的方法，由于沒有考慮到式（2）、式（3）中差值間的關(guān)系，無法有效地消除誤差。因此，須對式（2）、式（3）分子中的三組差值分別進行修正，既要保證每一組差值在數(shù)值上正確，又要保證相應(yīng)差值間的正負關(guān)系正確。綜合上述分析，本文提出了一種新的誤差消除方法——六分法。為避免誤差消除過程中各差值的相互影響，確保誤差修正的獨立性，在計算fx和fy前，首先需要將兩式的分子拆分為6個子式：式（5）—式（10），之后逐個子式進行計算與修正，使進行組合時差值的值與正負關(guān)系均正確。

1.3.2 差值修正標(biāo)準(zhǔn)的建立 計算獲得式（5）—式（10）的6對差值后，需要對其建立相應(yīng)標(biāo)準(zhǔn)進行修正，修正標(biāo)準(zhǔn)包含兩部分：

（1）差值間正負關(guān)系的修正標(biāo)準(zhǔn)。在鄰域坡向矩陣中，式（5）—式（7）對應(yīng)的坡向值相減順序均為“右減左”（X 方向），式（8）—式（10）均為“下減上”（Y方向），具有內(nèi)部的一致性，且在最后對fx和fy進行結(jié)合時，會對其分別取平方。因此，可以對X，Y方向上的差值采用相同的修正標(biāo)準(zhǔn)，保證它們的正負關(guān)系正確，且對兩方向上的差值進行結(jié)合時，不需考慮fx和fy之間的正負關(guān)系。

如圖3所示，坡向值的取值范圍為［0，360），并以正北方向為0，沿順時針方向增加。由于坡向為矢量，具有方向性，因此，計算SOA時，需采用向量方法計算坡向值的差值，而差值間的正負關(guān)系也必須利用向量的特征進行推導(dǎo)。將坡向按圖3中的規(guī)則以向量的形式表示后，式（5）—式（10）中的差值即為被減向量順時針轉(zhuǎn)到減向量方向時，需經(jīng)過的絕對值最小的角度α｜min｜（在正確情況下），該角度可以為負值；反之，在差值錯誤的情況下，其結(jié)果所反映的則是絕對值最大的角度α｜max｜，導(dǎo)致錯誤的原因則是坡向值差值的取值區(qū)間［0，360）與向量夾角的取值區(qū)間［0，180］不同，且坡向值不能連續(xù)地表示坡向的周期性，即無法連續(xù)地由360過渡到0。由于α｜min｜與α｜max｜均表示被減向量沿順時針方向轉(zhuǎn)過的角度，因此，α｜min｜·α｜max｜＜0，且｜α｜min｜｜＋｜α｜max｜｜＝360。

圖3 坡向相減示意Fig.3 Subtraction of aspects

本文以α｜min｜作為提取式（5）—式（10）中差值的標(biāo)準(zhǔn)（也可以逆時針方向建立相應(yīng)的標(biāo)準(zhǔn)）。圖3表示了在北坡的兩種情況下，方向向量A1和A2間差值的正負性。其中，圖3a中｜A1｜∈（0，90），｜A2｜∈（270，360），因此，｜A2｜－｜A1｜＞180。但因A1沿順時針方向轉(zhuǎn)到A2所經(jīng)過的α｜min｜∈（－180，0），對該差值進行修正時，需要變換其值的正負性。同理，圖3b中｜A1｜與｜A2｜的差值也需要進行正負性的變換，使得式（5）—式（7）與式（8）—式（10）中差值間的正負關(guān)系正確。

（2）差值絕對值的修正標(biāo)準(zhǔn)。根據(jù)誤差產(chǎn)生的原理，差值的絕對值必須在［0，180］區(qū)間內(nèi)。由于差值只能取α｜min｜與α｜max｜之一，其中｜α｜min｜｜∈［0，180］，｜α｜max｜｜∈（180，360），且｜α｜min｜｜＋｜α｜max｜｜＝360，因此，對于所有絕對值大于180的差值，需要以360－｜α｜max｜｜的形式將其轉(zhuǎn)換為｜α｜min｜｜。

綜上，由于α｜min｜·α｜max｜＜0 且｜α｜min｜｜＝360－｜α｜max｜｜，對差值的修正可根據(jù)式（11）完成。

其中，diff和diff＿r分別表示修正前、后的差值。

1.3.3 差值的結(jié)合——提取SOA值 根據(jù)式（11）完成對6對差值的修正后，可以得到絕對值與正負關(guān)系均正確的結(jié)果：f′xi，f′yi（i＝1，2，3）；再由式（2）、式（3）得到式（12）、式（13），將f′xi，f′yi分別帶入式（12）、式（13），即可得到結(jié)合后的正確的fx，fy。最后將fx，fy帶入式（1），即可得到最終的SOA提取結(jié)果。圖4所示的流程圖表示了六分法對SOA進行提取與修正的整個過程。

圖4 六分法流程Fig.4 Architecture of the senary division method

2 實驗結(jié)果與分析

由于采用低空間分辨率DEM提取地形因子時受空間尺度不確定性的影響［14］，很多細部的坡向變化無法表達，得到的僅為概括性的坡向結(jié)果，誤差已經(jīng)很大［14－17］，因此，通過該坡向數(shù)據(jù)難以獲得高精度的SOA。據(jù)此，實驗采用5m分辨率的DEM數(shù)據(jù)為基本信息源，并基于ArcEngine開發(fā)提取算法。

通過對兩種方法的結(jié)果求差，即可得到正反DEM方法結(jié)果的誤差分布情況。進一步提取差值圖的絕對值后，將其中絕對值大于1、大于5、大于10及大于30的點位提取出來（圖5），以直觀地顯示不同級別差值的具體分布狀況（黑色為高于閾值的點）。表1給出了不同閾值d所對應(yīng)區(qū)域的柵格數(shù)占總柵格數(shù)的百分比P（d）及區(qū)域內(nèi)部差值絕對值的均值E。

通過表1可以看出，正反DEM法的結(jié)果中仍然存在大量誤差。在實驗樣區(qū)內(nèi)，最大誤差值約為70.7，誤差區(qū)域柵格均值約為5.87。此外，由六分法得出的更加精確的SOA矩陣的平均SOA值約為58.0，相比正反DEM法的57.8，消除了其修正結(jié)果中占SOA總值約0.34%的誤差。

表1 六分法與正反DEM法的對比Table 1 Comparison between senary division and traditional method

圖5 六分法與正反DEM法的差值Fig.5 Differences between the results of senary division method and traditional method

圖6中的誤差點取的是新舊方法提取結(jié)果中差值的絕對值大于1的點位，可以看出，正反DEM法無法消除的誤差主要分布在山脊或山谷處，即坡向變率相對較大的位置。其原因主要是在這些位置，鄰域坡向矩陣中的差值通常較大，差值間關(guān)系更為復(fù)雜，正反DEM方法難以消除這些誤差。

SOA值較高的區(qū)域往往是數(shù)字地形分析中更值得關(guān)注的部分，在圖6中也可看出，正反DEM法的結(jié)果中誤差較大的點位又趨向分布在這些高值區(qū)域。實驗樣區(qū)的統(tǒng)計結(jié)果表明，在正反DEM法與六分法提取結(jié)果中差值大于5的區(qū)域內(nèi)，80.5%以上的SOA值都大于83且具有13.1以上的平均差值，其柵格數(shù)占樣區(qū)內(nèi)SOA值大于83的柵格總數(shù)的6.8%以上。鑒于較高的SOA值在反映很多地形要素特征時的重要性，表2給出了按不同SOA區(qū)間分類的情況下，六分法對正反DEM法結(jié)果的改進效果，包括在某區(qū)間［a，b）內(nèi)時，錯誤柵格數(shù)n占區(qū)間總柵格數(shù)的百分比P（n），區(qū)間內(nèi)的柵格數(shù)N占柵格總數(shù)的百分比P（N），誤差的均值Ew，及區(qū)間內(nèi)誤差總值w占全局誤差總值的百分比P（w），這里的誤差均指誤差的絕對值。

表2 六分法在不同SOA區(qū)間內(nèi)的改進效果Table 2 Performances of senary division method in different intervals of SOA

從表2可以看出，正反DEM法提取結(jié)果中存在的誤差分布特性與圖6中一致，集中分布在SOA值較高的區(qū)域，特別是當(dāng)SOA大于80時，其區(qū)域內(nèi)包含的誤差值約占誤差總值的86.98%，在［60，90）這一區(qū)間內(nèi)的累積值則達到94.82%。此外，在SOA大于80的區(qū)域內(nèi)，有12.7%以上的SOA值存在誤差，且其誤差絕對值的均值達到5.6以上，因此，六分法在SOA的提取、分析與應(yīng)用中有重要作用。

3 結(jié)語

本文通過對SOA誤差的產(chǎn)生原理與消除方法的研究，首先，對SOA誤差的成因進行了完善：由于坡向值［0，360）無法以連續(xù)的形式表達坡向的周期性，即在360到0時出現(xiàn)數(shù)值上的躍變（跳躍間斷點），使得在計算坡向值的差值時，會得到大于180°的錯誤向量夾角［12］。其次，由于正反DEM法更多地關(guān)注誤差產(chǎn)生的機理，而沒有進一步考慮與提取算法結(jié)合時會產(chǎn)生的問題，使得該模型難以全面消除誤差，且可能生成新的錯誤，給基于SOA的研究結(jié)果帶來不確定的影響。本文根據(jù)對SOA誤差成因及正反DEM法缺陷的深入分析，提出了新的SOA誤差修正模型，通過對坡度提取算法的拆分以及相應(yīng)誤差修正標(biāo)準(zhǔn)的建立，實現(xiàn)了對SOA誤差更為全面、理想的消除，也使得SOA可以更為廣泛、正確的應(yīng)用于未來的數(shù)字地形分析研究中。今后的研究將針對多地形因子對地貌特征表達時的相互影響及關(guān)系展開，以進一步完善地形因子的相關(guān)研究。

［1］ TANG G A.A Research on the Accuracy of Digital Elevation Models［M］.Beijing，New York：Science Press，2000.

［2］ 陳楠，林宗堅，湯國安，等.數(shù)字高程模型的空間信息不確定性分析［J］.測繪通報，2005（17）：14－17.

［3］ 王盛萍，張志強，張建軍，等.黃土殘塬溝壑區(qū)流域次生植被物種分布的地形響應(yīng)［J］.生態(tài)學(xué)報，2010，30（22）：6102－6112.

［4］ 周訪濱，劉學(xué)軍.基于柵格DEM自動劃分微觀地貌形態(tài)的研究［J］.武漢理工大學(xué)學(xué)報（信息與管理工程版），2008，30（2）：173－175.

［5］ 張勇.黃土高原地面坡譜研究［D］.西北大學(xué)，2003.55－75.

［6］ 何麗麗.基于不同比例尺和分辨率DEM的數(shù)字地形分析［D］.西南大學(xué)，2007.18－20.

［7］ 韓海輝.基于SRTM－DEM的青藏高原地貌特征分析［D］.蘭州大學(xué)，2009.27－29.

［8］ 周毅，湯國安，王春，等.基于DEM增強黃土典型地貌表達效果的方法研究［J］.測繪通報，2009（11）：34－36.

［9］ 劉洋，楊晏立，何政偉，等.分水線、匯水線的多分辨率多閾值提取分析［J］.地理空間信息，2010，8（1）：41－44.

［10］ 陳婷，周汝良，朱大運，等.基于DEM的2種提取地形特征線算法對比研究［J］.林業(yè)調(diào)查規(guī)劃，2011，36（6）：1－4.

［11］ 湯國安，李發(fā)源，劉學(xué)軍.數(shù)字高程模型教程［M］.北京，科學(xué)出版社，2010.145－149，157－158.

［12］ R.柯朗，F(xiàn).約翰.微積分和數(shù)學(xué)分析引論（第一卷，第二分冊）［M］.北京：科學(xué)出版社，2001.434.

［13］ 劉學(xué)軍，龔健雅，周啟鳴，等.基于DEM坡度坡向算法精度的分析研究［J］.測繪學(xué)報，2004，33（3）：258－263.

［14］ 湯國安，趙牡丹，李天文，等.DEM提取黃土高原地面坡度的不確定性［J］.地理學(xué)報，2003，58（6）：824－830.

［15］ 陳楠，湯國安，劉詠梅，等.基于不同比例尺的DEM地形信息比較［J］.西北大學(xué)學(xué)報，2003，33（2）：237－240.

［16］ 陳楠，林宗堅，李成名，等.1∶10000及1∶50000比例尺DEM信息容量的比較——以陜北韭園溝流域為例［J］.測繪科學(xué)，2004，39（3）：39－41.

［17］ 湯國安，陳楠，劉詠梅，等.黃土丘陵溝壑區(qū)1∶1萬及1∶5萬比例尺DEM地形信息容量對比［J］.水土保持通報，2001，21（2）：34－36.