基于多源数据的山前水系提取实验:以济南市及邻区为例

王泽利, 许洪泰, 肖泽后, 倪永进

(1.山东科技大学地球科学与工程学院, 青岛 266400; 2.山东省地震工程研究院, 济南 250000)

作为构造地貌的重要标志,河流已成为新构造运动研究中的一种主要研究对象[1-4]。河流的流量、弯曲度、分级、流域面积、纵剖面形态等特征对区域的构造隆升、断层活动等控制因素具有敏感的反应[5-6]。因此,快速、准确地提取河流信息,并通过对提取获取的水系形态特征与断裂构造之间的关系进行分析,来讨论其对构造运动的响应,进而为研究区域构造稳定性提供地貌上的证据至关重要。

随着空间数据采集方法的进步,基于数字高程模型(digital elevation model,DEM)的河流信息提取已在造山带地区受到广泛的研究与应用。但在沉积平原区,由于地貌平缓和人为影响干预较大的原因,该类方法应用受到极大限制。文献[7-16]分别借助雷达影像和光学卫星影像开展了对水系的提取的研究工作,形成了如阈值法、滤波法、分类器法等提取方法,但以上多关注于小面积范围的水系信息提取。值得指出的是,上述方法在各自针对的单一地质地貌区内都可以达成可接受的方法,但针对包含山地和平原的多相复杂地貌,均难以独立实现良好的水系提取效果。

济南市及邻区位于鲁西隆起区与济阳凹陷的边界部位,水系和断裂构造均较发育。拟针对该区域,选用DEM数据和遥感光学影像数据,尝试进行多源数据河流提取实验,以期探索山前平原及山间地带河流自动提取的方法。并通过分析提取出的水系在断裂附近的形态特征,探讨河流对构造活动的响应。

1 研究区概况

如图1所示,选取山东省中部的济南市及邻区作为研究区,范围为36°N～37.5°N,116.5°E～118°E。济南及邻区地处鲁中南低山丘陵与鲁西北冲积平原交界地带,构成了山前地带河流提取实验的良好实验地区。区内南部为泰山山脉,徂徕山山脉及两者所围限的莱芜盆地,有牟汶河、赢汶河等大汶河主要支流流经。北部多为平原地区,区内分布北东向的黄河,徒骇河等河流。如图2、图3所示,DEM高程图与坡度图表现出相应特征,区内海拔自北西向南东递增,最低点仅3 m,最高点为泰山山脉主峰玉皇顶,海拔1 545 m,与山前平原相对高差1 500 m以上,高程图(图2)显示的1 517 m应是由于DEM数据分辨率所限而导致的误差;坡度增加趋势与海拔相同,研究区北部基本上都是较低的区域,高坡度地区大都集中在研究区的南部。

图2 高程图Fig.2 Elevation map

图3 坡度图Fig.3 Slope map

2 数据来源与研究方法

2.1 数据来源

结合数字高程数据与遥感影像数据,尝试提取山前地带的河流水系。数字高程模型最早于20世纪50年代提出,在地质灾害,土地利用,景观格局等方面具有重大作用[17]。本次实验所用DEM数据于地理空间数据云中下载,名称为SRTMDEM(shuttle radar topography mission digital elevation model) 90 m分辨率原始高程数据,中心经度为117.5°,中心维度为37.5°,条带号60,行编号5,数据标识为srtm_60_05。

Landsat8卫星是美国陆地卫星计划(Landsat Program)的第8颗卫星,由美国国家航空航天局(National Aeronautics and Space Administration,NASA)与美国地质调查局(United States Geological Survey,USGS)合作开发、轨道科学公司建造。该卫星于2013年2月11号在加利福尼亚范登堡空军基地由Atlas-V火箭搭载发射成功,携带陆地成像仪(operational land imager,OLI)和热红外传感器(thermal infrared sensor,TIRS)。OLI陆地成像仪包括9个波段,空间分辨率为30 m,其中包括一个15 m的全色波段,成像宽幅为185 km×185 km;热红外传感器TIRS包括2个单独的热红外波段,分辨率100 m[18]。本次实验所用的Landsat8亦是于地理空间数据云中下载,云量低于5%,成像时间为2019年9月27日,条带号均为122,行编号分别为34、35,数据标识分别为LC81220342019270LGN00、LC81220352019270LGN00。

2.2 提取方法

理想情况下,地表水流向与地形坡向一致,并倾向于沿局部最低海拔和最大坡度形成地表径流,基于这一认识,发展出利用DEM提取线状河流的方法。在ArcGIS10.2平台的水文分析模块中,对下载的DEM数据进行镶嵌、填洼等预处理,然后计算流向、流量,确定流量阈值分析河网,对河流进行分级和矢量化工作。

