牛亮,楊勤科

(西北大學(xué)城市與環(huán)境學(xué)院,西安710127)

研究表明,DEM及其基礎(chǔ)上提取的地形要素因受到水平分辨率的強烈影響,表現(xiàn)出明顯的尺度效應(yīng)。認(rèn)識理解尺度效應(yīng),是有效利用DEM數(shù)據(jù)資源,準(zhǔn)確提取適應(yīng)與應(yīng)用尺度和目的地形數(shù)據(jù)的基礎(chǔ)[1]?；诙喾N比例尺地形圖生成的多種分辨率的DEM(原生DEM)為尺度效應(yīng)研究提供了基礎(chǔ),但由于不同比例尺地形圖之間并不具備完全相同的定位控制基礎(chǔ),限制了研究。為此,研究者通過對高分辨率DEM進(jìn)行變換方法生成多種分辨率(或者柵格)的DEM 數(shù)據(jù)。常見方法包括重采樣[2-4]、小波變換[5-6]、空間濾波[1]等。在這些方法中,如何將變換結(jié)果與原生DEM的分辨率建立聯(lián)系,也就是說如何確認(rèn)變換結(jié)果的分辨率,是一個有待研究解決的問題,其中直方圖的相似性比較是一個有效方法。

基于1∶1萬數(shù)字地形圖數(shù)據(jù)插值建立多種柵格尺寸的 DEM,通過直方圖相似度比較找到與1∶5萬、1∶10萬、1∶25萬原生 DEM 等分辨率的DEM,為坡度的尺度效應(yīng)分析研究和中低分辨率DEM坡度的變換提供依據(jù)。

1 數(shù)據(jù)基礎(chǔ)和研究方法

1.1 研究區(qū)域與數(shù)據(jù)基礎(chǔ)

本研究工作是在黃土高原中部延河流域的二級支流縣南溝流域內(nèi)進(jìn)行的。研究區(qū)為典型的黃土丘陵,地面起伏較大。基于10m分辨率DEM的坡度圖統(tǒng)計表明,坡度大于25°的地面約占總土地面積的50%。研究所用基礎(chǔ)數(shù)據(jù)為:(1)1∶1萬、1∶5萬、1∶10萬和1∶25萬四種不同比例尺的數(shù)字化地形圖,等高距分別為5 m、20 m、40m和100 m,與之對應(yīng)的原生DEM,分辨率應(yīng)分別為5 m、10 m、25 m 和50 m[7-9],標(biāo)記為 DEM5a 、DEM 10a 、DEM 25a、DEM50a 。

1.2 研究方法

本研究基于1∶1萬數(shù)字化地形圖,生成多種柵格尺寸的DEM,然后通過與原生DEM 的對比,使分辨率得到匹配,方法如下。

1.2.1 多種柵格DEM生成 為了弄清地形形態(tài)及其與坡面流水的關(guān)系,本研究采用 ANUDEM[10-11]軟件對1∶1萬數(shù)字化地形圖進(jìn)行空間插值,生成柵格尺寸為10m 、20 m、25 m、30m、40 m、50 m、60 m和70 m的 DEM,分別標(biāo)記為DEM10b、DEM 20b、DEM25b、DEM30b、DEM40b、 DEM 50b、DEM 60b和DEM70b。所用參數(shù)設(shè)置為第二粗糙系數(shù)0.5,迭代次數(shù)40。

1.2.2 DEM相似性比較參數(shù)計算 通常比較兩個DEM之間的相似性和差異性是在比較圖像統(tǒng)計特征基礎(chǔ)上,通過觀察來實現(xiàn)的[1]。本文引入直方圖相似度指數(shù)和逼真度指數(shù),以便客觀進(jìn)行比較。

(1)統(tǒng)計特征:將DEM 看成簡單的矩陣,利用經(jīng)典統(tǒng)計學(xué)統(tǒng)計其特征變量,包括:最小值、最大值、平均值、標(biāo)準(zhǔn)差、眾數(shù)等,并計算信息熵[12]。

(2)直方圖相似度:直方圖是數(shù)字圖像的重要特征之一,反映圖像(或其中一部分區(qū)域)的灰度的總體分布信息。在圖像變換中,變換前后2幅圖像的相似程度通常用直方圖相似性來度量。直方圖相似度指數(shù)包括:

①距離直方圖相似度指數(shù):直方圖相似性度量通常是基于2幅直方圖相同窗口之間的距離或經(jīng)過變換后的一些距離標(biāo)準(zhǔn),如L1距離(曼哈頓距離)、L2(歐幾里得距離)距離、χ2距離等[13],這些方法側(cè)重于衡量兩幅圖像的差異,得到的距離越小表示兩幅圖像越相似。其中L1、L2和 χ2距離的計算公式為:

