大气程辐射值的确定方法和装置与流程

本申请涉及大气环境遥感监测技术，特别涉及一种大气程辐射值的确定方法和装置。

背景技术：

在利用遥感监测技术评估大气颗粒物污染程度时，常常需要先确定表征大气气溶胶消光特性的关键参数：气溶胶光学厚度(Aerosol Optical Depth，AOD)。其中，气溶胶光学厚度具体可以是用于描述大气气溶胶消光特性的光学参数，也是用于推算气溶胶含量、研究气溶胶气候效应等的关键因子。而上述气溶胶光学厚度通常又具体与大气程辐射值有关。其中，大气程辐射值具体可以是指太阳辐射在大气传输过程中经大气分子、气溶胶、冰晶等粒子散射后直接到达传感器的辐射量，是一种用于表征大气质量状况的天空背景辐射遥感分量。因此，如何准确地计算大气程辐射值是本领域中一个重要的问题。

目前，为了能够确定大气程辐射值，大多是粗略地从遥感影像中提取局地均质性像元聚类点集，再根据局地均质性像元聚类点集建立对应的二维光谱特征空间，进而确定出大气程辐射值。现有的方法具体实施时，由于所获取的局地均质性像元聚类点集往往只能出粗略地表征一些植被、湿润土壤等暗像元的地表情况，其中，上述局地均质性像元聚类点集中常常会存在许多噪点等干扰；此外，获取的上述局地均质性像元聚类点集还会存在遗漏一些暗像元地表信息的情况。因此，直接利用局地均质性像元聚类点集确定大气程辐射值的准确度误差较大。综上可知，现有方法具体实施时，往往存在所确定的大气程辐射值误差大、准确度低的技术问题。

针对上述问题，目前尚未提出有效的解决方案。

技术实现要素：

本申请实施例提供了一种大气程辐射值的确定方法和装置，以解决现有方法中存在的确定大气程辐射值误差大、准确度低的技术问题，达到了精确确定大气程辐射值的技术效果。

本申请实施例提供了一种大气程辐射值的确定方法，包括：

获取目标区域的遥感影像；

从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据；

对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集；

根据所述有效暗像元集，确定目标区域的大气程辐射值。

在一个实施方式中，所述从所述遥感影像的像元中划分出暗像元和非暗像元，包括：

获取像元的绿光波段数据、第一短波红外波段数据和第二短波红外波段数据；

根据所述像元的绿光波段数据和第一短波红外波段数据，确定像元的归一化积雪指数；

根据所述像元的第二短波红外波段数据和归一化积雪指数，确定所述像元是否为暗像元。

在一个实施方式中，根据所述像元的第二短波红外波段数据和归一化积雪指数，确定所述像元是否为暗像元，包括：

将满足以下规则的像元确定为暗像元：

NDSI≤0且0.01≤B7≤0.25

其中，NDSI为像元的归一化积雪指数，B7为像元的第二短波红外波段数据。

在一个实施方式中，对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集，包括：

根据所述第一影像数据，建立二维光谱特征空间；其中，所述二维光谱特征空间以像元的第二短波红外波段数据为横轴，以像元的蓝光波段数据的反射率或者红光波段数据的反射率为纵轴；

通过坐标分区法，对所述二维光谱特征空间中的暗像元进行第一次筛选，得到第一暗像元集；

对所述第一暗像元集进行第二次筛选，以去除所述第一暗像元集中的噪声散点，得到第二暗像元集；

通过多次迭代拟合，对所述第二暗像元集进行第三次筛选，以剔除离散的异常点，得到第三暗像元集，并将所述第三暗像元集作为所述有效暗像元集。

在一个实施方式中，通过坐标分区法，对所述二维光谱特征空间中的暗像元进行第一次筛选，得到第一暗像元集，包括：

根据像元的横轴坐标数值，将所述二维光谱特征空间中的暗像元分别划分为多个子区间；

分别提取子区间中纵轴坐标数值最小的暗像元，作为第一暗像元集中的元素，以获取所述第一暗像元集。

在一个实施方式中，对所述第一暗像元集进行第二次筛选，以去除所述第一暗像元集中的噪声散点，得到第二暗像元集，包括：

获取所述第一暗像元集中暗像元的第二短波红外波段数据的反射率；

根据所述暗像元的第二短波红外波段数据的反射率，将所述第一暗像元集中的暗像元划分为多个子集；

分别确定子集的最小二乘相关系数、子集的一元线性拟合方程的斜率和子集的一元线性拟合方程的截距；

根据所述子集的最小二乘相关系数、所述子集的一元线性拟合方程的斜率和所述子集的一元线性拟合方程的截距，从所述多个子集中筛选出符合要求的子集作为所述第二暗像元集。

在一个实施方式中，通过多次迭代拟合，对所述第二暗像元集进行第三次筛选，以剔除离散的异常点，得到第三暗像元集，包括：

通过多次迭代拟合，将所述第二暗像元集中与所述第二暗像元集的一元线性拟合方程的垂直距离大于阈值距离的暗像元剔除。

在一个实施方式中，在获取目标区域的遥感影像后，所述方法还包括：

对所述遥感影像进行云像元检测，确定所述遥感影像中是否存在云像元；

在所述遥感影像中存在云像元的情况下，对所述遥感影像中的云像元进行掩膜处理，得到处理后的影像数据；

相应的，所述从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据，包括：

从所述处理后的影像数据的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据。

在一个实施方式中，在确定所述遥感影像中存在云像元的情况下，在确定目标区域的大气程辐射值后，所述方法还包括：

对所述第三暗像元集进行区间估计，以确定基于云像元邻近效应产生的大气程辐射值增量；

利用所述大气程辐射值增量对所述大气程辐射值进行校正。

本申请实施例还提供了一种大气程辐射值的确定装置，包括：

获取模块，用于获取目标区域的遥感影像；

划分模块，用于从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据；

筛选模块，用于对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集；

确定模块，用于根据所述有效暗像元集，确定目标区域的大气程辐射值。

在本申请实施例中，通过先划分出目标区域的遥感影像中的暗像元，再对遥感影像中的暗像元进行多次针对性的筛选，以去掉多种干扰因素的影响，得到表征效果较为准确的有效暗像元集，进而利用有效暗像元集确定大气程辐射值，从而解决了现有方法中存在的确定大气程辐射值误差大、准确度低的技术问题，达到了精确确定大气程辐射值的技术效果。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是根据本申请实施方式提供的大气程辐射值的确定方法的处理流程图；

图2是根据本申请实施方式提供的大气程辐射值的确定装置的组成结构图；

图3是根据本申请实施方式提供的用于确定大气程辐射值的电子设备示意图；

图4是在一个场景示例中应用本申请实施方式提供大气程辐射值的确定方法和装置确定气溶胶光学厚度的流程示意图；

图5是在一个场景示例中应用本申请实施方式提供大气程辐射值的确定方法和装置获得的二维光谱特征空间的示意图(其中，上图为以蓝光波段数据的反射率作为纵轴的二维光谱特征空间，下图为以红光波段数据的反射率作为纵轴的二维光谱特征空间)。

