基于SAR时间序列的建设用地扩展监测——以常州市新北区为例

图1 研究区范围示意图

Fig.1 Location of study area

1.2 数据源及预处理

本研究所使用的数据是从Sentinel-1A数据分发网站(https: //scihub.copernicus.eu)上下载,共计102幅Sentinel-1A影像,所有影像均为VV极化方式、升轨,可以完全覆盖研究区,数据列表如表1所示。SAR影像主要预处理步骤包括SAR影像之间的配准、空间滤波、地理编码和辐射定标等操作。其中配准、多视、空间滤波、地理编码和辐射定标是在SARscape5.1软件中进行,所有数据在使用之前均进行了预处理。

表1 数据列表

Tab.1 Data list

年份	升降轨	极化方式	获取时间	影像数量/个
2015年	A	VV	2015/07/01—2015/12/16	7
2016年	A	VV	2016/01/09—2016/12/22	18
2017年	A	VV	2017/01/03—2017/12/29	30
2018年	A	VV	2018/01/10—2018/12/24	30
2019年	A	VV	2019/01/05—2019/07/04	17

2 研究方法

本文基于在SAR时间序列中,建设用地扩展存在2种结构的时间序列,提出了一种时间序列自适应分段变化监测算法,算法流程如图2所示。算法分为时间序列分段、决策树分类以及众数滤波3个步骤: 首先,对SAR时间序列进行极值点约束时间序列重要点提取; 然后,分割成三分段或者五分段,接着用分段时间序列平均值作为特征,使用决策树提取发生建设用地扩展的像元; 最后,使用众数滤波的方法对结果进行后处理,填补空洞。

图2

图2 算法流程

Fig.2 Flowchart of algorithm

2.1 时间序列自适应分段

2.1.1 极值点约束时间序列重要点提取

极值点约束时间序列重要点提取是指在时间序列分段重要点提取的过程中考虑极值点的影响,将时间序列分段的非首尾重要点限定为时间序列极值点。极值点定义为局部窗口内的最值,可以是局部窗口内的最小值或是最大值,如图3所示,当窗口宽度为3时,时间序列的极值点提取结果。

图3

图3 极值点提取结果

Fig.3 Result of extreme points

时间序列的重要点可以定义为所在区域误差最大,距离区域端点最远的点。如图4所示,点到区域的端点可以采用3种距离度量方法,图中a为点m和p之间的距离,b为点n和p之间的距离,a+b为欧氏距离,c为垂直距离,d为正交距离。周大镯等^[21]已经证明采用垂直距离和正交距离产生的重要点一致,而采用欧氏距离产生的时间序列重要点不一致。

图4

图4 距离度量方式

Fig.4 Type of distance measure

本文采用垂直距离作为距离的量度,下面以“Z”形结构时间序列为示例介绍利用垂直距离来获得时间序列重要点过程: ①以窗口宽度为3提取时间序列极值点,初始化重要点集合(IPs)为空集,将时间序列的起点和终点加入重要点集合,如图5(a)所示,蓝点是重要点,红点是极值点; ②使用线性拟合对时间序列起点与终点之间的数据进行线性拟合,如果相关系数 $R^{2}$ 小于阈值,说明拟合效果较差,则需要继续分割,分割方法是连接起点与终点,遍历时间序列中起点与终点之间的极值点,计算极值点到连接起点与终点线段的垂直距离,选择垂直距离最大(MaxLength)的极值点作为重要点,并将该点加入重要点集合,如图5(b)所示; ③当重要点集合超过3个,每对相邻的重要点分别作为时间序列子段的起点与终点,如图5(c)所示,分别计算子段的垂直距离最大值(MaxLength-1和MaxLength-2),并将子段重要点加入重要点集合,重复步骤②过程,直至所有分段的相关系数 $R^{2}$ 均大于指定阈值或者时间序列子段内没有极值点; ④分段结束后,获得所有重要点,如图5(d)所示,时间序列被分割成6段(a—f)。

图5-1

图5-1 重要点提取过程

Fig.5-1 Process of important points extraction

图5-2

图5-2 重要点提取过程

Fig.5-2 Process of important points extraction

2.1.2 时间序列自适应分段

针对建设用地扩展过程中存在“Z”形结构和“V”形结构时间序列,本文使用分段平均值差的平方来确定时间序列是被分割成“Z”形或者“V”形时间序列,具体过程如下: ①假设该时间序列是“Z”形结构时间序列,被分割成3段,分割示意图如图6(a)所示,分割原则是以上升幅度最大的时间序列分段2作为时间序列的变化分段,变化分段两侧的时间序列作为平稳分段1和平稳分段3,计算分段1的平均值μ₁和分段3的平均值μ₃,再计算均值差的平方,即

(1)

t_{1} = (μ_{1} - μ_{3})^{2}

。

图6

图6 时间序列分段

Fig.6 Time series segmentation

②假设该时间序列是“V”形结构时间序列,被分割成5段,分割示意图如图6(b)所示,分割原则是以上升幅度最大和下降幅度最大的时间序列分段2和分段4作为变化分段,并且上升分段位于下降分段之前,计算分段1的平均值μ₁、分段3的平均值μ₃和分段5的平均值μ₅,再计算均值差的平方,即

(2)

t_{2} = (μ_{3} {- \frac{μ_{1}}{2} - \frac{μ_{5}}{2})}^{2}

。

③比较 $t_{1} 和 t_{2}$ 的大小,如果 $t_{1}$ > $t_{2}$ ,时间序列被分割成3段,反之,时间序列被分割成5段。

2.2 决策树提取目标变化像元

2.2.1 时间序列分段特征提取

完成时间序列分段之后,需要对平稳期时间序列以及整个时间序列进行量化,对于“Z”和“V”形时间序列,分别提取平稳分段1和平稳分段3的平均值μ₁和μ₃以及整个时间序列标准差s。2种结构时间序列的特征值如图7所示。

图7

图7 时间序列特征值

Fig.7 Eigenvalues of time series

2.2.2 构建决策树

决策树易于理解和实现,能够直接体现数据的特点。因此研究采用决策树来识别发生建设用地转变的像元,本文首先通过方差对变化像元进行粗选,过滤掉部分非变化像元,接下来过滤掉在变化前是建设用地的像元,保留变化前是耕地和水体的部分,最后通过阈值确定变化后是建设用地的部分,这样就提取出了发生建设用地扩展的像元。

样本点的分布和决策树如图8所示。决策树输入是待检测的时间序列标准差 $s$ 、分段1的平均值 $μ_{1}$ 、分段3的平均值 $μ_{3}$ ,决策树中阈值通过选择样本确定,在研究区内选择没有发生变化的建设用地、耕地和水体,首先利用样本点的位置获取所有样本点对应的SAR时间序列,分别计算时间序列的平均值 $μ_{建设用地}, μ_{耕地}, μ_{水体}$ 以及每个时间序列标准差的平均值 $s_{建设用地}, s_{耕地}, s_{水体}$ ,以 $s_{建设用地}, s_{耕地}, s_{水体}$ 中的最大值作为 $s_{1}$ , $(μ_{耕地} + μ_{水体}) / 2$ 作为 $n_{2}$ , $(μ_{耕地} + μ_{建设用地}) / 2$ 作为 $n_{3}$ ,阈值的选取如图8(c)所示。

图8

图8 样本选择以及决策树构建

Fig.8 Sample selection and the establishment of decision tree

2.3 众数滤波

对2.2节中识别出的发生建设用地扩展的像元,采用5×5窗口众数滤波的方法对提取结果进行空间滤波,对于中心像元,如果其邻域像元中,变化像元的数目大于非变化像元,则该中心像元也是变化像元,通过众数滤波可以达到去除孤立的变化像元,填补连片变化像元之间的缝隙和平滑变化像元边界的效果。

3 结果与分析

3.1 建设用地扩展识别结果

采用完整率和正确率2个统计指标来评定结果的精度,计算公式分别为:

(3)

COM = \frac{TP}{TP + FN}

(4)

COR = \frac{TP}{TP + FP}

式中: COM为完整率; COR为正确率; TP为在识别出发生建设用地转变的区域,实际上同样发生建设用地转变的像元个数; FN为在识别出未发生建设用地转变的区域,实际上发生了建设用地转变的像元个数; FP为在识别出的发生建设用地转变的区域,实际上未发生建设用地转变的像元个数。

变化检测结果如图9所示,由于缺少真实地表覆盖变化图,本文采用随机选点的方式进行精度评价。首先是在监测出的区域内随机设置500个样本点,利用Google Earth高分影像逐点验证,其中正确检出448个样本点,错误检出52个样本点; 其次在研究区内随机设置1 000个样本点,统计得到有20个变化样本没有被监测出。计算得到正确率为89.60%,完整率为92.73%。

图9

图9 变化检测结果

Fig.9 Results of change detection

3.2 对比分析

为了验证本文方法的实用性,与时间序列DTW相似性分析方法监测结果进行对比,DTW相似性方法首先计算待分类的时间序列与参考时间序列之间的相似性,通过设置阈值提取变化的区域。图10是选择了DTW方法监测结果中的2个区域进行分析,区域1与区域2分别是通过“V”形标准时间序列和“Z”形标准时间序列所提取出的变化区域,区域1实际上发生的是耕地到坑塘水面的转变,这部分像元的时间序列前期与耕地类似,后期与水体类似,因此这部分“V”形时间序列与耕地转变为水体时间序列相似,这就导致了部分耕地转变为水体的像元被误提; 区域2是通过“Z”形标准时间序列提取出的区域,提取结果较好,形态较完整。因此DTW方法对“Z”形时间序列的识别效果较好,对“V”形时间序列识别效果较差。

图10

图10 DTW方法检测结果

Fig.10 Change detection results of DTW

为了说明“V”形结构时间序列的重要性,研究统计了DTW方法和本文方法监测建设用地扩展结果中“Z”和“V”形结构数目,统计结果见表2,其中,本文方法共监测出13 949个像元,DTW方法共监测出11 925个像元,2种方法监测结果表明: 本文方法优于DTW方法,主要原因在于2种监测方法识别出的“Z”形结构数目类似,差距较小,但是识别出的“V”形结构数目差距较大,本文方法监测出“V”形结构数目4 719个,比DTW方法多2 173个,因此利用SAR时间序列研究建设用地扩展时,“V”形结构的时间序列是必须要考虑的。

表2 时间序列结构数目对比

Tab.2 Comparision of time series structure number(个)

方法	“Z”形结构数目	“V”形结构数目	合计
本文方法	9 230	4 719	13 949
DTW方法	9 379	2 546	11 925

同样采用正确率和完整率来评价建设用地扩展监测结果精度,在DTW方法监测结果内随机设置500个样本点,在监测结果外随机设置1 000个样本点,并使用Google Earth影像逐点验证,精度验证结果如表3所示,与DTW方法比较,本文所提出的顾及“Z”和“V”形结构的建设用地扩展监测方法可以不受时间序列变化时间和结构的影响,较好地提取出发生建设用地扩展的区域,正确率提高了1.80百分点,完整率提高了1.27百分点。

表3 DTW与本文方法精度对比

Tab.3 Accuracy of this method is compared with that of DTW

方法	正确/个	错误/个	漏分/个	正确率/%	完整率/%
本文方法	445	55	20	89.60	92.73
DTW方法	439	61	41	87.80	91.46

3.3 建设用地扩展时空特征

以时间序列分割节点作为变化节点,统计结果如表4所示,建设用地扩展时空分布如图11(a)所示。从表4中可以看出,自从2015年以来,新北区建设用地累计扩展557.96 hm²,平均每年扩展139 hm²,其中常州市新北区近2 a的建设用地扩展以2017年和2018年为主,分别增加建设用地209.68 hm²和210.82 hm²。为了揭示常州市新北区建设用地空间扩展方向,以新北区行政区域中心为圆心,绘制可覆盖研究区的圆,并将其等分成8个夹角与面积相等的扇形区域,统计各个扇形区域内建设用地扩展的面积,如图11(b)所示,可以看出新北区建设用地扩展具有方向性差异,城市扩展主要向南和东南方向扩展,西北和西方向几乎没有扩展,这主要是由于南部靠近常州市的主城区,经济发达,受主城区辐射较强。从变化类形来看,如图11(c)所示,建设用地扩展占用耕地552.91 hm²,占用水体5.05 hm²,占用耕地面积明显高于占用水体的面积,说明建设用地扩展以占用耕地为主,并且2016—2018年间占用耕地的面积逐渐增加,新北区的建设用地扩展仍然是以外延扩张为主。

表4 建设用地历年变化转移矩阵

Tab.4 Transfer matrix of construction land changes over the years(hm²)

土地利用类型	建设用地					合计
土地利用类型	2015年	2016年	2017年	2018年	2019年	合计
耕地	16.89	107.05	207.99	209.55	11.43	552.91
水体	0.62	1.42	1.69	1.27	0.04	5.05
合计	17.52	108.47	209.68	210.82	11.47	557.96

图11

图11 建设用地扩展

Fig.11 Expansion of construction land

3.4 方法的适用性和局限性

本文提出的时间序列分段变化检测方法,通过时间序列自适应分段,自动三分段或者五分段,获得稳定分段1和稳定分段3,通过决策树识别稳定分段1和稳定分段3的类别来判断像元是否发生了建设用地扩展。从本文实验结果上看,本文提出的方法既可以提取出发生建设用地扩展的像元也可以识别出建设用地扩展发生的时间点,是一种较为简便的建设用地扩展监测方法。与已有的方法相比,一方面,只需要选择未发生建设用地扩展的样本,不需要选择每种变化类形、变化时点的样本; 另一方面,本文方法充分考虑了时间序列的结构特征,可以同时针对“Z”和“V”形时间序列,提高了变化检测精度。因此本文方法操作简便,能够适用更一般的建设用地扩展监测。

本文提出的建设用地扩展方法,仅仅针对建设用地,没有考虑到其他的变化类形,比如耕地与水体之间的相互转变,本文方法并不能直接适用。但是,如果进一步扩展决策树,使得决策树可以识别出分段1和分段3的类别分别为耕地或水体,并且分段1和分段3的类别不同,说明本文方法是具有识别出其他变化类别的潜力。

4 结论

本文选取了常州市新北区作为研究对象,研究了2015年7月—2019年7月之间的建设用地扩展情况,顾及到变化时间序列存在“Z”和“V”形结构,采用自适应时间序列分段方法实现了变化区域的提取,最后与DTW方法进行了对比分析,其主要结论如下:

1)本文注意到在SAR时间序列中,建设用地扩展过程中存在“Z”和“V”形结构2种情况,提出了基于时间序列自适应分段的变化检测方法,效果较好,正确率为89.60%,完整率为92.73%。

