DEM 分辨率對提取東北地區(qū)坡度的影響研究1

劉佳軒 高 瑞 師黎靜

1）中國地震局工程力學(xué)研究所，地震工程與工程振動重點(diǎn)實(shí)驗(yàn)室，哈爾濱 154200 2）黑龍江科技大學(xué)，建筑工程學(xué)院，哈爾濱 154200

引言

相比鉆孔數(shù)據(jù)受經(jīng)費(fèi)、環(huán)境及鉆孔技術(shù)等的限制，基于地形坡度的大區(qū)域場地分類方法得到大量應(yīng)用。目前國內(nèi)外研究多基于Wald 等（2006）建立的坡度?VS30模型進(jìn)行場地類別劃分（史大成，2009；Thompson 等，2014；吳效勇，2019），其中坡度來源主要是基于30″分辨率數(shù)字高程數(shù)據(jù)（Digital Elevation Model，DEM）的提取。然而，提取的坡度受DEM 分辨率、算法及數(shù)據(jù)類型（網(wǎng)格點(diǎn)類型）等因素影響，且分辨率是最主要的影響參數(shù)。國內(nèi)外大量研究表明（Walker 等，1999；Mukherjee 等，2013；Grohmann，2015；劉利峰，2018；師動等，2020），提取的坡度隨著DEM 分辨率的降低而減小，即坡度均值與標(biāo)準(zhǔn)差明顯降低，且不同地區(qū)分辨率的影響不同，這對建立坡度-剪切波速、坡度-覆蓋層厚度等關(guān)系模型具有重要影響，因此應(yīng)選擇合適分辨率的DEM 進(jìn)行數(shù)字地形分析（Thornley，1976）。李昕蕾（2020）通過對新疆地區(qū)不同分辨率DEM 坡度的提取，發(fā)現(xiàn)利用精度高或大區(qū)域小比例尺數(shù)據(jù)進(jìn)行場地分類的效果較差，而30″分辨率效果較好；喬之軒（2019）研究發(fā)現(xiàn)，高分辨率DEM 提取的坡度具有更精細(xì)的局部變化，但有時會掩蓋區(qū)域坡度的整體變化趨勢和特征。受計(jì)算機(jī)能力與硬件顯示等條件限制，高分辨率DEM往往不能得到良好應(yīng)用，有時需采用與研究比例尺相適應(yīng)的分辨率。部分研究者在提取不同分辨率下的坡度時考慮地貌條件，發(fā)現(xiàn)不同地貌下的坡度尺度具有明顯差異性（Wang 等，2012；劉曉等，2017；土祥等，2018）。目前，對我國東北地區(qū)考慮地貌條件下的坡度研究相對較少，且不同研究區(qū)域可能搜集到不同分辨率DEM，需對提取的坡度進(jìn)行統(tǒng)一，王英等（2019）、劉飛等（2019）以分形理論為基礎(chǔ)，建立不同分辨率下DEM 坡度轉(zhuǎn)換模型，分別以涇河區(qū)域和四川丘陵地區(qū)某流域?yàn)檠芯繉ο?，研究成果不適用于我國東北地區(qū)。

基于上述問題，本文以美國地質(zhì)勘探局（United States Geological Survey，USGS）發(fā)布的我國東北地區(qū)1″、3″、30″分辨率DEM 數(shù)據(jù)為基礎(chǔ)，分別提取研究區(qū)域內(nèi)公里網(wǎng)格點(diǎn)的百分比坡度，分析不同分辨率下坡度差異原因。區(qū)分平原、丘陵和山地后進(jìn)行坡度-頻率分布曲線統(tǒng)計(jì)，分析不同地貌單元下的坡度隨分辨率變化特征。通過線性函數(shù)模型、多項(xiàng)式函數(shù)模型及冪函數(shù)模型，擬合不同分辨率坡度轉(zhuǎn)換關(guān)系，并對比擬合優(yōu)度R2和坡度分級平均相對差值Δmean，選取不同DEM 分辨率得到坡度最佳擬合公式。

1 研究區(qū)域及數(shù)據(jù)

1.1 研究區(qū)域

