基于不同算法的剖面曲率提取研究

1、第15頁（共19頁）基于不同算法的剖面曲率提取研究引言地形因子是為定量表達地貌形態(tài)特征而設(shè)定的具有一定數(shù)學(xué)意義的數(shù)學(xué)參數(shù)或指標(biāo)1。剖面曲率（Profile Curvature）KV 是地形因子之一，基于提取算法的坡面因子分，剖面曲率屬于二階坡面因子。在地理信息系統(tǒng)（GIS）中，剖面曲率的提取算法有三種：湯國安提出的SOS/SOA算法、四次表面模型算法和基于空間矢量模型的差分計算算法。數(shù)字高程模型（DEM）作為地理信息系統(tǒng)( GIS)空間數(shù)據(jù)庫的核心數(shù)據(jù)庫之一，是進行二維地形空間模擬和地形圖生產(chǎn)的基礎(chǔ)，DEM 能派生出各種地形因子，如坡度、坡向、平面曲率、剖面曲率、匯水面積等1。其中，剖面曲率是

2、最重要的地形因子之一，地面曲率影響著，也是制約生產(chǎn)力空間布局及自然地理研究的重要因子 。盡管利用GIS 平臺在DEM 上提取坡度的算法已經(jīng)成熟，但是由于算法會影響坡度的測算，所以需要比較不同算法對提取坡度的影響。與坡度、坡向等基本地形參數(shù)一樣，剖面曲率在數(shù)學(xué)定義上是明確的和唯一的，具有確定的表達式。盡管其理論定義是明確的，然而DEM是地形曲面的微分模擬，算法設(shè)計必然存在各種各樣的假設(shè)，不同假設(shè)和前提導(dǎo)致不同的剖面曲率計算模型和結(jié)果，這雖然對地形特征的可視化和地形分類的影響不大，但對以數(shù)值計算為主的地學(xué)分析模型的影響卻是非常顯著的2。剖面曲率的含義剖面曲率是地面曲率在垂直方向上的分量，是對地面坡

3、度的沿最大坡降方向的地面高程變化率的度量，是研究地形表面的重要因子之一。剖面曲率因子提取算法的基本原理為：將離散的高程數(shù)值擬為一個連續(xù)的曲面，基于微分幾何的思想，模擬曲面上每一個點所處的平行于水平面的曲線，利用曲線曲率的求算方法推導(dǎo)得出各個曲率因子的公式，進而求算出每一點的曲率值1。剖面曲率KV 具有重要的理論和實踐意義。KV 決定地表及土壤中物質(zhì)移動的相對速度，KV 0表明移動加速，KV 0表明移動減速，反映了侵蝕或沉積的程度。KV 通?？梢越沂境龉燃咕€呈突起狀或階梯狀，還可以用在很多領(lǐng)域，如在大尺度區(qū)域中模擬土壤濕度、土層厚度、土壤侵蝕、山崩石頭的散布、植被等；也可用于不同尺度的地貌研究，

4、研究區(qū)域尺度和次大陸尺度地質(zhì)構(gòu)造、地貌和土壤的空間相關(guān)規(guī)律；還用于識別顯示斷層、地質(zhì)斷陷線等3。笛卡爾坐標(biāo)系中，x、y 表示空間位置，z = f( x，y) 表示地形曲面函數(shù)，假定函數(shù)具有單值性、連續(xù)性和二階可導(dǎo)性，則剖面曲率的計算公式為：Kv=-p2r+2pqs+q2tp2+q2(1+p2+q2)3/2 （1）其中，p=zx ，是x方向高程變化率；q=zy ，是y方向高程變化率；r=2zx2 ，是x方向高程變化率的變化率；s=2zxy ，是x方向高程變化率在y方向的變化率；t=2zy2 ，是y方向高程變化率的變化率剖面曲率提取算法一般都是根據(jù)P.A.Borrough提出的窗口微分法遍歷整個圖