2)DTW方法对“Z”形时间序列提取结果较好,但对“V”形时间序列的提取结果较差,原因在于DTW方法容易将部分耕地转变为水体的像元与耕地转变为暗亮度建设用地的像元发生混淆。

3)新北区建设用地扩展以外延扩张为主,区域性差异显著。新北区南部与常州市主城区相邻,受主城区经济发展影响,建设用地扩展速度较快; 北部距离主城区较远,受主城区影响较弱,建设用地扩展不如南部。

本文通过对Sentinel-1A常州新北地区2015—2019年的时间序列变化监测分析,注意到建设用地扩展存在2种结构,并取得了较好的建设用地监测结果,但本文所提出的方法仅局限于耕地与水体转变为建设用地,在今后的研究中应该注意其他地类的转变情况。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Cui

X G

, Fang

C L

, Liu

H M

, et al.

Assessing sustainability of urbanization by a coordinated development index for an urbanization-resources-environment complex system:A case study of Jing-Jin-Ji region,China

[J]. Ecological Indicators, 2019,96:383-391.

DOI:10.1016/j.ecolind.2018.09.009 URL [本文引用: 1]

[2]

Savage

A M

, Hackett

, Guenard

, et al.

Fine-scale heterogeneity across Manhattan’s urban habitat mosaic is associated with variation in ant composition and richness

