CN102565778A

CN102565778A - 一种自动提取伪不变特征的遥感图像相对辐射校正方法

Info

Publication number: CN102565778A
Application number: CN2011104108929A
Authority: CN
Inventors: 胡昌苗; 唐娉; 唐亮
Original assignee: Institute of Remote Sensing Applications of CAS
Current assignee: Institute of Remote Sensing Applications of CAS
Priority date: 2011-12-12
Filing date: 2011-12-12
Publication date: 2012-07-11
Anticipated expiration: 2031-12-12
Also published as: CN102565778B

Abstract

本发明针对波长在可见光与近红外范围内的多光谱航天遥感数据提供一种自动提取伪不变特征点的遥感图像相对辐射校正技术。该技术主要包括以下步骤：针对DN值的目标图像，找到匹配的地表反射率的参考图像；对两图像进行预处理；利用典型相关分析技术寻找典型相关的点集；在典型相关点集中筛选伪不变特征点；用伪不变特征点拟合线性关系，并利用该关系对目标图像进行线性相对辐射校正处理。该算法直接将目标图像DN值转换到地表反射率，具有处理流程简单、无需人工交互、处理速度快的特点，算法稳定，可用于相同或不同传感器之间的相对辐射校正。尤其适用于辅助信息缺失或定标精度不高、不适合绝对辐射校正的遥感数据的辐射处理。

Description

一种自动提取伪不变特征的遥感图像相对辐射校正方法

技术领域

本发明涉及遥感图像辐射校正技术，具体的说，涉及一种自动提取伪不变特征的遥感图像相对辐射校正技术。

背景技术

遥感图像的辐射校正一直是遥感数据处理的主要难点之一。近些年来，随着定量遥感技术迅速发展，特别是利用多传感器、多时相遥感数据进行土地利用和土地覆盖变化监测、全球资源环境分析、气候变化监测等的广泛应用，使得遥感图像辐射校正方法的研究越来越受到重视。目前遥感图像的辐射校正方法可以分为两类，绝对辐射校正(Absolute Radiometric Correction)和相对辐射校正(Relative RadiometricCorrection)。

绝对辐射校正是基于大气辐射传输物理过程的，将遥感图像的DN(digital number)值转换为地表反射率或地表反射辐亮度的方法。研究者们发展了很多算法来消除太阳辐照度、大气和传感器差异带来的影响，例如：6S、MODTRAN等辐射传输模型。该类方法较复杂，需要利用一系列参数(例如，卫星过境时的地物反射率、大气能见度、太阳天顶角和卫星传感器标定参数等)。获得这些参数的代价昂贵，又不易实现，尤其对于历史存档数据，获取这些参数更加困难。

相对辐射校正是一种常用方法，它不需要遥感图像获取时的大气状况等参数，而是利用图像像元灰度值，建立多时相遥感图像各个波段之间的变换关系，对遥感图像进行归一化处理。根据变换关系的不同，相对辐射校正又可分为线性相对辐射校正与非线性相对辐射校正两类。线性相对辐射校正对地物的光谱特性保持较好，较为常用，这类方法假设：目标图像的像元与参考图像的像元是线性相关的。这个假设意味着在两图像时间间隔内采样点的光谱反射率特征没有发生变化。为使线性相对辐射校正中线性关系成立，一般假设目标图像与参考图像像元值有相同的物理单位，比如都是DN值或者反射率，这点限制了线性相对辐射校正的适用领域，使之无法用于不同物理量之间的转换。相对辐射校正可以提高数据的光谱一致性，而对物理量绝对精度缺乏验证。

基于伪不变特征(Pseudo-Invariant Features，PIFs)的线性相对辐射辐射校正是近年来的研究热点。PIFs是指遥感图像中那些光谱特性在短时间内变化不大的点，比如不透水层、路面、屋顶和停车场等。PIFs物理含义清晰，是线性辐射校正建立线性关系的依据。目前的PIFs提取算法有很多，存在的主要问题包括人机交互性强，数据通用性差，鲁棒性不高等。实际工程应用需要一个自动的、稳定的PIFs提取算法。