5、幅進行計算的。即是在33的DEM柵格窗口中進行的，窗口在DEM數(shù)據(jù)矩陣中連續(xù)移動后而能完成整幅圖的計算工作4。Z1Z2Z3Z4Z5Z6Z7Z8Z9圖1 33地形窗口(Z1Z9分別代表各個柵格的高程值)剖面曲率的計算方法GIS中剖面曲率的計算大多是基于DEM數(shù)據(jù)的，由離散的數(shù)據(jù)計算剖面曲率有不同算法，在考慮到算法的通用性和使用性的情況下，這里選擇了ArcGIS采用的湯國安提出的SOS/SOA算法、四次表面模型算法和基于空間矢量模型的差分計算算法1三種常用的剖面曲率算法進行分析。湯國安提出的SOS算法剖面曲率為地面高程變化的二次導(dǎo)數(shù)，即坡度在垂直方向的最大變化率；Tang已經(jīng)證明5，在Arc/In

6、fo等GIS軟件中，地面剖面曲率的數(shù)字矩陣可以直接通過對DEM數(shù)據(jù)的求取地面坡度的坡度（Slope of Slope,SOS）而獲得。SOS就是坡度的變率，是在所提取的坡度值的基礎(chǔ)上對地面每一點再求算一次坡度，坡度是地面高程變化率的求解，因此，坡度變率表示了地球表面高程相對于水平面變化的二階導(dǎo)數(shù)3。這里的“坡度”并非是坡度，而是按照坡度的算法實現(xiàn)了對坡度數(shù)據(jù)陣列中其變化率的量化提取。所獲得的剖面曲率值在0,90內(nèi)，并非是真正的曲率值，但卻真實的地反映了坡度和等高線的變化率，可提取山谷、山脊等地形因子，而且因其算法簡單易于實現(xiàn)，所以也稱剖面曲率，有些教材中稱之為地面變率因子3。3.2 四次表面模

7、型Zeverbergen 和 Thorne(1987)提出了一種計算算法，構(gòu)造一個含有9個參數(shù)的四次方程（x、y的最高次數(shù)為2）6，見公式（2）.函數(shù)的所有參數(shù)由33窗口（如圖1）內(nèi)的單元值唯一確定，函數(shù)經(jīng)過每個數(shù)據(jù)點，包括中心單元。系數(shù)A、B、C、I在計算曲率時未用到，可以忽略，系數(shù)DH由公式（3）（7）計算，則中心格網(wǎng)的剖面曲率的表達式變?yōu)楸砻婧瘮?shù)系數(shù)組成的式子，見公式（8）3。Z = Ax2y2 + Bx2y + Cxy2 + Dx2 + Ey2 + Fxy + Gx + Hy + I. (2)D =（Z4 + Z6）/2-Z5/L2 . (3)E =（Z2 + Z8）/2-Z5/L2

8、. (4)F =(- Z1 + Z3 +Z7 Z9)/4L2. (5)G =(- Z4 + Z6)/2L. (6)H =(Z2 Z8)/2L . (7)剖面曲率：KV = 2(DG2 + EH2 + FGH)/(G2 + H2). (8)3.3 基于空間矢量模型的差分計算 在DEM數(shù)據(jù)的基礎(chǔ)上，從微分幾何的思想出發(fā)，求出剖面曲率公式(1)中的各個導(dǎo)數(shù)pt,即可求出各曲率值，差分計算常用的是三階反距離平均權(quán)差分。對每個柵格點都確定一個33分析窗口，如圖2（1），ai分別代表各個柵格的高程值，L表示柵格大小，p、q可由以下公式求出3。p=zx=c+2f+i-(a+2d+g)8L （9）p=zx=a

9、+2b+c-(g+2h+i)8L (10)abcaxbxcxaybxcxdefdxexfxdxexfxghigxhxixgxhxix （1） （2） （3）圖2 提取地面曲率的分析窗口 求出所有 p、q 值，則組成 p 值柵格圖（ 見圖1（ 2 ） ） 和 q 值柵格圖（ 見圖1（ 3 ） ） 。 以此數(shù)據(jù)為基礎(chǔ)，進行差分運算，可求出r、s、t8,見公式( 10 ) 、（ 11 ） 和（ 12 ） 。r=2zx2=cx+2fx+ix-(ax+2dx+gx)8L （11）s=2zxy=ax+2bx+cx-(gx+2hx+ix)8L （12）t=2zy2=ay+2by+cy-(gy+2hy+iy)