[J]. Insect Conservation and Diversity, 2015,8(3):216-228.

DOI:10.1111/icad.2015.8.issue-3 URL [本文引用: 1]

[3]

Fang

C L

, Liu

H M

, Li

G D

, et al.

Estimating the impact of urbanization on air quality in China using spatial regression models

[J]. Sustainability, 2015,7(11):15570-15592.

DOI:10.3390/su71115570 URL [本文引用: 1]

[4]

Feng

T T

, Zhao

J Q

Review and comparison:Building extraction methods using high-resolution images

[C]// Proceedings of the 2009 Second International Symposium on Information Science and Engineering. IEEE, 2009.

[5]

Soh

M B C

Crime and urbanization:Revisited Malaysian case

[J]. Procedia-Social and Behavioral Sciences, 2012,42:291-299.

DOI:10.1016/j.sbspro.2012.04.193 URL [本文引用: 1]

[6]

Wang

, Chen

, Zhu

, et al.

Multi-level classification based on trajectory features of time series for monitoring impervious surface expansions

[J]. Remote Sensing, 2019,11(6):640.

[7]

刘警鉴, 李洪忠, 华璀, 等.

基于Sentinel-1A数据的临高县早稻面积提取

[J]. 国土资源遥感, 2020,32(1):191-199.doi: 10.6046/gtzyyg.2020.01.26.

