EN中文

Remote Sensing for Natural Resources >

2025 , Vol. 37 >Issue 1: 142 - 151

DOI: https://doi.org/10.6046/zrzyyg.2023284

Application of Sentinel-1/2 imagery in identifying hills cutting and land reclaiming in the Beishan Mountain of Lanzhou

  • NIU Quanfu , 1, 2, 3 ,
  • LEI Jiaojiao 1 ,
  • LIU Bo 1 ,
  • WANG Hao 1 ,
  • ZHANG Ruizhen 1 ,
  • WANG Gang 1
Expand
  • 1. School of Civil Engineering, Lanzhou University of Technology, Lanzhou 730050, China
  • 2. Emergency Mapping Engineering Research Center of Gansu Province, Lanzhou 730050, China
  • 3. Academician Expert Workstation of Gansu Dayu Jiuzhou Space Information Technology Co., Ltd., Lanzhou 730050, China

Received date: 2023-09-14

  Revised date: 2023-12-19

  Online published: 2026-06-03

Fold

Abstract

Hills cutting and land reclaiming (HCLR) serves as the most direct way for cities that are constrained by terrain to overcome land scarcity and facilitate urban spatial expansion. Obtaining the range of HCLR quickly and accurately using remote sensing technology is significant for assessing the regional ecosystem and new city development planning scientifically. Based on the Google Earth Engine (GEE) cloud computing platform for remote sensing, as well as multi-temporal data from Sentinel-1 SAR and Sentinel-2 multispectral imager (MSI), this study determined the spatiotemporal distribution of the 2017—2022 HCLR range in the study area using multi-temporal change monitoring and a random forest algorithm. Using a combination of Sentinel-1 ascending and descending images and based on noise filtering and multi-temporal image synthesis, this study calculated the difference in the backscattering intensity before and after HCLR. Then, the excavation range was determined using a threshold determined using the percentile threshold method combined with sample data. The results demonstrate that this method exhibited high operability, with an overall classification accuracy and Kappa coefficient of 85% and 0.83, respectively. By monitoring multi-temporal changes in the VH polarization band of Sentinel-1 and using a combination of Sentinel-1 ascending and descending images, this study acquired the spatial distribution of excavation areas within the study area from 2017 to 2022. Before 2019, the excavation areas were primarily concentrated in Jiuzhou Development Zone, Country Garden, and Poly Lingxiu Mountain. After 2020, new excavation areas, such as Liujiagou and Shuiyuan Station, emerged, with the scope and intensity of excavation gradually increasing. By combining the spectral, texture, and topographic features of SAR and optical data and based on feature optimization combined with a random forest algorithm, this study determined the spatial distribution of yearly HCLR from 2017 to 2022. Before 2018, the HCLR scale was small, with an area of 2.655 km2. After 2019, the scale increased each year, especially in 2021, when the area reached 12.607 km2, accounting for 34.56% of the total land reclamation area during the monitoring period. In 2022, the reclamation area obtained through further excavation on previously reclaimed land was estimated at only 2.686 km2 due to the increasing slope and excavation volume. The method developed in this study for monitoring excavation areas and extracting land reclaiming ranges enables effective monitoring and extraction of the HCLR range.

Key words： mountain cutting and land formation; Sentinel imagery; multi-temporal; change monitoring; random forests

Cite this article

NIU Quanfu , LEI Jiaojiao , LIU Bo , WANG Hao , ZHANG Ruizhen , WANG Gang . Application of Sentinel-1/2 imagery in identifying hills cutting and land reclaiming in the Beishan Mountain of Lanzhou[J]. Remote Sensing for Natural Resources, 2025 , 37(1) : 142 -151 . DOI: 10.6046/zrzyyg.2023284

0 引言