10、8L （13）由不同算法提取的剖面曲率由任何一種算法得出的剖面曲率，我們都要評估其精確性和制圖效果，并以此來作為我們求算地形的剖面曲率的依據(jù)。因此，一個科學(xué)、合理、準(zhǔn)確的剖面曲率分析是目前研究剖面曲率的首要任務(wù)9。在此，以西峽縣的某一地區(qū)為研究區(qū)，我將依據(jù)三種剖面曲率算法所作出的圖像和數(shù)據(jù)進行比較分析，得到柱狀圖和各統(tǒng)計量，為找到不同高程域的適合剖面曲率提取算法提供依據(jù)。4.1 研究區(qū)概況研究區(qū)地處伏牛山腹地，豫、鄂、陜交界地區(qū)，東臨南陽盆地，西嵌秦嶺，南襟荊襄平原，北負(fù)八百里伏牛山。境內(nèi)地形復(fù)雜，北部是海拔高、坡度大的中低山地，南部是鸛河谷地，兩側(cè)是起伏大的低山丘嶺。海拔181米，自然坡降為

11、33%。研究區(qū)內(nèi)的高程分布如表1所示。由表1知，分布在120405之間的高程數(shù)值最多，為3197380個，約占總數(shù)值個數(shù)的27.9%，分布在405690之間的高程數(shù)值約占總數(shù)值個數(shù)的22.5%，分布在690991之間的高程數(shù)值約占總數(shù)值個數(shù)的21.5%，分布在9911316之間的高程數(shù)值約占總數(shù)值個數(shù)的17.8%，而分布在13162171之間的高程數(shù)值最少，為1186783個，約占總數(shù)值個數(shù)的10.4%。并且根據(jù)圖3中所示可判斷該研究區(qū)的地形呈階梯性抬高，總體上是北高南低。表1 原DEM的高程值分布范圍圖3 原DEM重分類圖圖4和圖5中的三維圖是在ArcScene中制作的，原DEM的高程變化不

12、是很明顯，僅從顏色的分類可以知道高程的分布，和地形的起伏。但用高程值的3倍顯示，則可以看出此研究區(qū)的地形有一定的起伏度，但不是很大。圖4 原DEM的三維顯示側(cè)視圖圖5 原DEM高程的3倍的三維顯示側(cè)視圖由湯國安提出的SOS算法提取的剖面曲率由圖6知，由SOS算法提取的剖面曲率的圖像的南部，有兩道狹長而傾斜的淺顏色部分，這部分的變率在0.0013778 547.339714169之間，說明這部分的剖面曲率變化比較小，地形相對平坦。圖6 湯國安的SOS算法剖面曲率4.3 由四次表面模型算法提取的剖面曲率由圖7知，由四次表面模型算法提取的剖面曲率的圖像相比圖6的相對區(qū)域，圖7的這兩個狹長部分的變率則

13、在-0.7680859470.646213027之間，剖面曲率的變化比圖3變化更小，而且其間的其他數(shù)值范圍的變率也較之增加。圖7 四次表面模型算法剖面曲率4.4 由基于空間矢量模型的差分計算算法提取的剖面曲率由圖8知，由基于空間矢量模型的差分計算算法提取的剖面曲率的圖像中除了空數(shù)值為零外，幾乎所有的變率數(shù)值都分布在-0.1462174060.000303681之間，只有零星的部分分布在0.0003036810.027690800之間。與前兩種算法相比，此算法的差別較大。圖8 基于空間矢量模型的差分計算算法剖面曲率試驗結(jié)果及分析5.1 不同算法提取剖面曲率的統(tǒng)計指標(biāo)分析由圖9可知，由湯國安的SO

