EN中文

地球信息科学学报 >

2021 , Vol. 23 >Issue 5: 903 - 917

DOI: https://doi.org/10.12082/dqxxkx.2021.200266

遥感科学与应用技术

高分遥感图像相对辐射校正中的伪不变地物自动提取和优化选择

  • 施海霞 , 1 ,
  • 韦玉春 2, 3, 4 ,
  • 徐晗泽宇 , 2, 3, 4, * ,
  • 周爽 2, 3, 4 ,
  • 程琪 2, 3, 4
展开
  • 1.江苏省土地勘测规划院,南京 210000
  • 2.南京师范大学地理科学学院,南京 210023
  • 3.南京师范大学虚拟地理环境教育部重点实验室,南京 210023
  • 4.江苏省地理信息资源开发与利用协同创新中心,南京 210023
*徐晗泽宇（1993— ）,男,甘肃嘉峪关人,博士生,研究方向为遥感图像分析。E-mail:

施海霞（1979— ）,女,江苏南京人,高级工程师,主要研究方向为土地利用及地理信息应用研究。E-mail:

收稿日期: 2020-05-26

  要求修回日期: 2020-08-10

  网络出版日期: 2021-07-25

基金资助

江苏省研究生科研创新计划项目(KYCX20_1179)

国家自然科学基金项目(41471283)

版权

版权所有，未经授权，不得转载、摘编本刊文章，不得使用本刊的版式设计。
收起

Automatic Extraction and Optimal Selection of the Pseudo Invariant Features for the Relative Radiometric Normalization in High-Resolution Remote Sensing Imagery

  • SHI Haixia , 1 ,
  • WEI Yuchun 2, 3, 4 ,
  • XU Hanzeyu , 2, 3, 4, * ,
  • ZHOU Shuang 2, 3, 4 ,
  • CHENG Qi 2, 3, 4
Expand
  • 1. Jiangsu Provincial Land Survey and Planning Institute, Nanjing 210000, China
  • 2. School of Geography, Nanjing Normal University, Nanjing 210023, China
  • 3. Key Laboratory of Virtual Geographic Environment (Nanjing Normal University), Ministry of Education, Nanjing 210023, China
  • 4. Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application, Nanjing 210023, China
*XU Hanzeyu, E-mail:

Received date: 2020-05-26

  Request revised date: 2020-08-10

  Online published: 2021-07-25

Supported by

The Postgraduate Research and Practice Innovation Program of Jiangsu Province(KYCX20_1179)

National Natural Science Foundation of China(41471283)

Copyright

Copyright reserved © 2021
Fold

摘要

相对辐射校正是遥感变化检测中重要的预处理过程,伪不变地物（Pseudo-Invariant Features,PIF）是多时相影像中相对不变的地物,是相对辐射校正中的重要依据。针对高分遥感图像变化检测中相对辐射校正的要求,本文提出了一个自动提取和优化选择PIF的流程和方法：首先计算两期图像的亮度、光谱特征和空间特征的变化向量,然后对各变化向量的像元值从低到高进行排序,经多数投票后提取PIF,最后使用“迭代线性回归—去除异常值”方法选择获得最终PIF。以2016年11月27日和2017年7月18日的2期“北京二号”高空间分辨率多光谱影像为例,选择地物占比不同的两个实验区对流程和方法进行了验证,并与多元变化检测和迭代加权多元变化检测的PIF提取方法进行了比较。使用两期WorldView-2影像和Landsat-8 OLI影像对方法的适用性进行了验证。结果表明：① 2个实验区提取的PIF精度分别为98.74%和98.71%,PIF像元合理分布于未变化区域、包括了影像中主要的地物类型;② 使用本文方法提取的PIF建立的相对辐射校正模型具有显著的线性拟合效果（p<0.000 1）;③ 本文方法考虑了图像亮度、光谱信息以及空间信息的差异,使用参数少,可操作性高;④ 与多元变化检测和迭代加权多元变化检测方法相比,本文方法提取的PIF更为合理,建立的辐射校正方程拟合效果更佳;⑤ 本文方法适用于具有相同波段设置的中、高空间分辨率光学遥感影像。

关键词： 相对辐射校正; 高空间分辨率遥感图像; 伪不变地物; 遥感变化检测; 图像特征; 北京二号; 遥感信息提取; 特征选择

本文引用格式

施海霞 , 韦玉春 , 徐晗泽宇 , 周爽 , 程琪 . 高分遥感图像相对辐射校正中的伪不变地物自动提取和优化选择[J]. 地球信息科学学报, 2021 , 23(5) : 903 -917 . DOI: 10.12082/dqxxkx.2021.200266

Abstract

Relative radiometric normalization is an important process in the remote sensing image change detection. Identifying the Pseudo-Invariant Features (PIF) with invariant or near-invariant radiometric reflectance over a certain period in multi-temporal images is a key to radiometric normalization. This paper proposed a novel method for automatic extraction and optimal selection of PIF. First, the change vectors including brightness, and spectral and spatial domains of bitemporal images were generated. Then, the pixels of each change vector were sorted from the lowest to the highest value, and the majority vote algorithm was used to extract the initial PIF. Finally, the PIF was selected by the iterative linear regression and outlier analysis. Taking two multi-spectral high-resolution Tripesat-2 images acquired on November 27, 2016 and July 18, 2017 as example, two typical regions with different land cover types were selected to test the proposed method. The proposed method was compared with Multivariate Alteration Detection (MAD) and Iteratively Reweighted Multivariate Alteration Detection (IR-MAD) methods. The applicability of the proposed method was also validated by two WorldView-2 and Landsat OLI images. The results show that: (1) the accuracy of the PIF that extracted within two regions were 98.74% and 98.71%, respectively. The extracted PIF was distributed in the unchanged areas and covered the main land cover types in the images; (2) the linear regression models of the relative radiometric normalization using the PIF extracted by the proposed method were significant (p<0.000 1); (3) the differences in image brightness, spectral domain, and spatial domain were taken into account in this method with less parameters and high operability; (4) compared with MAD and IR-MAD, the proposed method showed a better performance in extraction precision and an significant linear regression model of the relative radiometric normalization; and (5) the proposed method was suitable for other medium- or high-resolution remote sensing images with same bands.