随着我国社会经济的快速发展,城镇化进程所带来的人地关系日趋紧张。因此,一些空间拓展受地形所限的城市不断将城市扩张的建设用地转移到荒山、荒沟等未利用地上,特别是地处典型河谷盆地的兰州市,黄河穿城而过,其“两山夹一沟”的地形地貌特点,使其可供城市拓展的建设用地十分有限[1]。为破解城市拓展与用地短缺的矛盾,对兰州北山的低丘缓坡荒沟等大量未利用地进行削山造地[2-3]。20世纪90年代,九州开发区实施削山造地建设; 2008年,沙井驿、安宁堡、西果园镇实施削山造地工程; 2012年启动的兰州新城建设中预计要推平700多座荒山,开发整理青白石区域低丘缓坡等未利用土地[4]。削山造地是将山丘的一部分或整体进行切割、削平,然后利用削平后的土地进行填筑,创造出新的平坦地面,通常用于城市建设、基础设施建设、土地开发或其他大型建设项目。削山造地区主要为黄土沟壑,生态环境极为脆弱,削山造地对生态环境的负面影响不容低估。因此,利用遥感技术对削山造地区进行动态监测,有效获取削山造地的开挖区和范围,是该区域生态环境影响评价和新城发展规划的基础工作。
自1989年起,针对削山造地城市空间拓展,国内外学者基于遥感技术从不同方面开展了研究工作。白建荣等[5-6]采用无人机遥感影像和人工矢量化等方法,按时间顺序依次提取兰州周边地区削山造地信息,获取了削山造地的规模和发展趋势; 何永刚等[7]利用多源数据和计算机技术构建Skyline系统,实现削山造地的三维监测; 赵元等[8]利用无人机倾斜摄影和土地利用等数据实现了削山造地精细化监测; 赵大卫等[9]利用多时相Landsat TM影像获取了多年新增的削山造地空间位置。上述研究多基于遥感影像和辅助数据,主要针对削山造地的平面位置进行监测,获取其空间分布和地块面积等。削山造地的规模,除了平面位置外,与其挖填方量相关的垂向监测也是重要的一环。高小龙等[10]使用高精度数字高程模型(digital elevation model,DEM)和数字正射影像(digital orthophoto map,DOM)构建三维地形信息来提取兰州削山造地范围。由于研究区削山造地监测区范围大,这类基于无人机摄影测量技术的高精度DEM数据需要频繁的航飞和三维建模来获取,耗时费力且经济消耗高,因而光学遥感数据(例如Landsat和Sentinel-2等)被广泛应用于削山造地等各类地类的分类提取[11⇓⇓⇓⇓-16]。然而,光学遥感影像易受数据获取的时效性、阴影和云覆盖等影响。合成孔径雷达(synthetic aperture Radar,SAR)数据有全天候、穿云透雾、相位干涉以及极化特点,尽管其受几何失真和噪声干扰等的影响[17],但随着深度学习算法出现,基于机器学习的SAR数据可较大程度地避免影响,已广泛应用于地表沉降、植被生长以及水体提取等领域[18⇓-20]。削山造地形成建设地类的过程可借助合成孔径雷达干涉测量(interferometric synthetic aperture Radar,InSAR),即测量2次数据采集之间的雷达相位变化,也可通过量化地表变形区来识别[21⇓-23]。当削山造地前后地表性质(如反射率、粗糙度、介电性质等)发生变化时,SAR数据的后向散射值和相干性也会发生变化[24]。造地区域在削山事件前地表性质的一致性很高,而削山后的变化会降低一致性,因而基于一致性的变化监测方法对削山造地开挖区识别提供了思路。目前在土地分类研究中,鲜有对高差变化的地类识别分类方法。因此,本文将基于SAR数据后向散射变化的地类识别方法应用于削山造地的开挖区识别,探讨削山造地这类在高差上变化的地类识别方法,为可靠识别削山造地开挖区以及丰富土地利用分类方法提供借鉴。
削山造地监测的重点为新增造地和开挖区,新增造地范围主要用于城市拓展空间,而开挖区为工程建设的适宜区。为此,本文以兰州北山为研究区,基于遥感云计算平台(Google Earth Engine,GEE)的Sentinel-1 SAR和Sentinel-2多光谱成像仪(multi spectral imager,MSI)数据,应用多时相变化监测和机器学习算法来获取削山和造地的范围,并基于样本点数据进行验证,为该区域生态环境影响评价和新城发展规划提供基础数据。

1 研究区概况与数据源

1.1 研究区概况

甘肃省省会兰州市,地处中国西北(E102°36’~104°34’,N35°35’~37°07’),总面积约1.3万km2,地形地貌独特,黄河自西南向东北穿城而过,城区位于南北两山之间,如串珠状分布于黄河两岸,南北较窄而东西狭长,是典型河谷型、条带状的内陆城市。随着城市经济的快速发展,其建设用地十分紧张。近年来,为扩张城市发展空间,兰州北山荒山沟壑实施了“削山造地”计划,旨在通过开发黄土丘陵、浅沟和山地来扩大城市发展空间,增加城市容积率和建设用地。

1.2 数据源及其预处理

1.2.1 Sentinel遥感数据获取