對于頻率直方圖,L1距離的取值范圍是[0,2],L2距離的取值范圍是,χ2距離的取值范圍是[0,1],值為0時最相似,值為最大值時最不相似。為了統(tǒng)一規(guī)范化相似性距離度量函數(shù)的取值范圍,將兩幅圖像的相似度定義在[0,1]之間,值越接近1則圖像越相似。對于累計頻率直方圖而言,統(tǒng)一規(guī)范化相似性距離度量函數(shù)的取值范圍比較困難,直接采用距離來判定相似性,距離越小越相似。

②相交直方圖相似度指數(shù)。另一種方法稱為直方圖的交[14],對于歸一化的直方圖,其具體計算公式如下:

直方圖的交運算就是取兩幅圖像的特征直方圖在各個數(shù)值級上的較小值,累加即表示圖像之間的相似程度。這種相似性實際上表示兩幅圖像在該直方圖特征意義下的公共部分。0≤D(X,Y)≤1,它的值越大,則兩幅圖像越相似,如果兩幅圖像完全一致 ,則 D(X,Y)=1。

(3)圖像逼真度指數(shù)。圖像逼真度,是描述待對比圖像是否存在較大誤差的一個重要技術(shù)指標(biāo)[15],本文根據(jù)公式(5)計算,以逼真度接近數(shù)值1的程度作為處理前后圖像的近似性度量[15],越接近于1,逼真度越高,反之越小。

式中:I——原生DEM 的坡度和曲率;I′——多柵格DEM 的坡度和曲率;Iij——原生DEM 的坡度和曲率的第i行j列柵格單元的值;I′ij——多柵格DEM的坡度和曲率第i行j列柵格單元值。

1.3 試驗步驟

本研究的試驗只要包括DEM 生成、統(tǒng)計表、直方圖相似度和逼真度計算、分析比較等。工作流程如圖1示。

(1)DEM生成。將經(jīng)過質(zhì)量檢查編輯好的1∶1萬、1∶5萬、1∶10萬、1∶25萬數(shù)字地形圖分別輸入到ANUDEM軟件中,所采用的主要參數(shù)是第二粗糙系數(shù)0.5,迭代次數(shù)40。1∶1萬選擇的柵格尺寸分別為 10m 、20m 、25m 、30 m 、40m 、50m 、60m和 70 m,生成 標(biāo)記為 DEM10b、DEM20b、DEM25b、DEM 30b、DEM 40b、DEM 50b、DEM 60b和 DEM 70b的多柵格DEM 。1∶5萬、1∶10萬和1∶25萬地形圖根據(jù)基于信息含量確定相應(yīng)分辨率分別為10 m、25 m和50 m,生成對應(yīng)分辨率 DEM,分別標(biāo)記為:DEM 10a 、DEM 25a 、DEM 50a 。

(2)坡度和曲率的生成及統(tǒng)計參數(shù)比較。在ARCGIS軟件中分別對所有的DEM求取坡度和曲率,并記錄各坡度和曲率的空間統(tǒng)計特征參數(shù),包括最小值、最大值、平均值、標(biāo)準(zhǔn)差、眾數(shù)和信息熵。通過比較各統(tǒng)計特征參數(shù),找到統(tǒng)計特征參數(shù)最為相似的兩個坡度或者曲率。

(3)坡度和曲率直方圖相似度計算。生成坡度和曲率做頻率和累計頻率圖。頻率和累計頻率圖實質(zhì)上也是直方圖的表現(xiàn)形式。通過直方圖相似度L1距離、L2距離、χ2距離和直方圖交四種方法分別計算基本比例尺1∶1萬地形圖生成DEM的坡度、曲率的頻率和累計頻率與 1∶5萬、1∶10萬和1∶25萬地形圖生成DEM的坡度、曲率的頻率和累計頻率的相似度,從而分別找到與1∶5萬、1∶10萬、1∶25萬比例尺DEM 的最為相似的1∶1萬比例尺插值DEM的柵格尺寸。

(4)坡度和曲率的直方圖逼真度計算。在A rcGIS軟件中,先使用重采樣工具對多柵格DEM重采樣生成目標(biāo)分辨率DEM的柵格大小,使用柵格計算器進(jìn)行柵格計算得出相應(yīng)參數(shù),帶入到公式(8)中得到逼真度值。

(5)分辨率匹配。根據(jù)統(tǒng)計參數(shù)、直方圖相似度、逼真度分別將多柵格DEM與多分辨率DEM進(jìn)行匹配,根據(jù)結(jié)果得出相關(guān)結(jié)論。

圖1 研究技術(shù)路線圖

2 結(jié)果與分析

2.1 坡度和曲率的統(tǒng)計特征參數(shù)

