一种图像地形辐射校正方法及装置与流程
本申请涉及图像辐射处理技术领域，更具体地说，尤其涉及一种图像地形辐射校正方法及装置。
背景技术：
随着遥感技术发展，高分辨率图像日益增多，相关应用越来越广泛。但随着分辨率的提高，由于地形起伏造成的图像辐射畸变表现的越来越明显，比如山体阴坡与阳坡的同一类地物在图像上表现出显著差异，对图像分类或者其他定量遥感应用带来很大困难。
因此，对图像进行地形辐射校正处理，是高分辨率图像在遥感技术中应用的迫切需求。
技术实现要素：
基于上述需求，本申请提出一种图像地形辐射校正方法及装置，能够实现对图像的地形辐射校正处理。
一种图像地形辐射校正方法，包括：
读取待校正图像的太阳天顶角和太阳方位角；
利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
根据所述坡度计算得到所述待校正图像的天空观测因子，以及根据所述太阳天顶角、所述坡度、所述坡向和所述太阳方位角计算得到太阳有效入射角；
提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子；
利用辐射传输模型计算得到所述待校正图像的大气相关辐射量；其中，所述大气相关辐射量包括大气上行透过率、下行透过率、地面直射辐照度和地面散射辐照度；
根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的反射率计算公式或相对反射率计算公式，计算得到所述待校正图像各像元的反射率。
可选的，在读取待校正图像的太阳天顶角和太阳方位角之前，所述方法还包括：
参照待校正图像的投影方式对数字高程模型进行重投影，使所述待校正图像的投影方式与所述数字高程模型的投影方式一致；
和/或，
根据数字高程模型的地理范围对待校正图像进行裁剪处理，使所述数字高程模型的地理范围覆盖所述待校正图像整个图像。
可选的，所述利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子，包括：
当数字高程模型与所述待校正图像的分辨率相同时，直接利用所述数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
当数字高程模型的分辨率高于所述待校正图像的分辨率时，先按照所述待校正图像的分辨率对所述数字高程模型进行重采样处理，然后利用重采样后的数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用重采样后的数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
当数字高程模型的分辨率低于所述待校正图像的分辨率时，在利用所述数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子后，按照所述待校正图像的分辨率对计算得到的坡度、坡向和参考遮蔽因子进行重采样处理。
可选的，所述提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子，包括：
利用直方图方法提取所述待校正图像中的阴影区，得到图像自身遮蔽因子；
利用所述参考遮蔽因子剔除所述图像自身遮蔽因子中的水体和云阴影，确定遮蔽因子。
可选的，所述根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的反射率计算公式或相对反射率计算公式，计算得到所述待校正图像各像元的反射率，包括：
当所述待校正图像具有辐射定标系数时，根据所述辐射定标系数，以及所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子和所述大气相关辐射量，利用预设的反射率计算公式，计算得到所述待校正图像的像元反射率；
当所述待校正图像不具备辐射定标系数时，根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的相对反射率计算公式，计算得到所述待校正图像的像元反射率。
可选的，所述预设的相对反射率计算公式为：
其中，dn(x,y)表示像元灰度值，dnba表示程辐射对应的像元值，τv表示大气上行透过率，b(x,y)表示遮蔽因子，edir表示地面直射辐照度，τs表示大气下行透过率，β表示太阳有效入射角，edif表示地面散射辐照度，θs表示太阳天顶角，vsky(x,y)表示天空观测因子，eiter表示邻近像元的多次反射辐照度。
可选的，所述方法还包括：
对所述待校正图像非阴影区域像元的反射率进行校正处理。
可选的，所述方法还包括：
对所述待校正图像阴影边界区域进行平滑处理。
一种图像地形辐射校正装置，包括：
数据读取单元，用于读取待校正图像的太阳天顶角和太阳方位角；
第一计算单元，用于利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
第二计算单元，用于根据所述坡度计算得到所述待校正图像的天空观测因子，以及根据所述太阳天顶角、所述坡度、所述坡向和所述太阳方位角计算得到太阳有效入射角；
第三计算单元，用于提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子；
第四计算单元，用于利用辐射传输模型计算得到所述待校正图像的大气相关辐射量；其中，所述大气相关辐射量包括大气上行透过率、下行透过率、地面直射辐照度和地面散射辐照度；
第五计算单元，用于根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的反射率计算公式或相对反射率计算公式，计算得到所述待校正图像各像元的反射率。
可选的，所述装置还包括：
预处理单元，用于参照待校正图像的投影方式对数字高程模型进行重投影，使所述待校正图像的投影方式与所述数字高程模型的投影方式一致；
和/或，
根据数字高程模型的地理范围对待校正图像进行裁剪处理，使所述数字高程模型的地理范围覆盖所述待校正图像整个图像。
可选的，所述第一计算单元利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子时，具体用于：
当数字高程模型与所述待校正图像的分辨率相同时，直接利用所述数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
当数字高程模型的分辨率高于所述待校正图像的分辨率时，先按照所述待校正图像的分辨率对所述数字高程模型进行重采样处理，然后利用重采样后的数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用重采样后的数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
当数字高程模型的分辨率低于所述待校正图像的分辨率时，在利用所述数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子后，按照所述待校正图像的分辨率对计算得到的坡度、坡向和参考遮蔽因子进行重采样处理。
可选的，所述第三计算单元提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子时，具体用于：
利用直方图方法提取所述待校正图像中的阴影区，得到图像自身遮蔽因子；
利用所述参考遮蔽因子剔除所述图像自身遮蔽因子中的水体和云阴影，确定遮蔽因子。
可选的，所述第五计算单元根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的反射率计算公式或相对反射率计算公式，计算得到所述待校正图像各像元的反射率时，具体用于：
当所述待校正图像具有辐射定标系数时，根据所述辐射定标系数，以及所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子和所述大气相关辐射量，利用预设的反射率计算公式，计算得到所述待校正图像的像元反射率；
当所述待校正图像不具备辐射定标系数时，根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的相对反射率计算公式，计算得到所述待校正图像的像元反射率。
可选的，所述预设的相对反射率计算公式为：
其中，dn(x,y)表示像元灰度值，dnba表示程辐射对应的像元值，τv表示大气上行透过率，b(x,y)表示遮蔽因子，edir表示地面直射辐照度，τs表示大气下行透过率，β表示太阳有效入射角，edif表示地面散射辐照度，θs表示太阳天顶角，vsky(x,y)表示天空观测因子，eiter表示邻近像元的多次反射辐照度。
可选的，所述装置还包括：
校正处理单元，用于对所述待校正图像非阴影区域像元的反射率进行校正处理。
可选的，所述装置还包括：
平滑处理单元，用于对所述待校正图像阴影边界区域进行平滑处理。
本申请提出的图像地形辐射校正方法，借助数字高程模型，依次获取像元反射率计算所需的各项参数，然后根据所获取的参数，利用预设的反射率计算公式或相对反射率计算公式计算得到图像各像元的反射率，达到图像地形辐射校正的目的，即该方法实现了对图像的地形辐射校正。
进一步的，由于本申请可以通过发射率计算公式计算得到像元的反射率，也可以通过相对反射率计算公式计算得到像元的相对反射率，因此本申请实施例技术方案对于待校正图像具有辐射定标系数和不具备辐射定标系数的情况，均能够对其进行地形辐射校正处理。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1是本申请实施例提供的一种图像地形辐射校正方法的流程示意图；
图2是本申请实施例提供的另一种图像地形辐射校正方法的流程示意图；
图3是本申请实施例提供的含有山区阴影的图像的示意图；
图4是本申请实施例提供的图像直方图示意图；
图5是本申请实施例提供的一种图像地形辐射校正装置的结构示意图。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
本申请实施例公开了一种图像地形辐射校正方法，参见图1所示，该方法包括：
s101、读取待校正图像的太阳天顶角和太阳方位角；
具体的，一般的遥感图像都会直接给出太阳天顶角θs和太阳方位角φs信息，因此，直接从上述的待校正图像中读取太阳天顶角和太阳方位角即可。
一般图像(例如landsattm或者etm+图像)只提供中心像元处的太阳天顶角θs和太阳方位角φs数值，全图所有像元的太阳天顶角θs和太阳方位角φs均采用该数值。对于宽幅图像，一般按照一定像元数间隔提供太阳天顶角和方位角，此时需要进行重采样得到每个像元处的θs和φs数值。例如sentinel-2数据每隔5km提供一组太阳天顶角和方位角数据(分别为22*22个像元)，对于10
s102、利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
示例性的，通过利用数字高程模型(digitalelevationmodel,dem)调用envi提供的工具topographic>topographicmodeling，或者提供的对应函数来计算得到待校正图像的坡度s、坡向a，以及结合上述太阳天顶角θs和太阳方位角φs，计算得到参考遮蔽因子b0(x,y)。其中，(x,y)表示像元坐标，后文各(x,y)的意义相同，均为像元坐标。
进一步的，本申请实施例还公开了，上述、利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子分三种情况处理：
当数字高程模型dem与所述待校正图像的分辨率相同时，直接利用所述数字高程模型dem计算得到所述待校正图像的地形坡度s和坡向a，以及利用所述数字高程模型dem结合所述太阳天顶角θs和太阳方位角φs，计算得到参考遮蔽因子b0(x,y)；
当数字高程模型dem的分辨率高于所述待校正图像的分辨率时，先按照所述待校正图像的分辨率对所述数字高程模型dem进行重采样处理，然后利用重采样后的数字高程模型dem计算得到所述待校正图像的地形坡度s和坡向a，以及利用重采样后的数字高程模型dem结合所述太阳天顶角θs和太阳方位角φs，计算得到参考遮蔽因子b0(x,y)；
当数字高程模型dem的分辨率低于所述待校正图像的分辨率时，在利用所述数字高程模型dem计算得到所述待校正图像的地形坡度s和坡向a，以及利用所述数字高程模型dem结合所述太阳天顶角θs和太阳方位角φs，计算得到参考遮蔽因子b0(x,y)后，按照所述待校正图像的分辨率对计算得到的坡度s、坡向a和参考遮蔽因子b0(x,y)进行重采样处理。
可以理解，本申请实施例技术方案在计算图像坡度、坡向和遮蔽因子时，充分考虑到数字高程模型的分辨率与图像分辨率是否匹配，对于不匹配的情况，按照图像分辨率对数字高程模型的分辨率进行重采样，或者对计算结果按照图像分辨率进行重采样，由此使得计算得到的参数为图像实际参数，解决了图像分辨率与数字高程模型分辨率不一致的问题。
s103、根据所述坡度计算得到所述待校正图像的天空观测因子，以及根据所述太阳天顶角、所述坡度、所述坡向和所述太阳方位角计算得到太阳有效入射角；
其中，上述天空观测因子可以根据上述的坡度s，通过如下三角法计算公式计算得到：
vsky(x,y)＝cos2(s(x,y)/2)
其中，vsky(x,y)表示像元的天空观测因子，s(x,y)表示像元坡度。
上述太阳有效入射角β，可以根据上述太阳天顶角θs、坡度s、坡向a和太阳方位角φs，通过如下计算公式计算得到：
cosβ＝cosθscoss+sinssinθscos(φs-a)
s104、提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子；
具体的，上述参考遮蔽因子是仅通过数字高程模型计算得到的遮蔽因子，并不是直接由待校正图像确定的遮蔽因子。为了使获取的遮蔽因子更准确，本申请实施例从待校正图像中直接提取遮蔽因子，然后利用上述的参考遮蔽因子对直接提取的遮蔽因子进行校正，得到最终确定的遮蔽因子。
上述待校正图像的阴影区，即为对待校正图像的阴影区和非阴影区进行区分识别，其识别结果即作为从待校正图像中直接提取的遮蔽因子，记为b1(x,y)。
本申请实施例对上述待校正图像的阴影区和非阴影区进行区分识别，得到上述的遮蔽因子b1(x,y)后，再利用上述的参考遮蔽因子b0(x,y)对b1(x,y)进行校正，剔除其中的误判因素等，得到最终确定的遮蔽因子b(x,y)。
s105、利用辐射传输模型计算得到所述待校正图像的大气相关辐射量；其中，所述大气相关辐射量包括大气上行透过率、下行透过率、地面直射辐照度和地面散射辐照度；
具体的，基于上述待校正图像自身(如暗像元算法、亮地表算法)或者其他辅助数据(如实测或者modis数据)估算图像气溶胶光学厚度、水汽含量，同时根据先验知识给出图像获取时的大气模型等辐射传输模型需要的输入。可利用大气辐射传输模型modtran、6s或者6sv计算以下大气参量：大气上行透过率τv、大气下行透过率τs、地面直射辐照度edir和地面散射辐照度edif。
s106、根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的反射率计算公式或相对反射率计算公式，计算得到所述待校正图像各像元的反射率。
具体的，通过执行上述步骤s101～s105分别得到上述的太阳天顶角θs、天空观测因子vsky(x,y)、太阳有效入射角β、遮蔽因子b(x,y)以及各项大气相关辐射量：大气上行透过率τv、大气下行透过率τs、地面直射辐照度edir和地面散射辐照度edif后，本申请实施例利用预设的反射率计算公式或相对反射率计算公式，计算得到上述待校正图像各像元的反射率。
其中，上述的反射率计算公式为：
其中，c0和c1为图像辐射定标系数，dn(x,y)表示像元灰度值，lba是图像程辐射对应的辐亮度，即图像有效像元区域最小值，τv表示大气上行透过率，b(x,y)表示遮蔽因子，edir表示地面直射辐照度，τs表示大气下行透过率，β表示太阳有效入射角，edif表示地面散射辐照度，θs表示太阳天顶角，vsky(x,y)表示天空观测因子，eiter表示邻近像元的多次反射辐照度。
邻近像元的多次反射辐照度eiter的计算公式如下：
式中：eg为水平表面接收的太阳的直射辐照度与漫散射辐照度之和；vsky为天空观测因子。
需要说明的是，由于上述的反射率计算公式需要图像的辐射定标系数，因此，当上述待校正图像具有辐射定标系数c0和c1时，可以根据该辐射定标系数c0和c1，以及上述的太阳天顶角θs、天空观测因子vsky(x,y)、太阳有效入射角β、遮蔽因子b(x,y)和各项大气相关辐射量：大气上行透过率τv、大气下行透过率τs、地面直射辐照度edir和地面散射辐照度edif，利用上述预设的反射率计算公式，计算得到待校正图像的像元反射率。
另一方面，当上述的待校正图像不具备辐射定标系数时，上述的反射率计算公式则无法应用。此时，本申请实施例设计相对反射率计算公式，用于计算待校正图像的相对反射率，该公式具体如下：
其中，ρt(x,y)表示像元的相对反射率，dn(x,y)表示像元灰度值，dnba表示程辐射对应的像元值，τv表示大气上行透过率，b(x,y)表示遮蔽因子，edir表示地面直射辐照度，τs表示大气下行透过率，β表示太阳有效入射角，edif表示地面散射辐照度，θs表示太阳天顶角，vsky(x,y)表示天空观测因子，eiter表示邻近像元的多次反射辐照度。
可见，上述的相对反射率计算公式不需要图像的辐射定标系数，因此，当上述的待校正图像不具备辐射定标系数时，根据上述步骤s101～s105中获取的太阳天顶角θs、天空观测因子vsky(x,y)、太阳有效入射角β、遮蔽因子b(x,y)和各项大气相关辐射量：大气上行透过率τv、大气下行透过率τs、地面直射辐照度edir和地面散射辐照度edif，利用上述的相对反射率计算公式，计算得到待校正图像的像元反射率，该像元反射率为像元的相对反射率。
需要说明的是，由于缺少辐射定标系数，无法求得真实地表反射率，因此这种情况下上述的邻近像元的多次反射辐照度eiter无法迭代计算得到。此时，本申请实施例设定周围像元的平均反射率为0.1，通过如下计算公式计算得到邻近像元多次反射辐照度：
eiter＝(edir+edif)ξ2/(1-ξ)
ξ＝0.1·(1-vsky(x,y))
通过上述介绍可见，本申请实施例提出的图像地形辐射校正方法，借助数字高程模型，依次获取像元反射率计算所需的各项参数，然后根据所获取的参数，利用预设的反射率计算公式或相对反射率计算公式计算得到图像各像元的反射率，达到图像地形辐射校正的目的，即该方法实现了对图像的地形辐射校正。
进一步的，由于本申请实施例可以通过发射率计算公式计算得到像元的反射率，也可以通过相对反射率计算公式计算得到像元的相对反射率，因此本申请实施例技术方案对于待校正图像具有辐射定标系数和不具备辐射定标系数的情况，均能够对其进行地形辐射校正处理。
并且，本申请实施例在获取地形辐射校正所需的各项参数时，充分考虑到数字高程模型和待校正图像的分辨率是否一致的情况，对数字高程模型或获取的参数按照待校正图像分辨率进行重采样，以及，在计算遮蔽因子时直接从图像中提取遮蔽因子，使得获取的参数切合图像分辨率，进而保证对图像的辐射校正更准确。
可选的，在本申请的另一个实施例中还公开了，在读取待校正图像的太阳天顶角和太阳方位角之前，所述方法还包括：
参照待校正图像的投影方式对数字高程模型进行重投影，使所述待校正图像的投影方式与所述数字高程模型的投影方式一致；
和/或，
根据数字高程模型的地理范围对待校正图像进行裁剪处理，使所述数字高程模型的地理范围覆盖所述待校正图像整个图像。
具体的，在从待校正图像中读取太阳天定角和太阳方位角信息之前，本申请实施例首先对数字高程模型和待校正图像进行一致性检查，主要可从两方面执行：
一方面，将待校正图像进行几何精校正或者正射校正，确保数字高程模型与待校正图像之间的配准精度，并且检查数字高程模型的投影方式与待校正图像的投影方式是否一致，如果不一致，则参照待校正图像的投影方式，对数字高程模型进行重投影，使待校正图像的投影方式与数字高程模型的投影方式一致。
另一方面，检查数字高程模型的地理范围是否覆盖整个待校正图像，如果没有覆盖整个待校正图像，则根据数字高程模型的地理范围对待校正图像进行裁剪处理，使数字高程模型的地理范围覆盖待校正图像整个图像。
需要说明的是，上述两方面的处理，可以根据实际情况选择执行或者不执行，当选择执行时可以单一执行，也可以全部执行。
可选的，参见图2所示，在本申请的另一个实施例中还公开了，所述提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子，包括：
s204、利用直方图方法提取所述待校正图像中的阴影区，得到图像自身遮蔽因子；
具体的，含有山区阴影的图像在直方图上表现为明显的双峰分布，例如对图3所示的含有山区阴影的图像，获取其直方图可以得到图4所示的结果，直方图分布情况明显呈双峰分布，两个峰值之间的谷值处的灰度值可作为区分阴影和非阴影的阈值，灰度值小于该阈值的区域即可判断为阴影区。
按照上述方法，本申请实施例获取上述的待校正图像的直方图，从直方图中查找两个峰值之间的谷值，确定待校正图像中的阴影区和非阴影区的划分阈值，然后利用该阈值识别待校正图像中的阴影区和非阴影区，并对识别的阴影区和非阴影区进行标记，例如阴影区标记为0，非阴影区标记为1，得到待校正图像的图像自身遮蔽因子，记为b1(x,y)。
s205、利用所述参考遮蔽因子剔除所述图像自身遮蔽因子中的水体和云阴影，确定遮蔽因子。
具体的，由于水体和云阴影在图像上的灰度值较低，在上述直方图判断时容易与山体遮蔽区混淆，被识别为图像阴影区，因此需要从上述的图像自身遮蔽因子中剔除水体和云阴影，得到最终的遮蔽因子。
本申请实施例以从待校正图像中提取的上述图像自身遮蔽因子b1(x,y)为主，结合上述的参考遮蔽因子b0(x,y)，去掉对于水体、云阴影的误判，最终确定遮蔽因子为b(x,y)。
本实施例中的步骤s201～s203、s206～s207分别对应图1所示的方法实施例中的步骤s101～s103、s105～s106，其具体内容请参见图1所示的方法实施例的内容，此处不再赘述。
可选的，在本申请的另一个实施例中还公开了，在计算得到待校正图像各像元的反射率后，所述方法还包括：
对所述待校正图像非阴影区域像元的反射率进行校正处理。
具体的，为了避免朗伯性假定情况下较大太阳天顶角入射时对非阴影区像元反射率的过校正，本申请实施例对上述待校正图像非阴影区域像元的相对反射率加入校正因子g进行校正，得到矫正后的反射率ρt_mm(x,y)为：
ρt_mm(x,y)＝ρt(x,y)·g
其中，校正因子g的计算如下：
根据在复杂地区校正的多数实验表明，g设置为0.2～0.25之间比较合适。βt根据太阳天顶角θs具有如下经验关系：
指数因子b根据波长和地表覆盖设置，对于非植被区域设置为b＝1/2；植被区但波段中心波长小于720nm设置为b＝3/4；植被区但中心波长大于等于720nm设置为b＝1/3。
可选的，在本申请的另一个实施例中还公开了，在计算得到待校正图像各像元的反射率后，所述方法还包括：
对所述待校正图像阴影边界区域进行平滑处理。
具体的，对于数字高程模型分辨率低于图像分辨率的情况，为了避免因阴影边界提取误差可能造成的过校正问题，本申请实施例对上述待校正图像阴影边界区域进行平滑处理。
上述平滑处理过程具体为：首先在上述待校正图像阴影边界区设定缓冲范围(例如左右5个像元)，然后对此范围内的像元反射率值和原始像元值的比值结果ρt(x,y)/dn(x,y)进行均值滤波处理得到mean_filter[ρt(x,y)/dn(x,y)]，最后用上述均值滤波处理结果再次乘以原始像元值mean_filter[ρt(x,y)/dn(x,y)]·dn(x,y)，完成对待校正图像阴影边界区域的平滑处理。
与上述的图像地形辐射校正方法相对应的，本申请另一实施例还公开了一种图像地形辐射校正装置，参见图5所示，该装置包括：
数据读取单元100，用于读取待校正图像的太阳天顶角和太阳方位角；
第一计算单元110，用于利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
第二计算单元120，用于根据所述坡度计算得到所述待校正图像的天空观测因子，以及根据所述太阳天顶角、所述坡度、所述坡向和所述太阳方位角计算得到太阳有效入射角；
第三计算单元130，用于提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子；
第四计算单元140，用于利用辐射传输模型计算得到所述待校正图像的大气相关辐射量；其中，所述大气相关辐射量包括大气上行透过率、下行透过率、地面直射辐照度和地面散射辐照度；
第五计算单元150，用于根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的反射率计算公式或相对反射率计算公式，计算得到所述待校正图像各像元的反射率。
可选的，在本申请的另一个实施例中还公开了，所述装置还包括：
预处理单元，用于参照待校正图像的投影方式对数字高程模型进行重投影，使所述待校正图像的投影方式与所述数字高程模型的投影方式一致；
和/或，
根据数字高程模型的地理范围对待校正图像进行裁剪处理，使所述数字高程模型的地理范围覆盖所述待校正图像整个图像。
其中，作为一种可选的实现方式，所述第一计算单元110利用数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子时，具体用于：
当数字高程模型与所述待校正图像的分辨率相同时，直接利用所述数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
当数字高程模型的分辨率高于所述待校正图像的分辨率时，先按照所述待校正图像的分辨率对所述数字高程模型进行重采样处理，然后利用重采样后的数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用重采样后的数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子；
当数字高程模型的分辨率低于所述待校正图像的分辨率时，在利用所述数字高程模型计算得到所述待校正图像的地形坡度和坡向，以及利用所述数字高程模型结合所述太阳天顶角和所述太阳方位角，计算得到参考遮蔽因子后，按照所述待校正图像的分辨率对计算得到的坡度、坡向和参考遮蔽因子进行重采样处理。
作为一种可选的实现方式，所述第三计算单元130提取所述待校正图像的阴影区，以及利用所述参考遮蔽因子对所述阴影区进行校正处理，确定遮蔽因子时，具体用于：
利用直方图方法提取所述待校正图像中的阴影区，得到图像自身遮蔽因子；
利用所述参考遮蔽因子剔除所述图像自身遮蔽因子中的水体和云阴影，确定遮蔽因子。
作为一种可选的实现方式，所述第五计算单元150根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的反射率计算公式或相对反射率计算公式，计算得到所述待校正图像各像元的反射率时，具体用于：
当所述待校正图像具有辐射定标系数时，根据所述辐射定标系数，以及所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子和所述大气相关辐射量，利用预设的反射率计算公式，计算得到所述待校正图像的像元反射率；
当所述待校正图像不具备辐射定标系数时，根据所述太阳天顶角、所述天空观测因子、所述太阳有效入射角、所述遮蔽因子以及所述大气相关辐射量，利用预设的相对反射率计算公式，计算得到所述待校正图像的像元反射率。
作为一种可选的实现方式，所述预设的相对反射率计算公式为：
其中，dn(x,y)表示像元灰度值，dnba表示程辐射对应的像元值，τv表示大气上行透过率，b(x,y)表示遮蔽因子，edir表示地面直射辐照度，τs表示大气下行透过率，β表示太阳有效入射角，edif表示地面散射辐照度，θs表示太阳天顶角，vsky(x,y)表示天空观测因子，eiter表示邻近像元的多次反射辐照度。
可选的，在本申请的另一个实施例中还公开了，所述装置还包括：
校正处理单元，用于对所述待校正图像非阴影区域像元的反射率进行校正处理。
可选的，在本申请的另一个实施例中还公开了，所述装置还包括：
平滑处理单元，用于对所述待校正图像阴影边界区域进行平滑处理。
具体的，上述图像地形辐射校正装置的各实施例中的各个单元的具体工作内容，请参见上述方法实施例的内容，此处不再赘述。
对于前述的各方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本申请并不受所描述的动作顺序的限制，因为依据本申请，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作和模块并不一定是本申请所必须的。
需要说明的是，本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。对于装置类实施例而言，由于其与方法实施例基本相似，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
本申请各实施例方法中的步骤可以根据实际需要进行顺序调整、合并和删减。
本申请各实施例种装置及终端中的模块和子模块可以根据实际需要进行合并、划分和删减。
本申请所提供的几个实施例中，应该理解到，所揭露的终端，装置和方法，可以通过其它的方式实现。例如，以上所描述的终端实施例仅仅是示意性的，例如，模块或子模块的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个子模块或模块可以结合或者可以集成到另一个模块，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或模块的间接耦合或通信连接，可以是电性，机械或其它的形式。
作为分离部件说明的模块或子模块可以是或者也可以不是物理上分开的，作为模块或子模块的部件可以是或者也可以不是物理模块或子模块，即可以位于一个地方，或者也可以分布到多个网络模块或子模块上。可以根据实际的需要选择其中的部分或者全部模块或子模块来实现本实施例方案的目的。
另外，在本申请各个实施例中的各功能模块或子模块可以集成在一个处理模块中，也可以是各个模块或子模块单独物理存在，也可以两个或两个以上模块或子模块集成在一个模块中。上述集成的模块或子模块既可以采用硬件的形式实现，也可以采用软件功能模块或子模块的形式实现。
专业人员还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件单元，或者二者的结合来实施。软件单元可以置于随机存储器(ram)、内存、只读存储器(rom)、电可编程rom、电可擦除可编程rom、寄存器、硬盘、可移动磁盘、cd-rom、或技术领域内所公知的任意其它形式的存储介质中。
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本申请。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。
- 还没有人留言评论。精彩留言会获得点赞！