本研究采用了Sentinel-1 SAR数据和Sentinel-2光学数据,下载自GEE (https://code.earthengine.google.com/)。其中,Sentinel-1由2颗极轨卫星组成,搭载C波段传感器,重返周期为6 d,包括条带模式(stripmap,SM)、超宽幅模式(extra wide swath,EW)、波谱模式(wave,WM)和干涉宽幅模式(interferometric wide swath,IW)4种成像模式[25]。本文主要以IW模式下的地距多视产品(ground range detected,GRD)为数据源,其2种极化方式为交叉极化(VH)和同极化(VV),含有升、降轨数据,以及地表后向散射系数信息,成像时间为2017—2022年。Sentinel-2A/2B为宽幅高分辨率多光谱成像卫星,该卫星拥有13个波段,包含4个空间分辨率为10 m的可见/近红外波段、6个空间分辨率为20 m的红边/短波红外波段和3个空间分辨率为60 m的波段,重返周期为5 d[26-27]。基于无云清晰、质量优良和无条带噪声原则,与Sentinel-1数据一致也选取了自2017—2022年,经过几何纠正和辐射校正的L1C级影像。

1.2.2 数据预处理

Sentinel-1地距多视影像(ground range detected,GRD)和Sentinel-2 MSI数据预处理均基于GEE云平台进行,基准面为WGS84,UTM投影。对于Sentinel-1数据预处理,首先,在GEE平台调用研究区的Sentinel-1数据,该数据经过去除热噪声、辐射校准和基于SRTM DEM的地形校准。由于研究区具有一定植被分布,而交叉极化VH对森林生物量结构较敏感[28],故本研究选择了IW模式、VH极化波段数据。同时,将所有Sentinel-1 GRD像元后向散射值均通过公式 10 × l g A计算,以分贝为单位表示,为进一步减少数据噪声,将后向散射值不高于-30 dB的像素删除。
对于Sentinel-2多光谱数据,调用filter函数筛选云量小于5%的影像,经过辐射校正、大气校正和几何纠正。利用QA60波段提供的云掩模信息,筛选有少量云影像的像元值,将大于阈值的像元值赋予空值,并掩模云影响[22]。对于GEE中影像镶嵌出现的影像重叠和不连续等问题,通过调用函数median()计算其中位数进行镶嵌[29],选择影像中Blue,Green,Red,NIR,SWIR1和SWIR2等6个常用波段进行影像合成,将Sentinel-2多光谱影像数据裁剪至研究区范围。

1.2.3 样本数据

基于实地调查,研究区土地覆盖类型主要包括削山造地类、建设用地、水域、森林、草地、耕地和裸地7类。本研究的样本数据主要采用GEE在线高分辨率影像并结合研究区Sentinel-2遥感影像,通过目视解译获得。各地物样本数据均匀分布,其中,森林为24个、建设用地为34个、削山造地为36个、水域为16个、草地为23个、裸地为34个、耕地为25个。在研究中,将70%的样本数据用于训练,30%的样本数据用于验证。

1.2.4 地形数据

本文采用的DEM为GEE遥感云计算平台下载的SRTM数据,该数据由美国国家航空航天局提供,空间分辨率为30 m。使用该DEM来计算坡度和曲率。

2 研究方法

2.1 削山造地开挖区识别

削山造地开挖区识别主要采用多时相的变化监测方法。双时相方法是从削山造地后SAR影像中减去削山造地前影像而生成差值影像,而多时相方法是在计算差值之前,将削山造地前后的多景影像分别生成合成影像,然后基于合成影像进行变化监测[30]。与双时相监测相比,多时相监测可减少噪声。本文采用多时相方法进行削山造地开挖区监测,首先收集研究区削山造地前和削山造地后2个时段的Sentinel-1 SAR数据,包括升轨和降轨影像,并对数据预处理。将数据分为升轨、降轨和升轨+降轨3种组合,为进一步消除噪声,对多期升轨、降轨数据集合取中值,升轨+降轨数据集合先取中值再求平均(若同时存在升轨和降轨数据则使用“组合升轨+降轨影像”,否则使用“升轨”或“降轨”)。其次,计算削山造地前后各期数据的强度变化值。强度变化值通常采用前后2期图像强度比值取对数的方法[31]。为计算简便,本文中采取求差值的等价方法,其公式为:
A r a t i o = A p r e - A p o s t,
式中: A r a t i o为后向散射变化值; A p r e为削山造地开挖前后向散射值; A p o s t为削山造地开挖后后向散射值。 A r a t i o值可正可负,正值表示地表开挖后SAR强度的增加。研究发现,并不是所有正值像元都代表地表开挖,为剔除这类误差,本文采用百分位阈值法[32],即将计算的n个像素 A r a t i o按升序排列,计算第90(10)百分位、95(5)百分位、99(1)百分位后向散射变化值,得到3组样本,与谷歌影像同期削山造地解译数据比较确定阈值。百分位值计算公式为:
A r a t i o 0 = ( 1 - α ) A r a t i o j + α A r a t i o j + 1,
j = [ p ( n + 1 ) ],
α = p ( n + 1 ) - j,
式中: j为按照升序排列后的序号; α为比例因子,用于插值计算特定百分位值; p为百分位值对应的概率;[]为数值取整。
研究区位于兰州市的北山,属于荒山浅沟等未利用地,为了掩模城区和平地,提取开挖范围,利用DEM制作坡度和曲率,并结合实地调查将阈值分别设定为5°和-0.005,技术路线如图1所示。
图1 削山造地开挖区提取技术路线

Fig.1 Technical route for extracting excavation areas from mountain cutting and land reclamation

2.2 削山造地范围识别

为有效提取研究区多年来削山造地范围,依据提取地类特点,经实地考察,选取了光谱特征、指数特征、地形特征、极化特征和纹理特征构建特征数据集(表1),结合随机森林算法(random forest,RF)提取削山造地范围。为进一步探讨不同特征组合对提取结果的影响,在参与特征变量选择上,设计了仅光谱、光谱+指数、光谱+纹理、光谱+地形、光谱+极化、全部特征以及RF算法优化后特征共7个方案。最优方案确定方法为: ①在样本点中,随机、独立产生5种不重复的随机数; ②对产生的样本,70%用于训练,30%用于检验,并计算其验证精度; ③将得到的5次验证精度取均值,作为最后的验证精度。
表1 特征变量

Tab.1 Characteristic variables

特征类型 特征变量 变量名称 变量描述/公式
光谱特征 波段 B2,B3,B4,B5,B8,B8a,B11 蓝光、绿光、红光、红边、近红外、窄近红外、短波红外
指数特征 比值植被指数
归一化植被指数
归一化水体指数
归一化建筑指数
改进土壤调节植被指数
增强植被指数
裸露土壤覆盖程度指数		 RVI(ratio vegetation index)
NDVI(normalized difference vegetation index)
NDWI(normalized difference water index)
NDBI(normalized difference built-up index)
MSAVI(modified soil-adujsted vegetation index)
EVI(enhancod vegetation index)
BSI(bare soil index)		 B 8 / B 4
( B 8 - B 4 ) / ( B 8 + B 4
( B 3 - B 8 ) / ( B 3 + B 8
( B 11 - B 8 ) / ( B 11 + B 8
[ ( 2 B 4 + 1 ) - 2 B 4 + 1 ) 2 - 8 ( B 8 B 4 ) ] / 2
2.5 ( B 8 - B 4 ) / ( B 7 + 6 B 4 - 7.5 B 2 + 1 )
[ ( B 4 + B 11 ) - ( B 4 + B 2 ) ] / [ ( B 4 + B 11 ) + ( B 8 + B 2 ) ]
地形特征 高程
坡度		 DEM
Slope
极化特征 VV极化后向散射系数
VH极化后向散射系数		 σ V V
σ V H
纹理特征 二阶矩
对比度
相关性
方差
逆差距
熵		 B 8 a s m , V V a s m , V H a s m
B 8 c o n , V V c o n , V H c o n
B 8 c o r r , V V c o r r , V H c o r r
B 8 v a r , V V v a r , V H v a r
B 8 i d m , V V i d m , V H i d m
B 8 e n t , V V e n t , V H e n t
RF算法是一种非线性和非参数分类器,通过集成学习的思想将多棵决策树集合成森林,再由相互独立的决策树投票表决出分类结果[33]。在构建RF模型时,选取了29个特征变量,为避免过多特征变量参与地物分类会降低分类精度的问题,利用袋外数据(out-of-bag,OOB)产生的OOB误差来计算不同特征变量重要性[34-35],选出最优分类特征变量。研究发现,当决策树个数不低于350时,所有特征OOB误差区域收敛,此时最大特征数被设置为总特征数的平方根。

2.3 后处理与精度分析

为减少分类误差,将提取结果采用形态学方法进行后处理,先去除小像素群的侵蚀和填充空隙的扩张运算,在此基础上,采用滤波去除小颗粒噪声和小孔。精度评价采用混淆矩阵计算,评价指标为总体精度(overall accuracy,OA)、Kappa系数、生产者精度(producer’s accuracy,PA)和用户精度(user’s accuracy,UA)[36]。其中,OA为正确分类样本点数与样本总数量之间的比值; Kappa系数为衡量地面值与预测值结果一致性的指标; PA为该类别的验证样本被正确验证的概率; UA为分类图上落在该类别的验证点被准确划分为该类别的概率。

3 结果与分析

3.1 削山造地开挖区多时相变化监测分析

在靠近水源站的地方(图2)自北向南绘制了横断面线以研究不同轨道数据的后向散射变化,横断面线长为920 m。其中,在2019年水源站还未进行开挖,到2021年已进行削山造地施工。以2019年和2021年为例,监测分析了升轨、降轨和升轨+降轨组合3种方案的削山造地开挖区(图3)。可以看出,无论是升轨、降轨还是升轨+降轨组合数据,都能明显地看出开挖后的后向散射强度比开挖前的后向散射强度要大(图3(a)—(c)),这主要是由于削山造地后的地形特征使得SAR数据的后向散射强度增大。从图3(d)可以看出红色、蓝色、绿色曲线交错,3种组合数据具有相似性,升轨+降轨数据曲线后向散射变化趋势平稳,当距离为0~200 m时升轨后向散射变化值比降轨大; 当距离为200~400 m时降轨后向散射变化值比升轨大; 升轨+降轨数据曲线位于两者中间,由此表明升轨+降轨数据监测可以更有效地反映削山造地开挖区的实际状况。
图2 开挖区空间分布(水源站)

Fig.2 Spatial distribution of excavation area (water source station)

图3 开挖前后的SAR后向散射强度对比

Fig.3 Comparison of SAR backscattering intensity before and after excavation

为确定3种方案的最优,对2021年削山造地开挖范围对比分析,首先基于2021年的数据选取随机点进行试验并计算精度。表2为基于Sentinel-1的VH极化波段采用升轨、降轨和升轨+降轨3种方案进行削山造地开挖区提取的精度分析。可以看出,采用升轨、降轨方案的总体精度分别为77%和73%,Kappa系数分别为0.72和0.67,而采用升轨+降轨组合方案时,精度较前2种方案都要高,其OA和Kappa系数分别为85%和0.83。由于升轨、降轨数据易受卫星轨道和拍摄侧视角的影响,因而采用单一升轨或降轨数据很难准确获取不同坡向的开挖区位置,采用升轨+降轨组合方案进行开挖区监测,可对地形影响区域相互补充,有效地提高了监测精度。
表2 2021年3种方案提取的开挖区统计精度

Tab.2 Statistical accuracy of excavation area extracted by three schemes in 2021

方案精度 升轨 降轨 升轨+降轨组合
OA/% 77 73 85
Kappa 0.72 0.67 0.83

3.2 2017—2022年削山造地开挖区提取

削山造地开挖区通常用于建设用地,填方区用于绿化等,因此对开挖区的监测显得十分重要。对于开挖区识别,本文采用Sentinel-1的VH极化波段的多时相变化监测。根据数据可用性和削山造地施工期选择开挖区的监测时段: 2017年以前,GEE平台可用的数据并不多,其升轨数据为2景,降轨为3景,故削山造地监测前时段确定为2015年1月—2017年2月,后时段为2018年1—2月(冬歇期); 2018年,监测前时段为2018年1—2月,后时段为2019年1—2月; 其他监测年时间段的数据选择参照2018年,各监测年所用Sentinel-1数据如表3所示。选用精度较高的升轨+降轨方案进行2017—2022年逐年开挖区遥感提取。
表3 获取SAR影像数

Tab.3 Obtaining the number of SAR images

监测年份 削山造地前时段/景 削山造地后时段/景
升轨 降轨 升轨 降轨
2017年 2 3 5 8
2018年 5 8 5 9
2019年 5 9 5 10
2020年 5 10 5 15
2021年 5 15 5 10
2022年 5 10 5 10
通过测试,使用95%下的百分位阈值法将后向散射变化值划分为削山造地区和非削山造地区,为消除因噪声产生的误差,对提取结果的单一像元或小于6个像素的区域进行剔除。具体做法为先计算像素的连通域,然后将小于某个阈值的区域用周围的像元值进行替换。图4为采用升轨+降轨组合数据,提取的研究区2017—2022年逐年开挖区空间分布。可见,开挖区主要分布在荒沟等未利用地的缓坡区域,由于2017年Sentinel-1可用数据较少、覆盖不全,仅提取出碧桂园等少量开挖范围,监测到开挖区面积较少。2018年和2019年开挖面积明显增大,重点围绕碧桂园、九洲开发区等区域进行。2020—2021年间,明显增加了刘家沟、水源站和保利领秀山等几个新开挖区,开挖强度达到高峰,特别是2021年,开挖区范围最大,面积达1.63 km2,占监测期间总开挖面积的20.4%。2022年则在现有开挖区基础上继续推进,开挖区面积相对较小。
图4 2017—2022年兰州北山削山造地开挖区空间分布

Fig.4 Spatial distribution of excavation areas in Lanzhou Beishan from 2017 to 2022

3.3 2017—2022年造地范围识别

为讨论本文遥感分类的精度,以2021年遥感数据为例,提取了在不同特征变量组合下的各地类的UA和PA(表4)。从表4中可以看出,在光谱特征基础上加入指数特征,草地精度变化不明显; 除草地外,其他地类提取精度,尤其是UA,均有明显提高。加入纹理后,削山造地类、建筑、草地和耕地有一定的改善,水体、森林和裸地表现并不明显。加入地形特征后,除了建筑,其他地类的提取精度均有一定的改善。加入极化特征后,除裸地外,其余地类均表现出明显的提取效果。在全部特征均参与分类的情况下,除建筑、耕地外的地类有一定的改善,但就数值上看提高并不明显,而经过特征优化后,除森林和耕地变化较小,其他各地类提取的UA和PA均具有明显的提高。为进一步探讨不同特征变量组合下的削山造地类提取精度,仅采用光谱特征进行削山造地范围提取时,其5次的平均总体精度为81%,Kappa系数为0.77; 在光谱特征基础上,统计依次增加指数特征、纹理特征、地形特征、极化特征以及全部特征,得其5次实验的平均OA分别为85%,78%,77%,82%和84%,Kappa系数分别为0.81,0.74,0.71,0.79,0.81和0.87。可见,不同特征变量参与,对削山造地范围的提取精度差异性明显,为避免过多的特征变量参与并不能明显提高削山造地范围识别精度的问题,利用RF算法对全部特征变量组合进行特征重要性分析,选取分类精度最高的一组作为最优特征组合,其中,由于研究区森林和耕地本身较少,因而其UA和PA提升有限,其平均总体精度为89%,Kappa系数为0.87。可见,特征优选后削山造地提取结果能有效保留分类最重要的特征变量,避免数据冗余,提高了分类精度。
表4 分类结果精度统计

Tab.4 Accuracy statistics of classification results for each category (%)

地类 评价
指标		 组合方案
光谱 光谱
+指数		 光谱
+纹理		 光谱
+地形		 光谱
+极化		 全部
特征		 优化
特征
削山造地 PA 84.2 89.0 85.0 80.0 79.1 86.5 92.1
UA 71.4 84.0 80.0 75.0 76.0 83.2 92.8
建筑 PA 85.0 80.0 81.1 83.1 75.3 76.4 83.7
UA 75.7 80.4 85.6 75.2 80.0 71.3 83.2
水体 PA 90.0 84.4 83.0 85.0 82.6 83.9 90.4
UA 83.2 87.4 82.3 89.0 86.6 85.4 90.9
草地 PA 70.4 70.0 72.0 72.8 76.5 75.0 81.3
UA 68.3 67.0 88.0 67.6 82.4 75.0 77.7
森林 PA 75.0 85.0 69.6 83.0 85.1 83.1 76.1
UA 75.0 85.2 68.0 78.7 82.2 82.2 77.2
裸地 PA 72.6 87.0 69.9 78.8 71.6 79.5 82.9
UA 69.7 72.0 77.9 76.2 66.3 74.8 81.4
耕地 PA 69.9 76.6 77.7 81.4 85.1 82.02 78.3
UA 71.3 88.2 73.6 74.4 81.5 78.7 81.1
图5为利用优化特征变量组合和RF算法提取的逐年削山造地范围空间分布。监测期间,造地的总体时空发展趋势呈现“东扩西展,南伸北拓”格局。统计削山造地发展趋势,在地形上,从荒沟缓坡逐渐向中等斜坡荒山地掘进; 在坡度上,主要分布于15°以内。在施工时间和规模上(图6),2017年前,研究区造地范围主要围绕九洲开发区进行,造地面积达8.163 km2,2018—2019年造地范围主要位于碧桂园和保利领秀山,其造地面积分别达2.655 km2和5.219 km2,2020—2021年造地规模最大; 在空间分布上呈现多区造地的局面,特别是2021年,新增了刘家沟、水源站和中央公园等几个施工区,造地面积达12.607 km2,占监测期间总造地面积的34.56%,2022年造地范围主要位于保利领秀山和水源站等区域,由于开挖区坡度逐渐增大,开挖难度增大,造地面积仅为2.686 km2,最终统计得出监测期间共造地36.494 km2。
图5 2017—2022年兰州北山削山造地时空分布

Fig.5 Spatio-temporal distribution of land-making in Lanzhou Beishan from 2017 to 2022

图6 2017—2022年兰州北山削山和造地面积

Fig.6 Area of mountain cutting and land reclamation in Lanzhou Beishan from 2017 to 2022

4 讨论与结论

4.1 讨论

本研究在获取削山和造地空间分布的同时,重点探讨削山造地前后Sentinel-1的差值影像对削山(开挖)区的识别能力。SAR数据地表沉降变化监测常用的方法有合成孔径雷达差分干涉测量技术(differential InSAR,D-InSAR)[37]、永久散射体合成孔径雷达干涉测量技术(persistent scatterer InSAR,PS-InSAR)[38]以及短基线集成干涉合成孔径雷达(small baseline subset InSAR,SBAS-InSAR)[39]等,其中,每一种监测方法具有其适用的环境和条件,且对参与计算的SAR数据数量也有要求,比如: PS-InSAR的SAR数据数量不小于5,而SBAS-InSAR的SAR数据数量不小于15。本文方法在削山造地范围识别中主要探讨了基于施工前后影像后向散射强度变化,通过计算监测前后SAR数据集平均后向散射强度的差值影像来获取削山造地开挖区范围。Mondini等[31]采用前后2期图像强度比值取对数的方法进行变化监测,本文的差值方法也较好地识别出了削山造地开挖区,以2021年削山造地为例,采用升轨+降轨影像提取的开挖区,其OA和Kappa系数分别为85%和0.83。同时,部分地方的SAR差值影像与实际削山造地无关,其主要原因是监测区的荒草地、灌木和森林覆盖等植被变化会干扰削山范围的准确识别。Jung等[22]通过增加前后监测期所用Sentinel-1图像数量来减弱植被变化引起的误差,本文的具体做法为延长削山造地前后的时段长度,例如,自2018年以来,将削山造地前时段从当年的1月份延长至2月底,后时段从次年的1月份延长至2月底,得到削山造地前和后的SAR时间序列数据集(表3),提取结果经样本数据验证,其OA和Kappa系数均较高(表4)。其次,为进一步提高监测精度,将前后监测期的SAR数据经辐射和地形校准,采用Lee滤波和取中值的方法,在有效去噪的同时,可较好体现SAR数据本身反映的变化趋势。同时,由于单一升轨或降轨数据易受卫星的轨道和拍摄侧视角的影响[40],很难准确获取不同坡向的削山造地范围,因而,本文采用升轨+降轨组合数据进行监测,有效地提高了监测精度。
相对于常规方法,本文方法使用灵活,不受算法、适宜环境和条件,以及SAR影像最少数量的限制,可实现削山造地和地质灾害这类地表发生变化情景的快速监测与提取。

4.2 结论

基于GEE遥感云计算平台,结合多时相Sentinel-1 SAR和Sentinel-2 MSI遥感数据,应用变化监测以及随机森林方法获取了兰州北山区的削山和造地范围空间分布。主要结论如下:
1)采用多时相变化监测进行研究区削山造地开挖区提取,通过组合升轨+降轨影像和合理选取削山造地前后时段来增加观测的影像数量,在进行噪声滤除和影像合成后,计算削山造地前后影像后向散射强度变化值,采用百分位阈值法结合谷歌影像数据确定阈值,进行研究区削山造地开挖区范围提取研究,得出其OA和Kappa系数分别达85%和0.83,开挖区提取精度可靠,可操作性强。
2)以Sentinel-1的VH极化波段为数据源,结合多时变化监测法,在对比研究升轨、降轨和升轨+降轨组合数据3种方案研究的基础上,提取了研究区2017—2022年逐年削山造地开挖区空间分布。由于2017年前Sentinel-1数据较少且有缺失,监测到开挖区面积较少。2018—2019年的开挖面积增大,2020—2021年明显增加了刘家沟、水源站和保利领秀山等几个新开挖区,特别是2021年开挖区范围最大,面积达1.63 km2,占监测期间总开挖面积的20.4%。2022年在原开挖基础上推进,由于开挖区的地形坡度逐渐增大,开挖难度增大,使得开挖区面积相对较小。
3)联合SAR和光学数据的光谱特征、纹理特征、极化和地形特征,在特征优化的基础上结合RF算法获取了2017—2022年逐年造地范围空间分布。发现在监测期间,研究区造地的总体时空发展趋势呈现“东扩西展,南伸北拓”的格局。在地形上主要集中于坡度小于15°以内,监测期间总造地面积达36.494 km2。其中,2017年前,研究区造地主要围绕九洲开发区,造地面积达8.163 km2; 2018—2019年主要位于碧桂园和保利领秀山,其造地面积分别达2.655 km2和5.219 km2; 2020—2021年造地规模较大,特别是2021年,新增了刘家沟、水源站和中央公园等开挖区,造地面积达12.607 km2,占监测期间总造地面积的34.56%; 2022年由于开挖区的开挖难度逐渐增大,其造地面积仅为2.686 km2。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
徐昔保, 杨桂山, 张建明. 兰州市城市土地利用优化研究[J]. 武汉大学学报(信息科学版), 2009, 34(7):878-881.

