基于DEM的京津冀地区地形起伏度分析

张 竞，杜 东，白耀楠，纪冬丽

(1.中国地质调查局 天津地质调查中心，天津 300170； 2.天津城建大学 市政与环境工程学院，天津 300384)

起伏度最佳分析窗口的大小因地形而异：国际地理学会地貌制图委员会在编制《欧洲1∶250万地貌图》时使用了16 km2的分析窗口[6]；陈志朋等[7]在《1∶400万中国及其邻区地貌图》的编制中使用了21 km2的分析窗口；王玲等[8]基于1∶25万DEM计算的新疆山地最佳分析窗口为2.56 km2；张伟等[9]基于SRTM DEM，认为黄土高原最佳分析窗口为3.88 km2，云贵山区为4.58 km2，青藏高原东缘为5.34 km2；胡最等[10]基于1∶25万DEM确定的湖南省最佳分析窗口为0.28 km2；王让虎等[11]基于ASTER GDEM提取的我国东北地区最佳分析窗口为2.62 km2；丁贤法[12]基于SRTM DEM得到的云南省富宁县最佳分析窗口为0.52 km2；郎玲玲等[13]基于1∶25万DEM得到福建低山丘陵区的最佳分析窗口为4.41 km2。

既然最佳分析窗口具有明显的空间差异性，当研究区幅员辽阔、地貌类型复杂时，就可能存在不止一个最佳分析窗口，如涂汉明等[4]指出中国存在5种不同规模的地形起伏度最佳分析窗口——2、10、16、20、22km2，适用于不同的地形地貌。但是，目前的大量研究都采用全局单值最佳分析窗口，其中不乏面积广大、地形复杂的地区[3,6,10-11,14-15]，此时的最佳分析窗口仅能够合理地表达一部分地形，对其他地形则失去了“最佳”属性。本研究选取地貌类型复杂的京津冀地区，针对不同地貌类型提取地形起伏度最佳分析窗口，分析地形对最佳分析窗口的影响，提出京津冀地区地形起伏度的合理表达。

1 研究区概况

京津冀地区位于华北平原北部，战略地位十分重要。研究区地理范围介于113°27′31″～119°51′33″E、36°02′47″～42°36′53″N，总面积21.54万km2，地质构造条件复杂，发育有平原、丘陵、盆地、山区、高原等多种地貌类型，北靠燕山山脉，南部平原展布，西倚太行山，东临渤海湾，西北部为坝上高原，丘陵主要分布于太行山东侧和燕山南侧，沿渤海海岸多滩涂、湿地，总体呈现出西北高东南低的地形特点。区内海拔2 000 m以上的山峰有10座，第一高峰位于张家口小五台山东台，海拔2 882 m[3]。

2 研究数据与研究方法

2.1 研究数据

本研究DEM数据使用90 m分辨率SRTM3 DEM，绝对高程精度为±16 m，绝对平面精度为±20 m[16]，数据获取自地理空间数据云网站。京津冀三地行政区划矢量数据获取自国家基础地理信息中心，比例尺为1∶25万。所有数据统一采用我国CGCS2000_GK_CM_117E投影，在该投影下DEM实际分辨率为84.8 m，使用双线性插值法重采样为90 m。利用矢量数据对DEM数据进行裁剪后得到京津冀地区DEM数据(图1)。

图1 京津冀地区DEM及研究样区

2.2 起伏度和最佳分析窗口提取原理

地形起伏度一般基于DEM数据采用窗口分析法求取，其表达式为

Ra=Hmax-Hmin

(1)

式中：Ra为以第a个栅格为中心的一定面积窗口内的地形起伏度；Hmax和Hmin分别为该窗口内的最大高程和最小高程。

起伏度概念的关键是最佳分析窗口，在该窗口下提取的地形起伏度可以满足山体完整性和区域普适性原则。山体完整性原则是指分析窗口起初只包含一部分山体，随着面积增加，窗口内高差迅速增大，当窗口恰好覆盖整个山体时，高差达到一较大值，之后随窗口面积增大高差增大的速度明显减小，那么这个分析窗口之内的高差就恰到好处地反映了该山体的地形起伏特征。区域普适性原则是指最佳分析窗口需要满足研究区域内高差最大山体的完整性，此时的分析窗口必然满足较小山体的完整性原则[1-2]。以往学者在提取最佳分析窗口时使用了不同的方法，如人工作图法[6,13,17]、最大高差法[13]、最大高差-面积比法[3,18]、模糊数学法[1]、均值变点法[8-9,11,14,19-26]等。赵斌滨等[27]对各种方法的精度进行了验证，认为人工作图法和均值变点法得到的最佳分析窗口对应的起伏度更符合实际情况，最大高差法和模糊数学法已不适用。由于人工作图法的主观性较强，近年来均值变点法成为最佳分析窗口计算中的主流方法，这种方法对恰有一个变点的检验最为有效[28]。本研究采用均值变点法计算地形起伏度的最佳分析窗口，计算原理如下。