Liu

J J

, Li

H Z

, Hua

, et al.

Extration of early paddy rice area in Lingao County based on Sentinel-1A data

[J]. Remote Sensing for Land and Resources, 2020,32(1):191-199.doi: 10.6046/gtzyyg.2020.01.26.

[8]

王世新, 田野, 周艺, 等.

基于后向散射模型的多极化SAR影像建筑物高度提取

[J]. 国土资源遥感, 2017,29(2):37-45.doi: 10.6046/gtzyyg.2017.02.06.

Wang

S X

, Tian

, Zhou

, et al.

Building heigh extration from multi-polarization SAR imagery based on backscattering model

[J]. Remote Sensing for Land and Resources, 2017,29(2):37-45.doi: 10.6046/gtzyyg.2017.02.06.

[9]

Ban

, Yousif

O A

Multitemporal spaceborne SAR data for urban change detection in China

[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012,5(4):1087-1094.

[10]

Yousif

, Ban

Improving urban change detection from multitemporal SAR images using PCA-NLM

[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013,51(4):2032-2041.

DOI:10.1109/TGRS.2013.2245900 URL [本文引用: 1]

Multitemporal synthetic aperture radar (SAR) images have been increasingly used in change detection studies. However, the presence of speckle is the main disadvantage of this type of data. To reduce speckle, many local adaptive filters have been developed. Although these filters are effective in reducing speckle in homogeneous areas, their use is often accompanied with the degradation of spatial details and fine structures. In this paper, we investigate a nonlocal means (NLM) denoising algorithm that combines local structures with a global averaging scheme in the context of change detection using multitemporal SAR images. First, the ratio image is logarithmically scaled to convert the multiplicative noise model to an additive model. A multidimensional change image is then constructed using image neighborhood feature vectors. Principle component analysis is then used to reduce the dimensionality of the neighborhood feature vectors. Recursive linear regression combined with fitting-accuracy assessment strategy is developed to determine the number of significant PC components to be retained for similarity weight computation. An intuitive method to estimate the unknown noise variance (necessary to run the NLM algorithm) based on the discarded PC components is also proposed. The efficiency of the method has been assessed using two different bitemporal SAR datasets acquired in Beijing and Shanghai, respectively. For comparison purposes, the algorithm is also tested against some of the most commonly used local adaptive filters. Qualitative and quantitative analyses of the algorithm have demonstrated the efficiency of the algorithm in recovering the noise-free change image while preserving the complex structures in urban areas.

[11]

易凤佳, 李仁东, 常变蓉, 等.

长株潭地区建设用地扩张遥感时空特征分析

[J]. 国土资源遥感, 2015,27(2):160-166.doi: 10.6046/gtzyyg.2015.02.25.

DOI:10.6046/gtzyyg.2015.02.25 URL [本文引用: 1]

2,前后2个时期的扩张面积分别占增加总量的57.30%和42.70%。对长株潭地区建设用地扩张进行时空特征分析,认为建设用地扩张引起区域内耕地和林地在数量和空间格局上的变化,城区周围9~12 km范围内耕地极易转入建设用地; 2个时期河流附近耕地流入建设用地随距离变化而呈现不同特征; 前5 a城区周围林地转入建设用地强度大,距离交通用地3 km附近的林地易转入建设用地。]]>