流向的确定是采用的D8算法,8个数字分别代表8个方向,以中心栅格为参照,计算与周围栅格的距离权落差,将结果最大的设定为流出位置,生成流向栅格数据。

流量阈值的确定工作是影响河网提取结果的关键因素,本实验采用河网密度法[19-22]来确定流量的具体阈值。河网密度指单位面积流域上的河流总长度,其计算公式为

(1)

式(1)中:D为河网密度;L为对应阈值下的河流总长度;S为研究区域面积。

利用ArcGIS生成不同阈值下的河网,统计河流长度并计算出河网密度;再拟合函数关系,绘制河网密度与阈值的图像;最后根据拐点位置找出河网密度随阈值变化趋于平缓的坐标,以该点的阈值进行河网提取可得到较好的准确度。依据此方法结果阈值提取河网并将其分级、矢量化后即可形成依据DEM数据提取的最终河流数据。

不同光学性质的地物对不同波段光线的反射率存在差异,构成了基于多波段遥感影像水体提取的理论基础,当前主要的方法包括归一化差异水体指数[9](normalized difference water index,NDWI)、改进的归一化差异水体指数[10](modified normalized difference water index,MNDWI)、Gauss归一化差异水体指数[23](Gauss normalized difference water index,GNDWI)、增强水体指数[24](enhanced water index,EWI)等,各自的特征和适用性具体如下。

(1)NDWI法。水体反射率从绿光到近红外波段逐渐递减,在近红外波段最弱;而植被的反射率则相反,在近红外波段反射率最强[9]。因此,本方法采用绿光和近红外波段构建模型,突出显示水体的同时抑制了植被等地物,特点是能够抑制植被信息,突出水体,适用于地势开阔平坦地区。

(2)MNDWI法。该方法以短波红外波段(short-wave infrared,SWIR)替换近红外波段(near infrared,NIR),显著增大了水体与建筑物的反差。该方法特点是能够较好地去除居民地和土壤的影响,有利于提取城市水体。

(3)GNDWI法。针对线状河流水体的精确提取,可较好保留河流完整性,受阴影影响较大,需要用到高程信息,实现过程比较麻烦。

(4)EWI法。能够抑制居民地、土壤和植被噪声,易受到阴影及浅滩的影响,适合半干旱地区水体提取。

考虑研究区面积较大,覆盖了山野、村落和中大型城市,既有陡峭山地,又包括开阔平原,选取NDWI法和MNDWI法。

NDWI的数学模型可表示为

(2)

MNDWI的数学模型可表示为

(3)

式中:Green为绿光波段反射率,对应Landsat8数据的第3波段;NIR为近红外波段反射率,对应Landsat8数据的第5波段;SWIR为短波红外波段反射率,对应landsat8数据的第6波段。

遥感影像处理在ENVI5.3平台上进行,在完成辐射定标、大气校正、图像镶嵌等预处理后,依据式(2)、式(3)分别计算获得NDWI和MNDWI图像。

3 提取结果与分析

3.1 DEM提取结果

经过预处理、填洼及流向分析后,分别选取100、200、400、600、…、3 600共19个集水流量阈值对研究区进行水文信息提取,计算出河网密度,如表1所示。绘制出阈值与河网密度的散点图,利用Origin对其进行曲线拟合,结果如图4所示,拟合后确定幂函数R2=0.999 7,此数值越接近一说明拟合效果越好,可以看出,河网密度开始随着阈值的增大而急剧减小,后下降趋势趋于平缓,对拟合曲线求二阶导数可知,阈值在约1 000时,其数值几乎接近于零,并且随着阈值增大其变化趋势也不明显,因此此法求得的适宜阈值为1 000,并以此结果作为流量阈值进行水文信息提取,然后完成河流分级及河网矢量化工作,即可获得区内基于DEM提取出的矢量化水系信息图,如图5所示。

表1 河网密度计算结果

图4 拟合曲线Fig.4 Fitting curves

图5 DEM提取图Fig.5 DEM extraction diagram

3.2 遥感数据提取结果

利用ENVI对遥感影像数据完成预处理后,按照上述比值模型公式,分别计算出NDWI、MNDWI的结果,如图6、图7所示,去除异常值后两图像的像元亮度值(digital number,DN)均在-1～1范围内,经多次试验调整后,确定NDWI阈值为0,MNDWI阈值为0.4,DN值大于阈值的像素判定为水体。

（2）经产母猪情期受胎率：采用人工授精为86.93%，采用自然交配为78.69%，差异显著（P<0.05）。

图7 MNDWI提取图Fig.7 MNDWI extraction diagram

从NDWI与MNDWI两者的结果来看,研究区內黄河的提取效果两者几乎相同,但对于徒骇河这种次级水系MNDWI提取结果的连续程度、完整程度要好于NDWI,此外对济南市区和山区的噪音处理方面MNDWI表现较为突出,结合上述两者的表现情况,认为MNDWI相较于NDWI体现出了更加丰富的水文信息,即MNDWI更加适用于研究区,所以将会用MNDWI方法的提取结果来参与接下来的研究步骤。