Xu X B, Yang G S, Zhang J M. Urban land use optimization in Lanzhou,China[J]. Geomatics and Information Science of Wuhan University, 2009, 34(7):878-881.

[2]
蒲小武, 王兰民, 吴志坚, 等. 兰州丘陵沟壑区挖方黄土高边坡面临的工程地质问题及稳定性分析[J]. 地震工程学报, 2016, 38(5):787-794.

Pu X W, Wang L M, Wu Z J, et al. Engineering geological problems of loess high excavation slope in loess hilly and gully region of Lanzhou and its stability analysis[J]. China Earthquake Engineering Journal, 2016, 38(5):787-794.

[3]
Niu Q, Bai J, Cheng W, et al. Mapping the dynamics of urban land creation from hilltop removing and gully filling Projects in the river-valley city of Lanzhou,China[J]. Journal of the Indian Society of Remote Sensing, 2022, 50(10):1813-1826.

[4]
Li Y, Li Y, Karácsonyi D, et al. Spatio-temporal pattern and driving forces of construction land change in a poverty-stricken county of China and implications for poverty-alleviation-oriented land use policies[J]. Land Use Policy, 2020,91:104267.

[5]
白建荣. 兰州北山区域削山造地监测技术研究[R]. 甘肃省,甘肃省地图院,2018-05-29.

