DEM坡向提取的向量場模板匹配方法

李潤超，袁林旺，李 碩，朱曉林，俞肇元

南京師范大學 虛擬地理環(huán)境教育部重點實驗室，江蘇 南京 210023

1 引 言

坡向在數學上定義為地表面一點切平面的法線矢量在水平面的投影與過該點的正北方向的夾角，坡向表征了地面點高程改變量的最大變化方向。目前計算坡向的數學計算模型主要利用DEM 數據采用二階或三 階差分 的方法［1－2］，均屬于直接基于標量場進行計算，受DEM精度誤差傳遞特性、地形復雜度等影響，多數坡向算法對DEM 高程數據誤差敏感，算法穩(wěn)定性不強［3－5］。DEM精度場數據可以描繪空間相關性規(guī)律，并可通過趨勢面擬合來揭示地形整體結構性的變化［6］。矢量場數據除包含表征地理現象影響強度大小的數值外，還同時包含了方向信息，因此，利用DEM數據構建矢量場來表征和計算坡向，更能真實地描述其本質特征，具有更好的整體結構性。

幾何代數（克利福德代數）是基于四元數和格拉斯曼擴張代數等發(fā)展的具有多維統(tǒng)一和維度無關特性的代數系統(tǒng)，已在計算機視覺和人工智能等領域得到廣泛應用。利用幾何代數的旋轉rotor算子，可以實現多維場數據局部特征的自適應模板匹配，進而可構建同時揭示數據整體和局部特征的自適應性通用模板，進而可構建求解場數據特征解析的算法［7－11］。本文通過構建正北方向的模板向量場和柵格高程梯度向量場，利用幾何代數中旋轉算子rotor來擬合兩個向量場進行坡向求解?；诟咚骨婺M數據對方法的正確性、窗口選擇、誤差分布進行了驗證，最終基于實際DEM數據進行了案例驗證與對比。

2 基于幾何代數的坡向表達方法

2.1 梯度場和模板向量場的幾何代數表達

在幾何代數中，空間可以直接被定義為向量集合間的運算，空間維數直接由運算法則確定。給定幾何代數空間Cln，空間格網場F中的任一向量a可定義為：a＝xae1＋yae2＋…＋waen，其中e1、e2、…、en為幾何代數空間Cln的基向量，xa、ya、…、wa為a在各基 向量 上的 投影大?。?2－13］。幾何代數中幾何積的結果為一個混合維度的向量，也叫做多重向量，多重向量是幾何代數的基本數學結構，其定義為

式中，A＝｛n，s，…，t｝。可定義以多重向量為基本元素的多重向量場函數f（M）為

上式將不同維度向量間計算直接統(tǒng)一于幾何積，實現了不同維度場的統(tǒng)一表達與計算。

對于柵格DEM數據，可以將其視為簡單的標量場數據（如圖1）。通過DEM的高程值來構建幾何代數中的多維向量（multivector），每個柵格點都對應一個multivector。坡向的高程變化梯度即高程變化率，可以分別求X方向和Y方向的高程變化梯度grad（x）（x，y）和grad（y）（x，y），利用文獻［14］介紹的中心梯度求解模板，X方向和Y方向的梯度求解公式為

式中，z（x，y）表示點（x，y）處的高程；d表示格網精度；grad（x）（x，y）表示點（x，y）處X方向的梯度；grad（y）（x，y）表示點（x，y）處Y方向的梯度；grad（x，y）表示點（x，y）的梯度。對于多重向量場f（x），上式可寫成如下形式

2.2 基于平行模板的坡向求解

對DEM的柵格點，可以通過構建以該點為中心和行列為奇數的窗口來表示目標梯度場。定義平行模板g（x，y）＝e1sinθ＋e2cosθ，其中角度θ用于限定模板朝向，圖2所示為窗口大小為3，朝向正北方向（θ＝0°）的平行模板。根據rotor的定義，多重向量場f（x，y）可寫成如下形式