3.3 适用地貌判定

比较基于DEM方法与MNDWI方法提取的水文信息,发现无论是何种方法,提取水文信息的准确程度在研究区中的不同区域是不同的,并且提取水文信息准确程度的高低与地貌类别具有极高的相关性,为了判定不同方法对不同地貌的适用程度,在研究区內选取了两处不同地貌的典型区域平原(M1)、山地(M2),如图8所示,并通过目视解译提取出M1、M2中的水文信息,通过两种方法分别在M1、M2中的提取结果与目视解译结果相对比来分析适用性,进而判定适用地貌,其准确率结果如表2所示。

表2 提取准确率比较

图8 不同地貌提取结果对比Fig.8 Comparison of extraction results of different landforms

如图8、表2所示,在平原(M1)中MNDWI方法效果较好,准确率可达90%以上,虽然由于遥感影像分辨率、人为建筑等因素的影响,有些细小水系没有提取出,但是没有像DEM方法这样出现大量与实际不符的伪河道,山区(M2)中DEM方法效果较好,准确率也可达90%以上,河流形态符合实际,而MNDWI方法只提取出了面状水系,故此准确率只有27%。这两种方法存在的局限性,导致它们均难以独自实现对此地貌的水系自动提取。而两种提取水文信息的方法在同一地貌表现出的互补特性,为将两者以地貌为依据进行综合应用提供了可能。

3.4 山前水系的综合提取

依据前述结果,在山地丘陵地貌区,利用DEM数据提取河流的准确率显著高于遥感影像自动提取,而在平原地区结果相反。因此在不同地貌区利用最优方法提取后再行合成可达到自动识别的最佳效果。为达成此目标,如何正确划分地貌区至关重要。

利用地形起伏度区分山地和平原是近十年来被广泛探讨的方法[25-26]。地形起伏度系指特定范围内的最大高度差,是描述区域地形特征的宏观指标,其特点在于具有尺度效应,即通过调节计算窗口的大小,可适用于不同面积的工作范围,利用均值变点法[27]来确定适合本区的计算窗口大小。

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

(4)

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

(3)对每个k(k≥ 2),将样本数列分为前后两段,即X1,X2,…,Xk-1和Xk,Xk+1,…,XN,分别计算前后两段数列的算术平均值Xk1、Xk2及总样本的算数平均值X。

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

(5)

式(5)中:Xt为某个具体分析窗口单位地势度T的对数值。

(6)

式(6)中:Sk为前后两段样本的离差平方和之和。

S为总样本的离差平方和,变点的存在会使S与Sk之间的差距增大,S-Sk的最大值对应的分析窗口,即为最佳分析窗口。

ArcGIS中“焦点分析”模块下计算并统计2×2～32×32,面积为3.24×104～829.44×104m2共31个递增矩形窗口,统计数据如表3所示,绘制出S-Sk随窗口大小变化的散点图并对其进行非线性曲线拟合,在窗口大小为11×11时,S-Sk具有最大值16.85,即适宜计算地形起伏度的窗口大小为0.98 km2。

表3 适宜计算窗口统计结果

选取大小为11×11网格计算出的区内平均起伏度数据,以30 m为界线将研究区地貌划分为平原与山地丘陵区,起伏度小于30 m的地区圈划为平原,否则圈划为山区[28],在山区应用DEM并保留MNDWI提取的湖泊等面状水系,在平原地区应用MNDWI水体指数,将两者基于ArcGIS10.2进行空间综合修正合并,得到图9所示的研究区内综合水系图。

从合成后的结果(图9)来看,山地丘陵区水系提取较完整,湖泊与河流空间位置吻合良好;平原地区受遥感影像分辨率、地表植被、建筑物和人工覆盖的影响,流量较大的河流(如黄河等)可以被近完整的提取,但流量较小的河流断续存在。两类地貌接触地带展现了类似的规律,大规模河流(如玉符河、北大沙河等)基本可以完整的连接,说明无论基于DEM亦或遥感影像自动提取,均正确识别了该河流。更多流量较小的河流被截断在山体丘陵区边界,即依据遥感影像自动提取未能识别该河流,为核实该类河流是否真实存在,对代表性区域进行了人工解译,结果发现该类河流真实存在,但由于是季节性河流或者被周边植被遮挡等原因,导致没有被遥感影像自动提取。全图未出现终止在山地丘陵区边界的平原河流,可以判断在转换带附近基于DEM数据提取方法获得的数据更丰富。值得关注的是,在转换带附近人工建筑对提取结果的干扰在遥感影像自动提取方法上更加明显。如济南市区,该市地表径流丰富,甚至形成了大明湖等较大面积的湖泊,但受人工建筑干扰的影响,流经该区的河流并未完整衔接。