Key words： relative radiometric normalization; high-resolution satellite images; pseudo-invariant features; remote sensing change detection; image feature; TripleSat-2; remote sensing information extraction; feature selection

1 引言

遥感变化检测使用同一区域的多时相遥感影像提取变化信息,分析地表的变化特征和过程[1]。相同地物在相近的成像条件下应具备较为一致的辐射水平,然而因传感器状态的差异、太阳辐射、大气环境的变化以及被观测地物的季节变化等,它们在不同时间的影像中并不具备相同的反射率数值[2],从而给影像的变化检测带来误差[3,4]。不同时相遥感影像之间的辐射特征差异是进行变化检测的基本依据,辐射校正通过减少同类地物的差异、突出不同类地物的区别,可以提升变化检测的精度[5,6]。
辐射校正方法包括绝对辐射校正（Absolute Radiometric Correction, ARC）和相对辐射归一化（Relative Radiometric Normalization,RRN）2种[4]。绝对辐射校正使用卫星同步观测数据,通过精确的遥感器校准、大气校正等将影像中像元值转换为地物的辐射亮度或反射率[7,8],实际应用中存在较多限制。相对辐射校正（即归一化）通过调整目标影像的像元值来实现目标影像与参考影像的近似匹配[4,7-8],不需要卫星同步观测数据,计算简便,应用更为广泛。相对辐射校正方法可分为基于分布和基于像元对的方法[9]。前者考虑整个影像的数据分布特征,通过线性拉伸的方法使得目标影像具有和参考影像相近的灰度分布,如直方图匹配、最大—最小归一化等[10],计算简便但容易造成原始光谱的变形[11]。后者使用伪不变地物（Pseudo-Invariant Features, PIF）作为辐射控制样本（Radiometric Control Set Samples）来建立目标影像和参考影像各个波段间的回归关系,然后使用回归方程进行相对辐射校正[12,13,14],该类方法简单易行并被认为具有较好的校正效果[10,11],但该类方法的关键在于PIF的准确选取。
PIF是多时相影像中随时间相对不变的地物,是建立影像相对辐射校正方程的基本依据[14]。针对PIF的提取问题,已有大量学者从理论和实践的角度进行了研究。其中,常用的方法是通过目视解译直接人工选取不同时相影像中相同的房屋屋顶、混凝土建筑等人造地物作为PIF[3,15],人工选取方法费时费力、主观性强。考虑到遥感信息自动提取的需要,自动化提取PIF成为研究的重点,出现了“暗-亮”地物识别[14,16]、回归方法[17,18]、变化向量分析方法[19]、图像自动匹配[20]、多规则处理[21,22]、基于多元变化检测（Multivariate Alteration Detection,MAD）[23,24]的提取方法等。
不同PIF提取方法的复杂性和自动化程度各有特点。如“暗-亮”地物识别方法是一种半自动PIF选取方法,使用缨帽变换计算绿度及亮度分量,通过阈值划分确定“暗-亮”地物作为PIF[16],不同传感器影像的缨帽变换系数差异和阈值的选取困难问题限制了该类方法的应用[22]。Elvidge等[17]提出的自动散点图控制回归（Automatic Scattergram-Controlled Regression,ASCR）方法,通过确定二维散点图中陆地和水体的聚集中心来建立初始回归线,然后通过设置阈值在初始回归线附近选择不变地物点,聚集中心和阈值的确定难以自动化。Canty等[24]将多元变化检测（MAD）方法应用于PIF的自动提取,通过进行2幅影像间的典型相关分析（Canonical Correlation Analysis）来提取不变像元[25],其优点是对图像的线性变换具有不变性、受大气和传感器差异影响较小[7,12-13],然而,MAD中协方差矩阵计算对于异常值的敏感性使其不适用于具有大量变化像元的影像[12,22]。基于MAD的PIF提取方法衍生的多种变形,如迭代加权MAD（Iteratively Reweighted Multivariate Alteration Detection,IR-MAD）[26,27]、多时相MAD（Multitemporal And Multivariate Alteration Detection,MMAD）[12]等,虽然具有较好的应用效果[7,12-13,27],但提取的PIF数量庞大,且未能充分利用遥感影像所包含的丰富的光谱和空间信息。多规则处理方法往往具有较强的专用性,但在其他同类型数据中的应用效果有待验证[21]。除此之外,上述方法在不同遥感影像中的应用还具有一定的性能差异[4,19]。因此,仍有必要构建精度高、适用性较强的PIF自动提取方法。
随着国产卫星遥感数据源的不断丰富,高效的多时相遥感变化信息提取对自动化相对辐射校正产生了更多需求。基于此,本文构建了一个PIF自动提取和选择的方法与流程。使用2期“北京二号”（TripleSat-2）高空间分辨率多光谱影像的两个典型区域进行方法验证,并与MAD和IR-MAD的PIF提取方法进行对比;使用WorldView-2和Landsat-8 OLI影像,对所提方法的适用性进行验证,期望为高空间分辨率遥感影像的变化检测提供方法支持和参考。

2 数据来源及预处理

2.1 遥感影像数据

