CN110632032A - 一种基于地表反射率库的沙尘暴监测方法 - Google Patents

Abstract

本发明公开了一种基于地表发射率库的沙尘暴监测方法，传统沙尘暴监测方法利用沙尘与其它地物光谱特性的差异，设定固定的阈值，实现沙尘信息识别。用最小值法合成月地表反射率数据库，通过6S模拟不同大气参数条件下晴空像元的最大表观反射率阈值，卫星遥感数据表光反射率大于模拟表观反射率，判别为沙尘像元，实现动态阈值法沙尘暴监测。通过验证结果分析比较，沙尘提取范围与真彩色影像沙尘区域一致，吻合较好，96.8%的地面观测站点沙尘天气与监测结果相同，OMI AI数据与沙尘监测结果基本一致。验证结果说明该方法有效降低了混合像元及大气环境等因素的影响，能够以较高精度实现沙尘暴监测。

Description

一种基于地表反射率库的沙尘暴监测方法

技术领域

本发明属于环境研究领域，更具体的说涉及一种基于地表反射率库的沙尘暴监测方法。

背景技术

沙尘暴(sand duststorm)是沙暴(sandstorm)和尘暴(duststorm)两者兼有的总称，是指强风把地面大量沙尘物质吹起并卷入空中，使空气特别混浊，水平能见度小于一千米的严重风沙天气现象。沙尘暴是一种区域性强灾害天气现象，对生态系统的破坏力极强，它能够加速土地荒漠化，对大气环境造成严重的污染，使空气质量显著下降，对人类健康、农业生产、交通运输、通讯等产生负面影响，还可以通过沙尘的沉降为海洋补充铁、氮等营养物质，改变海洋生态系统和生物化学循环，间接影响区域和全球的气候变化。沙尘暴研究手段主要有遥感监测和地面定点观测，由于沙尘暴主要发生在沙漠及其临近的干旱与半干旱地区，地面测站分布密度比较稀疏，因此用地面观测站网数据进行研究具有很大的局限性，而遥感数据的多源性、动态性、现时性和准确性，不仅弥补了地面观测数据时空分辨率的不足，同时还可与地面观测数据进行精度上的相互验证，互为补充，被广泛地应用于沙尘暴的监测、跟踪、预报和灾情评估等研究。

当前遥感监测技术主要使用卫星的可见光、近红外和热红外数据进行沙尘暴监测。Ackerman利用沙尘在热红外波段11μm和12μm处的亮温差为负的特性监测沙尘暴；胡秀清等根据沙尘粒子能导致地表发射到空间的红外信号产生衰减，用卫星观测到实时目标亮温减去地表背景亮温，构建红外差值沙尘指数IDDI (Infrared Difference Dust Index)监测沙尘暴；郭妮等使用MODIS数据的0.46μm 和2.1μm构建归一化沙尘判别指数NDSI>0，同时利用3.75μm和8.55μm构建另一沙尘判别指数DSI>33K识别沙尘暴；罗敬宁等给出了基于1.6μm通道e指数和红外分裂窗通道比值指数的沙尘强度指数(DDI)，DDI>1作为沙尘判识依据，值越大表示沙尘强度越强；Di等人利用可见光0.65μm和短波红外波段1.63μm 的表观反射率、红外波段3.9μm和10.8μm的亮温和大气气溶胶光学厚度构建增强沙尘指数EDI(Enhanced Dust Index)，EDI>0识别为沙尘像元。

以上遥感监测方法通过沙尘光谱特征分析，利用沙尘粒子在可见光、近红外和热红外等不同波段的光谱特征差异，选择一个或多个波段，采用不同的算法，给定合适阈值，从不同地物中提取沙尘信息，从而实现对沙尘暴的监测与预报，这种给定固定阈值的遥感监测方法称为固定阈值法。由于固定阈值法提取沙尘信息方便快捷，且具有一定的精度和稳定性，被广泛应用于沙尘暴遥感监测业务应用中。然而，沙尘暴遥感监测精度依赖于沙尘颗粒矿物的光谱特性、沙尘经过的下垫面特征、沙尘密度、云的类型和密度、传感器特性、沙尘层高度和沙尘监测技术等，传统固定阈值法由于受混合像元、复杂地表构成和大气环境背景等复杂因素的影响，在一些区域和条件下监测并不理想，存在薄沙尘区域监测困难，扬沙天气监测范围小，高反射率地表易产生误判，沙尘与云混淆等诸多问题。

