用于大观测角传感器的交叉辐射定标方法及系统与流程
本发明涉及一种针对大观测角传感器的基于数据同化的交叉辐射定标方法及系统，属于遥感影像处理技术领域。
背景技术：
随着遥感应用的深入，现在国内外主流的遥感技术已从定性分析发展到了定量研究，而辐射定标是定量遥感的基础和前提。辐射定标是将传感器记录的无量纲的数字量化值(dn值)转换成具有实际物理意义的大气顶层辐射亮度或反射率。它的主要作用体现在：(1)通过动态监测，校正传感器的性能衰变，修正大气、光照条件和环境变化等对测量结果的影响，保证测量结果真实反映目标和背景的特征；(2)建立不同平台、不同传感器的辐射量基准，使不同时间和空间尺度下获得的多源观测数据可以相互对比、转换和融合。
卫星发射后，传感器性能随着太空环境的变化而变化，必须对其进行在轨辐射定标，在轨辐射定标分为星上定标和替代定标，星上定标主要依赖于星上定标设备，具有较高的定标精度，但随着时间推移，定标设备会发生衰减，而且部分传感器没有配备星上定标设备，因此必须探寻星上定标的有效补充，即进行替代定标。目前常用的替代定标方法有场地定标和交叉定标，如我国的环境系列、风云系列卫星等,使用的定标数据来自于在敦煌绝对福射校正场及青海湖校正场进行的场地定标实验。然而场地定标方法需要地面和卫星同步观测,容易受到大气状况、太阳角度等各种自然条件的限制，而且需要耗费大量的人力、财力、物力；更重要的是，它不能实现对历史数据的定标。因此，无需地面观测数据的交叉定标方法被提出并得到广泛研究和应用。交叉定标的基本过程是,通过参考传感器与待定标传感器同时观测同一目标,用参考传感器获取的信息对待定标传感器进行定标。与场地定标相比,交叉定标不需要投入大量的人力、仪器及资金等进行地面测量,对实验场的要求也低于场地定标，还可以对历史数据进行定标。
交叉定标过程中必须消除参考传感器和待定标传感器之间差异，这些差异主要体现在光谱响应、大气辐射传输以及地表二向性反射分布三个方面，其中大气辐射传输和地表二向性反射分布都与卫星的观测角度显著相关。卫星的观测角就是成像时卫星与目标地物的角度关系，有卫星天顶角和卫星方位角，卫星天顶角是卫星和目标的连线与天顶方向的夹角，如果卫星垂直拍摄地面那么天顶角就是0°。卫星方位角是卫星和目标的连线在地面的投影与当地子午线的夹角。卫星观测角度的差异主要体现在卫星天顶角，普通卫星(大多是太阳同步轨道)一般近垂直拍摄，卫星天顶角一般接近0°，比如landsat-8oli的卫星天顶角范围为±7°。大观测角卫星，比如地球静止轨道卫星gf-4pms的卫星天顶角的范围为-65°-70°，观测地点距离星下点越远，角度越大。普通太阳同步轨道卫星的卫星方位角浮动范围较小，大观测角卫星，如静止轨道卫星的卫星方位角浮动范围相对较大，例如在敦煌定标场，landsat-8oli的卫星方位角为98°，而gf-4pms的卫星方位角为161°。可以作为参考的定标精度较高的传感器绝大多数观测角度都很小，那么待定标传感器的观测角度越大，他们之间的大气辐射传输和地表二向性反射分布差异就会越显著。因此，对于大观测角的卫星，由于不考虑大气辐射传输或者地表二向性反射分布差异性，传统交叉定标方法无法得到准确的辐射定标系数。
目前本领域急需一种针对大观测角传感器的交叉定标方法，消除参考传感器与待定标传感器之间的大气辐射传输和地表二向性反射分布差异，从而得到准确的交叉定标结果。
技术实现要素：
针对上述问题，本发明提出了一种针对大观测角传感器的基于数据同化的交叉辐射定标技术方案，同时考虑两传感器之间光谱响应、大气辐射传输以及地表二向性反射分布函数(brdf)的差异，将交叉定标看作是一个最佳逼近的问题，无需通过复杂模型计算brdf，而是建立一个brdf调整系数，与定标系数一起作为待同化参数，利用同化算法通过循环迭代的方式找到最优的辐射定标系数和brdf调整系数。
本发明的技术方案提供一种用于大观测角传感器的交叉辐射定标方法，包括以下步骤，
步骤1，配准待定标影像和参考影像，对待定标影像和参考影像的每一组对应波段分别进行随机采样得到对应的定标点；
步骤2，针对每个波段，分别对每一对定标点进行如下的交叉定标操作，
步骤2.1，基于参考影像的定标点的dn值计算大气表观反射率，并通过大气校正得到地表反射率；
步骤2.2，对参考影像的地表反射率进行匹配找到最接近的光谱曲线，计算参考影像与待定标影像之间的光谱调整系数sbafi，下标i用于表示第i波段；
步骤2.3，建立brdf调整系数ci，所述brdf调整系数ci为待定标影像第i波段的地表二向性反射分布函数与参考影像第i波段的地表二向性反射分布函数的比值；
步骤2.4，初始化参数，包括设定待定标影像的辐射定标系数gaini和offseti的初始值、有效范围和变化步长，以及brdf调整系数ci的初始值、有效范围和变化步长，作为同化算法的输入参数，其中gaini为定标增益系数，offseti为绝对定标偏移系数；然后以迭代的方式进行如下处理，
步骤2.4.1，基于待定标影像的dn值计算大气表观反射率，并通过大气校正得到地表反射率ρ(target,i)；
步骤2.4.2，基于待定标影像的地表反射率ρ(target,i)，利用步骤2.3所建立brdf调整系数ci的当前值和步骤2.2所得光谱调整系数sbafi通过下式模拟得到参考影像的地表反射率ρsimulated(reference,i)；
ρsimulated(reference,i)＝ρ(target,i)/ci/sbafi
步骤2.4.3，比较参考影像的模拟地表反射率ρsimulated(reference,i)与真实反演的地表反射率ρ(reference,i)之间的差异，若差异比小于预设阈值，便得到最佳的辐射定标系数，结束迭代；反之，利用同化算法优化更新辐射定标系数gaini和offseti以及brdf调整系数ci，然后回到步骤2.4.1；
步骤3，针对每个波段，将一个波段的交叉定标结果中偏离平均值大于预设阈值的点删除，然后取平均值作为最终的交叉定标结果，获得最终的辐射定标系数gaini和offseti。
而且，所述计算大气表观反射率，采用以下公式，
li＝gaini·dni+offseti
其中，i为波段序号，li为转换后辐亮度，dni为该波段的数字量化值，gaini为定标增益系数，offseti为绝对定标偏移系数；
ρtoa(i)＝π·li·d2/(ei·cosθ)
其中，ρtoa(i)为i波段的大气表观反射率，θ为太阳高度角，d是天文单位的日地距离；ei是大气层顶的平均太阳光谱辐照度。
而且，计算参考影像与待定标影像之间的光谱调整系数sbafi，采用以下公式，
其中，sbafi为某一地物第i波段的光谱匹配系数，λ表示波长，ρ(λ)为目标地物的连续光谱反射率，f(λ)为连续的大气层外太阳光谱辐照度，starget,i(λ)和sreference,i(λ)分别为待定标传感器和参考传感器第i波段的波谱响应函数，atarget,i，btarget,i和areference,i，breference,i分别为待定标传感器和参考传感器第i波段的波长范围下限和上限。
而且，参考影像的模拟地表反射率ρsimulated(reference,i)与真实反演的地表反射率ρ(reference,i)之间的差异比计算方式为，
而且，同化算法采用sce-ua算法、遗传算法、模拟退火算法、变分算法、kalman滤波或贝叶斯算法。
本发明相应提供一种用于大观测角传感器的交叉辐射定标系统，包括以下模块，
配准模块，用于配准待定标影像和参考影像，对待定标影像和参考影像的每一组对应波段分别进行随机采样得到对应的定标点；
交叉定标模块，用于针对每个波段，分别对每一对定标点进行如下的交叉定标操作，
步骤s1，基于参考影像的定标点的dn值计算大气表观反射率，并通过大气校正得到地表反射率；
步骤s2，对参考影像的地表反射率进行匹配找到最接近的光谱曲线，计算参考影像与待定标影像之间的光谱调整系数sbafi，下标i用于表示第i波段；
步骤s3，建立brdf调整系数ci，所述brdf调整系数ci为待定标影像第i波段的地表二向性反射分布函数与参考影像第i波段的地表二向性反射分布函数的比值；
步骤s4，初始化参数，包括设定待定标影像的辐射定标系数gaini和offseti的初始值、有效范围和变化步长，以及brdf调整系数ci的初始值、有效范围和变化步长，作为同化算法的输入参数，其中gaini为定标增益系数，offseti为绝对定标偏移系数；然后以迭代的方式进行如下处理，
步骤s4.1，基于待定标影像的dn值计算大气表观反射率，并通过大气校正得到地表反射率ρ(target,i)；
步骤s4.2，基于待定标影像的地表反射率ρ(target,i)，利用步骤s3所建立brdf调整系数ci的当前值和步骤s2所得光谱调整系数sbafi通过下式模拟得到参考影像的地表反射率ρsimulated(reference,i)；
ρsimulated(reference,i)＝ρ(target,i)/ci/sbafi
步骤s4.3，比较参考影像的模拟地表反射率ρsimulated(reference,i)与真实反演的地表反射率ρ(reference,i)之间的差异，若差异比小于预设阈值，便得到最佳的辐射定标系数，结束迭代；反之，利用同化算法优化更新辐射定标系数gaini和offseti以及brdf调整系数ci，然后回到步骤s4.1；
输出模块，用于针对每个波段，将一个波段的交叉定标结果中偏离平均值大于预设阈值的点删除，然后取平均值作为最终的交叉定标结果，获得最终的辐射定标系数gaini和offseti。
而且，所述计算大气表观反射率，采用以下公式，
li＝gaini·dni+offseti
其中，i为波段序号，li为转换后辐亮度，dni为该波段的数字量化值，gaini为定标增益系数，offseti为绝对定标偏移系数；
ρtoa(i)＝π·li·d2/(ei·cosθ)
其中，ρtoa(i)为i波段的大气表观反射率，θ为太阳高度角，d是天文单位的日地距离；ei是大气层顶的平均太阳光谱辐照度。
而且，计算参考影像与待定标影像之间的光谱调整系数sbafi，采用以下公式，
其中，sbafi为某一地物第i波段的光谱匹配系数，λ表示波长，ρ(λ)为目标地物的连续光谱反射率，f(λ)为连续的大气层外太阳光谱辐照度，starget,i(λ)和sreference,i(λ)分别为待定标传感器和参考传感器第i波段的波谱响应函数，atarget,i，btarget,i和areference,i，breference,i分别为待定标传感器和参考传感器第i波段的波长范围下限和上限。
而且，参考影像的模拟地表反射率ρsimulated(reference,i)与真实反演的地表反射率ρ(reference,i)之间的差异比计算方式为，
而且，同化算法采用sce-ua算法、遗传算法、模拟退火算法、变分算法、kalman滤波或贝叶斯算法。
本发明的特点是：
(1)本发明针对大观测角引起的交叉定标问题，消除待定标传感器和参考传感器在光谱响应、大气辐射传输以及brdf三个方面的差异性，使得交叉定标结果准确性更高；
(2)本发明建立了待定标传感器和参考传感器之间brdf的调整系数以消除两者之间的二向性反射差异，采用待定标传感器与参考传感器brdf的比值，利用同化算法以迭代优化的方式求解，无需引入其他辅助数据通过复杂模型计算brdf，计算方法简单稳定，效果可控。最后剔除误差较大的定标结果，取均值作为最终的辐射定标系数，可操作性强，而且可排除引入辅助数据计算brdf产生的误差以及异常定标点对结果的影响。
(3)本发明以迭代优化的方法得到辐射定标系数，操作简单，每对定标点都可以得到一组结果，可排除偶然误差，结果可控。
因此，本发明具有重要的实用意义和市场价值。
附图说明
图1为本发明的流程图。
具体实施方式
以下结合附图和实施例详细说明本发明的技术方案。
针对普通传感器传统方法无需进行大气校正和brdf校正，针对大观测角传感器必须解决大气和brdf的影响(因为二者和观测角度显著相关)，所以本发明针对大观测角传感器提出6s模型大气校正和brdf调整系数。
本发明实施例提供一种针对大观测角传感器的基于数据同化的交叉辐射定标方法，流程如图1所示，具体包括以下步骤：
(1)配准待定标影像和参考影像，对待定标影像和参考影像的每一组对应波段分别进行随机采样得到对应的定标点：
实施例中，对待定标影像和参考影像对进行几何配准并选取感兴趣区，以参考影像为基准，采用sift匹配方法选取数量足够且分布均匀的控制点对相同区域的参考影像进行几何配准，误差控制在1个像素以内。具体实施时，控制点应当均匀分布全局，优选地选取20个以上。在参考影像上随机选取100000个点，以每个点为中心，计算m×n(具体窗口大小可由本领域技术人员预设，和空间分辨率有关)窗口内所有像素标准差与均值的比值(cv)，如果该值小于1％，且位于待定标影像中相同位置的p×q(和参考影像的窗口大小一致，和空间分辨率有关)窗口的cv值也小于1％，则待定标影像和参考影像中的相应窗口被选为一对定标点，取窗口的均值作为定标点的dn值。
设参考影像的分辨率记为分辨率1，设待定标影像的分辨率记为分辨率2，存在以下关系：
(p×q×分辨率1×分辨率1)约等于(m×n×分辨率2×分辨率2)
(2)针对每个波段，分别对每一对定标点进行如下的交叉定标操作：
(2.1)基于参考影像的定标点的dn值计算大气表观反射率，并通过大气校正得到地表反射率：
本发明优选使用6s大气校正模型进行校正。6s大气校正模型建立在辐射传输理论的基础之上，描述了大气如何影响辐射在太阳-地表-传感器之间的传输，相比于其他大气校正模型，其应用范围广，不受研究区特点及目标类型的影响。6s模型进行大气校正所需的参数包括大气条件——大气模型，气溶胶类型，气溶胶光学厚度/大气能见度；观测几何——太阳天顶角，太阳方位角，卫星天顶角，卫星方位角，卫星高度；波普响应——波谱响应函数，具体实施时6s模型需要的气溶胶参数可从modis的aot产品数据获取。
实施例中，基于参考影像定标点的dn值dnreference,i利用以下的公式(1)和公式(2)计算得到大气表观反射率ρtoa(reference,i)，并通过6s模型大气校正得到地表反射率ρ(reference,i)；
li＝gaini·dni+offseti(1)
式中，i为波段序号，li为转换后辐亮度，单位为w·m-2·sr-1·μm-1，dni为该波段的数字量化值(dn值)，无量纲，gaini为定标增益系数，offseti为绝对定标偏移系数，参考卫星各个波段辐射定标系数gaini和offseti是已知的，待定标的目标卫星的相应系数待求解。然后可以通过如下公式将辐射亮度值转换为大气表观反射率：
ρtoa(i)＝π·li·d2/(ei·cosθ)(2)
其中，ρtoa(i)为i波段的大气表观反射率，无量纲，θ为太阳高度角，d是天文单位的日地距离。ei是大气层顶的平均太阳光谱辐照度，单位为w·m-2。
本步骤中，将基于参考影像定标点的dn值dnreference,i代入公式(1)作为dni计算，根据公式(2)得到相应的ρtoa(i)，即基于参考影像的大气表观反射率ρtoa(reference,i)。
(2.2)将参考影像的地表反射率与已有的地物光谱库(例如usgs光谱库)中的光谱曲线进行匹配找到最接近的光谱曲线，进而根据公式(3)计算参考影像与待定标影像之间的光谱调整系数sbafi；
其中，sbafi为某一地物第i波段的光谱匹配系数，λ表示波长，ρ(λ)为目标地物的连续光谱反射率，f(λ)为连续的大气层外太阳光谱辐照度，starget,i(λ)和sreference,i(λ)分别为待定标传感器和参考传感器第i波段的波谱响应函数，atarget,i，btarget,i和areference,i，breference,i分别为待定标传感器和参考传感器第i波段的波长范围下限和上限。
(2.3)建立brdf调整系数ci，消除待定标传感器与参考传感器之间成像时的角度的差异，模拟得到与参考影像相同观测角和相同太阳光入射角时待定标影像的地表反射率，ci的表达式如下：
其中，brdf是二向性反射分布函数，
θs-target和
(2.4)初始化参数，包括设定待定标影像的辐射定标系数gaini和offseti的初始值、有效范围和变化步长，以及brdf调整系数ci的初始值、有效范围和变化步长，作为同化算法的输入参数；具体实施时，初始化参数可由本领域技术人员根据具体情况设定，例如对gf4-pms卫星，这三个参数的初始值分别为0.1，5，1，有效范围分别为[0,0.5]，[0,20]，[0,1]，变化步长分别为0.001，0.1，0.01。
然后以迭代的方式进行如下处理：
利用辐射定标系数反演待定标影像定标点的地表反射率，进而通过brdf和光谱调整得到参考影像的模拟地表反射率，比较参考影像的模拟地表反射率与真实反演的地表反射率之间的差异，若满足收敛条件，便得到最佳的辐射定标系数；反之，利用同化算法优化更新辐射定标系数和brdf调整系数。不断重复该处理直到参考影像的模拟地表反射率与真实反演的地表反射率之间的差异小于设定阈值，便得到最优的辐射定标系数。
本发明中，利用同化算法求解brdf调整系数ci，通过迭代不断优化，直到满足收敛条件，便得到最优的brdf调整系数，无需通过复杂模型求解。
实施例中，具体实现包括以迭代方式进行以下流程：
(2.4.1)基于待定标影像的dn值dntarget,i，利用当前的辐射定标系数gaini和offseti取值，通过公式(1)和公式(2)计算大气表观反射率ρtoa(target,i)，进而通过6s模型大气校正得到地表反射率ρ(target,i)；
本步骤中，将dntarget,i代入公式(1)作为dni计算，根据公式(2)得到相应的ρtoa(i)，即基于待定标影像的大气表观反射率ρtoa(target,i)。第一次执行本步骤时，辐射定标系数gaini和offseti取值为初始值，后续采用上一次迭代优化后的结果。
(2.4.2)基于待定标影像的地表反射率ρ(target,i)利用步骤2.3所建立brdf调整系数ci的当前值和步骤2.2所得光谱调整系数sbafi通过公式(5)模拟得到参考影像的地表反射率ρsimulated(reference,i)；
ρsimulated(reference,i)＝ρ(target,i)/ci/sbafi(5)
其中ci是brdf调整系数，是待定标传感器与参考传感器brdf的比值。
(2.4.3)比较参考影像的模拟地表反射率ρsimulated(reference,i)与真实反演的地表反射率ρ(reference,i)之间的差异，若差异比小于预设阈值(具体实施时本领域技术人员可预设取值，例如0.01)，便得到最佳的辐射定标系数，结束迭代；反之，利用同化算法优化更新辐射定标系数gaini和offseti，以及brdf调整系数ci，然后回到步骤(2.4.1)。
差异比计算如下：
所用到的同化方法可以选择常见方法，对于具体采用哪一种方法不做限制，例如sce-ua算法、遗传算法、模拟退火算法、变分算法、kalman滤波、贝叶斯算法等数据同化方法。
利用同化算法求解brdf调整系数ci，通过迭代不断优化，直到满足收敛条件，便得到最优的brdf调整系数，无需通过复杂模型求解；
利用同化算法计算辐射定标系数gaini和offseti，通过迭代过程不断优化更新该系数值，直到参考影像模拟地表反射率与反演的地表反射率之间的差异最小。
执行同化算法时，通过参数的范围和变化步长设置，在这个有效范围内以变化步长为间隔生成大量的随机值。多个参数的情况下，随机值会随机组合。首先利用初始值得到模拟结果，将模拟结果与真实结果或者参考进行比较，不符合阈值条件时，同化算法会根据其规则找到更优的参数值，并淘汰大量随机值中的部分值，并将剩下的随机值按一定规则排列或组合成多个子集。不断重复这个过程，直到满足收敛条件，那么此时找到的参数值就是最佳解。
(3)当对于每个波段的所有的定标点进行上述步骤(2)的操作完成后，将一个波段的交叉定标结果中偏离平均值较远的点删除，然后取平均值作为最终的交叉定标结果，获得最终的辐射定标系数gaini和offseti。实施例中，去掉超过增益系数gaini平均值±10％的点，计算剩余点的平均值，如果剩余点中还有超过平均值±10％的点，将其剔除，不断重复此过程，直到没有超过平均值±10％的点，则将剩余点的均值作为该波段的增益系数。偏移系数offseti也通过同样的方法获取。
不同于传统方法利用所有定标点拟合得到一组辐射定标系数，本发明提供的针对方法每一对定标点都可以得到一组辐射定标系数，可以剔除误差较大的点，利用定标结果较为一致的定标点得到最终的辐射定标系数。因此该方法可以排除异常的定标点引起的误差。
具体实施时，以上流程可采用计算机软件技术实现自动运行。本发明还提出采用模块化方式提供相应系统。
本发明实施例相应提供一种用于大观测角传感器的交叉辐射定标系统，包括以下模块，
配准模块，用于配准待定标影像和参考影像，对待定标影像和参考影像的每一组对应波段分别进行随机采样得到对应的定标点；
交叉定标模块，用于针对每个波段，分别对每一对定标点进行如下的交叉定标操作，
步骤s1，基于参考影像的定标点的dn值计算大气表观反射率，并通过大气校正得到地表反射率；
步骤s2，对参考影像的地表反射率进行匹配找到最接近的光谱曲线，计算参考影像与待定标影像之间的光谱调整系数sbafi，下标i用于表示第i波段；
步骤s3，建立brdf调整系数ci，所述brdf调整系数ci为待定标影像第i波段的地表二向性反射分布函数与参考影像第i波段的地表二向性反射分布函数的比值；
步骤s4，初始化参数，包括设定待定标影像的辐射定标系数gaini和offseti的初始值、有效范围和变化步长，以及brdf调整系数ci的初始值、有效范围和变化步长，作为同化算法的输入参数，其中gaini为定标增益系数，offseti为绝对定标偏移系数；然后以迭代的方式进行如下处理，
步骤s4.1，基于待定标影像的dn值计算大气表观反射率，并通过大气校正得到地表反射率ρ(target,i)；
步骤s4.2，基于待定标影像的地表反射率ρ(target,i)，利用步骤s3所建立brdf调整系数ci的当前值和步骤s2所得光谱调整系数sbafi通过下式模拟得到参考影像的地表反射率ρsimulated(reference,i)；
ρsimulated(reference,i)＝ρ(target,i)/ci/sbafi
步骤s4.3，比较参考影像的模拟地表反射率ρsimulated(reference,i)与真实反演的地表反射率ρ(reference,i)之间的差异，若差异比小于预设阈值，便得到最佳的辐射定标系数，结束迭代；反之，利用同化算法优化更新辐射定标系数gaini和offseti以及brdf调整系数ci，然后回到步骤s4.1；
输出模块，用于针对每个波段，将一个波段的交叉定标结果中偏离平均值大于预设阈值的点删除，然后取平均值作为最终的交叉定标结果，获得最终的辐射定标系数gaini和offseti。
各模块具体实现可参见相应步骤，本发明不予赘述。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
- 还没有人留言评论。精彩留言会获得点赞！