式中，R表示由模板g（x，y）變換為f（x，y）的rotor；ε（x，y）為殘差項，由rotor的定義可知，當選取g（x，y）指向正北方向時，R即表示由正北方向變換到f（x，y）的旋轉，則坡向θ的求解可轉化為對rotor的求解。

圖1 高程標量場和正北方向模板Fig.1 The scalar field of elevation and the northern template

3 基于幾何代數的坡向求解與提取

3.1 基于矩陣分解SVD方法的rotor反解

幾何代數中的旋轉算子rotor具有可合成性、保形性和自反性等特性，可實現任意幾何體旋轉變 換 的 多 維 統(tǒng) 一 表 達［15－16］。 給 定 任 意 multivector對象和rotorR，其旋轉變換為

由于rotor的整體保形性，可將上述公式推廣至向量場，給定向量場F＝ ｛a1，a2，…，an｝，其旋轉變換為RFR－1＝ （Ra1R－1）＋ （Ra2R－1）＋…＋（RanR－1），R作用于向量ai的旋轉包括兩個要素：基平面和角度，基平面的對偶向量即ai繞過的旋轉軸。將R寫為指數形式，可顯式求得其所內蘊的旋轉軸和旋轉角度

式中，I為旋轉平面；φ為旋轉角度；p、q為旋轉平面上任意兩個角度為φ的vector；I也可用兩向量的外積p∧q來表示。

定義n維幾何代數空間Cln，rotorR可用幾何代數方程RAR－1＝B表示，也可記為RA＝BR。在該形式下等同于特征多重向量方程RA＝λR，λ為標量。如果存在一系列的多重向量JB：＝｛J∈Cln：JB＝BJ｝滿足（JBR）A－B（JBR）＝0，則JB為rotor的解空間，且JB的維度是2n。在三維歐式空間中，rotor是一個標量和一個二維向量的和，定義兩個單位向量a、b∈Cl3。需要求解R使其滿足RaR－1＝b，由于rotor可以寫作R＝eUabθ／2，Uab是bivector表示的單位平面的二重矢量，θ是旋轉角度。通過求解式中Ra－bR＝0的R，可得結果集JR：＝｛Rab，bRab，UbRab，IRab｝，其中I是Cl3空間中的偽標量并且Ub＝bI－1，解集JR中最優(yōu)解即為所求。

對標準向量場F′＝ ｛a1，a2，…，an｝和目標向量場F＝ ｛b1，b2，…，bn｝，可根據F＝RF′R－1來求解rotorR。構建向量場矩陣F，求解rotorR的最優(yōu)解可以轉化成對矩陣F的SVD分解，求得酉矩陣U、V，并通過公式R＝VUT求解最佳擬合rotor［11］。矩陣F的SVD分解為

3.2 基于幾何代數解析坡向角度

在幾何代數空間中，rotor算子不僅實現了任意對象的旋轉變化，其表達中也包含了變換的角度與旋轉平面等信息。據式（8）所示的rotor的指數形式表達，有如下關系

式中，θ和i2分別為rotor的旋轉角度和旋轉平面，再利用取維度算子〈〉n，旋轉角度θ可通過下式求得

解析出角度θ的余弦值，反解三角函數cosθ／2得到θ的角度值，即為坡向角度。將求得角度按照8個方向對柵格坡向角度進行分類賦值（如表1），即可得各位置點坡向值。

表1 坡向值分類Tab.1 Classification of slope value

4 算法設計與試驗驗證

4.1 基于幾何代數坡向提取算法流程

根據上述分析，基于幾何代數坡向提取算法的基本步驟如圖2。

圖2 幾何代數提取坡向算法流程Fig.2 The algorithm flow of slope extracting based on geometric algebra

（1）根據窗口大小n和公式（5），求解n×n窗口的每個柵格的梯度grad（x，y），構建梯度向量場Fi。

（2）構建同樣窗口大小的正北方向的向量模板M，初始化相關參數。