14、S算法提取的剖面曲率的數(shù)值偏正態(tài)分布，主要集中在2.335129512.303019之間，雖然數(shù)值不代表真正的曲率值，但能夠表示曲率的變化率。圖中的數(shù)值多處于圖像的正區(qū)，表明用此算法提取的該研究區(qū)的剖面曲率偏大，且地形多為凸地形，呈階梯性抬高地表。圖9 由SOS算法得出的剖面曲率柱狀圖由圖10知，由四次表面模型算法提取的剖面曲率的數(shù)值是正規(guī)的正態(tài)分布圖，且正區(qū)、負(fù)區(qū)基本持平，大多集中在-1.0182720.793010之間，表明用此算法提取的該研究區(qū)的地形剖面凸地形和凹地形分布均勻。圖10 由四次表面模型算法得出的剖面曲率柱狀圖由圖11知，由基于空間矢量的差分計算算法提取的剖面曲率的部分?jǐn)?shù)值發(fā)

15、生異常，導(dǎo)致數(shù)據(jù)分析失真，但從分布的數(shù)值來看，數(shù)值主要集中在-0.003119700.00030368之間，數(shù)值多分布在負(fù)區(qū)，表明用此算法提取的該研究區(qū)地形剖面凹地形多與凸地形。圖11 由基于空間矢量模型的差分計算算法得出的柱狀圖由表2知，湯國安的SOS算法的最大值與最小值之差最大，說明其值域幅度最大；基于空間矢量模型的差分計算算法的最大值與最小值之差最小，說明其值域幅度最小。而從標(biāo)準(zhǔn)差來看，基于空間矢量模型的差分計算算法最小，準(zhǔn)確度最高；四次表面模型算法的準(zhǔn)確度居中；湯國安的SOS算法標(biāo)準(zhǔn)差最大，準(zhǔn)確度最差；用標(biāo)準(zhǔn)差除以平均值，得出其各自的變異系數(shù)，SOS算法的變異系數(shù)為0.66459639

16、74616885，四次表面模型算法的變異系數(shù)為-32.43107466561735，基于空間矢量模型的差分計算算法的變異系數(shù)為-1.172042558940668。表2 三種算法的統(tǒng)計量比較算法最小值最大值和平均值標(biāo)準(zhǔn)差SOS0.00137785445.64225006102631668.98.9639350395.957398934四次表面模型-13.29473412.46571159-369645.7671-0.0322851681.047042694基于空間矢量模型的差分計算-0.1462174060.027690800-13322.68346-0.0015469370.00181307

17、65.2 基于三種算法提取的剖面曲率的優(yōu)缺點、誤差及制圖效果分析對運用不同算法提取的剖面曲率數(shù)據(jù)做統(tǒng)計分析,發(fā)現(xiàn)同一研究區(qū)不同算法所獲得的剖面曲率最大數(shù)值、最小數(shù)值、標(biāo)準(zhǔn)差是不同的。并且，通過對各算法所提取出的剖面曲率做重分類，更顯示出其差異。由圖12可知，由湯國安算法提取的剖面曲率數(shù)值分布在0.0013783.581054之間的數(shù)值比較集中，且數(shù)值最小，并與前面圖3、圖4、圖5相對應(yīng)，表明用此算法提取這一部分區(qū)域的剖面曲率是適合的，與原地形相比，在研究區(qū)南部的低平地區(qū)，制圖效果明顯，說明此算法在提取這種簡單地形有一定的優(yōu)勢。圖12 SOS算法提取的剖面曲率重分類由圖13可知，在圖12的相對應(yīng)

18、區(qū)域，用此算法提取的剖面曲率數(shù)值大約在-0.9701290.242128之間，數(shù)值偏小，且中間夾雜著其它類別的數(shù)據(jù)，并與前面圖3、圖4、圖5所示有一定的偏差，表明此算法不太適合用于該研究區(qū)的剖面曲率提取，制圖效果仍有一定的明顯性。圖13 四次表面模型算法的剖面曲率重分類由圖14可知，由基于空間矢量差分計算算法提取的剖面曲率的圖像在圖12、圖13的對應(yīng)區(qū)域，剖面曲率數(shù)值分布在-0.0010660.000988之間，中間也有其它類別的數(shù)值，比圖12顯示的數(shù)值稍小，表明制圖效果的地形起伏不太大，但也不平坦。圖14 基于空間矢量的差分計算算法的剖面曲率重分類6 結(jié)論通過在不考慮DEM來源誤差的情況下，