Bai J R. Research on Monitoring Technology of Land-cutting and Land-making in Beishan Area,Lanzhou[R]. Gansu Province,Gansu Map Institute,2018-05-29.

[6]
白建荣, 王桂钢, 杨腊梅. 兰州周边地区削山造地遥感监测及生态效应分析[J]. 测绘地理信息, 2020, 45(4):26-28.

Bai J R, Wang G G, Yang L M. Remote sensing monitoring and ecological effect assessment of land creation in Lanzhou City[J]. Journal of Geomatics, 2020, 45(4):26-28.

[7]
何永刚, 巨擘. 基于Skyline的兰州削山造地监测三维系统设计与实现[J]. 矿山测量, 2018, 46(1):44-48,108.

He Y G, Ju B. Design and implementation of land creation monitoring 3D system in Lanzhou based on Skyline Mine Surveying[J]. 2018, 46(1):44-48,108.

[8]
赵元, 王华栋. 无人机倾斜摄影在削山造地快速监测中的应用[J]. 矿山测量, 2018, 46(3):27-30.

Zhao Y, Wang H D. Application of UAV oblique photography in rapid monitoring in shaping mountain for land creation[J]. Mine Surveying, 2018, 46(3):27-30.

[9]
赵大卫, 贾永红, 白建荣, 等. 兰州北山区削山造地扬尘浓度预测方法及应用[J]. 武汉大学学报(信息科学版), 2021, 46(7):1106-1113.