2.1.1 方法验证影像
使用南京市2期无云雾覆盖的“北京二号”（TripleSat-2）多光谱影像的大气表观反射率（Top of Atmosphere Reflectance, TOA）数据进行方法验证。影像的空间分辨率为3.2 m,包括蓝光、绿光、红光和近红外4个波段（表1）,获取日期分别为2016年11月27日和2017年7月18日。选取位于南京市浦口区的两个典型区域R1和R2作为实验区,大小分别为 800行×600列,600行×800列。实验区位置及假彩色合成影像如图1所示,主要地物类型为建设用地、植被、水体和裸地,不同地物类型的占比具有差异。使用PCI Geomatica软件对两期图像进行正射校正和图像匹配,利用该软件的“快速傅里叶变换相位匹配”模块选择初始控制点,然后经人工选择确认最终控制点完成图像匹配,匹配后两期实验区图像的地表平均误差低于0.5个像元。
表1 “北京二号”影像的主要参数

Tab. 1 Main parameters of the TripleSat-2

类型 参数
卫星轨道信息 太阳同步轨道（SSO）
轨道高度：651 km
升交点地方时（LTAN）：10 : 30
侧摆能力 整星侧摆 ± 45°
影像分辨率/m 全色：0.8;多光谱：3.2
幅宽/km 24
信噪比 > 100:1
MTF 全色：> 0.1;多光谱：> 0.2
量化值/bit 10
成像模式 多景模式
沿轨立体
跨轨立体
条带模式（约4000 km）
区域成像（约40 km×40 km）
重复观测周期/d 1
波段范围/nm 蓝：440~510
绿：510~590
红：600~670
全色：450~650
近红外：760~910
运行寿命/a 7
图1 实验区位置及“北京二号”假彩色合成影像

Fig. 1 The location and the TripleSat-2 images of the study areas

2.1.2 方法适用性验证影像
使用南京市2期无云雾覆盖的WorldView-2高分辨率影像和Landsat-8 OLI中分辨率影像进行方法的适用性验证。WorldView-2影像为原始DN值数据（图2（a）—图2（b））,大小为400行×400列,空间分辨率为1.8 m,获取日期分别为2013年9月16日和2015年7月29日,使用最值归一化方法对每期WorldView-2影像进行归一化处理。Landsat-8 OLI（Operational Land Imager）影像为TOA数据（图2（c）—图2（d））,大小为400行×400列,空间分辨率为30 m,获取日期分别为2017年10月9日和2018年10月28日。图像匹配处理后的2期影像的地表平均误差均低于0.5个像元。
图2 方法适用性验证的假彩色合成影像

Fig. 2 The WorldView-2 and Landsat-8 OLI images for method applicability validation

2.2 地面真值数据

地面真值数据为经过目视解译和对比勘误的地面不变区域分布数据,用来进行PIF的空间分布准确性评价,包括用于方法验证和方法适用性验证的真值数据。

3 研究方法

本文方法的流程如图3所示。以2016年“北京二号”影像作为目标影像（T1）,以2017年影像作为参考影像（T2）。
图3 本文方法流程

Fig. 3 The flowchart of the proposed method

变化的地物在2期图像中应该具有较大差异,“伪不变地物”则应具有较小差异或无差异。从该基本点出发,基于图像特征、变化向量和迭代线性回归方法提取和选择PIF。首先,计算2期图像的亮度、光谱特征和空间特征的变化向量,得到3类特征的差异;然后,根据差异大小排序,通过决策融合提取PIF_initial;最后,使用“迭代线性回归—去除异常值”方法确定最终PIF_select。
通过精度评价和方法对比,对本文方法进行评价。通过方法适用性验证,对本文方法在其他光学遥感影像的应用效果进行评估。

3.1 特征提取及变化向量计算

3.1.1 图像亮度
遥感图像中的地物变化会导致图像亮度变化。使用T1、T2 图像在HSI（色相Hue、饱和度Saturation、强度Intensity）彩色空间中的强度成分 I 、HSV（色相Hue、饱和度Saturation,明度Value）彩色空间中的明度成分 V 作为图像亮度的表达,计算两期图像亮度的变化向量 D I 和 D V 。
选择真彩色和标准假彩色2种典型彩色合成图像作为彩色空间变换的输入,其中,真彩色合成（R,G,B）的对应波段为红光波段B3,绿光波段B2,蓝光波段B1,标准假彩色合成（R,G,B）的对应波段为近红外波段B4,红光波段B3,绿光波段B2。
I 1 = 1 3 R 1 + G 1 + B 1
I 2 = 1 3 R 2 + G 2 + B 2
D I = I 1 - I 2
式中： I 1 、 I 2 分别为T1、T2的强度; D I 为T1、T2强度的变化向量。
V 1 = max R 1 , G 1 , B 1
V 2 = max R 2 , G 2 , B 2
D V = V 1 - V 2
式中： V 1 、 V 2 分别为T1、T2的明度; D V 为T1、T2明度的变化向量。共产生4个变化向量。
3.1.2 光谱特征
（1）TOA
使用T1、T2各波段的光谱变化向量作为图像的TOA差异 D n 。
X 1 = X 11 , X 12 , X 13 , X 14
X 2 = X 21 , X 22 , X 23 , X 24
D n = X 1 - X 2
式中： X 1 、 X 2 分别为T1、T2各波段的TOA值; D n 为T1、T2的各波段光谱的变化向量,其中 n = 1,2 , 3,4 。对于4个波段的高分图像,共产生4个变化向量。
（2）光谱指数
归一化植被指数（Normalized Difference Vegetation Index,NDVI）[28]是进行植被变化信息检测的重要光谱指数。归一化水体指数（Normalized Difference Water Index, NDWI）[29]基于绿光波段和近红外波段的差异,能够有效增强水体信息、突出建筑用地等亮目标信息[22]。对于2期图像中植被、水体及建筑物的变化差异,计算T1、T2的NDVI和NDWI,分别以2期图像NDVI、NDWI的变化向量作为图像的光谱指数差异。
NDVI = ρ NIR - ρ Red ρ NIR + ρ Red
NDWI = ρ Green - ρ NIR ρ Green + ρ NIR
D NDVI = NDV I 1 - NDV I 2
D NDWI = NDW I 1 - NDW I 2
式中： D NDVI 、 D NDWI 分别为T1、T2的NDVI和NDWI的变化向量。共产生2个变化向量。
3.1.3 空间特征
二维Gabor滤波器在提取目标的局部空间和频率信息方面具有良好的特性,被广泛用于无监督变化检测中图像空间信息及上下文背景信息的提取和噪声去除[30,31]。常用的偶对称Gabor滤波器表示为[32]：
h x , y = 1 2 π σ u σ v exp - 1 2 u 2 σ u 2 + v 2 σ v 2 cos ωu
u = x cos θ + y sin θ
v = - x sin θ + y cos θ
式中： θ 为Gabor滤波器的方向; σ u 和 σ v 分别为高斯包络在 u 轴和 v 轴上的标准差（ u 轴平行于 θ 、 v 垂直于 θ ）; ω 为频率调制参数。
基于彩色合成图像,分别计算2期图像在HSV彩色空间中的明度成分作为G1、G2,使用大小为 M × N 的Gabor滤波器 F 对G1、G2进行滤波,得到滤波后的Gabor空间特征 G ˜ （式（17））,最终计算得到空间特征的变化向量 D G 。
G ˜ x , y = ∑ m = - M 2 M 2 ∑ n = - N 2 N 2 G x + m , y + n F m , n
D G = G ˜ 1 - G ˜ 2
式中： M = 3 , N = 3 ; G ˜ 为灰度图像G的滤波后结果; G ˜ 1 、 G ˜ 2 分别为T1、T2的空间特征; D G 为T1、T2的空间特征的变化向量。选择真彩色和标准假彩色2种典型彩色合成图像作为彩色空间变换的输入,共产生2个变化向量。
上述处理共产生12个变化向量。