（3）利用正北模板向量場M和梯度向量場Fi基于rotor擬合逼近，建立關系式，Fi＝RMR－1，基于SVD反解rotor，對該rotor解析，獲取rotor的類型及特征參數，通過對角度值θi的分類，得到n×n窗口中心柵格點的坡向。

（4）移動窗口，迭代數據，反解并解析rotor，得到θ并存儲在柵格中，直至所有數據都被處理得到全部柵格的坡向。

4.2 模擬試驗與驗證

傳統(tǒng)的DEM坡向計算均是基于高程變化率的窗口卷積運算，其計算模型為

式中，fx、fy分別是東西和方南北方向的高程變化率。根據求解fx和fy的不同而演化成不同方法?，F有算法中三階差分算法有平滑和過濾作用，避免極端數據影響，誤差相對較小，可靠性較高［17－18］，被 ArcGIS等行業(yè)軟件普遍采用。本文選取了三階反距離平方權差分法進行算法對比。

為考察算法的正確性，試驗應用模擬高斯地形合成曲面，并離散為特定尺度的DEM曲面，高斯合成曲面參數方程的定義為

式中，A、B、C為地勢起伏參數；m、n為范圍控制參數?；跀祵W曲面建立的DEM數據由于消除了格網點數據誤差，同時可根據Z的表達式直接求解出x和y方向上的偏導數fx和fy，結合公式（12）求得坡向真值，即可支撐算法對比與精度驗證［19－20］。

取A＝3，B＝10，C＝1／3，m＝500，n＝500，標準差為1.314 6，均值為0.627 1，格網分辨率為10×10構建地形曲面（圖3a）。在3×3、5×5、7×7窗口下，利用rotor算子擬合和基于ArcGIS的三階反距離平方權差分進行坡向計算，兩種算法三種窗口的計算結果與真值坡向均一致，表明算法計算結果的正確性。

為驗證算法的穩(wěn)定性，對模擬DEM數據分別添加0.02、0.04、0.06、0.08、0.1倍標準差的高斯噪聲，對加噪聲DEM數據同樣在3×3、5×5、7×7窗口下應用兩種方法分別計算坡向，與真值坡向的比較分析結果如圖3（b）所示。結果表明，隨噪聲增加，三階反距離平方權差分法的精度變低，表明當地形較破碎時，算法的穩(wěn)定性降低，噪聲干擾的影響加大，與真值坡向的一致性也變差。而幾何代數算法受到的干擾相對要小，且窗口越大受到噪聲的干擾越小。表明本文算法的抗噪性較反距離平方法要好。

圖3 高斯地形合成曲面和坡向求解結果對比Fig.3 The topography surface composited by Gaussian function and the compare of slope solving results

為驗證方法的普適性、抗噪能力以及不同窗口大小下坡向精度的分布特征，對添加數據0.02倍標準差噪音的每組DEM數據，分別提取不同大小窗口中各曲面的山脊線和山谷線，坡向計算誤差點的空間分布（圖4）。對于不同的DEM數據，各種窗口的算法對地形復雜度不敏感，正確率大致相當。在誤差點的空間分布上，誤差點隨著窗口增大而減少，并且總體構型和山脊線和山谷線吻合。這可能與山脊線和山谷線是兩個不同坡向地形面的交線有關?？傮w上，窗口越大，正確率越高，保形性越好，但窗口過大容易導致邊界數據缺失?；诓煌翱诖笮〉亩啻文M試驗結果顯示，當原始采樣數據精度較高，需要細節(jié)的坡向分布時，3×3窗口即可獲得較好的效果。當數據中存在明顯噪音時，選取5×5或7×7窗口可獲得很高精度。

圖4 不同場景下錯誤點位的空間分布Fig.4 The distribution of slope error computed by different algorithms

4.3 DEM坡向求解案例驗證