19、比較三種算法在研究區(qū)提取剖面曲率圖，并統(tǒng)計各剖面曲率的最大值、最小值、平均值和標(biāo)準(zhǔn)差，以及對各算法提取的剖面曲率做重分類，進行綜合分析。從三種算法提取的剖面曲率重分類圖像來看，曲率數(shù)值值域最大的是湯國安的SOS算法，曲率數(shù)值值域最小的是基于空間矢量的差分計算算法，但是這兩種算法都能凸顯出該研究區(qū)部分區(qū)域的剖面曲率的特征。由圖3的對原高程圖像進行的重分類可知，該研究區(qū)的大致地形是由南部的最低處向兩邊逐級抬升，有三個最高區(qū)域，兩個稍高區(qū)域，兩個低谷區(qū)域和一個平坦區(qū)域。通過比較圖像特征較明顯的區(qū)域，即平坦區(qū)域，認(rèn)為在該研究區(qū)提取剖面曲率選用湯國安的SOS算法、四次表面模型算法較為合理，而湯國安的SO

20、S算法在計算原理與方法上比四次表面模型算法更精細(xì)、更有特征性、邏輯性更強。雖然湯國安的SOS算法所計算的是坡度的坡度，不是真正的剖面曲率，但這種算法反映了剖面的變化率，且簡單易操作，更能準(zhǔn)確的突出研究區(qū)的特征。而四次表面模型算法所提取的剖面曲率在一定程度上體現(xiàn)了研究區(qū)的地形特征，具有極大的參考性。因此，在該研究區(qū)進行的基于不同算法剖面曲率的統(tǒng)計量分析時使用這兩種算法是科學(xué)的，由表2知，通過數(shù)據(jù)對比，該研究區(qū)的凹地形地區(qū)多于凸地形地區(qū)，而那部分平坦地區(qū)類似于盆地地形。本文對算法的效率問題在提取剖面曲率時沒有考慮，且僅分析了3種算法的成圖數(shù)據(jù)統(tǒng)計量，DEM的精度、DEM的比例尺、DEM的分辨率、數(shù)

21、據(jù)的組織方式等都對地形的剖面曲率提取具有重要影響，因此，對于使用不同算法提取剖面曲率的研究，有待考慮更多的因素，發(fā)現(xiàn)更加普遍的規(guī)律。參 考 文 獻1 湯國安，李發(fā)源，劉學(xué)軍等.數(shù)字高程模型 M . 北京：科學(xué)出版社，2011:142157.2 李天文，劉學(xué)軍，陳正江，湯國安，李軍鋒等.規(guī)則格網(wǎng)DEM坡度坡向算法的比較分析 J . 干旱區(qū)地理，2004, 27( 3) : 398 - 403. 3 范紅艷，戚鵬程，黃天勇等.關(guān)于GIS中三種地面曲率的探討 J .南陽師范學(xué)院學(xué)報，2011, 10( 6):66-69.4 劉學(xué)軍，龔健雅，周啟鳴，湯國安等.基于DEM 坡度坡向算法精度的分析研究 J

22、 .測繪學(xué)報，2004, 33( 3) : 258 - 263.5 陳楠，王欽敏，湯國安,趙元偉等.6種坡度提取算法的應(yīng)用范圍分析以在黃土丘陵溝壑區(qū)的研究為例 J .測繪信息工程，2006, 31( 4) : 20 - 21.6 Tang Guo-an A Research on the Accuracy of Digital Elevation ModelsM Beijing: Science Press Beijing， 20007 Zeverbergen L W，Thorne C R Quantitative Analysis of Land Surface TopographyJ Ea

23、rth Surface Processes and Landforms， 1987， 12( 1) : 47 568 劉學(xué)軍，王葉飛，曹志東，湯國安等.基于DEM的地形曲率計算模型誤差分析 J . 測繪學(xué)報，2006, 31( 5) : 50 - 53.9 陳楠，王欽敏，湯國安等.自DEM由不同算法提取坡度的對比分析 J .測繪工程，2006，15( 1) : 6 - 9. profile curvature extraction based on different algorithm researchAbstract:Without consideration of DEM error, this article through