3.2 PIF的自动提取

PIF的自动提取步骤为：
（1）按照3.1的定义提取变化向量,包括图像亮度 D I 、 D V 、光谱特征 D n 、 D NDVI 、 D NDWI 和空间特征 D G ;
（2）将变化向量按像元值从小到大进行排序;
（3）取排序后变化向量的前30%像元作为具有较小区分度的 D I ' 、 D V ' 、 D n ' 、 D NDVI ' 、 D NDWI ' 和 D G ' ,并转换为二值图像;
（4）使用多数投票算法（Majority vote algorithm）[33]对12个二值图像进行决策级别融合,将所有二值图像加和大于等于6的像元作为初步提取的PIF_initial。

3.3 PIF的优化选择

初步提取的PIF_initial内仍然可能存在较大的差异。为了提高PIF的准确性,构建了“迭代拟合—去除异常值”的优化选择方法。
使用PIF_initial得到2期图像相同波段的初始辐射控制样本。① 对于每个波段,使用最小二乘法求解一元线性回归方程 y = α ˆ + x β ˆ 及其上下95%置信区间（图4）;② 计算各点到拟合直线的垂直距离 d ,计算最大垂直距离 d max ,删除距离大于0.2倍 d max 的点并再次拟合;③重复步骤①、②,当占当前样点总数95%的点分布于上下95%置信区间内,且最大距离 d max 小于1.5倍平均距离 d mean 时停止迭代。若不满足该条件,则以0.1的步长增加 d mean 的倍数,至 d max < ( 1.5 + 0.1 n ) d mean 时终止,其中n为增加步长的个数。
图4 2期图像中指定波段的PIF的线性拟合和异常点

Fig. 4 The linear fit and outlier of the PIF for a specific band of bitemporal images

经上述处理得到PIF_B1、PIF_B2、PIF_B3、PIF_B4,使用求交运算取其共同部分作为最终的PIF_select。
2期影像相同波段的相对辐射校正线性回归直线方程表示为：
y = α ˆ + x β ˆ
式中： x 、 y 分别为影像T1、T2中各波段PIF的TOA值; α ˆ 、 β ˆ 分别为截距和斜率。使用最小二乘法求解,得到截距 α ˆ 和斜率 β ˆ 分别为：
β ˆ = ∑ i = 1 n x i - x ̅ y i - y ̅ ∑ i = 1 n x i - x ̅ 2 α ˆ = y - x β ˆ
式中： n 为PIF样本个数; x ̅ 、 y ̅ 分别为T1、T2中各波段PIF的平均TOA值。各样点距离拟合直线的垂直距离 d 表示为：
d = - y A + β ˆ x A + α ˆ 1 2 + β ˆ 2

3.4 精度评价

从空间分布准确性和拟合效果2个方面对PIF的精度进评价。
（1）空间分布准确性评价：首先,根据影像显示,对PIF所在位置的地物类别进行判别;然后,使用地面不变区域真值数据计算PIF的提取精度 P ,公式如下：
P = N a N
式中： N 为PIF总点数; N a 为分布于不变区域内的PIF点数。
（2）拟合效果评价：在提取精度大于85%的前提下,使用PIF构建相对辐射校正模型,即一元线性回归方程,使用回归方程的决定系数R2和均方根误差 RMSE（Root Mean Squared Error,式（23））,对拟合效果进行定量评价。
RMSE ( X , h ) = 1 m ∑ i = 1 m h x i - y i 2

3.5 方法对比

多元变化检测（MAD）[23]和迭代加权MAD（IR-MAD）[26]是常用的PIF提取方法[12]。分别使用MAD和IR-MAD对实验区的影像进行计算,选取包含主要变化信息的第1、2分量,通过最大类间方差法[34]选取分割阈值划分变化与未变化区域,以未变化区域作为提取的PIF,与本文方法最终提取的PIF进行对比。

3.6 方法适用性验证