(1)依次计算N个递增分析窗口下单位面积上的平均起伏度，即单位地势度T。

1.2.1 对照组 对患者每个月开展1次肺健康知识宣教，内容包括COPD相关知识、肺康复理念，指导戒烟，给予氧疗、药物治疗，共治疗20周。

Ti=ti/si(i=1，2，…，N)

(2)

式中：Ti为第i个分析窗口下的单位地势度；ti为平均起伏度；si为分析窗口面积。

(2)对数列T取对数lnT得到非线性数列样本Xk(k=1，2，…，N)。

(4)计算统计量，公式为

(3)

(4)

式中：S为总样本的离差平方和；Sk为前后两段样本的离差平方和之差。

变点的存在会使S和Sk的差距增大，S-Sk的最大值对应的分析窗口大小即为最佳分析窗口。

2.3 试验方案

根据DEM特征，选取7个包含各类地貌的样本:样本①～③为矩形样本，分别为山区小区、山区平原混合小区、平原小区；样本④～⑥为京津冀地区依地形划分的3个地貌大区，分别为高原大区、山区大区和平原大区；样本⑦为京津冀全区(图1)。对每个样本进行分析窗口递增的地形起伏度计算，窗口使用圆形邻域。以往研究中最佳分析窗口一般不超过25 km2，本研究将最大窗口控制在35 km2以内，分析窗口设置见表1。起伏度提取过程在ArcGIS中进行，利用Modelbuilder将Focal statistics工具(用于统计最大值、最小值)、Raster Calculator工具(用于计算高差)和Band Collection Statistics工具(用于将起伏度的统计数据输出)串联为工作流实现起伏度快速提取，在Excel中利用VBA编程实现均值变点法中S-Sk值的计算。

表1 地形起伏度分析窗口的设置

注：圆形窗口计算面积时半径栅格数需加0.5。

3 结果与分析

3.1 起伏度与分析窗口面积的统计特征

以每个分析窗口下计算的平均起伏度作为因变量，递增的分析窗口作为自变量绘制散点图(图2)。

图2 平均起伏度-分析窗口面积散点图

如图2所示，平均起伏度随分析窗口面积的增大而增加，变化曲线呈对数函数或幂函数特征。基于这两种函数进行回归分析得到拟合方程(表2)。由表2知，对数函数和幂函数均有较好的拟合效果，R2均大于0.96。相比之下，使用幂函数的拟合效果优于对数函数，且受地貌的影响较小，例如当使用对数函数拟合时，不同样本的R2在0.964 8～0.985 6之间变化，而使用幂函数拟合时，R2值更高且更集中，变化范围为0.983 8～0.996 9。当分析窗口面积一定时，不同样本的平均起伏度存在差异，这种差异在分析窗口面积较小时并不明显，随着分析窗口面积递增，差异逐渐放大，例如分析窗口从0.16 km2增加到33.88 km2时，7个样本平均起伏度的最大差值由85.54 m增加到578.43 m。

表2 平均起伏度与分析窗口面积的对数和幂函数拟合方程

3.2 最佳分析窗口

图2中的每条曲线存在一个增大速度由快变慢的点(非数学拐点)，该点对应的窗口即为最佳分析窗口。根据公式(3)和(4)用均值变点法计算得到各样本的系列S-Sk值，其与分析窗口面积的对应关系见图3。图3显示，各样本S-Sk值均呈单波峰曲线形态，峰值出现的位置较接近，说明各样本的最佳分析窗口面积差别不大。但是，统计分析显示7个样本仍表现出分区性：样本③和⑥属一个区，其他5个样本属另一个区，两区在S-Sk的均值、标准差、偏度、峰度等多个指标上均有明显差异，这种分区性同样体现在最佳分析窗口面积上，即样本③和⑥的最佳分析窗口面积为4.64km2，其他5个样本为5.35km2(表3)。

图3 S-Sk与分析窗口面积散点图

样本均值标准差偏度峰度最佳分析窗口面积(km2)①山区小区16.476.47-0.59-0.945.35②山区平原混合小区15.426.08-0.58-0.955.35④高原小区14.425.65-0.59-0.945.35⑤山区大区15.796.18-0.60-0.925.35⑦京津冀全区15.576.10-0.59-0.935.35③平原小区21.238.46-0.56-0.984.64⑥平原大区18.517.41-0.54-1.004.64

3.3 地形地貌对最佳分析窗口的影响