坡度和曲率統(tǒng)計特征參數(shù)見表1和表2。從表1可以看出 DEM10b的坡度統(tǒng)計特征與目標(biāo)數(shù)據(jù)DEM 10a上提取的坡度特征更為接近,DEM 25b的坡度統(tǒng)計特征與目標(biāo)數(shù)據(jù)DEM 25a上提取的坡度特征更為接近。從表2中可以看到DEM 10b的曲率統(tǒng)計特征與目標(biāo)數(shù)據(jù)DEM10a上提取的曲率特征更為接近,DEM25b的曲率統(tǒng)計特征與目標(biāo)數(shù)據(jù)DEM25a上提取的曲率特征更為接近。綜合坡度和曲率的統(tǒng)計特征,可以認(rèn)為DEM 10b、DEM25b的地形起伏反映能力分別相當(dāng)于DEM10a、DEM25a。從表1中找不到與DEM 50a的坡度統(tǒng)計特征相近的由1∶1萬地形圖生成的DEM,從表2中可以看到DEM60b與目標(biāo)數(shù)據(jù)DEM 50a的曲率統(tǒng)計特征較為相近。

表1 坡度的統(tǒng)計特征

表2 曲率的統(tǒng)計特征

2.2 坡度和曲率的相似度計算

坡度和曲率的頻率直方圖和累計頻率直方圖如下:以與DEM 10a的相似度比較為例(見圖2)。

對坡度和曲率的頻率和累計頻率分別使用直方圖L1距離、L2距離、χ2距離和直方圖交4種方法求取相關(guān)直方圖的相似度,與 DEM10a、DEM25a和DEM50a比較的相似度計算結(jié)果如表3-5所示(加黑數(shù)字為最大相似度值):

圖2 部分DEM坡度和曲率頻率分布

表3 與DEM 10a比較坡度、曲率直方圖相似度

表4 與DEM 25a比較坡度、曲率直方圖相似度

表5 與DEM 50a比較坡度、曲率直方圖相似度

從表3中可以看出DEM 10b與DEM 10a的坡度頻率、坡度累計頻率、曲率頻率和曲率累計頻率的所有相似度值都最高,因此可以認(rèn)為DEM10b與DEM 10a最為相似。從表4中可以看出 DEM25b與DEM25a的坡度頻率和坡度累計頻率的各相似度值最高,DEM 20b和DEM25b與DEM 25a的曲率頻率和曲率累計頻率的相似度值接近并最為相似,綜合認(rèn)為DEM25b與DEM 25a最為相似。從表5中可以看出隨著DEM柵格尺寸不斷變大,與DEM 50a的坡度頻率和累計頻率的相似度變大,是一個增函數(shù),但是相似度是不是一直增大,或者有個峰值,還需生成更大柵格尺寸的DEM 進(jìn)行驗證,再生成 DEM80b,DEM90b和 DEM100b進(jìn)行驗證,見表6所示,結(jié)果顯示,DEM 80b的頻率各相似度值最高,而DEM 90b的累計頻率相似度值最高。同時也可以從表5看出,DEM60b與DEM50a的曲率頻率、曲率累計頻率的相似度最高,綜合來看,坡度和曲率的相似度不統(tǒng)一。

2.3 坡度和曲率的逼真度計算

采用公式(8),計算出的坡度和曲率的逼真度如表7-8所示:

表6 與DEM 50a比較坡度、曲率直方圖相似度

表7 坡度的逼真度值

表8 曲率的逼真度值

從表7中可以看出DEM 10b與DEM 10a,DEM 20b與DEM 25a,DEM 60b與DEM50a的逼真度值最高,都接近于1。從表8中得出的結(jié)果是DEM 10b與DEM 10a,DEM 20b與DEM25a,DEM 40b與 DEM 50a的逼真度值較高,但是曲率逼真度值總體較坡度的逼真度值小。

3 結(jié)論與討論

3.1 結(jié)論

在變換前后DEM統(tǒng)計特征比較的基礎(chǔ)上,通過引入直方圖相似度和逼真度指數(shù)兩個指標(biāo),可以對分辨率匹配提供依據(jù)。

(1)初步分析表明,基于1∶1萬地形圖插值生成的 10 m、20 m或 25 m柵格的 DEM,與基于1∶5萬、1∶10萬地形圖生成的原生DEM(分辨率為10 m、25 m)相似。從統(tǒng)計特征參數(shù)、直方圖相似度指數(shù)和逼真度指數(shù)這三種方法得不到統(tǒng)一的與1∶25萬原生分辨率DEM相似的基于1∶1萬地形圖插值柵格大小的DEM,并且各方法得到的相似度值并不高。也可以看出隨著比例尺變小,原生DEM與1∶1萬地形圖插值生成的柵格DEM的相似度變小。

