本发明涉及遥感影像技术领域，具体涉及一种多源遥感影像辐射归一化方法。
背景技术：
遥感影像获取受传感器本身、光照、大气、地形等因素的影响，导致不同影像上相同地物的光谱特征存在很大差异。因此，在利用多源或多时相遥感影像进行变化检测或地物信息提取之前，需要对影像进行辐射归一化处理，控制和减少由于光照条件、大气效应、传感器响应等差异造成的地表景观的“伪变化”，保留真实的地表变化信息。
辐射归一化分为绝对辐射归一化和相对辐射归一化2种，后者由于不需要大气同步观测资料和计算相对简便而得到广泛应用。现有的相对辐射归一化方法大体上可分为2类：基于分布的相对归一化方法和基于像元对的相对辐射归一化方法。一般来说，基于分布的相对归一化如直方图匹配、平均值-标准偏差归一化法等，通过对影像的线性拉伸使2个影像的灰度值具有相似的灰度分布，具有计算简单的优点，但该方法容易造成原始光谱特征的扭曲，不利于后续的应用。基于像元对的方法从2个影像的重叠区域内选取伪不变特征点(PseudoInvariantFeatures，PIFs)，以PIFs的灰度变化作为辐射变化量度，建立影像间灰度的回归关系对目标影像进行处理，该类方法可以得到2个影像间较为准确的灰度映射关系，因而在准确选取PIFs的前提下可以得到更好的校正效果。目前，众多学者对基于PIFs的相对辐射归一化方法展开了研究，如自动散点控制回归、改进的自动散点控制回归、多元变化检测变换法以及迭代重新加权多元变化检测变换法、迭代加权最小二乘回归法等，且这些方法在一定条件下可取得较好的归一化效果。
然而，目前的辐射归一化方法研究和应用中，大多采用相同传感器的不同时相遥感影像，对多源、多时相传感器数据研究较少。即使对不同传感器数据进行辐射归一化时，也很少考虑传感器自身的辐射差异，或将这种差异性与光照、大气条件等外界因子看作是一个综合因素，与影像灰度值存在线性函数关系，对所有地类造成同等的影响。对于定量化遥感分析和应用而言，这种忽略不同传感器间的辐射差异或对其做全局线性假设有可能带来较大误差。吴荣华等研究了光谱响应函数对高精度交叉定标的影响，研究表明光谱响应差异是影响定标结果的重要因素，利用多源遥感数据做动态变化分析时必须考虑不同传感器光谱响应的影响，结果还表明，不同传感器光谱响应差异的影响还依赖于下垫面的类型。
随着国产遥感数据源的不断丰富及遥感应用的拓展，尤其是基于遥感时间序列分析应用的广泛开展，对多源、多时相遥感数据的协同利用以及遥感信息的定量化提取提出了日益迫切的需求。对于具有高观测频度的多源时间序列影像而言，在考虑多源影像间的光谱、辐射和几何分辨率差异的同时，有效提高辐射归一化的精度和自动化水平是满足其应用需求的关键所在；目前的辐射归一化方法大多集中于对同源和有限时相的数据进行研究，且自动化程度较低，鲜有针对多源时序数据辐射归一化的方法研究。鉴此，提出了基于分类的传感器辐射校正与基于NDVI差值直方图和类别约束相结合的相对辐射归一化方法，以期达到对多源时序影像的高精度、半自动化处理，为多源时间序列影像的协同利用提供方法借鉴。
技术实现要素：
针对以上问题，本发明提供了一种多源遥感影像辐射归一化方法，对传感器间的辐射偏差进行有效纠正，并在整体上获得比传统方法更好的辐射归一化精度；同时，多源时序影像的辐射校正结果也表明，方法能够有效地消除时序影像间的辐射特征波动，使植被等地类的季相变化信息得到更准确地表达，为多源时序影像的协同利用提供了借鉴方法，可以有效解决背景技术中的问题。
为了实现上述目的，本发明采用的技术方案如下：一种多源遥感影像辐射归一化方法，该方法将多源遥感影像的相对辐射归一化分为传感器辐射校正与针对光照等外部因素的辐射归一化两个过程，包括如下步骤：
S1、基于晴空影像，采用分类回归的方式获取传感器辐射校正系数；
S2、利用样本传递再分类的方法实现多源影像的半自动分类和传感器辐射偏差校正；
S3、基于NDVI差值直方图和类别约束的PIFs自动选取方法，实现影像的相对辐射归一化。
根据上述技术方案，所述步骤S1中，传感器光谱归一化系数获取方法步骤为：
(1)对影像进行辐射定标，将DN值转换为辐亮度，使不同影像像元值具有相同的量纲水平，消除传感器间的量化级数差异对拟合精度的影响；
(2)对重叠区影像进行植被、居民地、裸地、水体的大类的监督分类，并分别在各地类中采用样本抽选的方法选取足够数量的样本点样点大小视像元相对大小而定：像元大小相同时，以单个像元值为抽样值；像元大小不一致时，样点取以像元大小的最小公倍数为半径的圆，以圆内的像元均值为抽样值，通过样本抽选可减轻传感器间因分辨率差异带来的尺度效应误差；
(3)最后，根据获取的样本点集，针对两个影像中的不同波段和类别建立线性回归方程，求取回归系数即光谱归一化系数。
根据上述技术方案，所述步骤S2中影像分类与传感器辐射校正方法步骤如下：
(1)对影像数据集进行辐射定标，使影像像元值的量级与光谱归一化参数匹配，定标结果作为后继处理的数据基础；
(2)从定标影像序列中选取数据质量最优的影像作为参考影像，其余则为待纠正影像，采用步骤S1中的分类体系对参考影像进行样本选取与最大似然分类，获得参考影像的分类结果；
(3)在当前期影像及其分类结果基础上，对下一期待纠正影像进行样本筛选与样本纯化；
(4)利用样本纯化后得到的全类别样本对待纠正影像进行最大似然分类，并结合传感器类型，对影像中各类别对应的像元集进行传感器光谱归一化校正。
根据上述技术方案，所述步骤(3)中样本筛选以当前期影像及其分类结果为基础，从影像重叠区中计算和获取各地类具有代表性的像元空间位置集，作为下一期待纠正影像的候选样本空间位置，其具体步骤如下：
(1)获取当前影像及待纠正影像的重叠区范围，在重叠区内以分类图像作为类别范围约束，对当前期影像逐类别计算平均光谱向量作为该类在n维光谱空间的类别中心；
(2)以欧式距离为度量，计算每一类中所有像元光谱向量到类别中心的距离，并通过排序算像元光谱向量到类别中心的距离，并通过排序算法获取距离类别中心最近的m个像元作为该类的典型像元样本；
(3)将所有类别像元样本对应的空间位置集传递给待纠正影像，作为候选样本的空间分布。
根据上述技术方案，所述步骤(3)中样本纯化是针对待纠正影像进行样本特征的重新计算，剔除候选样本中由于地物类型变化及云影遮盖造成的噪声像元，提高样本纯度，其具体步骤如下：
(1)结合待纠正影像逐类读取候选样本位置的像元光谱向量，在此基础上求算每一类的平均光谱向量；
(2)逐类求解类内各像元的光谱向量与平均光谱向量的方差，并对类内的所有像元方差求方差平均值；
(3)针对每一地类取阈值，将各类中方差大于阈值的像元视为异质像元予以剔除，最终获得所有地类的纯样本像元集合。
根据上述技术方案，所述步骤S3基于NDVI差值直方图和类别约束的PIFs自动选取方法，构建待纠正影像与参考影像中各波段的线性回归方程，实现对待纠正影像的辐射归一化校正，NDVI及其差值的计算如式(1)、(2)所示：
ΔNDVI＝NDVIr-NDVIt (2)
式中：R_nir和R_red分别为影像近红外与红外波段反射率值；NDVIr和NDVIt分别为参考影像和待纠正影像的NDVI图像；
辐射稳定点位于ΔNDVI直方图的均值μ附近，而受噪声干扰的不稳定点位于分布图两侧。将位于μ±cσ范围内的点作为辐射稳定的PIFs，其中，σ为ΔNDVI的标准差，c为确定稳定点区间的常量，取值1，以PIFs为样本点，针对每一波段建立如式(3)的线性回归方程，根据最小二乘原理解算出每一波段的最优系数ki、bi，并对待纠正影像进行线性回归校正：
pri＝ki×pti+bi(i＝1，…，n) (3)
式中：pri和pti分别表示参考影像和待纠正影像的第i波段；ki和bi为拟合系数；n为波段数。
本发明的有益效果：
本发明提出了基于分类的传感器辐射校正与基于NDVI差值直方图和类别约束的PIFs自动选取相结合的相对辐射归一化方法，并利用多源、多时相影像进行了验证，结果表明该方法对传感器自身的辐射差异以及因光照等外界因素造成的差异具有良好的校正效果，可为多源时间序列影像的协同利用提供有益借鉴；还具有以下优点：
(1)采用样本传递的方式，针对待纠正影像本身进行样本特征计算和优化，有利于提高时序影像的分类精度和自动化水平；
(2)基于NDVI差值直方图与地类约束结合的PIFs选取和优化方法，避免了常规PIFs自动选取方法中可能包含植被、水体等辐射特征易随时间发生变化的地物类型，提高了回归方程的拟合精确度。
附图说明
图1为一种多源遥感影像辐射归一化方法的流程图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
实施例：
如图1所示，一种多源遥感影像辐射归一化方法，该方法将多源遥感影像的相对辐射归一化分为传感器辐射校正与针对光照等外部因素的辐射归一化两个过程；
步骤S1、在遥感分类中，样本质量与分类精度密切相关，为保证样本信息的高“保真度”，一般从待分类影像中进行样本选取。当同一样本应用于多景影像时，由于分辨率、光照和时相的差异容易造成同物异谱、同谱异物的现象，给分类带来较大不确定性。结合时序影像对同一区域连续观测的特点；提出了基于样本传递的半自动分类方法，将前期分类获得的类别空间位置作为下一期分类的候选样本位置，并针对新一期影像重新进行样本特征的计算、优化及再分类。因此，对于多源时序影像，只需对参考影像进行一次人工选样的监督分类便可实现全数据集的自动分类过程。基于样本传递的影像分类与传感器辐射校正过程如下：
(1)对影像进行辐射定标，将DN值转换为辐亮度，使不同影像像元值具有相同的量纲水平，消除传感器间的量化级数差异对拟合精度的影响；
(2)对重叠区影像进行植被、居民地、裸地、水体的大类的监督分类，并分别在各地类中采用样本抽选的方法选取足够数量的样本点样点大小视像元相对大小而定：像元大小相同时，以单个像元值为抽样值；像元大小不一致时，样点取以像元大小的最小公倍数为半径的圆，以圆内的像元均值为抽样值，通过样本抽选可减轻传感器间因分辨率差异带来的尺度效应误差；
(3)最后，根据获取的样本点集，针对两个影像中的不同波段和类别建立线性回归方程，求取回归系数即光谱归一化系数。
通过样本传递的自动分类，有效地提高了多源时序影像分类及整体辐射归一化过程的处理效率。其中，样本筛选和样本纯化结合了前、后期影像的各自特点进行样本的精化，对降低计算复杂度和提高再分类精度起关键作用。
步骤S2、利用样本传递再分类的方法实现多源影像的半自动分类和传感器辐射偏差校正；
影像分类与传感器辐射校正方法步骤如下：
(1)对影像数据集进行辐射定标，使影像像元值的量级与光谱归一化参数匹配，定标结果作为后继处理的数据基础；
(2)从定标影像序列中选取数据质量最优的影像作为参考影像，其余则为待纠正影像，采用步骤S1中的分类体系对参考影像进行样本选取与最大似然分类，获得参考影像的分类结果；
(3)在当前期影像及其分类结果基础上，对下一期待纠正影像进行样本筛选与样本纯化；
样本筛选以当前期影像及其分类结果为基础，从影像重叠区中计算和获取各地类具有代表性的像元空间位置集，作为下一期待纠正影像的候选样本空间位置，一般而言，同类地物具有相似的光谱特征，在n个影像波段构成的n维光谱空间呈集中分布，离类别中心越近的像元具有更高的代表性，其具体步骤如下：
(1)获取当前影像及待纠正影像的重叠区范围，在重叠区内以分类图像作为类别范围约束，对当前期影像逐类别计算平均光谱向量作为该类在n维光谱空间的类别中心；
(2)以欧式距离为度量，计算每一类中所有像元光谱向量到类别中心的距离，并通过排序算像元光谱向量到类别中心的距离，并通过排序算法获取距离类别中心最近的m个像元作为该类的典型像元样本；
(3)将所有类别像元样本对应的空间位置集传递给待纠正影像，作为候选样本的空间分布。
根据上述技术方案，所述步骤(3)中样本纯化是针对待纠正影像进行样本特征的重新计算，剔除候选样本中由于地物类型变化及云影遮盖造成的噪声像元，提高样本纯度，其具体步骤如下：
(1)结合待纠正影像逐类读取候选样本位置的像元光谱向量，在此基础上求算每一类的平均光谱向量；
(2)逐类求解类内各像元的光谱向量与平均光谱向量的方差，并对类内的所有像元方差求方差平均值；
(3)针对每一地类取阈值，将各类中方差大于阈值的像元视为异质像元予以剔除，最终获得所有地类的纯样本像元集合。
(4)利用样本纯化后得到的全类别样本对待纠正影像进行最大似然分类，并结合传感器类型，对影像中各类别对应的像元集进行传感器光谱归一化校正。
S3、基于NDVI差值直方图和类别约束的PIFs自动选取方法，实现影像的相对辐射归一化。基于NDVI差值直方图和类别约束的PIFs自动选取方法，构建待纠正影像与参考影像中各波段的线性回归方程，实现对待纠正影像的辐射归一化校正，NDVI及其差值的计算如式(1)、(2)所示：
ΔNDVI＝NDVIr-NDVIt (2)
式中：R_nir和R_red分别为影像近红外与红外波段反射率值；NDVIr和NDVIt分别为参考影像和待纠正影像的NDVI图像；
辐射稳定点位于ΔNDVI直方图的均值μ附近，而受噪声干扰的不稳定点位于分布图两侧。将位于μ±cσ范围内的点作为辐射稳定的PIFs，其中，σ为ΔNDVI的标准差，c为确定稳定点区间的常量，取值1，以PIFs为样本点，针对每一波段建立如式(3)的线性回归方程，根据最小二乘原理解算出每一波段的最优系数ki、bi，并对待纠正影像进行线性回归校正：
pri＝ki×pti+bi(i＝1，…，n) (3)
式中：pri和pti分别表示参考影像和待纠正影像的第i波段；ki和bi为拟合系数；n为波段数。
基于上述，本发明的优点在于，本发明提出了基于分类的传感器辐射校正与基于NDVI差值直方图和类别约束的PIFs自动选取相结合的相对辐射归一化方法，并利用多源、多时相影像进行了验证，结果表明该方法对传感器自身的辐射差异以及因光照等外界因素造成的差异具有良好的校正效果，可为多源时间序列影像的协同利用提供有益借鉴；还具有以下优点：采用样本传递的方式，针对待纠正影像本身进行样本特征计算和优化，有利于提高时序影像的分类精度和自动化水平；基于NDVI差值直方图与地类约束结合的PIFs选取和优化方法，避免了常规PIFs自动选取方法中可能包含植被、水体等辐射特征易随时间发生变化的地物类型，提高了回归方程的拟合精确度。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。