為驗證算法在不同窗口下的穩(wěn)定性和適用性，選取不同窗口的坡向求解模板，對黃土高原地區(qū)丘陵、梁、峁等典型地貌類型5m×5m分辨率的DEM數據進行統(tǒng)計和對比分析。選用方向統(tǒng)計量計算工具包CircStat［21］分別計算坡向數據的均值、標準差、峰度和偏度。計算結果如圖5和表2所示。各地貌類型的坡向統(tǒng)計值均存在躍變過程（圖中虛線），且不同統(tǒng)計量發(fā)生躍變的窗口大小大體一致?？紤]到小窗口對地形保真度較好，而大窗口去噪能力較強，該躍變大致對應于坡向計算的最優(yōu)窗口。黃土丘陵、梁、峁、塬和風沙地貌的最優(yōu)窗口大小分別是13、11、11、11和17。

圖5 真實地形數據坡向結果Fig.5 The slope results of real terrain data

圖6 不同窗口坡向統(tǒng)計量結果Fig.6 The statistical result of terrains slope with different window size

5 結 論

坡向作為向量的方向，常用算法多將其視為標量來計算，基于向量場擬合方法對求解實際地形坡向更具有現實意義。本文通過幾何代數模板匹配方法構建了用于地形坡向計算的標準模板與地形變化率向量場模板，用rotor旋轉算子進行向量場擬合，反解rotor進行角度解析獲取地表坡向值。結果顯示，該算法較反比例平方權算法有更好的整體性，且抗噪性好?；?種典型地形的實例分析顯示，黃土丘陵、梁、峁、塬和風沙地貌的最優(yōu)窗口大小分別是13、11、11、11、17。黃土梁、黃土塬地貌一致性較強，所需計算窗口較小，而風沙和峁由于細節(jié)結構相對較多，選取更大的窗口更有助于獲得穩(wěn)定的結果。而黃土丘陵的窗口則介于上述兩類之間。由于窗口越大，計算復雜度增加，計算效率趨于降低，建立面向大范圍坡向計算的并行算法并進行經驗參數率定是本文未來的研究方向。

［1］ AI Tinghua，LIU Yaolin，HUANG Yafeng.The Hierarchical Watershed Partitioning and Generalization of River Network ［J］.Acta Geodaetica et Cartographica Sinica，2007，36（2）：231－236，243.（艾廷華，劉耀林，黃亞鋒.河網匯水區(qū)域的層次化剖分與地圖綜合［J］.測繪學報，2007，36（2）：231－236，243.）

［2］ CHEN Nan，WANG Qinmin，TANG Guoan.Comparative Analysis of the Algorithms of Aspect Deriving from DEM：Take the Study of Loess Hill and Gully Area as the Case［J］Remote Sensing Information，2007，（1）：70－75.（陳楠，王欽敏，湯國安.基于DEM的坡向提取算法對比分析：以黃土丘陵溝壑區(qū)的研究為例［J］.遙感信息，2007，（1）：70－75.）

［3］ LIU Xuejun，BIAN Lu，LU Huaxing，et al.The Accuracy Assessment on Slope Algorithms with DEM Error Spatial Autocorrelation［J］.Acta Geodaetica et Cartographica Sinica，2008，37（2）：200－206.（劉學軍，卞璐，盧興華，等.顧及DEM誤差自相關的坡度計算模型精度分析［J］.測繪學報，2008，37（2）：200－206.）

［4］ LU Huanxing，LIU Xuejun，JIN Bei.Stochastic Process Based Accuracy Field Model for Grid DEM ［J］.Acta Geodaetica Et Cartographica Sinica，2012，41（2）：273－277.（盧華興，劉學軍，晉蓓.基于隨機過程的格網DEM精度場模型［J］.測繪學報，2012，41（2）：273－277.）

［5］ LI Tianwen ，LIU Xuejun，TANG Guoan.Influence of Terrain Complexity on Slope and Aspect［J］.Journal of Mountain Research，2004，22（3）：272－277.（李天文，劉學軍，湯國安.地形復雜度對坡度坡向的影響［J］.山地學報.2004，22（3）：272－277.）