Zhao D W, Jia Y H, Bai J R, et al. Prediction method and application of dust from land creation in Lanzhou northern mountain area[J]. Geomatics and Information Science of Wuhan University, 2021, 46(7):1106-1113.

[10]
高小龙, 白建荣, 王桂钢. 基于航空遥感数据的地形变化分析[J]. 地理空间信息, 2020, 18(3):79-81,86,7.

Gao X L, Bai J R, Wang G G. Analysis of the terrain changes based on airborne remote sensing data[J]. Geospatial Information, 2020, 18(3):79-81,86,7.

[11]
Pauleit S, Ennos R, Golding Y. Modeling the environmental impacts of urban land use and land cover change—a study in Merseyside,UK[J]. Landscape and Urban Planning, 2005, 71(2/3/4):295-310.

[12]
Sharma J, Prasad R, Mishra V N, et al. Land use and land cover classification of multispectral Landsat8 satellite imagery using discrete wavelet transform[J]. The International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2018,XLII-5:703-706.

[13]
周珂, 杨永清, 张俨娜, 等. 光学遥感影像土地利用分类方法综述[J]. 科学技术与工程, 2021, 21(32):13603-13613.

Zhou K, Yang Y Q, Zhang Y N, et al. Review of land use classification methods based on optical remote sensing images[J]. Science Technology and Engineering, 2021, 21(32):13603-13613.