具体实施方式

为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。

考虑到现有方法，具体实施时，通常是基于表观反射率标准差，剔除标准差大于0.02的异质性像元，获取局地均质性像元聚类点集，再利用上述局地均质性像元聚类点集确定相应的大气程辐射值。由于没有仔细地针对暗像元的具体特征，对暗像元进行细致的筛选，导致所获取的局地均质性像元聚类点集往往只能出粗略地表征一些植被、湿润土壤等暗像元的地表情况。例如，上述局地均质性像元聚类点集中常常会存在许多噪点等干扰；又例如，获取的上述局地均质性像元聚类点集还会存在遗漏一些暗像元地表信息的情况。因此，现有方法具体实施时往往会存在确定大气程辐射值误差大、准确度低的技术问题。针对产生上述技术问题的根本原因，本申请考虑可以充分考虑暗像元本身的性质特点以及遥感影像中其他的干扰因素，通过多次针对性的筛选处理，从目标区域的遥感影像中获取表征效果较为准确的有效暗像元集，利用有效暗像元集确定大气程辐射值，从而解决现有方法中存在的确定大气程辐射值误差大、准确度低的技术问题，达到精确确定大气程辐射值的技术效果。

基于上述思考思路，本申请实施例提供了一种大气程辐射值的确定方法。具体请参阅图1所示的根据本申请实施方式提供的大气程辐射值的确定方法的处理流程图。本申请实施例提供的大气程辐射值的确定方法，具体实施时，可以包括以下步骤。

S11：获取目标区域的遥感影像。

在一个实施方式中，上述遥感影像具体可以多种类型的遥感影像中的至少一种：LandSat卫星TM、ETM、OLI影像数据、Terra和Aqua卫星MODIS影像数据等。其中，上述TM、ETM、OLI影像数据具体可以为通过LANDSAT系列卫星(一种陆地卫星)采集的遥感数据，MODIS影像数据具体可以为通过Terra和Aqua卫星(一种对地观测卫星)采集的遥感数据。其中，LandSat卫星TM、ETM、OLI影像数据、Terra和Aqua卫星MODIS影像数据通常为中等分辨率数据。当然，需要说明的是，上述所列举的多种类型的遥感影像数据只是为了更好地说明本申请实施方式，具体实施时，也可以根据具体情况和实施要求，引入除上述所列举的遥感影像数据类型以外的其他遥感影像数据作为目标区域的遥感影像。对此，本申请不作限定。

在本实施方式中，所涉及的目标区域的遥感影像主要是Landsat影像数据。具体的，主要包括的是Landsat影像数据中Landsat7影像数据和Landsat8影像数据。其中，Landsat运行在近极地太阳同步圆形轨道上，Landsat 5、7、8轨道高度均为705km，扫描幅宽185km，单颗卫星重访周期16天，降交点地方时10点左右。基于运行轨道设计特点，集成Landsat 7和8的数据应用，在赤道和高纬度地区的重访周期可缩短到8天以内。Landsat 5的有效载荷为专题制图仪(TM)和多光谱成像仪(MSS)，TM被动感应地表反射的太阳辐射和发射的热辐射，有7个波段的感应器，覆盖了从可见光到红外的不同波长范围，除热红外波段外均为30米分辨率。Landsat 7卫星装备有增强型专题制图仪(ETM+)，与TM传感器相比，ETM+增加分辨率为15米的全色波段(PAN波段)，热红外波段分为低增益和高增益数据，其分辨率从120米提高到60米，因此有更高的准确性。总体上，Landsat 8的设计和特征与Landsat 7基本相同，这使得Landsat 8数据可以和前期的Landsat数据保持很高的一致性和可比性。Landsat 8卫星有陆地成像仪(OLI)和热红外传感器(TIRS)两个传感器。除了保持原有Landsat 7卫星的基本特点外，Landsat 8在原蓝光波段之外新增了1个深蓝波段(分辨率30米)，在原近红外波段与短波红外波段间新增了1个卷云波段(分辨率30米)，在原有热红外波段的光谱范围内设置了2个热红外波段(分辨率100米)。除此之外，Landsat 8还收窄了近红外波段和全色波段的光谱范围，将辐射分辨率从8bit提高到12bit。

S12：从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据。

在本实施方式中，上述像元也称像素点或像元点，即影像单元(picture element)，具体可以是组成数字化影像的最小单元。在遥感数据采集，例如扫描成像时，它是传感器对地面景物进行扫描采样的最小单元。

在本实施方式中，上述暗像元具体可以是指反射率或者辐射亮度相对较小的黑暗像元，其中，该类像元的地物反射率在蓝红可见光和短波红外通道之间呈线性关系。具体的，可以是指遥感影像中用于表征植被、湿润土壤和地形阴影等的像元。上述非暗像元具体可以是指除上述暗像元以外的其他像元。具体的，可以是遥感影像中用表征冰雪、水体、云阴影、薄云及部分建筑用地等的像元。

在本实施方式中，考虑到后续需要利用暗像元的相关数据确定大气程辐射值，因此，需要对上述遥感影像的像元进行分类，划分出暗像元和非暗像元，并对非暗像元进行掩膜处理，得到剔除非暗像元之后的遥感影像，即所述第一影像数据。因此，上述第一影像数据具体是暗像元被掩盖的目标区域的遥感影像。

在本实施方式中，上述掩膜处理具体可以是用选定的图像、图形或物体，对待处理的图像(全部或局部)进行遮挡，来控制图像处理的区域或处理过程。用其中，于覆盖的特定图像或物体通常可以称为掩模或模板。具体的，在光学图像处理中，掩模可以是胶片、滤光片等。在数字图像处理中，掩模可以为二维矩阵数组，有时也可以为多值图像。

在本实施方式中，上述对所述非暗像元进行掩膜处理具体可以是指对所述遥感影像中的非暗像元进行掩盖。

在一个实施方式中，综合考虑目标区域的遥感影像中干扰因素，发现遥感影像中以表征冰雪和水体的像元为代表的非暗像元的干扰相对较大，因此考虑主要通过区分以表征冰雪和水体的像元为代表的非暗像元来划分出暗像元和非暗像元。具体实施时，可以按照以下方法对上述遥感影像的像元进行分类，划分出暗像元和非暗像元：

S1：获取像元的绿光波段数据(B2)、第一短波红外波段数据(B5)和第二短波红外波段数据(B7)；

S2：根据所述像元的绿光波段数据和第一短波红外波段数据，确定像元的归一化积雪指数；

S3：根据所述像元的第二短波红外波段数据和归一化积雪指数，确定所述像元是否为暗像元。