3.5 新构造活动响应

当代构造地貌学研究普遍认识到,随着构造运动的进行,河流会持续记录地貌演化过程,因此对河流地貌分析可反映出长时间尺度的构造活动信息[29-31]。在均衡状态下,河流的空间形态和纵剖面表现相对平滑;当河流流经的断层两盘发生足以影响地表形态的相对运动时,这一均衡状态被打破,进入一种瞬时不均衡状态,即在纵剖面上形成裂点或在地表形成流向变化,这种瞬态地貌会随着时间发生空间迁移,如裂点的溯源侵蚀[32]和弯曲度的变化[33-34],并逐渐消失已达到新的均衡状态。

研究区南部为泰山山脉、徂徕山山脉及由二者围限的莱芜盆地。泰山基底为古老的太古宇泰山岩群,受太古宙-元古宙多次岩浆活动、褶皱构造的影响,变质变形作用明显,主要由斜长角闪岩组成,其次为黑云变粒岩、角闪变粒岩等,其原岩为超基性、基性火山岩和火山凝灰岩建造,后因泰山大幅抬升,部分盖层风化剥蚀,使得基底重新出露;盖层为古生界寒武系-奥陶系灰岩、页岩,主要分布于北侧张夏崮山一带,从老到新出露馒头组、毛庄组、徐庄组、张夏组、崮山组、长山组、凤山组,与基底呈角度不整合接触。受泰山山前断裂等限制,所述山脉普遍表现为南断北超的单斜断块山系[35-36]。

研究区北部为大面积的冲洪积平原。黄河、徒骇河由该部分北东向经过,经黄河等冲、洪积形成了巨厚的第四系沉积,可达60～310 m,自底向顶可分为平原组、黑土湖组、寒亭组、黄河组等岩石地层。该部分地势平缓,起伏较小,最低海拔约为5 m。

研究区以广饶-齐河断裂为界,北部为华北东部盆地断坳区,区内断裂以近东西向为主,如临邑断裂、兴隆断裂等,第四纪以来不活动或仅在早、中更新世发生过活动;南部鲁西断块隆起区内北西向莱芜断裂、泰山西麓断裂、长清断裂等也多为早、中更新世活动,但北西向的铜冶店-孙祖断裂和近东西向的泰山山前断裂等被认为是晚更新世活动断裂[37-39]。

将区内的地层与第四纪断裂资料完成矢量化后叠加到综合水系图中,得到图10所示的综合水系地质图,通过对比研究区的断层和水系分布,分析河流对活动构造的响应:从水系与断裂水平空间分布结果来看,黄河在长清断裂附近、徒骇河在桑梓店断裂附近流向发生了显著变化,但这种变化与走滑断裂引起的河流变形现象存在较明显差异,后者多表现为短距离的快速曲折,断层影响带两侧流向整体保持一致。黄河在广饶-齐河断裂、千佛山断裂附近,徒骇河在夏口断裂、林南断裂附近、马颊河在宁南断裂附近、玉符河在千佛山断裂附近等,水系未发生明显的流向、弯曲程度变化。泰山山前断裂、莱芜断裂本身位于两种地貌区的边界,因此有石汶河、赢汶河河流等终止于该断裂,此种河流变化的起因更多的是基于DEM和遥感影像提取水系效果的差异。但断裂与地貌分区边界吻合良好这一现象本身也可说明其对地貌可能存在较大的控制作用。

图10 综合水系地质图Fig.10 Geological map of integrated drainage system

在断裂与河道交汇处利用GIS软件提取出断裂两侧河道的高程数据,发现断裂两侧河道高程变化不明显,表明河道的现在垂直形态没有表现出受断层影响后的特征。

从总体的河流水平及垂直形态的变化来看,研究区范围内识别的河流目前未有受断层影响的显著表现,至少说明研究区作为新构造运动的弱活动区,断层走滑及垂直运动造成的地貌影响已充分衰减,河流已经达到了平衡状态。

4 结论

以济南及邻区作为研究区,以SRTMDEM 90 m分辨率原始高程数据和Landsat 8 OLI_TIRS卫星影像作为源数据,通过对比不同方法(DEM、NDWI、MNDWI)的特点及在各地貌区的应用效果,尝试了一种能够自动提取上述地区水系的方法。结果显示,在区内的山地丘陵区和平原区,所对应的最优水系提取方法分别为基于DEM和MNDWI的提取方法,两分区的分界线可通过地形起伏度进行确定。通过对比研究区的断层和水系分布,认为研究区内的早中更新世断裂对河流造成的影响已经衰减,并且河流已达到了新的平衡状态。