发明内容

本发明地表反射率数据库支持的动态阈值监测沙尘暴的方法，以期降低混合像元及大气环境等因素的影响，提高沙尘暴监测精度和稳定性。在气溶胶、大气校正、云检测等领域进行了深入的研究，较好地解决了困扰遥感监测的混合像元和复杂大气环境影响的问题。

为了实现上述目的，本发明是通过以下技术方案实现的：所述的基于地表反射率库的沙尘暴监测方法包含以下步骤：

步骤1.构建地表反射率数据库；

步骤2.通过地表反射率得到表观反射率；

步骤3.确定地表反射率和模拟表观反射率的函数关系；

步骤4.求出像元对应的模拟表观反射率，再与遥感影像实际像元值比较，识别沙尘像元。

优选的，所述的步骤1.构建地表反射率数据库；具体方式如下：

从地球观测系统的数据和信息系统(EOSDIS)下载8天合成的MOD09A1地表反射率数据，每个月可有四幅MOD09A1数据产品影像，用最小值合成技术，将四幅影像中地表反射率最低的像元作为这个月的像元值，公式如下：

I(i,j)＝Min(I₁(i,j),I₂(i,j),I₃(i,j),I₄(i,j))

式中I为合成影像，I₁、I₂、I₃、I₄分别为一个月中四幅MOD09A1数据产品影像，i、j分别为一幅影像的行和列。

优选的，所述的步骤2.通过地表反射率得到表观反射率；具体方式如下：

由大气辐射传输方程可知，假设在气溶胶模式、大气和观测几何等条件已知的前提下，可以通过地表反射率得到表观反射率，所述的大气辐射传输方程为：

其中，ρ_a为大气固有反射率(程辐射)。由公式(3)可以看出表光反射率是大气光学厚度和地表反射率的函数。因此，在已知表观反射率ρ^*和地表反射率ρ的情况下，即可获得大气信息的ρ_a、s和T(τ_a,μ_s)T(τ_a,μ_v)。同样，如果在地表反射率ρ、ρ_a、s和T(τ_a,μ_s)T(τ_a,μ_v)已知的情况下，可以计算出相应的表观反射率ρ^*。

优选的，所述的步骤3.并确定地表反射率和模拟表观反射率的函数关系；具体方式如下：分别对MODIS数据的1、3、6和7四个波段进行模拟，得到每个像元地表反射率在晴空大气条件下的不同太阳和卫星天顶角下的表观反射率，模拟表观反射率取各种条件下的最大值，并确定地表反射率和模拟表观反射率的函数关系，

沙尘监测动态阈值法方程

。

优选的，所述的步骤4.求出像元对应的模拟表观反射率，再与遥感影像实际像元值比较，识别沙尘像元；具体方式如下：

沙尘判别满足以下条件:

R＝R_B∪R_R∪R_SWRI1∪R_SWIR2

式中

为不同波段的模拟表观反射率，ρ_i为已知地表反射率，α为太阳天顶角，β为卫星天顶角，

为表观反射率，R_i每个波段监测的结果，R为沙尘监测结果。如果影像数据像元的表观反射率

大于模拟表观反射率即R_i为正数，则像元判定为沙尘像元。

优选的，所述的步骤4，沙尘监测动态阈值法识别沙尘暴时，大气中的云被误判为沙尘，因此还需把可能的云像元分离。

优选的，采用归一化处理的方法，构建归一化沙尘判别指数NDDI(NormalizedDifference Dust Index)：

式中

和

分别为MODIS数据的第7波段和第3波段的表观反射率， NDDI为归一化的沙尘指数，如果NDDI>0，则该像元为沙尘，反之为云像元。

本发明有益效果：