显然，地形差异是出现2种最佳分析窗口的原因。张伟等[9]曾基于SRTMDEM提取我国平原、高原、丘陵、山地4种地貌的最佳分析窗口，结果分别为4.83、4.97、4.79、4.83km2，即山地的分析窗口小于等于高原和平原。此结论与常理恰好相反，也就是说，最佳分析窗口的大小并不能直观地衡量一个地区的宏观地貌特征。事实上，根据前人对最佳分析窗口需满足的山体完整性和区域普适性原则[1]的阐述，最佳分析窗口的大小取决于一个区域内的最大山体，而与小起伏度地区无关，即便后者占据更大面积。因此，张伟等的研究中平原和高原样本很可能存在局部大起伏山体。

为了量化最佳分析窗口对起伏度的表达，以每个样本在其最佳分析窗口下的最大起伏度Rm作为该分析窗口能够表达的最大高差。7个样本的Rm值见表4，其中：样本③平原小区最小，为69m；样本⑥虽同为平原区，但检查DEM数据发现该样本在北京顺义区大孙各庄镇北侧发育小面积孤山，因而其Rm值远大于样本③，达到392m；高原样本④Rm值为445m；样本①、②、⑤、⑦均覆盖一定山区，因而其Rm为高值；样本⑦京津冀全区包含了区内所有的山体，其Rm值必然是最大的。

表4 各样本的Rm值

由于DEM和分析窗口增大的步长均为分辨率较粗的离散数据，最佳分析窗口面积并没有随Rm单调增大，而是呈阶梯状增大(图4)，这种阶梯特征曾在涂汉明等[1]的研究中被证实。考虑到7个样本中，样本③平原小区起伏度非常小，不可能存在Rm值大于样本⑦的地貌区，因此在SRTMDEM数据下京津冀地区仅存在4.64和5.35km2两个最佳分析窗口(这2个窗口半径栅格数为连续值13和14)。此外，样本⑥和样本④位于阶梯抬升处且两者Rm值较接近，因此可以大致将Rm值400m作为选择2种最佳分析窗口的临界值，即4.64km2窗口可以表达400m以内的高差，5.35km2窗口更适合于表达400m以上的高差。

图4 最佳分析窗口面积与Rm的关系

3.4 京津冀地区的地形起伏度

在上文分析的基础上，采用4.64和5.35km2双窗口方案提取京津冀地区的地形起伏度。首先将4.64km2窗口下提取的起伏度中高差小于400m的部分保留，而将高差在400m以上的区域采用5.35km2分析窗口重新提取起伏度，将两部分合并后，依据《中国1∶1 000 000地貌图制图规范》得到京津冀地区地形起伏度分级图，并与5.35km2单值窗口提取方案进行对比(图5)。图5显示，双窗口方案和单窗口方案整体结果较接近，但从细节看，单窗口方案会将平原区起伏度夸大，例如在河北燕郊镇周边的平坦地区，单窗口方案识别出7个起伏度大于30m的微起伏栅格，在现有的多个地貌形态划分方案下，这些栅格均被划分为丘陵地貌[4]，这是不符合实际情况的，相比之下，双窗口方案仅识别出2个起伏度大于30m的栅格(图6)。在整个研究区内，两种方案得到的平坦区域(起伏度小于30m)面积总计相差1 126km2。可见，双窗口方案兼顾了京津冀地区半山半平原的地势特点，因而更具有优势。

(a)使用单窗口提取 (b)使用双窗口提取

(a)使用单窗口提取 (b)使用双窗口提取

京津冀地区的地形起伏度在0～1 145m之间，以平坦区、中起伏区和小起伏区为主，其中平坦区占比为43.99%，中起伏区和小起伏区占比分别为30.46%和14.12%，大起伏区占比为5.02%，极大起伏区极少(表5)。宏观上各种地貌大致以东北—西南向展布，而在西北—东南向则表现出明显的分带性，即中间为中、大起伏山区，东南部为平坦的华北平原，西北部为微、小起伏和平坦的坝上高原。

4 结 论

表5 京津冀地区各级地形起伏度面积

京津冀地区幅员辽阔、 地貌类型复杂， 全局单值最佳分析窗口在该地区地形起伏度的提取中具有局限性。本研究计算了研究区内不同地貌样本的最佳分析窗口，分析了地形地貌对最佳分析窗口的影响，在此基础上提取了京津冀地区的地形起伏度。

京津冀地区地形起伏度随分析窗口面积增大而增大，对两者对应关系的拟合效果幂函数优于对数函数，前者各样本R2均大于0.98。区内最佳分析窗口随地形高差的增大而呈阶梯状增大，存在4.64和5.35km2两个最佳分析窗口，前者可以表达400m以内的高差，后者更适合于表达400m以上的高差。利用4.64和5.35km2双窗口方案提取的地形起伏度优于单窗口方案，两种方案在中、大起伏区结果一致，但在平坦地区前者改善了后者对起伏度的夸大，更好地兼顾了研究区半山半平原的地形特征。分析结果显示，京津冀地区的地形起伏度在0～1 145m之间，以平坦、中起伏和小起伏为主。