[14]
Borràs J, Delegido J, Pezzola A, et al. Land use classification from Sentinel-2 imagery[J]. Revista De Teledetección, 2017(48):55.

[15]
Jiao L, Liu Y. Analyzing the shape characteristics of land use classes in remote sensing imagery[J]. ISPRS Annals of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2012,I-7:135-140.

[16]
张磊, 宫兆宁, 王启为, 等. Sentinel-2影像多特征优选的黄河三角洲湿地信息提取[J]. 遥感学报, 2019, 23(2):313-326.

Zhang L, Gong Z N, Wang Q W, et al. Wetland mapping of Yellow River Delta wetlands based on multi-feature optimization of Sentinel-2 images[J]. Journal of Remote Sensing, 2019, 23(2):313-326.

[17]
Colesanti C, Wasowski J. Investigating landslides with space-borne synthetic aperture Radar (SAR) interferometry[J]. Engineering Geology, 2006, 88(3/4):173-199.

[18]
罗卿莉, 崔峰志, 魏钜杰, 等. SAR影像变化检测的前景特征流形排序法[J]. 测绘学报, 2022, 51(11):2365-2378.

DOI

Luo Q L, Cui F Z, Wei J J, et al. Foreground feature manifold ranking method for SAR image change detection[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(11):2365-2378.

DOI

[19]
Huang H, Chen Y, Clinton N, et al. Mapping major land cover dynamics in Beijing using all Landsat images in Google Earth Engine[J]. Remote Sensing of Environment, 2017,202:166-176.

[20]
罗时雨, 童玲, 陈彦. 全极化SAR图像的山地低矮植被区域土壤含水量估计[J]. 遥感学报, 2017, 21(6):907-916.

Luo S Y, Tong L, Chen Y. Soil moisture estimation for mountainous areas covered with low vegetation using fully polarimetric SAR images[J]. Journal of Remote Sensing, 2017, 21(6):907-916.

[21]
张昊, 高小红, 史飞飞, 等. 基于Sentinel-2 MSI与Sentinel-1 SAR相结合的黄土高原西部撂荒地提取——以青海民和县为例[J]. 自然资源遥感, 2022, 34(4):144-154.doi:10.6046/zrzyyg.2021385.

Zhang H, Gao X H, Shi F F, et al. Sentinel-2 MSI and Sentinel-1 SAR based information extraction of abandoned land in the western Loess Plateau:A case study of Minhe County in Qinghai[J] Remote Sensing for Natural Resources, 2022, 34(4):144-154.doi:10.6046/zrzyyg.2021385.

[22]
Jung J, Yun S H. Evaluation of coherent and incoherent landslide detection methods based on synthetic aperture radar for rapid response:A case study for the 2018 Hokkaido landslides[J]. Remote Sensing, 2020, 12(2):265.

[23]
Mondini A C, Guzzetti F, Reichenbach P, et al. Semi-automatic recognition and mapping of rainfall induced shallow landslides using optical satellite images[J]. Remote Sensing of Environment, 2011, 115(7):1743-1757.

[24]
Handwerger A L, Jones S Y, Huang M H, et al. Rapid landslide identification using synthetic aperture radar amplitude change detection on the Google Earth Engine[J]. Natural Hazards and Earth System Sciences Discussions,2020:1-24.

[25]
Van Tricht K, Gobin A, Gilliams S, et al. Synergistic use of radar Sentinel-1 and optical Sentinel-2 imagery for crop mapping:A case study for Belgium[J]. Remote Sensing, 2018, 10(10):1642.

[26]
Carrasco L, O’Neil A W, Morton R D, et al. Evaluating combinations of temporally aggregated Sentinel-1,Sentinel-2 and Landsat 8 for land cover mapping with Google Earth Engine[J]. Remote Sensing, 2019, 11(3):288.

[27]
宋宏利, 雷海梅, 尚明. 基于Sentinel 2A/B时序数据的黑龙港流域主要农作物分类[J]. 江苏农业学报, 2021, 37(1):83-92.

Song H L, Lei H M, Shang M. Crop classification based on Sentinel 2A/B time series data in Heilonggang River basin[J]. Jiangsu Journal of Agricultural Sciences, 2021, 37(1):83-92.

[28]
Le Toan T, Beaudoin A, Riom J, et al. Relating forest biomass to SAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 1992, 30(2):403-411.

[29]
Svoboda J, Štych P, Laštovička J, et al. Random forest classification of land use,land-use change and forestry (LULUCF) using Sentinel-2 data:A case study of Czechia[J]. Remote Sensing,2022,14(5):1189.

[30]
Lindsay E, Frauenfelder R, Rüther D, et al. Multi-temporal satellite image composites in google earth engine for improved landslide visibility:A case study of a glacial landscape[J]. Remote Sensing, 2022, 14(10):2301.

[31]
Mondini A C, Santangelo M, Rocchetti M, et al. Sentinel-1 SAR amplitude imagery for rapid landslide detection[J]., Remote Sensing, 2019, 11(7):760.

[32]
赵怡卿, 李玫洁, 吴瑕. 包头地区极端温度天气事件的百分位阈值计算与检验[J]. 农业灾害研究, 2023, 13(3):86-88.

Zhao Y Q, Li M J, Wu X. Percentile threshold calculation and test of extreme temperature weather events in Baotou area[J]. Journal of Agricultural Catastrophology, 2023, 13(3):86-88.

[33]
Breiman L. Random forests[J]. Machine Learning, 2001, 45(1):5-32.

[34]
帅爽, 张志, 张天, 等. 特征优化结合随机森林算法的干旱区植被高光谱遥感分类方法[J]. 农业工程学报, 2023, 39(9):287-293.

Shuai S, Zhang Z, Zhang T, et al. Hyperspectral image classification method for dryland vegetation by combining feature optimization and random forest algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering, 2023, 39(9):287-293.

[35]
魏鑫, 任雨, 陈曦东, 等. 基于时序Sentinel-2数据水体动态变化监测——以河南省为例[J/OL]. 自然资源遥感, 2023:1-11.(2023-06-19)[2023-09-11].https://kns.cnki.net/kcms/detail/10.1759.P.20230619.1013.018.html.

Wei X, Ren Y, Chen X D, ed al. Monitoring dynamic changes of water bodies based on time-series Sentinel-2 data:A case study in Henan Province[J/OL] Remote Sensing for Natural Resources, 2023:1-11.(2023-06-19) [2023-09-11].https://kns.cnki.net/kcms/detail/10.1759.P.20230619.1013.018.html.

[36]
Stehman S V. Selecting and interpreting measures of thematic classification accuracy[J]. Remote Sensing of Environment, 1997, 62(1):77-89.

[37]
Berardino P, Fornaro G, Lanari R, et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11):2375-2383.

[38]
Du W B, Ji W Q, Xu L J, et al. Deformation time series and driving-force analysis of glaciers in the eastern Tienshan Mountains using the SBAS InSAR method[J]. International Journal of Environmental Research and Public Health, 2020, 17(8):2836.

[39]
Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry[C]// IEEE 1999 International Geoscience and Remote Sensing Symposium. IEEE,2000:1528-1530.

[40]
Adriano B, Yokoya N, Miura H, et al. A semiautomatic pixel-object method for detecting landslides using multitemporal ALOS-2 intensity images[J]. Remote Sensing, 2020, 12(3):561.

Options
Outlines

/

〈 〉