选取南京市2期无云雾覆盖的WorldView-2高分辨率影像和Landsat-8 OLI中分辨率影像,对本文方法在其他光学遥感影像的应用效果进行评估。其中,Landsat-8 OLI影像仅使用蓝光、绿光、红光和近红外波段。

4 结果与分析

4.1 本文方法的PIF提取结果

4.1.1 PIF提取精度和分布
本文方法在2个实验区提取的PIF如图5所示。对比人工目视解译的不变区域,提取精度 P 1 、 P 2 分别为98.74%和98.71%。实验区R1共提取39 559个PIF像元,对应地物类型分别为建设用地（19 559个）、植被（17 899个）和水体（2101个）;实验区R2共提取25 779个PIF像元,对应地物类型分别为建设用地（6023个）、植被（18 462个）、水体（1208个）和裸地（86个）。PIF涵盖了各影像范围内的主要地物类别,满足人工选择PIF准则中对于暗、亮地物选取以及土地利用/覆盖类型分布的要求[10,35]。
图5 本文方法的PIF分布

Fig. 5 The PIF extracted by the proposed method

4.1.2 相对辐射校正方程的性能
使用PIF_initial和PIF_select构建的两期图像相同波段辐射控制样本散点图如图6和图7所示。在原始图像中,各波段样点分布散乱,由于两期图像的辐射亮度差异,各地物类别和波段的一致性较差（图6（a）—图6（d）、图7（a）—图7（d））。对于PIF_initial,除B4波段外,图中样点均呈纺锤状分布于拟合直线两侧,波段间的辐亮度差异被进一步缩小（图6（e）—图6（h）、图7（e）—图7（h））。对于PIF_select,图中样点沿拟合直线两侧均匀分布,聚集分布于TOA值的低值和中值区域,高值区域分布较少（图6（i）—图6（l）、图7（i）—图7（l））。由于本文方法使用自适应的距离值作为迭代拟合过程的终止条件,PIF_select的样点具有不同的主轴宽度及分布。
图6 实验区R1中各波段辐射控制样本点散点图（本文方法）

注：横纵坐标分别为目标影像T1和参考影像T2的像元值。

Fig. 6 The scatter plot of radiometric control set samples in region R1 (the proposed method)

图7 实验区R2中各波段辐射控制样本点散点图（本文方法）

注：横纵坐标分别为目标影像T1和参考影像T2的像元值。

Fig. 7 The scatter plot of radiometric control set samples in region R2 (the proposed method)

从原始图像到PIF_select的决定系数R2和均方根误差RMSE（表2）可知,实验区R1中,B1-B4的R2分别由0.4958、0.4140、0.3818、0.1067提升至0.9540、0.9624、0.9720、0.9267;RMSE分别由0.0131、0.0163、0.0227、0.0419降低至0.0037、0.0037、0.0041、0.0093。实验区R2中,B1—B4的R2分别由0.5033、0.3924、0.3334、0.1205提升至0.9562、0.9529、0.9737和0.7884;RMSE分别由0.0119、0.0162、0.0262、0.0398降低至0.0026、0.0032、0.0032、0.0140。基于初步提取的PIF_initial, B1-B3波段校正模型的R2均大于0.80,已具有一定的线性相关性;基于优化选择的PIF_select,各波段校正模型具有高显著性水平（p<0.0001）。PIF的优化选择提高了相对辐射校正方程的性能,其中,对B4波段的相对辐射校正方程的改进最为明显。
表2 基于PIF的各波段相对辐射校正方程（本文方法）

Tab. 2 The relative radiation normalization equation of each band based on PIF (the proposed method)

波段 R1 R2
一元线性回归方程 R2 RMSE 一元线性回归方程 R2 RMSE
原始图像 B1 y=0.6978x+0.0685 0.4958 0.0131 y=0.7654x+0.0530 0.5033 0.0119
B2 y=0.6133x+0.0774 0.4140 0.0163 y=0.6533x+0.0670 0.3924 0.0162
B3 y=0.6212x+0.0718 0.3818 0.0227 y=0.6976x+0.0597 0.3334 0.0262
B4 y=0.3412x+0.1846 0.1067 0.0419 y=0.3629x+0.1941 0.1205 0.0398
PIF_initial B1 y=0.7740x+0.0499 0.8917 0.0066 y=0.8233x+0.0377 0.8719 0.0065
B2 y=0.7542x+0.0487 0.8715 0.0079 y=0.7685x+0.0419 0.8572 0.0077
B3 y=0.7688x+0.0403 0.8378 0.0115 y=0.8027x+0.0312 0.8194 0.0118
B4 y=0.6207x+0.1148 0.3251 0.0370 y=0.5447x+0.1404 0.2572 0.0369
PIF_select B1 y=0.7748x+0.0507 0.9540 0.0037 y=0.8242x+0.0374 0.9562 0.0026
B2 y=0.7630x+0.0479 0.9624 0.0037 y=0.7760x+0.0405 0.9529 0.0032
B3 y=0.7724x+0.0415 0.9720 0.0041 y=0.8073x+0.0304 0.9737 0.0032
B4 y=0.7547x+0.0820 0.9267 0.0093 y=0.6577x+0.1155 0.7884 0.0140
总体上,最终的PIF对各波段的相对辐射校正方程的拟合效果均有改进。由于实验区R2中的植被占比较高,且2期图像的获取的日期差异较大,因此,对于不同实验区影像,B4波段的辐射校正方程的拟合效果具有一定差异。
4.1.3 PIF空间分布的合理性分析
本文方法提取的PIF（图5）能够集中于“严格不变”的像元部分,空间分布更为合理。
图8为典型区域 PIF的放大显示。对比显示表明,两期实验影像存在成像条件差异及地物光谱的变化差异。图8（a）—图8（b）和图8（g）—图8（h）分别为具有较大占地面积的高速入口和厂房,由于光照条件及成像角度差异,大型建筑的顶棚和部分道路具有较为直观的亮度差异。图8（c）—图8（d）和图8（i）—图8（j）分别为近城区绿地和郊区的植被区域,由于植被生长的季节差异,造成视觉上较为一致的植被区域具有一定的辐射差异。图8（e）—图8（f）和图8（k）—图8（l）区域分别为河渠和湖泊,由于水体性质差异、水体中悬浮物、叶绿素和有机溶解物质等对水体的光学特征的影响,使得水体范围内像元并不具有完全一致的辐射水平。仅使用目视甄别及一般PIF筛选方法,难以对上述差异进行有效的识别。
图8 实验区典型区域PIF的放大显示