采用MOD09A1地表反射率8天合成产品，用最小值法合成技术构建全球地表反射率数据库，通过6S模拟不同条件下的地表反射率和表观反射率的关系，实现了每个像元的动态阈值识别。沙尘暴遥感监测动态阈值法运用于四次沙尘暴事件，通过目视解译沙尘提取范围与真彩色影像沙尘区域一致，96.8％的MICAPS 地面观测站点沙尘天气与监测结果相同，OMI AI数据与沙尘监测结果基本一致，验证结果表明动态阈值法监测沙尘暴具有较高的监测精度，有效降低了混合像元及大气环境等因素的影响。

附图说明

图1，不同地物类型在MODIS(1-7波段)光谱响应曲线图；

图2，不同太阳和卫星天顶角下的地表反射率和表观反射率模拟结果图；

图3，图3为沙尘暴监测结果与MODIS原始影像、MICAPS数据和OMI AI产品数据的比较；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，下面结合附图与实例对本发明作进一步详细说明，但所举实例不作为对本发明的限定。

从图1中可以看出，沙尘(厚沙尘和薄沙尘)和沙漠(戈壁)的光谱曲线一致，反射率值随着波长的增加而增大，在蓝光波段0.469μm处最低，到了1.64μm 处达到最大值，在2.13μm处略有下降。云和雪(冰)的反射率值在可见光和近红外区较高，通常在0.50以上，但到了短波红外区呈直线下降的趋势，在2.13μm处达到最低值。植被光谱曲线变化不明显，其各波段的反射率值在0.16以下。

通过各地物类型的光谱曲线变化特征分析，沙尘的反射率值在可见光、近红外波段，一般介于云和地表之间，高于反射率在0.16以下的植被，低于反射率通常在0.5以上云和雪(冰)；在短波红外波段，沙尘光谱曲线变化和云、雪(冰) 正好相反，沙尘在此区域反射率值较高，通常在0.30上，而云、雪(冰)和植被正好相反；沙漠(戈壁)反射率值均明显低于厚沙尘，但与薄沙尘区界限不明显。

根据沙尘与不同地物在不同波段的光谱曲线的差异，选取反射率值差异大的波段，设定固定的阈值，将沙尘从云、冰雪、植被等地物中分离，实现沙尘暴的监测。

传统固定阈值法确定阈值的前提是从纯像元(pure pixel)中提取信息，但实际的遥感影像中很少是由单一均匀的地物组成，大都是由几种不同的地物类型混合而成。遥感影像中的单个像元的光谱信息不单单是某一种地物的光谱特征，而是反映了几种地物光谱的混合信息，观测的混合像元(mixed pixel)的反射率为像元内地物的光谱和它们面积的函数，公式如下：

式中P像元的反射率，N地物类型的总数，e_i为像元中每个地物的反射率， c_i为e_i在像元中所占的百分比，n不确定因子(随机噪声)。

从公式(4)中可以看出，像元值由地物类型的总数、地物的反射率及所占的百分比和不缺定因子决定。遥感影像中的像元可能是由沙尘和沙漠组成，可能是沙尘、云和戈壁的组合，也可能是冰雪与植被构成等，每个像元的反射率值代表这个像元不同性质的地物辐射能量的总和。混合像元是遥感影像普遍存在的问题，不同区域的不同地物构成的像元的反射率值各不相同，使用统一的固定阈值监测沙尘会产生大量的漏判和误判现象。混合像元问题使传统的固定阈值法沙尘暴遥感监测变得更加困难。

MODIS(Moderate Resolution Imaging Spectroradiometer)数据共36个波段，光谱范围为0.4～14.44μm，涵盖可见光、近红外和中远红外波段，空间分辨率分别为250m、500m、1km，带宽2330km，能够实现一到两天对地球观测覆盖一次。 MODIS数据具有高空间分辨率、高时间分辨率、高光谱分辨率等特点，可以为沙尘暴监测提供丰富的信息。根据图1不同地物类型的光谱曲线特征分析，选择沙尘反射率值与其它地物类型有较大差异的波段，在可见光、近红外和短波红外区选取1、3、4、6和7共五个波段进行沙尘暴监测研究，表1为选择的MODIS数据1-7波段的详细参数。

Table1 MODIS数据1-7波段参数

构建地表反射率数据库