(2)關(guān)于相似度指數(shù),從本文的試驗結(jié)果看,四種直方圖相似度的計算方法對計算DEM統(tǒng)計特征的直方圖相似度同樣適用,取得基本相同的結(jié)果,但是L1距離、L2距離、χ2距離相比直方圖的交方法的優(yōu)勢在于可以對累計直方圖進(jìn)行計算。L1距離所反映的相似度值的幅度變化較大,所以更加直觀,因此L1距離是4種方法中最適合的相似度指數(shù)計算方法。

3.2 討論

(1)由于直方圖相似度和逼真度計算只是基于統(tǒng)計學(xué)原理,而DEM具有空間格局屬性,因此在直方圖相似度和逼真度計算基礎(chǔ)上,還必須充分考慮DEM的格局和地形信息含量等,才能形成全面準(zhǔn)確的認(rèn)識;

(2)隨著比例尺的變小,原生DEM 與1∶1萬地形圖插值生成的柵格DEM的相似度變小。可能的原因是,DEM地形分析的尺度推繹要考慮DEM格網(wǎng)單元異質(zhì)性和尺度推繹范圍,由于尺度推繹非線性的根本原因,尺度推繹不能無限度的進(jìn)行[16],隨著從高分辨率到低分辨率的轉(zhuǎn)換過程,尺度轉(zhuǎn)換的有效性在降低。

(3)利用1∶1萬地形圖插值,生成與1∶25萬地形圖原生DEM分辨率對應(yīng)的DEM的有效性還需進(jìn)一步探討。

[1] 楊勤科,Jupp D B L,郭偉玲,等.基于濾波方法的DEM尺度變換方法研究[J].水土保持通報,2008,28(6):58-62.

[2] 郝振純,池宸星,王玲,等,DEM空間分辨率的初步分析[J].地地球科學(xué)進(jìn)展,2005,20(5):499-504.

[3] Wo lck D M,M cCabe G J.Differences in topographic characteristics computed from 100-m and 1000-m resolution digital elevation model data[J].hyd ro logical p rocesses,2000,14(2):987-1002.

[4] 吳險峰,劉昌明,王中根.柵格DEM的水平分辨率對流域特征的影響分析[J].自然資源學(xué)報,2003,18(2):168-173.

[5] 吳凡,祝國瑞.基于小波分析的地貌多尺度表達(dá)與自動綜合[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2001,26(2):170-176.

[6] 于浩,楊勤科,張曉萍,等.基于小波多尺度分析的DEM 數(shù)據(jù)綜合研究[J].測繪科學(xué),2008,33(3):93-95.

[7] 張彩霞,楊勤科,段建軍.高分辨率數(shù)字高程模型的構(gòu)建方法[J].水利學(xué)報,2006,37(8):1009-1013.

[8] 師維娟,楊勤科,趙東波,等.中分辨率水文地貌關(guān)系正確DEM 建立方法研究:以黃土丘陵區(qū)為例[J].西北農(nóng)林科技大學(xué)學(xué)報:自然科學(xué)版,2007,35(2):143-148.

[9] Yang QK,McVicar T R,Van Niel TG,etal.Imp roving terrain rep resentation of a digital elevationmodel by reducing source data errors and optim ising interpolation algorithm parameters:an example in the Loess Plateau[J].China International Journal of App lied Earth Observation and Geoinformation,2007,9(3):235-246.

[10] Hutchinson M F.ANUDEM Version 5.1.A.N.U.Centre for Resource and Environmental Studies[M].Canberra,2004.

[11] 楊勤科,M cvicar T R,李領(lǐng)濤,等.ANUDEM-專業(yè)化數(shù)字高程模型插值算法及其特點[J].干旱地區(qū)農(nóng)業(yè)研究[J].2006,24(3):36-41.

[12] 朱偉,王東華,周曉光.基于信息熵的DEM最佳分辨率確定方法研究[J].遙感信息,2008(5):79-82.

[13] 馬煜,顧曉東,汪源源.基于平均窗口平移的直方圖相似性度量[J].應(yīng)用科學(xué)學(xué)報,2008,26(1):28-33.

[14] 林麗惠.一種改進(jìn)的顏色直方圖相似性度量算法[J].武夷學(xué)院學(xué)報,2009,28(2):58-61.

[15] 尹忠海,簡劍鋒,周利華,等.現(xiàn)有影像逼真度定義的缺陷及其改進(jìn)[J].西安電子科技大學(xué)學(xué)報:自然科學(xué)版[J].2004,31(6):833-836,869.

[16] 劉學(xué)軍,盧華興,仁政,等.論DEM 地形分析中的尺度問題[J].地理研究,2007,26(3):433-442.