在本实施方式中，需要补充的是，基于上述遥感影像本身的性质特点，通常可以根据遥感影像获取遥感影像中各个像元的7个波段数据。具体的，上述7个波段数据具体可以包括：蓝光波段数据(B1)，该波段位于水体衰减系数最小的部位，对水体的穿透力最大，可以用于判别水深，研究浅海水下地形、水体浑浊度等，进行水系及浅海水域制图；绿光波段数据(B2)，该波段位于绿色植物的反射峰附近，对健康茂盛植物反射敏感，可以用于识别植物类别和评价植物生产力，对水体具有一定的穿透力，可反映出水下地形、沙洲、沿岸沙坝等特征；红光波段数据(B3)，该波段位于叶绿素的主要吸收带，可用于区分植物类型、覆盖度、判断植物生长状况等，此外该波段对裸露地表、植被、岩性、地层、构造、地貌、水文等特征均可提供丰富的植物信息；近红外波段数据(B4)，该波段位于植物的高反射区，反映了大量的植物信息，多用于植物的识别、分类，同时它也位于水体的强吸收区，用于勾绘水体边界，识别与水有关的地质构造、地貌等；第一短波红外波段数据(B5)，该波段位于两个水体吸收带之间，对植物和土壤水分含量敏感，从而提高了区分作物的能力，此外，在该波段上雪比云的反射率低，两者易于区分，B5的信息量大，应用率较高；热红外波段数据(B6)，该波段对地物热量辐射敏感，根据辐射热差异可用于作物与森林区分、水体、岩石等地表特征识别；第二短波红外波段数据(B7)，波长比B5大，是专为地质调查追加的波段，该波段对岩石、特定矿物反应敏感，用于区分主要岩石类型、岩石水热蚀变，探测与交代岩石有关的粘土矿物等；B8为全色波段(Pan)，该波段为Landsat-7新增波段，它覆盖的光谱范围较广，空间分辨率较其他波段高，因而多用于获取地面的几何特征。其中，需要补充的是上述第一短波红外波段数据(B5)和第二短波红外波段数据(B7)都是第一短波红外波段数据，但两者对应的波长范围不同。具体的，第一短波红外波段数据(B5)对应的波长范围具体为1.57-1.65μm，第二短波红外波段数据(B7)对应的波长范围具体为2.11-2.29μm。

在本实施方式中，上述归一化积雪指数(normalized difference snow index)具体是用于区分冰雪的积雪指示参数。

在一个实施方式中，上述根据所述像元的绿光波段数据(B2)和第一短波红外波段数据(B5)，确定像元的归一化积雪指数具体实施时，可以包括按照以下公式确定像元的归一化积雪指数：

上式中，NDSI为像元的归一化积雪指数，B2为像元的绿光波段数据，B5为像元的第一短波红外波段数据。

在一个实施方式中，考虑到除冰雪对暗像元的影响外，还存在水体的影响，为了综合区分出冰雪、水体对应的非暗像元，具体实施时，可以综合利用像元的归一化冰雪指数和像元的第二短波红外波段数据(B7)作为指示参数，确定像元是否为暗像元。具体实施时，可以将满足以下规则的像元确定为暗像元：

NDSI≤0且0.01≤B7≤0.25

上式中，NDSI为像元的归一化积雪指数，B7为像元的第二短波红外波段数据。

在一个实施方式中，在获得了剔除非暗像元的第一影像数据后，为了后续能更加精确地对第一影像数据中的暗像元进行筛选，具体实施时，可以将所述多个第一影像数据划分为多个图像区块；相应的，所述对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集，具体可以包括：对所述第一影像数据中的多个图像区块中的暗像元分别进行针对性的多次筛选。在本实施方式中，具体实施时，可以按每个图像区块为101×101的像元范围，即以101×101的像元范围(对应3km×3km的区域范围)作为一个图像区块，将所述第一影像划分成多个不同的图像区块，以便后续可以针对各个图像区块分别建立精细的二维光谱特征空间，以便可以更加准确地进行对应的暗像元筛选，进而确定出更加精确的大气程辐射值。

S13：对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集。

在本实施方式中，具体实施时，可以通过基于暗像元本身的特征性质，对第一影像中的暗像元进行多次针对性的筛选，以排除其他干扰因素的，得到表征效果较好的有效暗像元集，以便后续可以更加精确地确定出大气程辐射值。

在一个实施方式中，上述对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集，具体实施时，可以包括以下内容：

S13-1：根据所述第一影像数据，建立二维光谱特征空间；其中，所述二维光谱特征空间以像元的第二短波红外波段数据为横轴，以像元的蓝光波段数据的反射率或者红光波段数据的反射率为纵轴。

在本实施方式中，上述二维光谱特征空间(spectrum feature space)，又称散点图，是一种基于二维遥感影像数据中的不同波段影像数据所构成的测度空间。具体的，上述二维光谱特征空间具体可以以像元的第二短波红外波段数据为横轴，可以以像元的蓝光波段数据的反射率为纵轴，也可以以红光波段数据的反射率为纵轴。其中，第一影像数据中的暗像元根据各个暗像元的第二短波红外波段数据、蓝光波段数据的反射率或者红光波段数据的反射率的具体数值，分布在上述二维光谱特征空间中。

在本实施方式中，需要说明的是，对于划分了为多个不同图像区块的第一影像数据，上述根据所述第一影像数据，建立二维光谱特征空间，具体可以包括：根据所述第一影像数据中的多个图像区块，分别建立针对所述多个图像区块中各个图像区块的二维光谱特征空间。相应的，后续对暗像元的多次筛选，具体可以是分别在各个图像区块的二维光谱特征空间中，对各个图像区块的暗像元进行针对性的多次筛选。从而可以提高筛选暗像元的精度，以便更加准确地确定出各个区块图像中大气程辐射值。

S13-2：通过坐标分区法，对所述二维光谱特征空间中的暗像元进行第一次筛选，得到第一暗像元集。

在一个实施方式中，上述通过坐标分区法，对所述二维光谱特征空间中的暗像元进行第一次筛选，可以较为有效地降低第一影像数据中的地表反射率波动产生的影响。具体实施时，上述步骤可以包括以下内容：

S1：根据像元的横轴坐标数值，将所述二维光谱特征空间中的暗像元分别划分为多个子区间；

S2：分别提取子区间中纵轴坐标数值最小的暗像元，作为第一暗像元集中的元素，以获取所述第一暗像元集。

在本实施方式中，具体实施时，可以在上述二维光谱特征空间中沿横坐标轴划分出多个子区间，例如，可以按照不同的横轴坐标值划分多个子区间。具体的，可以将横坐标值相同纵坐标值不同的多个像元点作为一个子区间。再从上述各个子区间的暗像元中分别提取纵坐标值最小的像元点组成集合，即可以得到上述第一暗像元集。上述筛选的依据具体是：在二维光谱特征空间中，对于足够数量的呈线性分布的暗像元的散点集合，其下边界(即同一横坐标值所对应的纵坐标值最小的点)往往可以代表最暗像元集合的线性分布，是确定基于暗像元的拟合方程的重要依据。而处于下边界以上的暗像元的散点所表征的像元相对较亮，有可能会掩盖掉微弱的气溶胶信号。因此，在本实施方式中，为了降低地表反射率波动的影响，在所建立的二维光谱空间中，通过基于横坐标划分多个子区间，再从各个子区间中提取纵坐标最小所对应的像元，即从各个子区间中分别获取蓝光波段数据或红光波段数据的反射率最小的像元点，构成初始的最暗像元点集，即上述第一暗像元集。