研究區(qū)域主要包括黑龍江省、吉林省、遼寧省及內(nèi)蒙古自治區(qū)部分地區(qū)，地理坐標(biāo)為116°E～135°E，41°N～53°N。東西向最大橫距約1 441 km，南北向最大縱距約1 656 km，面積約1 450 000 km2。東北地區(qū)海拔高度為?275～2 692 m，整個地勢大致分為3 環(huán)，外圍由黑龍江、烏蘇里江、鴨綠江等流域構(gòu)成，整體地勢較低；中部為由大興安嶺、小興安嶺及長白山構(gòu)成的山地與丘陵，整體地勢較高；內(nèi)部為由松花江、嫩江等流域構(gòu)成的東北平原，整體地勢較低。

1.2 DEM 高程數(shù)據(jù)

研究數(shù)據(jù)來源于USGS 官網(wǎng)發(fā)布的DEM，原始數(shù)據(jù)采樣間隔分別為地球等角坐標(biāo)系的1″（約30 m）、3″（約90 m）和30″（約1 km），圖1 所示為3 種分辨率下高程數(shù)據(jù)模型，不同分辨率下DEM 表現(xiàn)出的整體地勢結(jié)構(gòu)基本相同，高分辨率DEM 對研究區(qū)域的描述更詳細(xì)，具有更精確的高程數(shù)據(jù)范圍及更精細(xì)的地形刻畫，隨著分辨率的降低，區(qū)域內(nèi)高程數(shù)據(jù)產(chǎn)生一定程度的過濾與平滑，使區(qū)域整體高程數(shù)據(jù)范圍減小，但整體地勢結(jié)構(gòu)更突出。

圖1 東北地區(qū)不同分辨率下高程數(shù)據(jù)模型圖Fig. 1 Elevation maps of 1 arc second, 3 arc second, and 30 arc second resolution in the Northeast, China

2 坡度提取

2.1 提取方法

目前最常用的坡度提取方法為擬合曲面法（Olaya，2009），計(jì)算式為：

式中，slope 表示坡度，本文使用廣泛應(yīng)用于大區(qū)域場地分類的百分比坡度；dz/dx，dz/dy分別表示計(jì)算點(diǎn)處2 個正交水平方向高程梯度變化率，無量綱，考慮復(fù)雜地貌坡度影響時，dz/dx，dz/dy優(yōu)選算法建議使用三階反距離權(quán)算法，即Horn 算法（Horn，1981）。

圖2 所示為計(jì)算高程梯度變化率dz/dx，dz/dy的3×3 像元窗口，Horn 算法中高程梯度變化率dz/dx，dz/dy計(jì)算如下：

式中，w為像元寬度，即分辨率。

以1″和30″分辨率DEM 坡度計(jì)算過程為例，說明不同分辨率下坡度提取差異。當(dāng)分辨率為1″時，式（2）和式（3）是對圖2 代表的約90 m×90 m 范圍內(nèi)高程數(shù)據(jù)點(diǎn)的計(jì)算，最終通過代表中間柵格數(shù)據(jù)的高程梯度變化率得到坡度。當(dāng)分辨率為30″時，式（2）和式（3）是對圖2 代表的約3 km×3 km 范圍內(nèi)高程數(shù)據(jù)點(diǎn)的計(jì)算，進(jìn)而計(jì)算得到坡度。雖然1″、30″分辨率下計(jì)算中心點(diǎn)可以相同，但坡度體現(xiàn)計(jì)算面積平均特性，隨著分辨率的增大，柵格計(jì)算面積變小，因此高程數(shù)據(jù)有所差異，導(dǎo)致不同分辨率下的坡度有時存在較大差距，且區(qū)域附近坡度變化越大，可能導(dǎo)致不同分辨率下坡度差異越大。

圖2 計(jì)算高程梯度變化率的像元窗口Fig. 2 3×3 unit schematic diagram of slope calculation point

2.2 提取結(jié)果