Fig. 8 The zoomed in PIF in the typical region of R1 and R2

本文方法不仅可以避免选取具有明显地表覆盖变化的像元,还能够有效去除相同地表覆盖下辐射差异明显的像元,从而将PIF集中于“严格不变”的像元部分。本文方法在PIF的自动提取和优化选择方面具有更好的性能。
通过可靠的PIF进行相对辐射校正,是提升相同地表覆盖物间辐射一致性的有效手段,有助于避免变化检测中无关变化的检出、从而有效地提高检测的精度[6]。

4.2 对比方法的PIF提取结果

使用MAD提取的2个实验区的PIF如图9（a）—图9（b）所示,提取精度分别为96.15%和92.34%。使用IR-MAD提取的2个实验区的PIF如图9（c）—图9（d）所示,提取精度分别为95.70%和88.54%。在PIF提取精度方面,对比方法低于本文方法,在空间分布方面,对比方法得到的PIF数量过多,涵盖了影像中未发生明显变化的各主要区域,代表性较差。
图9 MAD和IR-MAD的PIF分布

Fig. 9 The PIF extracted by MAD and IR-MAD

使用PIF_MAD和PIF_IR-MAD构建的两期图像相同波段辐射控制样本散点图如图10和图11所示。图中,B1—B3波段均呈菱形分布于拟合直线两侧,B4波段均呈不规则梯形分布,存在大量明显远离拟合直线的异常点。
图10 实验区R1中各波段辐射控制样本点散点图（MAD和IR-MAD）

注：横纵坐标分别为目标影像T1和参考影像T2的像元值。

Fig. 10 The scatter plot of radiometric control set samples in region R1 (MAD and IR-MAD)

图11 实验区R2中各波段辐射控制样本点散点图（MAD和IR-MAD）

注：横纵坐标分别为目标影像T1和参考影像T2的像元值。

Fig. 11 The scatter plot of radiometric control set samples in region R2 (MAD and IR-MAD)

使用PIF_MAD和PIF_IR-MAD建立的相对辐射校正方程及其决定系数R2和均方根误差 RMSE如表3所示。基于PIF_MAD建立的校正方程中,各波段的R2均小于0.70,实验区R2各波段的拟合效果略好于实验区R1,实验区R1中的B4波段具有最差的拟合效果（R2仅0.1169）;基于PIF_IR-MAD建立的校正方程中,各波段的R2均小于0.70, B4波段的拟合效果较MAD略有提升。
表3 基于PIF的各波段相对辐射校正方程（MAD和IR-MAD）

Tab. 3 The relative radiation normalization equation of each band based on PIF (MAD and IR-MAD)

波段 R1 R2
一元线性回归方程 R2 RMSE 一元线性回归方程 R2 RMSE
PIF_MAD


 B1 y=0.6823x+0.0682 0.6434 0.0090 y=0.7945x+0.0460 0.6680 0.0084
B2 y=0.5894x+0.0769 0.5655 0.0108 y=0.6883x+0.0582 0.6028 0.0106
B3 y=0.6092x+0.0676 0.5554 0.0148 y=0.7523x+0.0454 0.6073 0.0156
B4 y=0.3837x+0.1764 0.1169 0.0427 y=0.4517x+0.1786 0.1575 0.0404
PIF_IR-MAD


 B1 y=0.7700x+0.0570 0.6549 0.0085 y=0.8699x+0.0373 0.6342 0.0082
B2 y=0.6665x+0.0693 0.5500 0.0109 y=0.7440x+0.0542 0.5060 0.0119
B3 y=0.6839x+0.06264 0.4901 0.0164 y=0.8381x+0.0412 0.4656 0.0198
B4 y=0.4750x+0.1638 0.1871 0.0376 y=0.5272x+0.1701 0.1990 0.0381
在PIF提取精度方面,与本文方法相比,MAD方法较为相近,IR-MAD相差较大。在建立的相对辐射校正方程中,本文方法具有更高的精度和显著性水平,而MAD和IR-MAD各波段的回归方程存在有较大差异,难以直接使用,表明本文的PIF选择方法是有效的。

4.3 方法适用性的验证结果

对于WorldView-2和Landsat-8 OLI影像,本文方法提取的PIF如图12所示。在2期WorldView-2和Landsat-8 OLI影像中,分别得到9640个PIF像元（图12（a））和13 307个PIF像元（图12（b））,提取精度分别为99.15%和99.78%,涵盖了影像中的各主要地物类别,满足PIF的选择要求。
图12 WorldView-2（a）和Landsat-8 OLI（b）影像的PIF分布

Fig. 12 The PIF of WorldView-2 (a) and Landsat-8 OLI (b) images

基于WorldView-2和Landsat-8 OLI影像提取的PIF构建的相同波段辐射控制样本散点图分别如图13（a）—图13（d）和图13（e）—图13（h）所示。通过自动提取和优化选择,各图中散点能够均匀分布于拟合直线两侧。
图13 WorldView-2和Landsat-8 OLI影像中各波段辐射控制样本点散点图

注：横纵坐标分别为目标影像T1和参考影像T2的像元值。

Fig. 13 The scatter plot of radiometric control set samples in WorldView-2and Landsat-8 OLI images