S13-3：对所述第一暗像元集进行第二次筛选，以去除所述第一暗像元集中的噪声散点，得到第二暗像元集。

在本实施方式中，具体实施时，可以根据暗像元的第二短波红外波段数据的反射率，将第一暗像元集划分为多个子集，进行第二次筛选，以去除上述第一暗像元集中噪声散点(像元)的干扰，得到精度相对第一暗像元集更高的第二暗像元集。

在一个实施方式中，上述对所述第一暗像元集进行第二次筛选，具体实施时，可以包括以下内容：

S1：获取所述第一暗像元集中暗像元的第二短波红外波段数据的反射率；

S2：根据所述暗像元的第二短波红外波段数据的反射率，将所述第一暗像元集中的暗像元划分为多个子集；

S3：分别确定子集的最小二乘相关系数、子集的一元线性拟合方程的斜率和子集的一元线性拟合方程的截距；

S4：根据所述子集的最小二乘相关系数、所述子集的一元线性拟合方程的斜率和所述子集的一元线性拟合方程的截距，从所述多个子集中筛选出符合要求的子集作为所述第二暗像元集。

在一个实施方式中，具体实施时，可以按照以下公式分别计算第一暗像元中多个子集的各个子集的最小二乘相关系数：

按照以下公式，分别确定子集的最小二乘相关系数：

上式中，rm为编号为m的子集的最小二乘相关系数，Xi为编号为m的子集编号为i的像元的第二短波红外波段数据的反射率，为编号为m的子集中像元的第二短波红外波段数据的反射率的平均值，Yi为编号为m的子集编号为i的像元的蓝光波段数据的反射率或者红光波段数据的反射率，为编号为m的子集中像元的蓝光波段数据的反射率平均值或者红光波段数据的反射率的平均值。

在一个实施方式中，上述符合要求的子集作具体可以为所述第二暗像元集为子集的一元线性拟合方程的斜率小于1，子集的一元线性拟合方程的截距大于0的子集中子集的最小二乘相关系数最大的子集。具体的，例如，可以从第二暗像元集中筛选出多个满足子集的一元线性拟合方程的斜率小于1的子集作为第一筛选集；再从第一筛选集中筛选出多个满足子集的一元线性拟合方程的截距大于0的子集作为第二筛选集；最后从第三筛选集中筛选出最小二乘相关系数最大的子集作为上述符合要求的子集。

在本实施方式中，考虑到噪声散点的影响，具体实施时，可以按照以下方式先将第一暗像元集划分为多个子集，再对多个子集中进行筛选，选出符合要求的子集作为上述第二暗像元集。具体的，例如，可以以八分位点将第一暗像元集中暗像元的第二短波红外波段数据的反射率划分为以下几个子集：0％～25％、0％～37.5％、0％～50％、0％～62.5％、0％～75％、0％～87.5％、0％～100％、12.5％～37.5％、12.5％～50％、12.5％～62.5％、12.5％～75％、12.5％～87.5％、12.5％～100％。分别根据各个子集中的暗像元，计算对应子集的最小二乘相关系数、子集的一元线性拟合方程的斜率和子集的一元线性拟合方程的截距。在同时满足子集的一元线性拟合方程的斜率小于1，子集的一元线性拟合方程的截距大于0的多个子集中筛选出子集的最小二乘相关系数最大的子集作为上述第二暗像元集。

S13-4：通过多次迭代拟合，对所述第二暗像元集进行第三次筛选，以剔除离散的异常点，得到第三暗像元集，并将所述第三暗像元集作为所述有效暗像元集。

在本实施方式中，具体实施时，可以通过多次迭代拟合得到多个拟合方程，对所述第二暗像元集进行第三次筛选，以剔除离散的异常点，得到相对第二暗像元集精度更高的第三暗像元集。

在一个实施方式中，上述通过多次迭代拟合，对所述第二暗像元集进行第三次筛选，以剔除离散的异常点，得到第三暗像元集，具体实施时，可以包括以下内容：通过多次迭代拟合，将所述第二暗像元集中与所述第二暗像元集的一元线性拟合方程的垂直距离大于阈值距离的暗像元剔除。其中，上述阈值距离具体可以根据所要求的大气程辐射值的精度要求设定。

在本实施方式中，上述第二暗像元集中通常可能还会保留有一些离散的异常点(像元)，例如传感器噪声、湿地、云阴影像元等。上述离散的异常点在后续确定大气程辐射值的过程中也会对所确定的大气程辐射值的准确度造成影响。具体的，例如，可以首先计算第二暗像元集的一元线性拟合方程和相关系数。然后在二维光谱特征空间中按照第二短波红外波段数据的反射率(即横坐标波段反射率)进行排序，获取第二暗像元集中像元总数的一半点集(例如，数量为m)作为迭代子集，并计算迭代子集的一元线性拟合方程和相关系数。比较两种拟合方程和相关系数，判断迭代子集的拟合方程是否与第二暗像元集的拟合方程是否一致，如果不一致，则将垂直偏差最大的点从第二暗像元集中去除。按照相同的比较方法，依次去除第二暗像元集中垂直偏差较大的像元点(即离散的异常点)。之后，得到新的第二暗像元集，再计算新的第二暗像元集的一元线性拟合方程和相关系数，此时，第二次迭代子集的像元数量可以设置为m+1，重复执行以上相同的迭代过程，进行多次的迭代拟合，以逐次去除第二暗像元集中的离散的异常点。直到新的迭代子集包含了新的第二暗像元集中全部的散点(即暗像元)后迭代拟合结束。最后，经过多次迭代拟合后得到的新的第二暗像元集中剩余的像元点便是通过多次针对性的筛选处理最终获取的表征效果较好、干扰因素的影响较少的暗像元点，即得到了上述的第三暗像元集。进而可以将上述第三暗像元集确定为有效暗像元集，以便后续用于确定大气程辐射值。

S14：根据所述有效暗像元集，确定目标区域的大气程辐射值。

在一个实施方式中，上述根据所述有效暗像元集，确定目标区域的大气程辐射值，具体实施时，可以包括以下内容：

S14-1：在所述二维光谱特征空间中，根据所述有效暗像元集中的暗像元，建立拟合回归线；

S14-2：将所述拟合回归线与所述二维光谱特征空间的纵轴的截距确定为所述大气程辐射值。

在本申请实施例中，相较于现有技术，由于充分考虑了暗像元自身的性质特征以及遥感影像中的干扰因素的影响，先划分出了目标区域的遥感影像中的暗像元，再对遥感影像中的暗像元进行多次针对性的筛选，以去掉多种干扰因素的影响，得到表征效果较为准确的有效暗像元集，并利用有效暗像元集确定大气程辐射值，从而解决了现有方法中存在的确定大气程辐射值误差大、准确度低的技术问题，达到了精确确定大气程辐射值的技术效果。

在一个实施方式中，考虑到所获取的目标区域的遥感影像中可能会存在云，而云在遥感影像中会对后续确定大气程辐射值造成干扰，因此，为了能够事先去除云的影响具体实施，在获取目标区域的遥感影像后，所述方法具体还可以包括以下内容：