F J

, Li

R D

, Chang

B R

, et al.

Spatial-temporal features of construction land expansion in Changzhutan (Changsha-Zhuzhou-Xiangtan) area based on remote sensing

[J]. Remote Sensing for Land and Resources, 2015,27(2):160-166.doi: 10.6046/gtzyyg.2015.02.25.

[12]

郭晶鹏, 赵映慧, 谌慧倩, 等.

黑龙江工业城市建设用地时空扩张遥感监测

[J]. 国土资源遥感, 2018,30(3):204-212. doi: 10.6046/gtzyyg.2018.03.28.

Guo

J P

, Zhao

Y H

, Chen

H Q

, et al.

Remote sensing monitoring for temporal and spatial expansion of construction land of industrial cities in Heilongjiang

[J]. Remote Sensing for Land and Resources, 2018,30(3):204-212.doi: 10.6046/gtzyyg.2018.03.28.

[13]

Z X

, Yeh

A G-O

, Li

, et al.

Monthly short-term detection of land development using RADARSAT-2 polarimetric SAR imagery

[J]. Remote Sensing of Environment, 2015,164:179-196.

DOI:10.1016/j.rse.2015.04.018 URL [本文引用: 1]

[14]

Walter

Object-based classification of remote sensing data for change detection

[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2004,58(3-4):225-38.

DOI:10.1016/j.isprsjprs.2003.09.007 URL [本文引用: 1]

[15]

, Jensen

J R

A change detection model based on neighborhood correlation image analysis and decision tree classification

[J]. Remote Sensing of Environment, 2005,99(3):326-340.

DOI:10.1016/j.rse.2005.09.008 URL [本文引用: 1]

[16]

Wang

, Chen

, Li

, et al.

Detecting spatio-temporal and typological changes in land use from Landsat image time series

[J]. Journal of Applied Remote Sensing, 2017,11(3):035006.

DOI:10.1117/1.JRS.11.035006 URL [本文引用: 1]

[17]

Wang

Y F

, Zhao

F F

, Cheng

, et al.

Framework for monitoring the conversion of cultivated land to construction land using SAR image time series

[J]. Remote Sensing Letters, 2015,6(10):794-803.

DOI:10.1080/2150704X.2015.1081304 URL [本文引用: 1]

[18]

, Wang

Y F

, Chen

P P

, et al.

Spatiotemporal fuzzy clustering strategy for urban expansion monitoring based on time series of pixel-level optical and SAR images

[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2017,10(5):1769-1779.

DOI:10.1109/JSTARS.2017.2657607 URL [本文引用: 1]

[19]

Feng

, Zhang

, Liao

Building change detection using coherent and incoherent features from multitemporal SAR images

[C]// Proceedings of the 2019 10th International Workshop on the Analysis of Multitemporal Remote Sensing Images (MultiTemp). IEEE, 2019.

[20]

车美琴, 阿里木·赛买提, 杜培军, 等.

利用旋转不变特征提取全极化SAR影像人工地物

[J]. 遥感学报, 2016,20(2):303-314.

Che

M Q

, Samat

, Du

P J

, et al.

Urban man-made target extraction from Quad-PolSAR imagery with roll invariant parameters

[J]. Journal of Remote Sensing, 2016,20(2):303-314.

[21]

周大镯, 李敏强.

基于序列重要点的时间序列分割

[J]. 计算机工程, 2008,34(23):14-16.

Zhou

D Z

, Li

M Q

Time series segmentation based on series importance point

[J]. Computer Engineering, 2008,34(23):14-16.