基于对应PIF建立的各传感器图像的辐射校正方程及其决定系数R2和均方根误差RMSE如表4所示。在WorldView-2影像的PIF建立的校正方程中,各波段的R2均大于0.84,其中,B3波段具有最高的R2（0.9112）,各波段的RMSE均低于0.01,其中,B1波段具有最低的RMSE（0.0050）;在Landsat-8 OLI影像的PIF建立的校正方程中,各波段的R2均大于0.98,各波段RMSE均小于0.0080。
表4 基于PIF的各波段相对辐射校正方程（方法适用性验证）

Tab. 4 The relative radiation normalization equation of each band based on PIF (method applicability validation)

波段 WorldView-2 Landsat-8 OLI
一元线性回归方程 R2 RMSE 一元线性回归方程 R2 RMSE
B1 y=0.6550x-0.0076 0.8672 0.0050 y=0.9660x+0.0033 0.9882 0.0017
B2 y=1.1694x-0.0186 0.8408 0.0096 y=0.9457x+0.0046 0.9899 0.0018
B3 y=0.6819x-0.0054 0.9112 0.0057 y=0.9534x+0.0054 0.9909 0.0024
B4 y=0.2154x+0.0845 0.8446 0.0135 y=0.9391x+0.0015 0.9826 0.0070
结果表明,本文方法能够适用于具有相同波段设置的中、高分辨率光学遥感影像,PIF具有较高的提取精度,使用PIF建立的相对辐射校正方程具有较好的拟合效果。

5 结论与讨论

面向高空间分辨率卫星图像相对辐射校正中的PIF提取,本文基于图像特征、变化向量和迭代线性回归方法,提出了一个自动提取和优化选择PIF的流程和方法。使用2期“北京二号”多光谱影像对该方法进行了验证,并与MAD和IR-MAD的PIF提取方法进行了对比,使用两期WorldView-2和Landsat-8 OLI影像对本文方法的适用性进行了验证。实验结果与分析表明：
（1）本文方法提取的PIF具有较高精度,能够合理分布于影像的未变化区域并涵盖影像范围内的主要地物类别,满足人工选择PIF准则的各项要求：本文方法在“北京二号”影像两个实验区的提取精度分别为98.74%和98.71%;
（2）使用本文方法提取的PIF建立的相对辐射校正模型具有显著的线性拟合效果：在使用初步提取的PIF_initial建立的相对辐射校正模型已具有较好的线性关系的基础上,通过优化选择,能够有效提高各波段相对辐射校正方程的性能,特别对B4波段的校正方程的改进最为明显。使用优化选择后的PIF_select建立的相对辐射校正模型具有高显著水平（p<0.000 1）;
（3）本文方法针对具有4个波段的高空间分辨率光学遥感影像设计,在特征差异计算过程中考虑了影像变化造成的图像亮度、光谱信息以及空间信息的差异,提取过程中仅涉及差异比例作为提取参数,迭代线性回归的选择过程中仅涉及距离判断和迭代终止参数,参数少、可操作性高;
（4）在使用相同影像的情况下,与多元变化检测（MAD）和迭代加权多元变化检测（IR-MAD）的PIF提取方法相比,本文方法提取的PIF代表性强、空间分布更为合理,在提取精度和所建立的相对辐射校正模型的线性拟合效果方面均具有优势;
（5）本文方法适用于具有相同波段设置的中、高分辨率光学遥感影像中的PIF提取。
今后的工作中,还需对多源光学遥感影像相对辐射校正中的PIF提取展开探究和论证。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
Singh A. Review Article digital change detection techniques using remotely-sensed data[J]. International Journal of Remote Sensing, 1989,10(6):989-1003.

DOI

[2]
Teillet P M. Image correction for radiometric effects in remote sensing[J]. International Journal of Remote Sensing, 1986,7(12):1637-1651.

DOI

[3]
郭丽峰, 高小红, 亢健, 等. 伪不变特征法在遥感影像归一化处理中的应用[J]. 遥感技术与应用, 2009,24(5):588-595.

[ Guo L F, Gao X H, Kang J, et al. Application of the pseudo-invariant feature in normalization process of the remote sensing images[J]. Remote Sensing Technology and Application, 2009,24(5):588-595. ]

[4]
Hong G, Zhang Y. A comparative study on radiometric normalization using high resolution satellite images[J]. International Journal of Remote Sensing, 2008,29(2):425-438.

DOI

[5]
李德仁. 利用遥感影像进行变化检测[J]. 武汉大学学报·信息科学版, 2003,28(S1):7-12.

[ Li D R. Change detection from remote sensing images[J]. Journal of Wuhan University, 2003,28(S1):7-12. ]