S1：对所述遥感影像进行云像元检测，确定所述遥感影像中是否存在云像元；

S2：在所述遥感影像中存在云像元的情况下，对所述遥感影像中的云像元进行掩膜处理，得到处理后的影像数据。

相应的，所述从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据，具体可以包括：从所述处理后的影像数据的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据。

在本实施方式中，上述云像元具体可以认为是遥感影像中用于表征云的像元。具体实施时，可以通过检测出目标区域的遥感影像中的云像元，先对所述遥感影像中的云像元进行掩膜处理，以消除云本身对遥感影像中暗像元的干扰。

在一个实施方式中，上述对所述遥感影像进行云像元检测，确定所述遥感影像中是否存在云像元，具体实施时，可以包括以下内容：

S1：根据所述遥感影像，获取遥感影像中多个像元中各个像元的绿光波段数据(B2)、红光波段数据(B3)、第一短波红外波段数据(B5)、热红外波段数据(B6)；

S2：将所述像元的红光波段数据、第一短波红外波段数据、热红外波段数据分别进行HSV色彩空间变换，得到各个像元的色调数据、饱和度数据、明度数据；

S3：根据像元的绿光波段数据、红光波段数据、第一短波红外波段数据和像元的色调数据、饱和度数据、明度数据，确定像元是否为云像元。

在本实施方式中，上述HSV(Hue、Saturation、Value)是一种对对RGB色彩空间中点的关系的表述。具体的，上述HSV模型通常用于计算机图形应用中。

在一个实施方式中，上述根据像元的绿光波段数据、红光波段数据、第一短波红外波段数据和像元的色调数据、饱和度数据、明度数据，确定像元是否为云像元，具体可以包括：将满足以下规则中的至少一种的像元确定为云像元：

0.07≤B3≤0.08，且S>0.7，V>0.6，30≤H≤90，或者，S>0.7，V>0.6，270≤H≤330；

或者，

B3>0.08，-0.25<(B2-B5)/(B2+B5)<0.7，B6<300K，(1-B5)*B6≥225，B5≥0.08，且30≤H≤90，S>0.7，V>0.6，或者，270≤H≤330，S>0.7，V>0.6；

或者，

B3>0.08，-0.25<(B2-B5)/(B2+B5)<0.7，B6<300K，(1-B5)*B6<225，且30≤H≤90，S>0.7，V>0.6，或者，270≤H≤330，S>0.7，V>0.6；

其中，B2为像元的绿光波段数据，B3为像元的红光波段数据，B5为像元的第一短波红外波段数据，B6为像元的热红外波段数据，H为像元的色调数据，S为像元的饱和度数据，V为像元的明度数据。

在一个实施方式中，在所述目标区域的遥感影像为OLI影像数据(即一种Landsat8的遥感影像数据)的情况下，所述对所述遥感影像进行云像元检测，确定所述遥感影像中的云像元，具体还可以包括以下内容：

S1：从所述目标区域的遥感影像中获取质量评估波段(Quality Assessment Band)；

S2：将所述质量评估波段转换为16位的二进制数据；

S3：从所述转换后的16位的二进制数据中提取置信度信息；

S4：根据所述置信度信息，从所述目标区域的遥感影像的像元中提取第一云像元；

S5：对所述第一云像元的红光波段数据、第一短波红外波段数据、热红外波段数据分别进行HSV色彩空间变换，得到第一云像元的色调数据、饱和度数据、明度数据；

S6：根据所述第一云像元的色调数据、饱和度数据、明度数据，从所述第一云像元中筛选出第二云像元，作为所述云像元。

在一个实施方式中，在确定所述遥感影像中存在云像元的情况下，考虑到虽然通过云掩膜处理可以掩盖掉云自身直接造成的干扰，但是由于云像元在去除之前已经由于自身的光学特征对邻近的像元形成了影响，即云邻近效应。这种影响实质上反应在了邻近云像元的非云像元上，这种云邻近效应也会对后续确定大气程辐射值造成干扰。因此，在确定目标区域的大气程辐射值后，所述方法具体还可以包括以下内容：

S1：对所述第三暗像元集进行区间估计，以确定基于云像元邻近效应产生的大气程辐射值增量；

S2：利用所述大气程辐射值增量对所述大气程辐射值进行校正。

在本实施方式中，由于考虑到了云像元的云邻近效应对大气程辐射值的影响，根据具体影响对所确定的大气程辐射值进行校正，从而可以得到更加准确的大气程辐射值。

在本实施方式中，需要补充的是，对于多云天气，在二维光谱空间中，由于云邻近效应往往会产生程辐射率增量使得卷积窗口中的暗像元集的拟合回归线较真实情况沿纵轴方向有所升高。因此，可以根据云影响下暗像元拟合回归线的分布变化特征，采用区间估计的处理方法来消除上述云邻近效应对的大气程辐射值的影响。

在本实施方式中，上述区间估计具体是一种统计学中的参数估计形式。具体的，可以通过从总体中抽取的样本，根据一定的正确度与精确度的要求，构造出适当的区间，以作为总体的分布参数(或参数的函数)的真值所在范围的估计。具体实施时，当成对的两个变量数据分布大体上呈直线趋势时，可以运用合适的参数估计方法，求出一元线性回归模型，然后根据自变量与因变量之间的关系，对于自变量x的一个给定值xp，根据估计的回归方程可以得到因变量y的个别值的估计区间，这一区间称为预测区间。再划定上述预测区间的两个数值分别作为预测区间的下限和预测区间的上限。在置信水平固定的情况下，通常样本量越多，预测区间越宽。在样本量相同的情况下，通常置信水平越高，预测区间越宽。在本实施方式中，具体可以将置信水平设置为99.9％，并取预测区间下限作为消减云影响之后的暗像元拟合回归线，继而可以确定出基于云像元邻近效应产生的大气程辐射值增量，对大气程辐射值进行校正，得到不受云邻近效应影响的，更为精确的大气程辐射值。

在一个实施方式中，在确定大气程辐射值后，为了能够准确地确定目标区域的气溶胶光学厚度(AOD)具体实施时，可以包括以下内容：将地表反射率设为0；根据所述大气程辐射值进行气溶胶光学厚度反演，以确定目标区域的气溶胶光学厚度。

在一个实施方式中，在确定出目标区域的气溶胶光学厚度后，为了能够对目标区域进行相应的大气颗粒物污染评价，以便后续可以更有针对性地对目标区域的大气进行相应保护，具体实施时，可以按照以下方式评价目标区域的大气污染情况：

S1：根据所述目标区域的气溶胶光学厚度，确定目标区域中的颗粒物含量；

S2：根据所述目标区域中的颗粒物含量，对目标区域进行大气污染评价。