发明内容

本发明的目的是针对遥感应用，提供一种自动的航天遥感图像相对辐射校正技术，特别是对于那些大气参数及传感器辅助信息缺失或传感器定标精度不高的遥感图像的辐射校正。本技术基于典型相关分析自动提取PIFs，并将相对辐射校正技术适用性扩展到不同物理量纲的数据。

本发明的基本思路为：针对波长在可见光与近红外范围内的多光谱航天遥感图像，对于DN值的目标图像，在满足传感器观测角度在正负十度之间，能见度大于30km的情况下，选择一个与其具有重叠区域的、已经转换到地表反射率的参考图像。先对两图像进行预处理，选出重叠区域对应的象素点，继而在这些象素点中采用典型相关分析的方法自动寻找其中典型相关点集，然后从这些典型相关点集中进一步筛选出PIFs，用最终筛选出的PIFs拟合出一个线性关系，最后用该线性关系对目标图像进行相对辐射校正处理。

所述的目标图像与参考图像可以来自同种传感器，比如Landsat TM/ETM+，也可以来自不同传感器，比如目标图像为Landsat TM/ETM+，参考图像为MODIS Terra MOD09GA产品。

本发明的技术方案提供的自动提取伪不变特征的遥感图像相对辐射校正方法，其特征在于包括以下实施步骤：

A根据目标图像信息，确定与其匹配的参考图像；

B对目标图像与参考图像进行预处理，包括坐标重投影、几何配准、计算地理重叠区域、各波段象素点对应，并使用掩膜技术排除极值点、过饱和点、云、云下阴影及水体的象素点；

C使用典型相关分析从获取的点集中提取典型相关点集；

D从提取的典型相关点集中筛选PIFs；

E使用提取的PIFs计算各波段的线性关系，用这些线性关系对目标图像各波段象素点进行线性变换，将DN值直接转换到地表反射率。

上述实施步骤的特征在于：

步骤A中所述根据目标图像信息，主要包括目标图像的传感器类型、成像大致时间、坐标范围与投影等信息，这些信息用于确定需要的参考图像。步骤A中确定与其匹配的参考图像，主要包括：要求参考图像为地表反射率数据，与目标图像的成像时间相近(可以是不同年份的，但月份要求在一个月以内，成像时刻在1小时以内)，具有相近的成像几何，传感器观测角差异在10°以内，同时要求目标图像能见度大于30km。

步骤B中所述的预处理主要包括坐标重投影、几何配准、计算重叠区域及波段对应等处理。对于不同传感器数据，考虑分辨率一致性处理。步骤B中使用掩膜技术，将一些过亮及过暗的点及一些饱和的点去除，这样可以去除图像中的云与阴影，或者针对某一具体传感器数据，设计专门的云与阴影检测算法。

步骤C中所述典型相关分析，主要包括具有仿射不变性的特征的典型相关分析算法，以便从目标图像DN值与参考图像地表反射率之间选出具有近似线性关系的典型相关点集。比如多元变化检测(MultivariateAlteration Detection，MAD)变换、迭代的多元变化检测算法(iteratively re-weighted modification of the MADtransformation，IR-MAD)。

步骤D中从提取的典型相关点集中筛选PIFs，筛选的步骤包括通过NDVI排除植被象素点，选择传感器观测角差异小的点，保留同质区域内的象素点，对于不同传感器数据，处理的步骤及采用的算法会有差异。

步骤E中计算线性关系，常用最小二乘法，也可使用正交回归法。步骤E的其它处理与一般的线性相对辐射校正相同。

本发明与现有技术相比有如下特点：该算法直接将目标图像DN值转换到地表反射率，具有处理流程简单、无需人工交互、运算速度快的特点，而且在算法的稳定性以及适用性上都具有优势。对于定标精度高、辅助信息完整的遥感图像建议采用传统的处理流程，而对于无法使用传统方式处理的遥感图像，本发明提供了一种简单的解决方案。该技术使用典型相关分析的办法自动提取典型相关点集，并从这些点中自动提取PIFs，自动化程度高、通用性强，不仅适用于同种传感器数据，对于满足要求的不同种传感器数据同样适用。