构建地表反射率数据库使用MOD09A1产品。MOD09A1是Terra卫星获取的 MODIS地表反射率8天合成产品，空间分辨率500m，涵盖可见光、近红外和短波红外共7个波段，是MODIS1B产品1-7波段数据经过大气、气溶胶及卷云校正后获得。从美国国家航空航天局NASA地球观测系统的数据和信息系统(EOSDIS) 下载8天合成的MOD09A1地表反射率数据产品，假设地表反射率在一定时期内不发生变化，每个月可有四幅MOD09A1数据产品影像，用最小值合成技术，将四幅影像中地表反射率最低的像元作为这个月的像元值，公式如下：

I(i,j)＝Min(I₁(i,j),I₂(i,j),I₃(i,j),I₄(i,j))

式中I为合成影像，I₁、I₂、I₃、I₄分别为一个月中四幅MOD09A1数据产品影像，i、j分别为一幅影像的行和列。

动态阈值法

由大气辐射传输方程可知，假设在气溶胶模式、大气和观测几何等条件已知的前提下，可以通过地表反射率得到表观反射率。本文将大气模式、气溶胶光学厚度、观测几何、气溶胶模式等参数输入6S(the Second Simulation of the Satellite Signal in theSolar Spectrum)模型，分别对MODIS数据的1、3、6和7四个波段进行模拟，得到每个像元地表反射率在晴空大气条件下的不同太阳和卫星天顶角下的表观反射率，模拟表观反射率取各种条件下的最大值，模拟结果见图2，并确定地表反射率和模拟表观反射率的函数关系(见方程6-10)。

沙尘暴监测动态阈值法根据大气辐射传输理论，按照地表反射率和模拟表观反射率函数关系，将从构建的地表发射率数据库中获得的像元地表反射率带入函数方程，求出像元对应的模拟表观反射率，再与遥感影像实际像元值比较，如果卫星实测表观反射率大于模拟表观反射率，像元被认为是沙尘像元，实现了每个像元的动态阈值识别。

沙尘监测动态阈值法方程6-10：

沙尘判别满足以下条件:

R＝R_B∪R_R∪R_SWRI1∪R_SWIR2 (10)

式中

为不同波段的模拟表观反射率，ρ_i为已知地表反射率，α为太阳天顶角，β为卫星天顶角，

为表观反射率，R_i每个波段监测的结果，R为沙尘监测结果。如果影像数据像元的表观反射率

大于模拟表观反射率即R_i为正数，则像元判定为沙尘像元。

云的识别

沙尘暴发生时，往往伴随云一块产生。沙尘监测动态阈值法识别沙尘暴时，大气中的云被误判为沙尘，因此还需把可能的云像元分离。从图1可以看出，云的反射率在可见光和近红外波段较高，在蓝光波段第3波段，云的反射率表现为高值区，在短波红外区第7波段云的反射率最低，利用云在第7波段和第3波段上的反射率差值为负的特性识别云。因此，采用归一化处理的方法，构建归一化沙尘判别指数NDDI(Normalized Difference DustIndex)：

式中

和

分别为MODIS数据的第7波段和第3波段的表观反射率， NDDI为归一化的沙尘指数，如果NDDI>0，则该像元为沙尘，反之为云像元。

图3为沙尘暴监测结果与MODIS原始影像、MICAPS数据和OMI AI产品数据的比较，其中图3-a为四次沙尘暴的MODIS真彩色合成图像(RGB:1、4、 3)，图3-b为沙尘暴监测结果，图3-c为MICAPS地基站点实测沙尘数据，图 3-d为OMI AI产品数据，相关数据的成像和获取时间如表2所示：

Table2沙尘暴原始影像及验证数据获取时间统计表

通过图3-a真彩色合成影像和图3-b沙尘监测结果目视解译比较，无论是厚沙尘还是薄沙尘都能够较为准确地识别，特别是在低植被覆盖的高反射率地表薄沙尘提取效果较好，薄沙尘形状和纹理清晰可见，减小了云和高亮地表等因素的干扰问题。四次沙尘暴沙尘提取范围与真彩色影像沙尘区域一致，吻合较好。96.8％的MICAPS地面观测站点沙尘天气数据(图3-c)与监测结果相同，OMI AI数据(图3-d)与沙尘监测结果基本一致，验证结果表明动态阈值法监测沙尘暴具有较高的监测精度。