从以上的描述中，可以看出，本申请实施例提供的大气程辐射的确定方法，由于充分考虑了暗像元自身的性质特征以及遥感影像中的干扰因素的影响，先划分出了目标区域的遥感影像中的暗像元，再对遥感影像中的暗像元进行多次针对性的筛选，以去掉多种干扰因素的影响，得到表征效果较为准确的有效暗像元集，并利用有效暗像元集确定大气程辐射值，从而解决了现有方法中存在的确定大气程辐射值误差大、准确度低的技术问题，达到了精确确定大气程辐射值的技术效果；又通过云像元检测，确定遥感影像中的云像元，对云像元进行云掩膜处理，消除了干扰，提高了所确定的大气程辐射值的精度；还通过确定基于云像元邻近效应产生的大气程辐射值增量，并利用上述大气程辐射值增量对所确定的目标区域的大气程辐射值进行针对性校正，进一步提高了大气程辐射值的精度。

基于同一发明构思，本发明实施例中还提供了一种大气程辐射值的确定装置，如下面的实施例所述。由于大气程辐射值的确定装置解决问题的原理与大气程辐射值的确定方法相似，因此大气程辐射值的确定装置的实施可以参见方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。请参阅图2，是本申请实施例提供的大气程辐射值的确定装置的一种组成结构图，该装置具体可以包括：获取模块21、划分模块22、筛选模块23和确定模块24，下面对该结构进行具体说明。

获取模块21，具体可以用于获取目标区域的遥感影像；

划分模块22，具体可以用于从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据；

筛选模块23，具体可以用于对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集；

确定模块24，具体可以用于根据所述有效暗像元集，确定目标区域的大气程辐射值。

在一个实施方式中，上述装置具体还可以包括评价模块，其中，上述评价模块具体实施时可以按照以下方式执行：将地表反射率设为0，根据所述大气程辐射值进行气溶胶光学厚度反演，以确定目标区域中的气溶胶光学厚度；根据所述目标区域中的气溶胶光学厚度，确定目标区域中的颗粒物含量；根据所述目标区域中的颗粒物含量，对目标区域进行大气污染评价。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

需要说明的是，上述实施方式阐明的系统、装置、模块或单元，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。为了描述的方便，在本说明书中，描述以上装置时以功能分为各种单元分别描述。当然，在实施本申请时可以把各单元的功能在同一个或多个软件和/或硬件中实现。

此外，在本说明书中，诸如第一和第二这样的形容词仅可以用于将一个元素或动作与另一元素或动作进行区分，而不必要求或暗示任何实际的这种关系或顺序。在环境允许的情况下，参照元素或部件或步骤(等)不应解释为局限于仅元素、部件、或步骤中的一个，而可以是元素、部件、或步骤中的一个或多个等。

从以上的描述中，可以看出，本申请实施例提供的大气程辐射值的确定装置，由于充分考虑了暗像元自身的性质特征以及遥感影像中的干扰因素的影响，通过划分模块先划分出了目标区域的遥感影像中的暗像元，再通过筛选模块对遥感影像中的暗像元进行多次针对性的筛选，以去掉多种干扰因素的影响，得到表征效果较为准确的有效暗像元集，并利用确定模块根据有效暗像元集确定大气程辐射值，从而解决了现有方法中存在的确定大气程辐射值误差大、准确度低的技术问题，达到了精确确定大气程辐射值的技术效果；又通过云像元检测模块，确定遥感影像中的云像元，通过云掩膜处理模块对云像元进行云掩膜处理，消除了干扰，提高了所确定的大气程辐射值的精度；还通过校正模块确定基于云像元邻近效应产生的大气程辐射值增量，并利用上述大气程辐射值增量对所确定的目标区域的大气程辐射值进行针对性校正，进一步提高了大气程辐射值的精度。

本申请实施方式还提供了一种电子设备，具体可以参阅图3所示的根据本申请实施方式提供的用于信息投放的电子设备示意图，所述电子设备具体可以包括输入设备31、处理器32、存储器33。其中，所述输入设备31具体实施时，可以用于输入目标区域的遥感影像。所述处理器32具体实施时，可以用于从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据；对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集；根据所述有效暗像元集，确定目标区域的大气程辐射值。所述存储器33具体可以用于存储所所输入的目标区域的遥感影像，以及处理器32运行过程中生成的第一影像数据、有效暗像元集等数据。

需要说明的是，对于上述电子设备而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处可以参见方法实施例的部分说明。本申请在此不作赘述。

在本实施方式中，所述输入设备具体可以是用户和计算机系统之间进行信息交换的主要装置之一。所述输入设备可以包括键盘、鼠标、摄像头、扫描仪、光笔、手写输入板、语音输入装置等；输入设备用于把原始数据和处理这些数的程序输入到计算机中。所述输入设备还可以获取接收其他模块、单元、设备传输过来的数据。所述处理器可以按任何适当的方式实现。例如，处理器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(Application Specific Integrated Circuit，ASIC)、可编程逻辑控制器和嵌入微控制器的形式等等。所述存储器具体可以是现代信息技术中用于保存信息的记忆设备。所述存储器包括很多层次，在数字系统中，只要能保存二进制数据的都可以是存储器；在集成电路中，一个没有实物形式的具有存储功能的电路也叫存储器，如RAM、FIFO等；在系统中，具有实物形式的存储设备也叫存储器，如内存条、TF卡等。

本说申请实施方式中还提供了一种计算机存储介质，所述计算机存储介质存储有计算机程序指令，在所述计算机程序指令被执行时实现：获取目标区域的遥感影像；从所述遥感影像的像元中划分出暗像元和非暗像元，并对所述非暗像元进行掩膜处理，得到第一影像数据；对所述第一影像数据中的暗像元进行多次筛选，得到有效暗像元集；根据所述有效暗像元集，确定目标区域的大气程辐射值。

在本实施方式中，上述存储介质包括但不限于随机存取存储器(Random Access Memory,RAM)、只读存储器(Read-Only Memory,ROM)、缓存(Cache)、硬盘(Hard Disk Drive,HDD)或者存储卡(Memory Card)。所述存储器可以用于存储计算机程序指令。网络通信单元可以是依照通信协议规定的标准设置的，用于进行网络连接通信的接口。

在本实施方式中，该计算机存储介质存储的程序指令具体实现的功能和效果，可以与其它实施方式对照解释，在此不再赘述。

在一个具体的实施场景示例中，应用本申请提供大气程辐射值的确定方法和装置确定某地区的大气程辐射值，具体实施过程，可以结合图4所示的在一个场景示例中应用本申请实施方式提供大气程辐射值的确定方法和装置确定气溶胶光学厚度的流程示意图，参阅以下内容执行。

S1：对获取的目标区域的遥感影像进行云像元检测，并对遥感影像中的云像元进行云掩膜处理。

在本实施方式中，具体实施时，可以按照以下规则，确定目标区域的遥感影像中的像元是否为云像元：

规则①.0.07≤B3≤0.08 and 30≤H≤90 and S>0.7 and V>0.6，为云像元；

规则②.0.07≤B3≤0.08 and 270≤H≤330 and S>0.7 and V>0.6，为云像元；

规则③.B3>0.08 and -0.25<(B2-B5)/(B2+B5)<0.7 and B6<300K and(1-B5)*B6≥225 and B5≥0.08 and 30≤H≤90 and S>0.7 and V>0.6，为云像元；