附图说明：

图1是定量化的遥感图像相对辐射校正流程图

图2是大气辐射传输示意图

图3是Landsat TM/ETM+数据与MODIS Terra MOD09GA产品象素对应示意图

图4是Landsat TM/ETM+数据与MODIS Terra MOD09GA产品拓扑关系示意图

图5是PIFs筛选流程图

图6是传统的绝对辐射校正处理流程示意图

图7是本发明的处理流程示意图

图8是Landsat TM/ETM+标准分幅数据示意图

具体实施方式：

本技术的思想是使用线性相对辐射校正直接将DN值转换到地表反射率，其核心假设为：在一定条件下PIFs在目标图像的DN值与在参考图像的地表反射率之间存在近似线性关系。该假设基于如下推导：

一般遥感图像辐射定标过程使用的是线性公式，比如(1)：

L_b＝Gain*DN_b+Bias (1)

其中Gain与Bias是定标系数。公式(1)将DN值转换到表观辐射亮度L_b。

由表观辐射亮度L_b转换到表观反射率ρ，计算采用公式(2)：

ρ = \frac{π * L_{b} * d^{2}}{{ESUN}_{λ} * \cos θ} - - - (2)

其中d为日地距离，ESUN_λ为太阳常数，θ为太阳天顶角。这些变量对于成像时间近似、地理覆盖范围不大的单景遥感图像而言近似是常量，此时L_b与ρ之间满足近似的线性关系。

由表观反射率ρ_j(λ)转换到地表反射率

是基于遥感辐射传输方程，见公式(3)：

ρ_{j} (λ) = A (λ) + \frac{B (λ)}{1 - S (λ) < ρ (λ) >} ρ_{j}^{o} (λ) + \frac{C (λ)}{1 - S (λ) < ρ (λ) >} < ρ (λ) > - - - (3)

其中λ为中心波长，<ρ(λ)>为大气平均反射率，A(λ)，B(λ)，C(λ)与S(λ)是传感器获取的不同组分，其在大气中的传输路径见图2。A(λ)为光线射入大气，经大气的传输与散射后进入传感器的部分；B(λ)为太阳→地表→传感器的路径传输；C(λ)为地表“临近效应”；S(λ)为大气半球反照率。

公式(3)在能见度相对高并且单景图像覆盖地理范围不大时，S(λ)<ρ(λ)>值很小，并且<ρ(λ)>值近似为常量，此时公式(3)退化为线性形式，见公式(4)：

ρ_{j} (λ) = A (λ) + B (λ) ρ_{j}^{o} (λ) + C (λ) < ρ (λ) > - - - (4)

综合公式(1)、(2)、(4)，三个线性关系的叠加还是线性的，故由DN值到地表反射率可以通过一个线性公式直接转换。

该线性关系成立的条件有：目标图与参考图成像几何近似；目标图能见度高(大于30km)，且气溶胶分布均匀；单景目标图覆盖地理范围不大，观测角度变化范围不大(比如-10～+10度)。这些条件相对于传统的辐射处理方式要宽松很多，不需要精确的定标系数、成像时的几何及大气参数，结果精度主要取决于参考图像。该线性关系对于大多数处于可见光与近红外范围的多光谱航空遥感图像都是适用的。

采用本发明实现遥感图像相对辐射校正的实施例如图1所示，现结合附图对其进行描述。

处理单元111寻找参考图像，根据目标图像的信息，在已有数据中寻找最佳匹配的参考图像。寻找的依据为：参考图像为地表反射率数据；参考图像与目标图像有重叠区域；两图像成像时间与几何近似。

处理单元112对目标图像与检索得到的参考图像进行数据预处理，预处理主要包括：坐标重投影、几何配准、计算地理重叠区域、各波段象素点对应、掩膜过滤异常点。

坐标重投影与几何配准属通常的遥感图像处理步骤。