圖3 所示為通過式（1）及DEM 數(shù)據(jù)計(jì)算得到的不同分辨率下百分比坡度圖，由圖3 可知，高分辨率下的坡度圖對地表地形的描述更精細(xì)，坡度特征表現(xiàn)層次更豐富；隨著DEM 分辨率的降低，坡度值域逐漸減小，部分地區(qū)坡度趨于0，整體地形結(jié)構(gòu)更突出，但缺少地形細(xì)節(jié)的表達(dá)。通過GIS 軟件建立與不同分辨率下坡度圖相同的漁網(wǎng)式公里網(wǎng)格點(diǎn)，并提取網(wǎng)格點(diǎn)對應(yīng)的坡度。

圖3 東北地區(qū)不同分辨率下坡度圖Fig. 3 Slope maps of 1 arc second, 3 arc second, and 30 arc second resolution in the Northeast, China

圖4 所示為東北地區(qū)不同DEM 分辨率下坡度統(tǒng)計(jì)曲線，由圖4（a）可知，高分辨率下坡度累計(jì)頻率曲線較低分辨率更平緩，即高分辨率下坡度整體相對偏高；對于坡度值域而言，高分辨率下坡度值域范圍較低分辨率廣，1″、3″、30″分辨率下坡度分別約為0.6、0.4、0.2 m/m。由圖4（b）可知，隨著分辨率的降低，坡度均值和標(biāo)準(zhǔn)差明顯減小，且在高分辨率下減小速率較快，在低分辨率下減小速率較慢，標(biāo)準(zhǔn)差的降低說明了坡度值域整體向均值附近靠近。圖4 基本表現(xiàn)出東北地區(qū)整體坡度隨著DEM 分辨率的降低而減小的現(xiàn)象。

圖4 不同DEM 分辨率下的坡度統(tǒng)計(jì)Fig. 4 Slope statistics of different DEM resolutions

3 不同地貌坡度統(tǒng)計(jì)分析

區(qū)分不同地貌條件，對研究區(qū)域網(wǎng)格采樣點(diǎn)坡度進(jìn)行統(tǒng)計(jì)分析，研究不同地貌下坡度隨分辨率變化特征，為擬合公式按照不同地貌條件提供理論基礎(chǔ)。

圖5 所示為我國東北地區(qū)地貌分布，平原、丘陵和山地地貌占比分別約為33.3%、12.2%、43.1%，其他地貌占比約為11.4%。通常，平原地貌較平坦或起伏較??；丘陵地貌由多個低矮山丘組成，具有一定起伏；山地地貌由眾多山脈區(qū)域組成，地面起伏較大。對于海拔高度而言，通常山地海拔較高，丘陵次之，平原較低。對于坡度而言，通常山地坡度較大、較陡，丘陵次之，平原坡度較小、較緩。

圖5 地貌分布圖Fig. 5 Landform distribution map

圖6 所示為不同分辨率下不同地貌單元坡度-頻率分布曲線。由圖6 可知，平原在坡度較小時頻率分布較大，山地在坡度較大時頻率分布相對較大；隨著坡度的增大，平原坡度-頻率分布曲線下降速率較山地快，而丘陵坡度-頻率分布曲線變化特征處于兩者中間；隨著分辨率的減小，平原坡度-頻率分布曲線峰值位置迅速向低坡度段移動，峰值區(qū)域變窄，坡度減小至0.05 m/m 以下，而山地坡度-頻率分布曲線峰值位置向低坡度移動速度相對較緩，峰值區(qū)域變窄程度相對較小，坡度減小至0.15 m/m 以下，丘陵坡度-頻率分布曲線變化特征處于平原與山地之間，坡度減小至0.1 m/m 以下。各地貌單元坡度減小程度不同，因此需考慮地貌條件，建立不同DEM 分辨率的坡度轉(zhuǎn)換關(guān)系。

圖6 不同分辨率下不同地貌單元坡度-頻率分布曲線Fig. 6 Slope-frequency distribution curves of different topography at different resolutions

4 轉(zhuǎn)換關(guān)系擬合回歸

4.1 擬合模型

本文采用常用函數(shù)模型對不同DEM 分辨率提取的坡度進(jìn)行擬合，各函數(shù)模型如下：

線性函數(shù)模型：