规则④.B3>0.08 and -0.25<(B2-B5)/(B2+B5)<0.7 and B6<300K and(1-B5)*B6≥225 and B5≥0.08 and 270≤H≤330 and S>0.7 and V>0.6，为云像元；

规则⑤.B3>0.08 and -0.25<(B2-B5)/(B2+B5)<0.7 and B6<300K and(1-B5)*B6<225 and 30≤H≤90 and S>0.7 and V>0.6，为云像元；

规则⑥.B3>0.08 and -0.25<(B2-B5)/(B2+B5)<0.7 and B6<300K and(1-B5)*B6<225 and 270≤H≤330 and S>0.7 and V>0.6，为云像元；

在本实施方式中，色调、饱和度和明度参数的获取需要对RGB的彩色影像进行HSV正变换。通过大量试验发现：在TM波段B3(对应R分量)、B5(对应G分量)、B6(对应B分量)假彩色合成后的影像中，云在RGB色彩空间中为紫色或黄色，即在HSV空间中色调位于270°～330°或30°～90°之间；而雪基本为红色，陆地为绿色，水体和阴影为黑色，云与其他地物色彩差别很大，从而可以很容易地进行区分。为了增加算法的鲁棒性，特别是增强云与阴影、裸露地表的区分，可以在色调阈值限定的基础上，又添加了饱和度和明度的阈值限定。具体的，例如可以以云与其他地物的自身光谱特性和光谱差异为理论基础，综合ACCA(Automatic Cloud Cover Assessment)算法和HSV色彩空间的应用，应用到Landsat 7 ETM+的遥感影像中可以很好地将云、雪、裸露地表、植被和水体区分开来，具有简单可行、客观性强、精度高和计算速度快等特点，并且能够适用于不同的下垫面和季节，云检测效果理想。实验证明，该算法同样适用于Landsat 5 TM遥感影像中的云检测。

S2：对云掩膜处理后的遥感影像进行暗像元检测，以从处理后的遥感影像中划分出非暗像元和暗像元，并对非暗像元进行掩膜处理。

在本实施方式中，通常蓝红可见光和短波红外通道之间的线性关系在冰雪、湿地、水体等地表类型上表现较差，因此，在云掩膜处理之后，还需要将这些非暗像元进行检测，以便进一步进行掩膜处理。

在本实施方式中，为了能够快速地确定水体和冰雪等非暗像元，综合上述不同非暗像元的具体特点，通过大量试验，并参考MODIS(中分辨率成像光谱仪)气溶胶业务化反演算法中暗像元的选择方法，以Landsat 8 OLI影像数据为例，具体可以按照以下规则确定像元是否为暗像元：NDSI≤0 and 0.01≤B7≤0.25，

其中，B7是OLI数据的短波红外波段(即第二短波红外波段数据)(2.11-2.29μm)的反射率。针对非云像元应用以上规则，结合掩膜处理，可以有效地剔除雪、水体、云阴影、薄云及部分建筑用地的影响。

S3：建立二维光谱特征空间。

在本实施方式中，基于二维光谱特征空间的大气程辐射值自动提取，具体实施时，首先需要将遥感影像转换到光谱特征空间中。以OLI Landsat 8影像数据为例，可以按照101×101卷积窗口(约3km×3km)分块统计窗口内所有暗像元的蓝(0.45-0.51μm)、红(0.64-0.67μm)、短波红外(2.11-2.29μm)波段表观反射率。进而可以根据各暗像元在横轴(短波红外)和纵轴(蓝或红)2个波段的反射率值，绘制其在二维光谱空间内的分布特征，具体可以参阅图5所示的在一个场景示例中应用本申请实施方式提供大气程辐射值的确定方法和装置获得的二维光谱特征空间的示意图(上图为以蓝光波段数据的(反观)反射率作为纵轴的二维光谱特征空间，下图为以红光波段数据的(反观)反射率作为纵轴的二维光谱特征空间)。图中所示的像元散点包含有本实施方式所涉及到的多种类型的像元散点，其中，初始暗像元即上述掩膜处理后划分出的暗像元。

S4：初始最暗像元集合提取(即通过第一次筛选获得第一暗像元集)。

在本实施方式中，通过上述的步骤S1至S3中一系列定量化的规则，能够区分出绝大部分的云、冰雪、水体及其他暗像元。但是，由于各规则阈值很难实现自适应设定，再加上TM、ETM+数据缺少卷云波段，这使得构建的二维光谱空间中，有可能保留了非暗像元如薄云、裸土、浅水等像元噪声散点造成的影响，导致获得的暗像元并不一定完全准确。又考虑到，在晴朗无云的天气下，暗像元检测结果能保持一定精度，但卷积窗口内同一类型地物地表反射率的波动，也会造成可见光和短波红外通道表观反射率之间相互关系的波动。

在本实施方式中，具体实施时，对于足够数量的呈线性分布的暗像元散点集合，其下边界往往可以代表最暗像元集合的线性分布，即横坐标所对应的纵坐标值最小的点。以此可以作为确定暗像元线的重要依据。相对的，处于下边界以上的散点所表征的地表情况的像元通常相对较亮，有可能掩盖微弱的气溶胶信号。为了降低地表反射率波动的影响，在构建的二维光谱空间中，可以通过横坐标分组的方法，将所有的像元点按其短波红外波段反射率归入不同的组中，在各组中分别获取蓝或红波段反射率最小的像元点，构成初始的最暗像元点集(即第一暗像元集)。

S5：通过自适应区间选取提取暗像元(即通过第二次筛选获得第二暗像元集)。

在本实施方式中，考虑到在二维光谱空间横坐标的全区间内，各组纵坐标最小值所对应的像元点并不一定都是暗像元，也有可能存在薄云、裸土、浅水等像元噪声散点。基于上述考虑，可以引入一种简单的自适应方法，具体的，划分横坐标区间，从而缩小初始暗像元点集，完成暗像元的初步筛选。

具体实施时，可以包括以下内容：首先，以八分位点对初始暗像元点集横轴波段反射率范围进行划分，然后分别确定初始暗像元点集的0％～25％、0％～37.5％、0％～50％、0％～62.5％、0％～75％、0％～87.5％、0％～100％、12.5％～37.5％、12.5％～50％、12.5％～62.5％、12.5％～75％、12.5％～87.5％、12.5％～100％子集，并计算各子集的最小二乘相关系数r、一元线性拟合方程的斜率和截距。其中，上述各子集的最小二乘相关系数可以按照以下公式确定：

式中，Xi和Yi分别为该像元点横坐标波段和纵坐标波段的反射率值，和分别为Xi和Yi的平均值，i为各子集内散点数量。

最后，可以选取相关系数最大并且拟合直线的斜率小于1，截距大于0的子集，作为暗像元点的有效子集(即第二暗像元集)，其横坐标区间也将确定为暗像元的横坐标初始区间。

S6：通过迭代筛选暗像元点(即通过第三次筛选获得第三暗像元集)。