计算地理重叠区域。对于单个目标图像，有可能找到多景参考图像，此时就需要进行重叠区域内数据的裁切与拼接整合。比如对于单景TM/ETM+，对应到MODIS的参考图像上可能会有4种情况，分别为单景TM/ETM+数据对应一景MODIS数据，单景TM/ETM+数据对应两景MODIS数据，单景TM/ETM+数据对应三景MODIS数据，单景TM/ETM+数据对应四景的MODIS数据，见图4。

各波段象素点对应。对于目标图像与参考图像为同种传感器数据，分辨率相同，当坐标重投影、几何配准精度足够高时，象素点本身就是对应的。对于不同传感器的遥感图像，分辨率差异超过3倍以上时，就需要插值。比如将Landsat TM/ETM+，30米分辨率的目标图像对应到MODIS Terra MOD09GA，500米分辨率的参考图像上，如图3所示。象素对应计算采用公式(5)：

D_{t} = Σ_{i = 1}^{n} (p_{t}^{i} {* S}_{t}^{i}) - - - (5)

公式(5)描述了第t个参考图像象素点覆盖范围内包含了n个目标图像的象素，目标图像象素值与所占面积比率乘积的和即为对应象素点的值。

掩膜过滤异常点。该步先通过统计目标图像与参考图像重叠区内的直方图，过滤掉过亮和过暗的象素点。然后检测参考图像中的云、云下阴影、水体并将其掩膜掉。具体的检测算法依据不同的传感器数据而不同。

处理单元113对掩膜过滤后得到的数据点集进行典型相关分析，进一步从这些点集中选出典型相关点集。所选用的典型相关分析算法应当具有不受总体的大气状况和传感器标定导致的线性变换的影响的特性。

典型相关分析算法为了遮蔽两个不同时相图像中的变化像素，首先形成两幅图像N个波段内像素值的线性组合。用随机向量X和Y分别表示目标图与参考图重叠区内筛选出的像素值。根据变换公式(6)：

U = a^{T} X = a_{1} X_{1} + a_{2} X_{2} + \cdot \cdot \cdot + a_{N} X_{N}

(6)

V＝b^TY＝b₁Y₁+b₂Y₂+…+b_NY_N

其中a_i与b_i为变换系数，为了最小化U与V之间的正相关，在服从约束Var(U)＝Var(V)＝1的前提下，使得：

Var(U-V)＝Var(U)+Var(V)-2cov(U，V)＝2(1-corr(U，V))→Maximum

最小化正相关系数corr(U，V)是一个标准的统计过程，即所谓的广义特征值问题。求出的Var(U-V)各个分量相互正交，并且是线性变换的不变量，这意味着这种变换对测量尺度和测量装置的增益与偏移不敏感，因而对图像数据X和Y测量尺度的一致性没有要求，可以为不同的物理量纲。

典型相关分析方法实质上是把图像X与Y之间的差异总信息分配到互不相关的一组变量上，以达到最大限度保持这一差异的总信息量不改变的情况下，检测出图像X与Y的差异，从差异较小的部分中便可以选出典型相关点集。

MAD变换是一种满足线性不变特性的典型相关分析算法，IR-MAD变换进一步提高了MAD算法的精度与稳定性。

目标图像DN值与参考图像地表反射率之间虽然存在近似的线性关系，但是两者在数值上的差异还是巨大的。DN值一般为8位Byte型保存，数值0～255，而地表反射率为0～1的值。为了减少数值上的差异，在进行典型相关分析前，我们先将DN值转换到0～1范围，采用的办法是先乘以400再除以10000。

处理单元114对典型相关分析后得到的数据点集进行进一步的筛选，筛选的结果作为最终的PIFs。筛选的细节流程见图5。筛选主要是植被象元的筛选、传感器观测角度的筛选与匀质区域象素点筛选。植被由于受时空因素影响，光谱变化剧烈，所以一般不将其选为PIFs，但对于重叠区几乎全为植被覆盖的情况，如果目标图与参考图时间与成像几何非常近似，植被光谱变化不大的情况下也是可以选为PIFs的。可以通过计算植被指数(NDVI)，并设定一个域值进行非植被象元的筛选。分类数据也可用于非植被象元的筛选。成像几何中传感器观测角度的影响最为显著，如果最后筛选的结果仍然包含大量的点，还可以依据观测角度进一步筛选。筛选的依据是看每个象素点在目标图与参考图上对应的观测角度是否相近，优先选择角度最相近的点，当最终选取的点数达到一定比例时停止筛选。匀质区域象素点的筛选主要是针对不同传感器数据，减小由点扩散函数的不同与地表“临近效应”所带来的影响。具体的筛选方法是使用一个3*3的模板，计算每个象素所在3*3模板内所有象素的最大值与最小值的差，如果差值小于某个域值，则认为当前象素点为匀质区域的点。域值根据具体数据与经验确定。

