一种基于分段统计的大气气溶胶反演方法与流程

本发明属于遥感信息
技术领域：
，更为具体地讲，涉及一种基于分段统计的大气气溶胶反演法。
背景技术：
：大气气溶胶是指大气与悬浮在其中的固体和液体微粒共同组成的多相体系，是指悬浮在地球大气中的具有一定稳定性的，沉降速度小的，尺度范围在微米到几十微米之间的分子团、液态或固态粒子所组成的混合物。气溶胶粒子主要来源于工业活动、生物燃烧等人为源以及沙尘、近海海洋粒子等自然源。研究表明，大气气溶胶不仅影响地球表层系统的辐射收支平衡，而且通过对太阳短波辐射的散射调节地气系统的反射率从而导致地气系统的降温过程，通过对太阳辐射的吸收作用导致升温过程，通过对云的形成微物理过程的影响，可以改变云的微物理性质，对大气化学过程以及生物地球化学循环也起着重要的作用。同时，气溶胶粒子中包含一些对人体有害的粒子，对人类的健康状况造成严重的影响，尤其是对于人口密集、工业集中的城市地区。所以，大气气溶胶的反演对于全球气候变化的研究和大气污染的监测和治理工作具有重大意义。目前国内外学者在气溶胶遥感反演方面取得了一定的成果，并且在新型传感器的设计中也越来越考虑到对气溶胶的监测。但是目前大多数的反演算法都是针对大尺度平均性质的气溶胶，而且具有很多局限性，例如，只适用于暗地表。目前还没有一种有效的适用于城市等亮地表地区的气溶胶反演算法。技术实现要素：本发明的目的在于克服现有技术的不足，提供一种基于分段统计的大气气溶胶反演法，采用分段统计的方式划分为两个部分，采用不同方式进行气溶胶厚度值反演，可以提高亮地表地区反演结果的精度和分辨率。为实现上述发明目的，一种基于分段统计的大气气溶胶反演方法，该方法包括以下步骤：S1：通过卫星获取所分析区域的多波段遥感影像，并进行辐射定标得到含有表观反射率的多波段遥感影像，所述多波段遥感影像中包括中红外1.6微米波段、中红外2.1微米波段，从所述多波段遥感影像中选取一个波段作为反演波段；S2：确定大气辐射传输模型；根据反演波段构建对应的气溶胶光学厚度查找表，查找表中包含每个气溶胶厚度值AOT所对应的用于计算表观反射率的参数；S3：对于多波段遥感影像中的所有像素，采用分段统计方法进行分段，分段方法包括以下步骤：S3.1：将中红外2.1微米波段的表观反射率范围按照预设间隔λ1划分为N个区间S3.2：根据多波段遥感影像中每个像元在中红外2.1微米波段的表观反射率，将各个像元划分到对应的区间，得到每个区间对应的像元集合xn；S3.3：对于步骤S3.2得到的每个像元集合xn，如果像元集合xn中的像元数量|xn|≥T1，T1表示预设阈值，则保留该像元集合，否则删除；记筛选得到的像元集合为xn′，n′＝1,2,…,N′，N′表示筛选得到的像元集合数量；S3.4：将中红外1.6微米波段的表观反射率范围按照预设间隔λ2划分为M个区间γm，m＝1,2,…,M；S3.5：对每个像元集合xn′，根据其每个像元在中红外1.6微米波段的表观反射率，将各个像元划分到对应的区间，得到该像元集合下每个区间γm对应的像元集合yn′,m；S3.6：如果像元集合yn′,m中的像元数量|yn′,m|≥T2，T2表示预设阈值，则保留该像元集合，否则删除；记筛选得到的像元集合为yk，k＝1,2,…,K，K表示筛选得到的像元集合数量；S3.7：对于步骤S3.6中筛选得到的K个像元集合yk，根据预设阈值T3进行分类，如果像元集合yk中的像元数量|yk|≥T3，T3表示预设阈值，则将该像元集合yk划入集合Ys中，否则划入集合Yt中；S4：记集合Ys中第p个像元集合为ys,p，p＝1,2,…,|Ys|，|Ys|表示集合Ys中的像元集合数量，对每个像元集合ys,p依次进行反演，具体方法包括以下步骤：S4.1：将反演波段的表观反射率范围按照预设间隔λ3划分为D个区间ωd，d＝1,2,…,D；对于像元集合ys,p，根据其每个像元在反演波段的表观反射率，将各个像元划分到对应的区间，得到该像元集合下每个区间ωd对应的像元集合zs,p,d；S4.2：从D个像元集合zs,p,d中搜索得到像元数量|zs,p,d|≥T4的像元集合，T4表示预设阈值，从搜索得到的像元集合中选择区间序号d最小的像元集合d*作为像元集合ys,p中的清洁区段；S4.3：令清洁区段在反演波段的表观反射率其对应的气溶胶厚度值AOT0为预设的清洁区段气溶胶厚度值；在气溶胶光学厚度查找表中查找得到对应的参数，根据表观反射率计算清洁区段对应的地表反射率；S4.4：将清洁区段对应的地表反射率作为整个像元集合ys,p的地表反射率，在步骤2确定的大气辐射传输模型框架下反演得到像元集合ys,p中每个像元的气溶胶厚度值；S5：对集合Yt进行反演，其具体方法为：S5.1：遍历集合Ys中的每个像元，搜索预设半径内的无值像元，令无值像元的气溶胶厚度值设置为该像元的气溶胶厚度值；记所有被赋值的像元集合为Ps；S5.2：对于集合Ps中每个像元，在气溶胶光学厚度查找表中查找得到所填充气溶胶厚度值对应的参数，然后根据表观反射率计算得到该像元的地表反射率；S5.3：按照预设的边长将多波段遥感影像划分为网格；S5.4：记集合Yt中第q个像元集合为yt,q，q＝1,2,…,|Yt|，|Yt|表示集合Yt中的像元集合数量；对于每个像元集合yt,q，在步骤S5.3划分出的每个网格中求取与Ps的交集，令像元集合yt,q在当前网格中所有像元的地表反射率等于该交集中所有像元的地表反射率均值；S5.5：在步骤2确定的大气辐射传输模型框架下对步骤S5.4被赋予地表反射率的每个像元进行气溶胶厚度值反演，得到该像元的气溶胶厚度值。进一步的，所述步骤2中采用6S辐射传输模型；表观反射率ρTOA的计算公式为：其中，ρ0表示大气层辐射反射率，ρs为地表反射率，S表示大气底层向下的半球反射率，Ts表示入射光线从大气顶层到达地表面的总的透过率，Tv表示向上进入卫星传感器视场方向的总透过率；构建的气溶胶光学厚度查找表包括：AOT、ρ0、S、Ts·Tv之间的关系。进一步的，所述步骤S3.1中间隔λ1的取值范围为0.002≤λ1≤0.01。进一步的，所述步骤S3.4中间隔λ2的取值范围为0.002≤λ2≤0.01。进一步的，所述步骤S4.2中阈值T4的取值范围为0.02|Ys|≤T4≤0.1|Ys|。进一步的，所述步骤S4.3中间隔λ3的取值范围为0.002≤λ3≤0.01。进一步的，所述步骤S4.4和步骤S5.5中气溶胶厚度值的反演方法为：从查找表中查找得到各个气溶胶厚度值AOT所对应的参数，根据得到的地表反射率计算每个气溶胶厚度值AOT对应的表观反射率ρi，其中i＝1,2,…,L，L表示查找表中气溶胶厚度值的个数；对于每个像元，从L个表观反射率ρi中搜索与其自身表观反射率最接近的表观反射率，其对应的气溶胶厚度值AOT即为该像元的气溶胶厚度值。本发明基于分段统计的大气气溶胶反演法，首先得到多波段遥感影像和反演波段对应的气溶胶光学厚度查找表，然后按照中红外2.1微米波段的表观反射率区间对像元进行划分和筛选，然后对得到的像元集合按照中红外1.6微米波段的表观反射率区间进行进一步划分和筛选，将最终得到的像元集合按照像元数量划分为两类，像元较多的划为一类，其他为另一类；将像元较多的一类作为基准部分进行反演，采用方法是先从像元集合中搜索出清洁区段，以清洁区段的地表反射率作为整个像元集合的地表反射率，反演得到气溶胶厚度值，然后以这些像元作为基准，对另一类进行反演。本发明在进行地表反射率的计算时采用了一种新方法，该方法在计算地表反射率时，是依据清洁像元和基准区段来确定地表反射率的，并不依赖于暗像元，因此只要满足完成分段统计即可进行反演，所以较传统的暗像元算法，本发明对亮地表地区也具有良好的反演效果。具有更广泛的适用性。附图说明图1是本发明基于分段统计的大气气溶胶反演法的具体实施方式流程图；图2是本实施例所采用的影像图；图3是本发明中像素分段的流程图；图4是本实施例中波段7的像元划分示例图；图5是地表反射率波谱曲线示例图；图6是对集合Ys中每个像元集合ys,p进行反演的流程图；图7是集合Yt进行反演的流程图；图8是图2所示影像的网格化结果；图9是采用本发明对图2所示多波段遥感图像的反演结果。具体实施方式下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。实施例图1是本发明基于分段统计的大气气溶胶反演法的具体实施方式流程图。如图1所示，本发明基于分段统计的大气气溶胶反演法包括以下步骤：S101：获取多波段遥感影像：通过卫星获取所分析区域的多波段遥感影像，并进行辐射定标得到含有表观反射率的多波段遥感影像，多波段遥感影像中包括中红外1.6微米波段、中红外2.1微米波段，从多波段遥感影像中选取一个波段作为反演波段。目前获取多波段遥感影像的卫星为LANDSAT卫星，在LANDSAT8OLI影像中，包含了9个波段，分别是Band1海蓝波段、Band2蓝绿波段、Band3绿波段、Band4红波段、Band5近红外、Band6中红外1.6微米波段、Band7中红外2.1微米波段、Band8全色、Band9卷云、Band10热红外、Band11热红外，本发明中需要使用的是Band6中红外和Band7中红外。本实施例中所选择的反演波段为Band1蓝色波段，这是因为LANDSAT8OLI影像中的Band1对气溶胶更敏感，可以提高反演的精度。本实施例中所采用的是北京地区2014年第103天的影像。图2是本实施例所采用的影像图。如图2所示，该影像图中包含了亮地表区域(图中颜色较浅区域)。S102：获取气溶胶光学厚度查找表：将步骤S101所获得的多波段遥感影像输入大气辐射传输模型，根据反演波段得到对应的气溶胶光学厚度查找表，查找表中包含每个气溶胶厚度值AOT所对应的用于计算表观反射率的参数。本实施例中，假设下垫面是朗伯体表面，其表观反射率ρTOA的计算公式为：其中，ρ0表示大气层辐射反射率，ρs为地表反射率，S表示大气底层向下的半球反射率，Ts表示入射光线从大气顶层到达地表面的总的透过率，Tv表示向上进入卫星传感器视场方向的总透过率。因此可知，本实施例中查找表中需要包含的参数包括大气程辐射反射率、大气半球反射率、入射光线从大气顶层到达地表面的总的透过率和向上进入卫星传感器视场方向的总透过率的乘积。目前行业内已经开发出多种大气辐射传输模型，本实施例中采用6S辐射传输模型。除了多波段遥感影像外，其他需要输入6S辐射传输模型的参数包括：多波段遥感影像所对应的日期、天顶角、方位角，气溶胶类型和大气模式，气溶胶厚度值的步长等等。本实施例中选取的大气模式是中纬度夏季，气溶胶类型是大陆型，所选择的反演波段为Band1海蓝波段。表1是本实施例中获取的查找表。AOTρ0STs·Tv0.010.090940.758230.1740.020.091780.753230.17572…………2.990.243860.074490.306720S103：像元分段统计：对于多波段遥感影像中的所有像素，采用分段统计方法进行分段。因为中红外1.6微米波段、Band7中红外2.1微米波段具有比较长的波长，所以它们的表观反射率可以近似等于地表反射率，因此本发明采用中红外1.6微米波段、中红外2.1微米波段的表观反射率来对像素进行分段。由于本实施例中采用的LANDSAT8OLI影像，因此中红外1.6微米波段为波段6，中红外2.1微米波段为波段7。图3是本发明中像素分段的流程图。如图3所示，本发明中像素分段的具体步骤为：S301：划分波段7表观反射率区间：将波段7的表观反射率范围按照预设间隔λ1划分为N个区间划分间隔λ1可以根据实际需要来确定，一般来说，间隔λ1的取值范围为0.002≤λ1≤0.01。S302：像元划分：根据多波段遥感影像中每个像元在波段7的表观反射率，将各个像元划分到对应的区间，得到每个区间对应的像元集合xn。记某个像元在波段7的表观反射率为ρ，那么其对应的区间序号即为表示向上取整。图4是本实施例中波段7的部分像元以0.005的间隔划分示例图。图4中，图(a)是间隔划分的整体图，图(b)是间隔划分的局部图。如图4所示，通过对波段7的表观反射率范围划分区间，从而可以对多波段遥感影像中的像元进行区段划分。S303：像元集合筛选：依次对步骤S302得到的N个像元集合xn进行判定，如果像元集合xn中的像元数量|xn|≥T1，T1表示预设阈值，则保留该像元集合，否则删除。记筛选得到的像元集合为xn′，n′＝1,2,…,N′，N′表示筛选得到的像元集合数量。阈值T1的数值也可以根据实际需要来设置，显然，当多波段遥感影像的分辨率越高，间隔λ1越大，每个区间内所包含的像元数量就越多，相应地阈值T1也可以设置得略高一些，反之阈值T1可以设置得略低一些。本实施例中T1＝100000。S304：划分波段6表观反射率区间：将波段6的表观反射率范围按照预设间隔λ2划分为M个区间γm，m＝1,2,…,M。同样地，划分间隔λ2可以根据实际需要来确定，一般来说，预设间隔λ1的取值范围为：0.002≤λ1≤0.01。S305：像元集合细分：对步骤S303筛选得到的每个像元集合xn′分别进行细分，细分方法为：对每个像元集合xn′，根据其每个像元在波段6的表观反射率，将各个像元划分到对应的区间，得到该像元集合下每个区间γm对应的像元集合yn′,m。显然，对像元集合进行细分后，共计可以得到N′×M个像元集合。S306：像元集合再筛选：对于步骤S305得到的N′×M个像元集合进行筛选。如果像元集合yn′,m中的像元数量|yn′,m|≥T2，T2表示预设阈值，则保留该像元集合，否则删除。记筛选得到的像元集合为yk，k＝1,2,…,K，K表示筛选得到的像元集合数量。与阈值T1一样地，阈值T2也是根据实际情况来设置的，由于像元集合yn′,m是细分后的像元集合，显然像元集合yn′,m中的像元数量要更少，因此T2也应当小于T1。本实施例中T2＝10000S307：像元集合分类：对于步骤S306中筛选得到的K个像元集合yk，根据预设阈值T3进行分类，如果像元集合yk中的像元数量|yk|≥T3，T3表示预设阈值，则将该像元集合yk划入集合Ys中，否则划入集合Yt中。相应地，阈值T3也是根据实际需要进行设置的。本实施例中，T3的取值范围为100000＜T3＜200000，本实例选择T3＝200000，显然，集合Ys和集合Yt包含了许多个小的像元集合。根据以上步骤就将多波段遥感影像中的像元按照波段6和波段7的表观反射率进行了分段。图5是地表反射率波谱曲线示例图。图5中展示了ENVI标准波谱库中一些地物的地表反射率波谱曲线，通过研究大量地物的波谱曲线可以发现同种地物或不同地物在不同的条件下具有不同的地表反射率曲线，它们可能在中红外1.6微米波段和中红外2.1微米波段同时有重合，但是这种情况很少，尤其是对于常见的地物。而且可以通过后续的两种对分段结果不同的处理方法来减小误差，同时也可以通过规定每个区段的最小像元数量来减小误差。因此，尽管本发明中所采用的分段方法存在误差，但是仍旧可以利用中红外1.6微米波段和中红外2.1微米波段将多波段遥感影像中的像元分割成许多区段。如果这些区段没有误差的存在，那么属于同一区段的像元应该具有近似相同的地表反射率曲线，它们在反演波段的地表反射率应该近似相等。S104：对集合Ys进行反演：本发明将集合Ys作为基准部分，采用一种基于清洁像元反演污染像元的方法进行处理。在气溶胶的反演中，地表反射率的获得至关重要。实验证明基准集合Ys中的每个像元区段存在很少的误差，这主要因为基准集合Ys中的每个像元集合都是同一种常见地物，例如，植被和裸土等，同时还可以通过设置一个大的阈值可以进一步排除一些误差。而且Ys中每个像元集合含有清洁像元的概率也比较大。记Ys中第p个像元集合为ys,p，p＝1,2,…,|Ys|，|Ys|表示集合Ys中的像元集合数量。那么根据分段方法的特点，可以得到其中，ρs,p表示每个像元集合ys,p在波段1的近似地表反射率，表示像元集合ys,p中清洁像元在反演波段(本实施例为波段1)的地表反射率。因此，只要求得每个像元集合中清洁像元的地表反射率即可得到像元集合中每个像元的地表反射率。因此，首先要得到每个像元集合中清洁像元的地表反射率。图6是对集合Ys中每个像元集合ys,p进行反演的流程图。如图6所示，像元集合ys,p进行反演的具体方法为：S601：像元集合细分：将波段1的表观反射率范围按照预设间隔λ3划分为D个区间ωd，d＝1,2,…,D。同样地，划分间隔λ3可以根据实际需要来确定，间隔λ3的取值范围为0.005≤λ3≤0.02。对于像元集合ys,p，根据其每个像元在波段1的表观反射率，将各个像元划分到对应的区间，得到该像元集合下每个区间ωd对应的像元集合zs,p,d。S602：搜索清洁区段：从D个像元集合zs,p,d中搜索得到像元数量|zs,p,d|≥T4的像元集合，T4表示预设阈值，从搜索得到的像元集合中选择区间序号d最小的像元集合d*作为像元集合ys,p中的清洁区段。阈值T4也是根据实际情况来设置的，一般来说，阈值T4的取值范围为0.02|Ys|≤T4≤0.1|Ys|，可以根据多波段遥感影像的污染范围进行选择，污染范围越大该值可以设置的越小，本实例选择T4＝0.03|Ys|。S603：计算地表反射率：令清洁区段在波段1的表观反射率其对应的气溶胶厚度值AOT0为预设的清洁区段气溶胶厚度值。显然，AOT0应当是一个极小值，是预先根据需要来设置的。在气溶胶光学厚度查找表中查找得到对应的参数，计算清洁区段对应的地表反射率。显然，本实施例中地表反射率是通过表观反射率计算公式反向求解得到的。S604：对集合ys,p进行反演：根据之前的分析可知，清洁区段对应的地表反射率就是整个像元集合ys,p的地表反射率，因此就能以此来反演得到像元集合ys,p中每个像元的气溶胶厚度值。反演方法可以根据需要来进行选择，本实施例中采用的反演方法为：从查找表中查找得到各个气溶胶厚度值AOT所对应的参数，计算每个气溶胶厚度值AOT对应的表观反射率ρi，其中i＝1,2,…,L，L表示查找表中气溶胶厚度值的个数。对于每个像元，从L个表观反射率ρi中搜索与其自身表观反射率最接近的表观反射率，其对应的气溶胶厚度值AOT即为该像元的气溶胶厚度值。S105：对集合Yt进行反演：由于步骤S104中对集合Ys中的像元进行了反演，得到了各个像元的气溶胶厚度值，对于集合Yt中的像元，以集合Ys中像元为基准，进行反演。图7是集合Yt进行反演的流程图。如图7所示，集合Yt进行反演的具体步骤为：S701：基准像元扩展：遍历集合Ys中的每个像元，搜索预设半径内的无值像元，令无值像元的气溶胶厚度值设置为该像元的气溶胶厚度值；记所有被填充赋值的像元集合为Ps。搜索半径r一般根据需要进行设置，搜索半径越小，所填充的气溶胶厚度值越准确，但是后续剩余的无值像元数量会增多，一般其取值范围为1≤r≤10。S702：计算集合Ps中像元的地表反射率：对于集合Ps中每个像元，在气溶胶光学厚度查找表中查找得到所填充气溶胶厚度值对应的参数，然后根据表观反射率计算得到该像元的地表反射率。S703：多波段遥感图像网格化：考虑到地物在地理空间上具有一定的相关性，按照预设的边长将多波段遥感影像划分为网格。预设边长可以根据多波段遥感图像的分辨率来设置。图8是图2所示影像的网格化结果。S704：获取地表反射率：记集合Yt中第q个像元集合为yt,q，q＝1,2,…,|Yt|，|Yt|表示集合Yt中的像元集合数量。对于每个像元集合yt,q，在步骤S703划分出的每个网格中求取与Ps的交集，令像元集合yt,q在当前网格中所有像元的地表反射率等于该交集中所有像元的地表反射率均值。由此，就可以得到集合Yt中像元的地表反射率。S705：气溶胶厚度值反演：对步骤S705被赋予地表反射率的每个像元进行气溶胶厚度值反演，得到该像元的气溶胶厚度值。S106：无值像元插值。由于本发明在步骤S103中对像元集合进行了筛选，因此集合Ys和Yt中并未包含多波段遥感图像的所有像元，在步骤S105中也难以实现对所有像元的反演。对于剩余的无值像元，采用插值方法补充其气溶胶厚度值。插值方法是一种图像处理领域的常用方法，其具体过程在此再赘述。图9是采用本发明对图2所示多波段遥感图像的反演结果。如图9所示的结果可知，采用本发明对亮地表地区也具有良好的反演效果，较传统的暗像元算法具有更广泛的适用性。为了验证其准确性，在http://aeronet.gsfc.nasa.gov/这个网站上下载了AERONET(AerosolRoboticNetwork)地面观测网络的数据。将本发明的反演结果和地面观测数据以及暗像元算法的反演结果进行比较，所选择的两个站点均位于城市中，即属于亮地表区域。表2是图2所示影像图本发明反演结果、暗像元算法反演结果与地面观测数据对比。站点BeijingBeijing-RADI观测数据1.3511.389暗像元法反演结果1.661.73本发明反演结果1.381.42暗像元法误差0.3090.341本发明误差0.0290.031表2如表2所示，本发明的反演结果与实际站点的观测数据更为接近，其误差更小。可见，对于亮地表地区，本发明具有良好的反演效果。因此与暗像元法相比，本发明的适用范围更广。尽管上面对本发明说明性的具体实施方式进行了描述，以便于本
技术领域：
的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本
技术领域：
的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。当前第1页1 2 3