在本实施方式中，具体实施时，对于有效子集中的像元点，还需进一步筛选，以剔除一些离散的异常值，比如传感器噪声、湿地、云阴影像元等。因此，可以对有效子集中的像元点进行循环迭代，每次循环都将垂直偏差最大的点从子集中去除。具体迭代过程可以是，首先计算有效子集的一元线性拟合方程和相关系数，然后按照横坐标波段反射率排序，获取有效子集中像元总数的一半点集(数量为m)作为迭代子集，并计算迭代子集的一元线性拟合方程。判断迭代子集拟合方程是否与有效子集一致，如果不一致，则将垂直偏差最大的点从子集中去除。按照相同的判断过程，去除剩余点集中垂直偏差最大点。之后，计算新的有效子集的一元线性拟合方程和相关系数，迭代子集像元数量变为m+1，重复执行以上相同的迭代过程，直到迭代子集包含了新的有效子集全部散点后循环结束。最后，有效子集中剩余的像元点便是自动算法最终获取的暗像元点即第三暗像元集，具体可以参阅图5，图中的筛选暗像元即为本实施方式获得的暗像元。然后可以通过最小二乘拟合，便可得到暗像元线的方程(即拟合回归线)，确定该暗像元线与纵轴的截距作为大气程辐射值。

S7：消减云影响(即对所确定的大气程辐射值进行校正)。

在本实施方式中，对于多云天气，在二维光谱空间中，由于3D云邻近效应会产生的程辐射率增量使得卷积窗口暗像元线较真实情况沿纵轴方向有所升高。根据云影响下暗像元线的分布变化特征，本实施方式提出可以采用预测区间估计的处理方法来消减多次散射产生的程辐射率误差影响，得到更加准确的暗像元，进而可以更加准确地确定出大气程辐射值。具体可以参阅图5，图中的预测区间估计即为消除云影响后得到暗像元。

具体的，在统计学中，上述区间估计是参数估计的一种形式。具体可以通过从总体中抽取的样本，根据一定的正确度与精确度的要求，构造出适当的区间，以作为总体的分布参数(或参数的函数)的真值所在范围的估计。当成对的2个变量数据分布大体上呈直线趋势时，运用合适的参数估计方法，求出一元线性回归模型，然后根据自变量与因变量之间的关系，对于自变量x的一个给定值xp，根据估计的回归方程可以得到因变量y的个别值的估计区间，这一区间称为预测区间。因变量个别值的方差的估计值为：

再在1-α置信水平下因变量个别值的预测区间为：

其中，sy为估计标准误差，具体可以表示为：

其中，Sind为因变量个别值的方差的估计值，sy为估计标准误差，xp或yp为一个给定值，为均值，i＝1，2，3...n，1-α为置信水平即总体参数值落在样本统计值某一区内的概率，t为用t分布理论来推论差异发生的概率即t检验统计量，和线性回归模型的截距和斜率，n为样本数量。

进而可以划定预测区间的两个数值分别称为预测下限和预测上限。通常在置信水平固定的情况下，样本量越多，预测区间会越宽。而在样本量相同的情况下，置信水平越高，预测区间会越宽。在本实施方式中所得到的卷积窗口尺度大大变小，并且暗像元线的提取经过了严格的、客观的、自动迭代筛选处理过程，使得最终保留的有效暗像元数量有限，因此，本算法中将置信水平设置为99.9％，并取预测区间下限作为消减云影响之后的暗像元线。

S8：AOD反演(即根据所确定的大气程辐射值确定气溶胶光学厚度)。

在本实施方式中，具体实施时可以将地表反射率设为零，利用6S辐射传输模型构建0.49μm或0.66μm程辐射率随气溶胶光学厚度(从0.0-2.0变化)和几何角度(由卫星观测几何和时间决定)变化的查找表。进而可以通过遍历气溶胶光学厚度和几何角度的组合，利用暗像元线拟合方程确定的截距，在查找表中检索与观测程辐射率最邻近的两个程辐射率值，然后根据这两个程辐射率对应的气溶胶光学厚度。再进行插值运算，实现云影响下卷积窗口中心像元的气溶胶光学厚度反演。最后通过移动卷积窗口，重复以上S3-S8的各步骤，直到完成整幅影像各像元的AOD反演。

在本实施方式中，需要补充的是上述气溶胶光学厚度反演的精度通常会受到云和暗像元检测阈值等的影响，导致可能造成不同地区影像局部气溶胶光学厚度反演异常值出现。通常在一定范围内大气程辐射和气溶胶分布相对稳定，因此利用相邻像元AOD的空间自相关性，可以通过空间一致性检验将异常值剔除，而后再采用9像元×9像元的距离加权平均的滤波方法进行平滑处理，从而进一步内插部分非暗像元点的监测值及抑制异常点，输出最终反演结果，即较为准确的气溶胶光学厚度。

通过上述场景示例，验证了本申请实施例提供的大气程辐射值的确定方法和装置，由于充分考虑了暗像元自身的性质特征以及遥感影像中的干扰因素的影响，通过先划分出了目标区域的遥感影像中的暗像元，再对遥感影像中的暗像元进行多次针对性的筛选，以去掉多种干扰因素的影响，得到表征效果较为准确的有效暗像元集，并利用有效暗像元集确定大气程辐射值，确实解决了现有方法中存在的确定大气程辐射值误差大、准确度低的技术问题，达到了精确确定大气程辐射值的技术效果。

尽管本申请内容中提到不同的具体实施例，但是，本申请并不局限于必须是行业标准或实施例所描述的情况等，某些行业标准或者使用自定义方式或实施例描述的实施基础上略加修改后的实施方案也可以实现上述实施例相同、等同或相近、或变形后可预料的实施效果。应用这些修改或变形后的数据获取、处理、输出、判断方式等的实施例，仍然可以属于本申请的可选实施方案范围之内。

虽然本申请提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的装置或客户端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境，甚至为分布式数据处理环境)。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备中还存在另外的相同或等同要素。

上述实施例阐明的装置或模块等，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本申请时可以把各模块的功能在同一个或多个软件和/或硬件中实现，也可以将实现同一功能的模块由多个子模块的组合实现等。以上所描述的装置实施例仅仅是示意性的，例如，所述模块的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个模块或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。

本领域技术人员也知道，除了以纯计算机可读程序代码方式实现控制器以外，完全可以通过将方法步骤进行逻辑编程来使得控制器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种控制器可以被认为是一种硬件部件，而对其内部包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。

本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构、类等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。

通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，移动终端，服务器，或者网络设备等)执行本申请各个实施例或者实施例的某些部分所述的方法。

本说明书中的各个实施例采用递进的方式描述，各个实施例之间相同或相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。本申请可用于众多通用或专用的计算机系统环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器系统、基于微处理器的系统、置顶盒、可编程的电子设备、网络PC、小型计算机、大型计算机、包括以上任何系统或设备的分布式计算环境等等。

虽然通过实施例描绘了本申请，本领域普通技术人员知道，本申请有许多变形和变化而不脱离本申请的精神，希望所附的实施方式包括这些变形和变化而不脱离本申请。