式中，Si表示分辨率為i（其中i=1″、3″、30″）時計(jì)算得到的百分比坡度（m/m）；aj、bj（j=1、2、3）、c1、c2、d1分別表示不同函數(shù)模型的回歸參數(shù)。為使擬合公式計(jì)算結(jié)果更合理，需考慮邊界效應(yīng)，即當(dāng)1″、3″分辨率下坡度趨于0 時，30″下坡度也應(yīng)趨于0。因此，當(dāng)Si=0 時，S30=0，即a1、d1、c2=0。

4.2 擬合回歸

基于上述函數(shù)模型，分別對不同地貌單元1″、3″分辨率下坡度與30″分辨率下坡度進(jìn)行擬合，坡度轉(zhuǎn)換關(guān)系參數(shù)如表1、2 所示。由表1、2 可知，平原回歸關(guān)系參數(shù)b最小，而山地回歸關(guān)系參數(shù)b最大，丘陵回歸關(guān)系參數(shù)b處于二者之間，基本符合隨著分辨率減小，平原坡度減小程度較大，而山地坡度減小程度相對較小的變化規(guī)律；多項(xiàng)式函數(shù)及冪函數(shù)模型擬合參數(shù)具有類似規(guī)律。

表1 不同地貌單元回歸模型擬合參數(shù)及評價(jià)參數(shù)（1″與30″分辨率下）Table 1 Fitting and evaluation parameters of regression models for different landforms (1 arc second and 30 arc second)

5 回歸方程評價(jià)方法

5.1 擬合優(yōu)度

表2 不同地貌單元回歸模型擬合參數(shù)及評價(jià)參數(shù)（3″與30″分辨率下）Table 2 Fitting and evaluation parameters of regression models for different landforms (3 arc second and 30 arc second)

擬合優(yōu)度R2計(jì)算如下：式中，SSR為回歸平方和，SST為總平方和，ESS為誤差平方和，TSS為總離差平方和。

由表1、2 可知，冪函數(shù)模型轉(zhuǎn)換的擬合優(yōu)度R2最小，即冪函數(shù)模型擬合效果最差；線性函數(shù)模型擬合優(yōu)度R2略小于多項(xiàng)式函數(shù)模型；整體來說，多項(xiàng)式函數(shù)模型擬合結(jié)果最優(yōu)。

5.2 坡度分級平均相對差值

高分辨率DEM 可獲得高分辨率坡度，進(jìn)而獲得高分辨率等效剪切波速或覆蓋層厚度，但場地分類主要使用30″分辨率DEM 提取的坡度，分類方法基本為建立坡度與等效剪切波速的分級對應(yīng)關(guān)系，因此還需考慮轉(zhuǎn)換前后30″坡度分級結(jié)果。將坡度按0.01 m/m 寬度范圍進(jìn)行分級統(tǒng)計(jì)，以30″分辨率DEM 提取的坡度分級作為參考值，分別計(jì)算不同擬合模型對高分辨率坡度轉(zhuǎn)換后的坡度分級與參考值的平均相對差值Δmean，平均相對差值越接近于0，說明該擬合公式向30″坡度轉(zhuǎn)化越好，平均相對差值越接近于1，說明該擬合公式向30″坡度轉(zhuǎn)化越差，平均相對差值計(jì)算如下：

式中，Pi,k表示在i（i=1″、3″）分辨率下轉(zhuǎn)換后第k組分段坡度，n表示坡度分級級數(shù)。

不同地貌單元不同分辨率下坡度轉(zhuǎn)換結(jié)果如圖7 所示，由圖7 可知，平原1″、3″與30″真值接近，基本將分布在0～0.2 m/m 的坡度轉(zhuǎn)換為0～0.03 m/m；丘陵1″向30″轉(zhuǎn)換時，線性函數(shù)模型轉(zhuǎn)換與多項(xiàng)式函數(shù)模型轉(zhuǎn)換結(jié)果較好，而3″向30″轉(zhuǎn)換時多項(xiàng)式函數(shù)模型轉(zhuǎn)換結(jié)果較好；山地線性函數(shù)模型轉(zhuǎn)換效果較好，而多項(xiàng)式函數(shù)模型與冪函數(shù)模型轉(zhuǎn)換效果較差。由表1、2 可知，線性函數(shù)模型分級平均相對差值Δmean最小，即僅考慮坡度分級時，線性函數(shù)擬合為最佳擬合。