上述仅以实施例来进一步说明本发明的技术内容，以便于读者更容易理解。本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其他实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，任何依本发明所做的技术延伸或再创造，均受本发明的保护。

Claims (7)

1.一种基于地表反射率库的沙尘暴监测方法，其特征在于：所述的基于地表反射率库的沙尘暴监测方法包含以下步骤：

步骤1.构建地表反射率数据库；

步骤2.通过地表反射率得到表观反射率；

步骤3.并确定地表反射率和模拟表观反射率的函数关系；

步骤4.求出像元对应的模拟表观反射率，再与遥感影像实际像元值比较，识别沙尘像元。

2.根据权利要求1所述的一种基于地表反射率库的沙尘暴监测方法，其特征在于：所述的步骤1.构建地表反射率数据库；具体方式如下：

从地球观测系统的数据和信息系统(EOSDIS)下载8天合成的MOD09A1地表反射率数据，每个月可有四幅MOD09A1数据产品影像，用最小值合成技术，将四幅影像中地表反射率最低的像元作为这个月的像元值，公式如下：

I(i,j)＝Min(I₁(i,j),I₂(i,j),I₃(i,j),I₄(i,j))

式中I为合成影像，I₁、I₂、I₃、I₄分别为一个月中四幅MOD09A1数据产品影像，i、j分别为一幅影像的行和列。

3.根据权利要求1所述的一种基于地表反射率库的沙尘暴监测方法，其特征在于：所述的步骤2.通过地表反射率得到表观反射率；具体方式如下：

由大气辐射传输方程可知，假设在气溶胶模式、大气和观测几何等条件已知的前提下，可以通过地表反射率得到表观反射率，所述的大气辐射传输方程为：

其中，ρ_a为大气固有反射率(程辐射)。由公式可以看出表光反射率是大气光学厚度和地表反射率的函数。因此，在已知表观反射率ρ^*和地表反射率ρ的情况下，即可获得大气信息的ρ_a、s和T(τ_a,μ_s)T(τ_a,μ_v)。同样，如果在地表反射率ρ、ρ_a、s和T(τ_a,μ_s)T(τ_a,μ_v)已知的情况下，可以计算出相应的表观反射率ρ^*。

4.根据权利要求1所述的一种基于地表反射率库的沙尘暴监测方法，其特征在于:所述的步骤3.并确定地表反射率和模拟表观反射率的函数关系；具体方式如下：分别对MODIS数据的1、3、6和7四个波段进行模拟，得到每个像元地表反射率在晴空大气条件下的不同太阳和卫星天顶角下的表观反射率，模拟表观反射率取各种条件下的最大值，并确定地表反射率和模拟表观反射率的函数关系，

沙尘监测动态阈值法方程

。

5.根据权利要求1所述的一种基于地表反射率库的沙尘暴监测方法，其特征在于：所述的步骤4.求出像元对应的模拟表观反射率，再与遥感影像实际像元值比较，识别沙尘像元；具体方式如下：

沙尘判别满足以下条件:

R＝R_B∪R_R∪R_SWRI1∪R_SWIR2

式中

为不同波段的模拟表观反射率，ρ_i为已知地表反射率，α为太阳天顶角，β为卫星天顶角，

为表观反射率，R_i每个波段监测的结果，R为沙尘监测结果。如果影像数据像元的表观反射率

大于模拟表观反射率

即R_i为正数，则像元判定为沙尘像元。

6.根据权利要求5所述的一种基于地表反射率库的沙尘暴监测方法，其特征在于：所述的步骤4，沙尘监测动态阈值法识别沙尘暴时，大气中的云被误判为沙尘，因此还需把可能的云像元分离。

7.根据权利要求6所述的一种基于地表反射率库的沙尘暴监测方法，其特征在于：采用归一化处理的方法，构建归一化沙尘判别指数NDDI(Normalized Difference DustIndex)：

式中和

分别为MODIS数据的第7波段和第3波段的表观反射率，NDDI为归一化的沙尘指数，如果NDDI>0，则该像元为沙尘，反之为云像元。

Images