[6]
Gómez C, White J C, Wulder M A. Optical remotely sensed time series data for land cover classification:A review[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016,116:55-72.

DOI

[7]
Zhang Y, Yu L, Sun M, et al. A mixed radiometric normalization method for mosaicking of high-resolution satellite imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017,55(5):2972-2984.

DOI

[8]
Garcia-Torres L, Caballero-Novella J J, Gómez-Candón D, et al. Semi-automatic normalization of multitemporal remote images based on vegetative pseudo-invariant features[J]. PLoS One, 2014,9(3):e91275.

DOI

[9]
Yu X, Zhan F B, Hu J, et al. Radiometric normalization for multitemporal and multispectral high resolution satellite images using ordinal conversion[C]. Proceedings of the 2010 18th International Conference on Geoinformatics. IEEE, 2010.

[10]
余晓敏, 邹勤. 多时相遥感影像辐射归一化方法综述[J]. 测绘与空间地理信息, 2012(6):8-12.

[ Yu X M, Zou Q. Methods of radiometric normalization for multi-temporal remote sensing images: a review[J]. Geomatics & Spatial Information Technology, 2012(6):8-12. ]

[11]
黄启厅, 覃泽林, 曾志康. 多源多时相遥感影像相对辐射归一化方法研究[J]. 地球信息科学学报, 2016,18(5):606-614.

DOI

[ Huang Q T, Qin Z L, Zeng Z K. A study on the relative radiometric normalization of multi-sources and multi-temporal remote sensing data[J]. Journal of Geo-information Science, 2016,18(5):606-614. ]

[12]
Lin B, Wang Z, Syariz M A, et al. Pseudoinvariant feature selection using multitemporal MAD for optical satellite images[J]. IEEE Geoscience and Remote Sensing Letters, 2019,16(9):1353-1357.

DOI

[13]
Byun Y, Han D. Relative radiometric normalization of bitemporal very high-resolution satellite images for flood change detection[J]. Journal of Applied Remote Sensing, 2018,12(2):026021.

[14]
Schott J R, Salvaggio C, Volchok W J. Radiometric scene normalization using pseudoinvariant features[J]. Remote Sensing of Environment, 1988,26(1):1-16.

DOI

[15]
Coppin P R, Bauer M E. Processing of multitemporal Landsat TM imagery to optimize extraction of forest cover change features[J]. IEEE Transactions on Geoscience and Remote Sensing, 1994,32(4):918-927.

DOI

[16]
Hall F G, Strebel D E, Nickeson J E, et al. Radiometric rectification: Toward a common radiometric response among multidate, multisensor images[J]. Remote Sensing of Environment, 1991,35(1):11-27.

DOI

[17]
Elvidge C, Yuan D, Weerackoon R D, et al. Relative radiometric normalization of Landsat Multispectral Scanner (MSS) data using an automatic scattergram-controlled regression[J]. Photogrammetric Engineering & Remote Sensing, 1995,61(10):1255-1260.

[18]
Zhong C, Xu Q, Li B. Relative radiometric normalization for multitemporal remote sensing images by hierarchical regression[J]. IEEE Geoscience and Remote Sensing Letters, 2016,13(2):217-221.

DOI

[19]
Sadeghi V, Ebadi H, Ahmadi F F. A new model for automatic normalization of multitemporal satellite images using Artificial Neural Network and mathematical methods[J]. Applied Mathematical Modelling, 2013,37(9):6437-6445.

DOI

[20]
Deng X, Wang C, Lei B. Automatic relative radiometric normalization algorithm based on pseudo-invariant neighborhood[C]. Proceedings of the 2008 Congress on Image and Signal Processing. IEEE, 2008.

[21]
De Carvalho A O, Guimarães F R, Silva c n, et al. Radiometric Normalization of temporal images combining automatic detection of pseudo-invariant features from the distance and similarity spectral measures, density scatterplot analysis, and robust regression[J]. Remote Sensing, 2013,5(6):2763-2794.

DOI

[22]
Zhou H Z, Liu S M, He J J, et al. A new model for the automatic relative radiometric normalization of multiple images with pseudo-invariant features[J]. International Journal of Remote Sensing, 2016,37(19):4554-4573.

DOI

[23]
Nielsen A A, Conradsen K, Simpson J J. Multivariate Alteration Detection (MAD) and MAF postprocessing in multispectral, bitemporal image data: New approaches to change detection studies[J]. Remote Sensing of Environment, 1998,64(1):1-19.

DOI

[24]
Canty M J, Nielsen A A, Schmidt M. Automatic radiometric normalization of multitemporal satellite imagery[J]. Remote Sensing of Environment, 2004,91(3):441-451.

DOI

[25]
白洋. 基于核典型相关分析的遥感图像辐射归一化研究[D]. 北京:中国科学院大学, 2018.

[ Bai Y. Radiometric normalization of remote sensing image based on kernel canonical correlation analysis[D]. Beijing: University of Chinese Academy of Scienes, 2018. ]

[26]
Nielsen AA. The regularized iteratively reweighted MAD method for change detection in multi- and hyperspectral data[J]. IEEE Transactions on Image Processing, 2007,16(2):463-478.

DOI

[27]
Canty M J, Nielsen A A. Automatic radiometric normalization of multitemporal satellite imagery with the iteratively re-weighted MAD transformation[J]. Remote Sensing of Environment, 2008,112(3):1025-1036.

DOI

[28]
Tucker C J. Red and photographic infrared linear combinations for monitoring vegetation[J]. Remote Sensing of Environment, 1979,8(2):127-150.

DOI

[29]
McFeeters S K. The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features[J]. International Journal of Remote Sensing, 1996,17(7):1425-1432.

DOI

[30]
Li Z X, Shi W Z, Zhang H, et al. Change detection based on gabor wavelet features for very high resolution remote sensing images[J]. IEEE Geoscience and Remote Sensing Letters, 2017,14(5):783-787.

DOI

[31]
Li H, Celik T, Longbotham N, et al. Gabor feature based unsupervised change detection of multitemporal SAR images based on two-level clustering[J]. IEEE Geoscience and Remote Sensing Letters, 2015,12(12):2458-2462.

DOI

[32]
陈小光, 封举富. Gabor滤波器的快速实现[J]. 自动化学报, 2007,33(5):456-461.

[ Chen X G, Feng J F. Fast Gabor filtering[J]. Acta Automatica Sinica, 2007,33(5):456-461. ]

[33]
Boyer R S, Moore J S. Automated reasoning: Essays in honor of Woody Bledsoe[M]. Dordrecht: Springer, 1991.

[34]
Otsu N. A threshold selection method from gray-level histograms[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1979,9(1):62-66.

DOI

[35]
Eckhardt D W, Verdin J P, Lyford G R. Automated update of an irrigated lands GIS using SPOT HRV imagery[J]. Photogrammetric Engineering & Remote Sensing, 1990,56(11):1515-1522.

Options
文章导航

/

〈 〉