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

技术领域

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

背景技术

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

绝对辐射校正是基于大气辐射传输物理过程的，将遥感图像的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值与参考图像地表反射率之间选出具有近似线性关系的典型相关点集。比如多元变化检测(Multivariate  Alteration Detection，MAD)变换、迭代的多元变化检测算法(iteratively re-weighted modification of the MAD  transformation，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)：

$\rho =\frac{\pi *L_b*d^2}{{\mathrm{ESUN}}_{\lambda }*\cos \theta }---\left(2\right)$

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

由表观反射率ρ_j(λ)转换到地表反射率是基于遥感辐射传输方程，见公式(3)：

${\rho }_j\left(\lambda \right)=A\left(\lambda \right)+\frac{B\left(\lambda \right)}{1-S\left(\lambda \right)<\rho \left(\lambda \right)>}{\rho }_j^o\left(\lambda \right)+\frac{C\left(\lambda \right)}{1-S\left(\lambda \right)<\rho \left(\lambda \right)>}<\rho \left(\lambda \right)>---\left(3\right)$

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

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

${\rho }_j\left(\lambda \right)=A\left(\lambda \right)+B\left(\lambda \right){\rho }_j^o\left(\lambda \right)+C\left(\lambda \right)<\rho \left(\lambda \right)>---\left(4\right)$

综合公式(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=\underset{i=1}{\overset{n}{\Sigma }}\left(p_t^i{*S}_t^i\right)---\left(5\right)$

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

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

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

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

$U=a^{T}X=a_1{X}_1+a_2{X}_2+···+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天，所以对于每景Landsat  TM/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值与地表反射率之间的近似线性关系。经过大量数据的验证，该发明的处理结果精度较 高，满足遥感应用对数据精度的要求。

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