［6］ BI Xiaoling，LI Xiaojuan，HU Zhuowei，et al.Analysis of the Influence of DEM Grid Size on the Accuracy of Topographical Factors［J］.Science of Surveying and Mapping，2012，37（6）：150－152.（畢曉玲，李小娟，胡卓瑋，等.DEM網格尺寸對地形因子精度的影響分析［J］.測繪科學，2012，37（6）：150－152.）

［7］ ZHANG Juqing，LIU Pingzhi.Combining Fitting Based on Robust Trend Surface and Orthogonal Multiquadrics with Application in DEM Fitting［J］.Acta Geodaetica et Cartographica Sinica，2008，37（4）：526－530.（張菊清，劉平芝.抗差趨勢面與正交多面函數結合擬合DEM數據［J］.測繪學報，2008，37（4）：526－530.）

［8］ GREEN A C，MARSHALL S，GREENHALGH D，et al.Design of Multi－mask Aperture Filters［J］.Signal Processing，2003，83（9）：1961－1971.

［9］ YANG J，LI S.Smoothness of Multivariate Refinable Functions with Infinitely Supported Masks［J］.Journal of Approximation Theory，2010，162（6）：1279－1293.

［10］ MA T，WANG S.Structural Classification and Stability of Divergence－free Vector Fields［J］.Physica D：Nonlinear Phenomena，2002，171（1－2）：107－126.

［11］ LUO Wen，YUAN Linwang，YU Zhaoyuan，et al.Adaptive Template Matching Method for Convergence and Divergence Characteristics of Multi－Dimensional Vector Fields［J］.Acta Electronica Sinica，2012，（9）：1739－1734.（羅文，袁林旺，俞肇元，等.多維向量場輻散輻合結構特征自適應匹配方法［J］.電子學報，2012，（9）：1739－1734.）

［12］ DORST L，FONTIJNE D，MANN S.Geometric Algebra for Computer Science［M］∥The Morgan Kaufmann Series in Computer Graphics.San Mateo：Morgan Kaufmann，2008.

［13］ PERWASS C，LASENBY J.A Unified Description of Multiple View Geometry［C］∥Proceedings of Geometric Computing with Clifford Algebra（G.Som－mer，ed.）.Berlin：Springer－Verlag，2001.

［14］ EBLING J，SCHEUERMANN G，CLIFFORD.Convolution and Pattern Matching on Vector Fields［C］∥ Proceedings of IEEE Visualization＇03，IEEE Computer Society.Los Alamitos：［s.n.］，2003：193－200.

［15］ DORAN C，LASENBY A.Geometric Algebra for Physicists［M］.Cambridge：Cambridge University Press，2003.

［16］ HESTENES D.New Foundations for Classical Mechanics［M］.New York：Kluwer，2002.

［17］ LI Tianwen，LIU Xuejun，CHENG Zhengjiang，et al.Study on the Accuracy and Algorithms for Calculating Slopes and Aspects Based on the Digital Elevation Model［J］.Arid Land Geography，2004，27（3）：398－403.（李天文，劉學軍，陳正江，等.規(guī)則格網DEM坡度坡向算法的比較分析［J］.干旱區(qū)地理，2004，27（3）：398－403.）

［18］ 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 Cartographic Sinica，2004，33（3）：258－263.（劉學軍，龔健雅，周啟鳴，等.基于DEM坡度坡向算法精度的分析研究［J］.測繪學報，2004，33（3）：258－263.）

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

［20］ JIA Yini，TANG Guoan，LIU Xuejun，et al.The Impact of Elevation Interpolation on the Accuracy of Gradient and Aspect from DEMs［J］.Journal of Geo－Information Science，2009，11（1）：36－41.（賈旖旎，湯國安，劉學軍.高程內插方法對DEM所提取坡度、坡向精度的影響［J］.地球信息科學學報，2009，11（1）：36－41.）

［21］ BERENS P.CircStat：A MATLAB Toolbox for Circular Statistics［J］.Journal of Statistical Software，2009，31（10）：1－21.