处理单元115使用最终筛选出的PIFs点集对各波段分别拟合出一个整体的线性关系。可以使用基于最小二乘原理的拟合、也可用正交回归的方式。

处理单元116线性相对辐射校正，使用单元115计算的线性关系对目标图像各波段进行变换。该步与传统的线性相对辐射校正处理的方式相同。

本发明实现的遥感数据DN值到地表反射率的处理流程比传统的绝对辐射校正流程要简化。传统的绝对辐射校正流程如图6，需要辐射定标与绝对大气校正两步，并且这两步的处理需要大量的参数，包括定标系数、传感器参数与大气参数等，并且基于辐射传输理论的绝对大气校正运算复杂，处理速度慢。本发明的处理流程如图7，处理流程只有相对辐射校正一步，输入的参数也少，无需定标参数及精确的大气参数。并且本发明的计算过程简单，处理速度快。

本发明对Landsat TM/ETM+数据的适用性分析如下：

Landsat轨道设计为与太阳同步的近极地圆形轨道，确保了北半球中纬度地区获得中等太阳高度角(25°-30°)的上午成像，卫星以同一地方时、同一方向通过同一地点，保证了观测条件的基本一致。分辨率为30米，传感器观测角度为±5°，近似垂直观测，幅宽为185km，单景数据地表近似平面。轨道是严格回归的，数据以WSR-2(USGS Worldwide Reference System-2)体系按照Path与Row索引进行标准分幅后，相邻的数据之间存在固定范围的重叠区域，如图8所示，利于相邻数据的比较。Landsat TM/ETM+数据的成像几何近似满足本发明对数据的要求，同时Landsat TM/ETM+数据的绝对辐射校正技术目前相对成熟，转换到地表反射率的参考数据容易获取，在选取参考数据时重点关注成像时间的近似即可。由于Landsat的重访周期为16天，Landsat TM/ETM+数据的参考图像的大气状况及地表情况与目标图像存在差异，进而影响最终结果的精度。

本发明对目标图像为Landsat TM/ETM+、参考图像为MODIS Terra MOD09GA产品的适用性分析如下：

在Terra卫星上的MODIS传感器数据在可见光与近红外的波段设置与Landsat TM/ETM+近似、并且Terra卫星的过境时间仅仅比Landsat晚30分钟，Terra卫星的重访周期为1天，所以对于每景LandsatTM/ETM+数据，容易找到与其对应的Terra数据，两数据成像时间差小于30分钟，大气状况与地表情况可以认为近似不变。MODIS数据的绝对辐射校正算法成熟，其地表反射率产品MOD09GA的精度被遥感领域所广泛认可。由于MODIS数据为500米分辨率，所以需要进行处理单元111中所述的各波段象素点对应处理。考虑到不同传感器点扩散函数以及地表临近效应的影响，需要进行处理单元114中所述的选择同质区域内的象素点。由于MODIS传感器观测角范围大(±50°)，需要进行处理单元114中所述的观测角度的筛选。

本发明的一个实例在PC平台上实现，目前支持的处理数据包括Landsat TM/ETM+之间的相对辐射校正；目标图像为Landsat TM/ETM+、参考图像为MODIS Terra MOD09GA产品的相对辐射校正。对于国产环境星CCD图像，由于成像范围过大，只能对整景图像中的局部区域进行处理，此时参考图像除了环境星CCD数据外，也可以选用Landsat TM/ETM+数据。在经过大量的数据处理与验证的基础上，确定了针对上述数据的相关参数，对上述数据的处理做到了无需人机交互的自动化。并且处理速度较传统的处理流程有较大的提高，单景Landsat TM/ETM+图像平均处理时间小于5分钟。算法鲁棒性较强，典型相关分析能较好地适应DN值与地表反射率之间的近似线性关系。经过大量数据的验证，该发明的处理结果精度较高，满足遥感应用对数据精度的要求。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。