圖7 不同地貌單元不同分辨率下坡度轉(zhuǎn)換結(jié)果Fig. 7 Conversion results of DEM slopes with different resolutions under different landforms

5.3 最佳擬合模型

由于分別考慮擬合優(yōu)度R2和分級平均相對差值Δmean，各模型在坡度轉(zhuǎn)換上均有優(yōu)劣，因此，以比例參數(shù)K綜合考慮R2、Δmeans，計(jì)算公式如下：

K值越大，說明擬合模型擬合效果越好。

根據(jù)表1、2 中的比例參數(shù)K，最終確定不同地貌單元不同分辨率下坡度轉(zhuǎn)換公式，如表3 所示。平原和山地最佳坡度擬合模型均為線性函數(shù)模型；丘陵1″與30″分辨率最佳坡度擬合模型為線性函數(shù)模型，3″與30″分辨率最佳坡度擬合模型為多項(xiàng)式函數(shù)模型。

表3 不同地貌單元不同分辨率下坡度轉(zhuǎn)換公式Table 3 The slope conversion formula of different resolution DEM for different landforms

6 結(jié)論

本文以我國東北地區(qū)1″、3″、30″分辨率DEM 數(shù)據(jù)為基礎(chǔ)，分別提取研究區(qū)域公里網(wǎng)格點(diǎn)百分比坡度，分析不同分辨率下坡度提取差異原因。區(qū)分平原、丘陵和山地進(jìn)行坡度-頻率分布曲線統(tǒng)計(jì)，分析不同地貌單元坡度變化特征。通過線性函數(shù)模型、多項(xiàng)式函數(shù)模型及冪函數(shù)模型擬合不同分辨率坡度轉(zhuǎn)換關(guān)系，并計(jì)算擬合公式擬合優(yōu)度R2及坡度分級平均相對差值Δmean，得出以下結(jié)論：

（1）不同分辨率提取的坡度不同，主要原因是坡度計(jì)算方法體現(xiàn)計(jì)算面積平均特性，隨著分辨率的增大，柵格計(jì)算面積變小，高程數(shù)據(jù)有所差異，導(dǎo)致不同分辨率下坡度存在較大差異。由于坡度計(jì)算平均特性，區(qū)域附近坡度變化越大，可能導(dǎo)致不同分辨率下獲得的坡度差異越大。

（2）我國東北地區(qū)坡度隨DEM 分辨率的降低而減小，高分辨率下坡度圖對地表地形的描述更精細(xì)，坡度特征表現(xiàn)層次更豐富，而低分辨率下整體地形結(jié)構(gòu)更突出，但缺少地形細(xì)節(jié)的表達(dá)。隨著分辨率的降低，我國東北地區(qū)坡度均值、方差及變化范圍減小，坡度-頻率分布曲線向低坡度段移動。

（3）隨著分辨率的降低，不同地貌單元坡度-頻率分布曲線變化趨勢相同，變化程度不同。各地貌單元坡度-頻率分布曲線均向低坡度方向移動，且峰值區(qū)域變窄，其中平原地貌變化趨勢性顯著，丘陵次之，山地最小。因此，建立坡度轉(zhuǎn)換關(guān)系時需區(qū)分不同地貌條件。

（4）給出不同地貌單元下東北地區(qū)不同分辨率坡度轉(zhuǎn)換模型，為利用30″分辨率的坡度進(jìn)行大區(qū)域場地分類提供轉(zhuǎn)換關(guān)系。平原和山地最佳坡度擬合模型均為線性函數(shù)模型；丘陵1″與30″分辨率最佳坡度擬合模型為線性函數(shù)模型，3″與30″分辨率最佳坡度擬合模型為多項(xiàng)式函數(shù)模型。

（5）本文給出的最佳坡度轉(zhuǎn)換公式主要應(yīng)用于工程抗震領(lǐng)域中的場地分類。