Claims

1.一种自动提取伪不变特征的遥感图像相对辐射校正方法，该方法针对波长在可见光到近红外范围内的多光谱航天遥感数据，直接将目标图像DN(digital number)值转换到地表反射率，具体包括以下步骤：

步骤1，根据目标图像，确定匹配的参考图像，要求参考图像为地表反射率数据，与目标图像的成像时间近似(可以是不同年的，但月份要求在一个月以内，成像时刻在1小时以内)，传感器观测角差异在10°以内，同时要求目标图像能见度大于30km；

步骤2，对目标图像与参考图像进行预处理，包括重投影，几何配准，计算地理重叠区域，各波段象素点对应，排除极值点、过饱和点、云、云下阴影及水体所在象素点；

步骤3，使用典型相关分析从步骤2获取的点集中提取典型相关点集；

步骤4，从步骤3提取的典型相关点集中筛选伪不变特征点集(Pseudo-Invariant Features，PIFs)，筛选的步骤包括通过植被指数(NDVI)排除植被象素点，按照传感器观测角差异大小对剩余点排序，保留角度相对一致的点，对于不同传感器数据，保留同质区域内的象素点；

步骤5，使用步骤4提取的PIFs计算各波段线性关系，并使用这些线性关系对目标图像各波段象素点进行线性变换，将DN值直接转换到地表反射率。

2.根据权利要求1中所述的技术方案，其特征在于：

用线性关系将DN值转换到地表反射率，前提是DN值转换到表观辐射亮度、再转换到表观反射率、最后转换到地表反射率，三步转换的过程线性关系都近似成立，步骤1中所述的“要求参考图像为地表反射率数据，与目标图像的成像时间近似(可以是不同年的，但月份要求在一个月以内，成像时刻在1小时以内)，传感器观测角差异在10°以内，同时要求目标图像能见度大于30km”仅是保证线性关系成立的一种情况，根据具体数据及对结果精度的要求不同，对目标图像及参考图像的要求可以变化。

3.根据权利要求1中所述的技术方案，其特征在于：

步骤2中对数据的预处理中，对不同传感器数据、分辨率差异超过3倍时，采用求加权乘积再求和的方式完成每个象素点的对应，采用公式为

即第t个参考数据象素点覆盖范围内包含了n个目标图像的象素，目标图像象素值与所占面积比率乘积的和即为对应象素点的值。

4.根据权利要求1中所述的技术方案，其特征在于：

步骤2中排除的象素点包括极值点、过饱和点、云、云下阴影及水体所在象素点，对不同源的数据，确定这些象素点的算法也不同，对具体的确定异常点、云、云下阴影、水体的方法并不做限定。

5.根据权利要求1中所述的技术方案，其特征在于：

选用具有线性不变性的典型相关分析算法提取典型相关点集。

6.根据权利要求1中所述的技术方案，其特征在于：

步骤4从典型相关点集中提取PIFs的步骤包括排除植被象素点、筛选观测角度差异小的象素点、选择同质区域内的象素点三步，对不同源的数据，这三步的处理会有差异，对具体的排除植被象素点与筛选观测角度差异小的象素点的处理算法并不做限定。

7.根据权利要求1中所述的技术方案，其特征在于：

步骤4对于不同传感器数据，保留位于同质区域的象素点，方法是通过一个3*3的窗口，计算每个特征不变点所在3*3窗口内的所有象素的最大值与最小值的差，当差值小于某一域值时将其划分为同质点，不同传感器数据域值不同，域值的确定依据统计经验。

8.根据权利要求1中所述的技术方案，其特征在于：

所述处理流程不仅适用于同种传感器数据，还可以用于满足条件的不同种传感器数据。