CN103913750B - 一种基于卫星云参数产品确定灰霾分布的方法 - Google Patents

Abstract

本发明公布了种基于卫星云参数产品确定灰霾分布的方法，该方法包括以下步骤：读取卫星数据中云量产品并和对应的判据进行比较，确定云或霾出现可能性较大的范围；利用卫星数据中云光学厚度信息和对应的判据进行比较，进一步提取霾分布的区域；利用卫星数据中云相位信息和对应的判据进行比较,进一步提取霾分布的区域；利用云顶压力数据和对应的判据进行比较，进一步提取霾分布的区域。本发明还设计确定所述判据的方法。本发明尤其以POLDER卫星数据P3L2TRGB为例说明实施的过程。该方法在不需要其他辅助数据前提下，能够快速稳定地实现POLDER卫星数据霾的检测，计算量小且适用范围广。

Description

一种基于卫星云参数产品确定灰霾分布的方法

技术领域

本发明涉及卫星遥感数据云检测技术领域，涉及一种利用云参数产品进行霾检测方法。

背景技术

改革开放以来，我国的城市空气污染日趋严重，工业烟尘、汽车尾气、生物燃烧和沙尘和自然粉尘混合使灰霾天气现象迅速增加。尤其是以京津冀、长三角和珠三角为代表的重点城市群。灰霾天气主要是由于以细粒子为主的气溶胶粒子造成的。这些气溶胶粒子不但对人体有害，而且能够影响地表辐射平衡，并通过充当云的凝结核改变云的含水量和生命周期等特性，从而对气候和降水造成影响。目前在对全球气候模拟的过程中，云和灰霾对气候变化的影响仍然有很多不确定性。一方面是云和霾本身的分布特征很难精确的确定，另一方面云和气溶胶影响着全球水汽和能量循环。因此，云和灰霾分布范围的区分是研究它们对气候变化贡献的前提和基础。

目前，遥感领域利用卫星影像识别云和霾分布的方法有：判据法、聚类分析法和人工神经网络法等。其中判据法应用最为广泛，该方法基于卫星影像中象元对应的亮温、亮温差、反射率与预设的判据进行对比，来判断该象元是否被云或者霾覆盖。这种方法实施方便，但是对云和霾的错分情况严重，而且单一的判据的适应性较差；为了提高判据算法的精度，有些学者引入了纹理特征约束，这种约束是基于卫星影像上灰霾的一般呈现平滑的纹理特征，而云层的表面细节变化较大。本质上是在设置判据的时候考虑到纹理特征。显然，当大范围表面平滑的云存在时，纹理特征将不能有效区分二者；另外，有些学者通过建立地面的地表反射率库，将库中地面象元对应的反射率值作为设置判据的标准，来判断象元是否收到云或者灰霾的覆盖。这一方法需要长时间序列的地表观测数据来建库，操作起来不方便。聚类分析法通过将影像分割成为一系列像元阵列，采用聚类的方法将象元归属到不同的地物类型中，然后确定各组象元的类别：云、霾或其他地物类型。当像元阵内各种类型象元观测值存在较小差别时，无法对象元进行有效的分组，这种方法对云和霾分布的识别误差较大。人工神经网络法则先对云和霾区域样本区域进行学习，获取其在不同波段反射率、亮温等观测结果的先验知识，借助这些先验知识来识别云和霾的分布。这种方法得到的识别结果很大程度上依赖于训练资料的选取。

卫星遥感一直重视对云和灰霾等重要大气现象的观测和研究，很多传感器如美国的中分辨率成像光谱仪(MODIS)、法国的地表反射率偏振和方向传感器(POLDER)等都有专门的产品(云参数)来描述云的光学和微物理特性，比如云量、云顶压力、云相态、云光学厚度等云参数。其中云量代表每个象元中被云覆盖的百分比，100％表示该象元完全被云所覆盖，相反，0则代表该象元无云覆盖；云顶压力为卫星观测到的云顶的大气压，单位为hPa；云相态反映卫星观测到的象元对应的类别，包括水云、冰云、冰水混合云、晴空和无法确定相态五类；云光学厚度则反演云的消光能力的强弱，光学厚度越大，消光能力越强。卫星遥感首先识别影像中云的分布，然后对确定有云的象元进行上述云产品的生产。灰霾覆盖的象元通常被错误地归属为云区域，然后按照流程反演出相应的云量、云顶压力、云相态结果。然而，云和灰霾的光学和微物理特性都存在明显的差异(云的粒子更大、在对流层不同高度都有分布、光学厚度可以达到几十甚至上百；灰霾主要有细颗粒物构成、主要分布在地表、消光能力远不如一些厚层云)，这些差别也将表现在反演得到的云产品中。

现有的霾识别算法一般基于卫星观测的反射率和亮温来构造判据，这种做法对于目前重污染天气的灰霾范围识别精度较差。很多卫星提供云量、云顶压力、云相态、云光学厚度等观测数据，这些数据包含灰霾的分布数据，但是并没有相关的技术采用这些数据来识别灰霾的范围。

发明内容

本发明提供一种基于卫星云参数产品的灰霾识别方法,该方法能够精确的识别重污染天气下灰霾的分布范围。

本发明基于卫星云参数产品的确定灰霾分布的方法包括如下步骤：

确定云参数识别灰霾的判据。针对待研究区域，下载历史云参数产品中卫星云量、云光学厚度、云相态频率和云顶压力等云参数产品组建样本集，对样本集中影像的云和灰霾区域进行区分，分别对云、灰霾、晴空区域的所述各个云参数结果进行采样统计，确定对应于所述各个云参数的识别判据，分别为T1、T2、T3、T4。

确定云和灰霾共存的范围。利用当期云参数产品中云量数据和对应的判据T1进行比较，确定云和灰霾都可能存在的区域；

利用云光学厚度数据识别灰霾的分布。在其它步骤选出的区域，利用当期云参数产品中云光学厚度数据和对应的判据T2进行比较，进一步提取霾分布的区域；

利用云相态数据识别灰霾分布。在其它步骤选出的区域，利用当期云参数产品中云相态数据和对应的判据T3进行比较，进一步提取霾分布的区域；

利用云顶压力数据识别灰霾分布。在其它步骤选出的区域，利用当期云参数产品中云顶压力数据和对应的判据T4进行比较，进一步提取霾分布的区域。

其中，确定判据T1、T2、T3、T4的方法包括如下：

统计灰霾或云覆盖区域的云量数据，构建“灰霾或云”频率随云量变化的拟合曲线；统计晴空区域的云量数据，构建晴空出现频率随云量变化的拟合曲线,两曲线交点的云量值作为云量判据T1。

统计灰霾区域的云光学厚度数据，构建灰霾频率随云光学厚度变化的拟合曲线；统计云区域的光学厚度数据，构建云出现频率随云光学厚度变化的拟合曲线，两曲线交点的云光学厚度值为作为云光学厚度判据T2。

统计灰霾区域各种云相态(水云、冰云、冰水混合云、无法确定类型、晴空)出现的频率，超过一定比例R的云相态类型作为识别判据T3。

统计灰霾区域的云顶压力数据，构建灰霾出现频率随云顶压力变化的拟合曲线；统计云区域的云顶压力数据，构建云出现频率随云顶压力变化的拟合曲线；两曲线交点对应的云顶压力值作为云顶压力判据T4。

本发明针对现有霾识别算法在重污染天气下识别精度较差的问题，第一次使用卫星观测的云量、云光学厚度、云相态和云顶压力数据来识别灰霾的范围。该方法在不需要其他辅助数据前提下，能够快速稳定的实现霾的检测，计算量小且适用范围广。为进一步研究霾的成因及分布，气溶胶的传输，云和气溶胶的相互作用等提供必要的前提。

附图说明

图1是基于卫星云参数产品确定灰霾分布方法流程图

图2是云量判据T1的确定方法示意图

图3是云光学厚度判据T2的确定方法示意图

图4是云相态判据T3的确定方法示意图

图5是云顶压力判据T4的确定方法示意图

具体实施方式

下面结合附图，说明本方法的实施方式。

如图1所示，是本发明基于卫星云参数产品确定灰霾分布方法的一个实施例。步骤如下：

S1、确定云参数识别灰霾的判据。根据卫星云量、云光学厚度、云相态频率和云顶压力数据确定识别判据，分别为T1、T2、T3、T4。

S2、确定云和灰霾共存的范围。利用云参数产品中云量数据，利用云量数据和对应的判据(T1)进行比较，确定云和灰霾都可能存在的区域。

S3、利用云光学厚度数据识别S2结果中灰霾的分布。针对S2选出的区域，利用云光学厚度数据和对应的判据T2进行比较，进一步提取霾分布的区域。

S4、利用云相态数据识别S3中的灰霾分布。针对S3获取的区域，利用云相态数据和对应的判据T3进行比较，进一步提取霾分布的区域。

S5、利用云顶压力数据识别S4中的灰霾分布。在S4提取的霾区域的基础上，利用云顶压力数据和对应的判据T4进行比较，进一步提取霾分布的区域。

以上步骤中S3、S4和S5的顺序可以互换，不影响识别结果。

其中，步骤S1进一步包括：

S1.1、针对待研究区域，筛选该区域范围内灰霾事件及这些事件对应的日期，下载这些日期的卫星云量、云光学厚度、云相态和云顶压力的历史数据，组建样本集。

S1.2、用人工目视解译的方法对样本集中影像的云和灰霾区域进行区分，分别对云、灰霾、晴空区域的上述云参数结果进行采样统计。

S1.3、统计灰霾或云覆盖区域的云量数据，构建“灰霾或云”频率随云量变化的拟合曲线；统计晴空区域的云量数据，构建晴空出现频率随云量变化的拟合曲线,两曲线交点的云量值作为云量判据T1。云量判据T1的确定方法如图2所示，虚线是根据晴空象元的云量频率分布拟合得到，实线是根据云或霾象元的云量频率分布拟合得到。一般情况下，晴空象元的云量较小，因此虚线在左边。二者的交点对应的云量值即为T1。

S1.4、统计灰霾区域的云光学厚度数据，构建灰霾频率随云光学厚度变化的拟合曲线；统计云区域的光学厚度数据，构建云出现频率随云光学厚度变化的拟合曲线，两曲线交点的云光学厚度值为作为云光学厚度判据T2。云光学厚度判据T2的确定方法如图3所示，虚线是根据霾象元的云光学厚度频率分布拟合得到，实线是根据云象元的云量频率分布拟合得到，一般情况下，霾象元的云光学厚度较小，因此虚线在左边。二者的交点对应的云光学厚度值即为T2。

S1.5、统计灰霾区域各种云相态出现的频率(水云、冰云、冰水混合云、无法确定类型、晴空)，超过一定比例R的云相态类型作为识别判据T3。云相态判据的确定方法如图4所示例如灰霾区域的超过所述比例R的云相态类型有水云、冰水混合云，那么T3即为水云、冰水混合云。表示为T3＝{水云，冰水混合云}。

S1.6、统计灰霾区域的云顶压力数据，构建灰霾出现频率随云顶压力变化的拟合曲线；统计云区域的云顶压力数据，构建云出现频率随云顶压力变化的拟合曲线；两曲线交点对应的云顶压力值作为云顶压力判据T4。云顶压力判据T4的确定方法如图5所示，虚线是根据云象元的云顶压力频率分布拟合得到，实线是根据霾象元的云顶压力频率分布拟合得到，一般情况下，云的位置较高，因此该区域的云顶压力也较小。二者的交点对应的云顶压力值即为T4。

其中，步骤S2进一步包括：

S2.1、针对已经确定霾识别判据的区域，下载当前卫星云参数产品，包括云量、云顶压力、云相态和云光学厚度。

S2.2、利用当前云参数产品中的图像高度和宽度，象元经纬度数据以及云量数据。重新生成与输入数据等高和等宽的影像，像素值均为1，作为霾识别的掩膜。

S2.3、利用逐个象元的云量结果和对应的判据T1进行比较，大于或等于T1的象元为云或灰霾，该象元在掩膜中对应的值保持不变，而没有通过判据判断的象元的掩膜值乘以0，作为步骤S2的提取结果。

其中，步骤S3进一步包括：

S3.1、利用遥感影像中的云光学厚度数据。

S3.2、对于前一步骤S2识别出的象元，将云光学厚度结果和对应的判据T2进行比较，大于或等于T2的象元为灰霾，其对应的掩膜中的象元值保持不变，而没有通过判据判断的象元的掩膜值乘以0，作为S3的提取结果。

其中，步骤S4进一步包括：

S4.1、利用卫星影像中的云相态数据，即每个象元的在以下几类中的归属：水云、冰云、冰水混合云、晴空或者无法确定相态。

S4.2、对于前一步骤S3识别出的象元，逐个判断其云相态类别，云参数产品中云相态类别与T3一致的象元为灰霾，其对应的掩膜值保持不变，而没有通过判据判断的象元的掩膜值乘以0，作为S4的提取结果。

其中，步骤S5进一步包括：

S5.1、利用卫星影像中的云顶压力数据。

S5.2、对于前一步骤S4识别出的象元，依次判断其云顶压力结果，云顶压力大于等于判据T4的象元为灰霾,其对应的掩膜值保持不变，而没有通过判据判断的象元的掩膜值乘以0，从而得到霾掩膜。

以上方法中所使用的云参数产品为卫星观测的云量、云顶压力、云光学厚度、云相态等数据，所有能够提供此类数据的卫星传感器都可以使用该方法识别灰霾范围。例如POLDER传感器，它搭载在法国PARASOL卫星上，其二级辐射收支产品中包含云量、云光学厚度、云相态和云顶压力数据，同时包括成像的几何数据和图像数据，分辨率为18km，格式为文本格式。

法国的地表反射率偏振和方向传感器(POLDER)包括以下波段：443nm，490nm，565nm，670nm，763nm，765nm，865nm，910nm，1020nm。其中490nm，670nm和860nm为偏振通道，其余为非偏振通道。除了偏振观测以外，POLDER具有多角度观测的优势，同一象元至多可以获得16个角度的观测值。POLDER一级产品中存放基本观测数据，如所有非偏振波段的反射率、偏振波段的偏振辐射强度等，一级产品数据的分辨率为6km。二级产品按照类型分为陆地产品、海洋产品和辐射收支产品三种。二级产品是在一级产品的基础上生产而来，其分辨率为18km。本发明中使用的云参数产品属于POLDER二级辐射收支产品，数据自带象元的经纬度值、云参数的校正系数等，数据格式为文本格式。

以通过POLDER卫星影像数据分析霾分布为例，步骤如下：

S1：选定研究区域，这里以雾霾现象常见的华北平原为例。筛选2012年1月至2013年1月期间华北平原灰霾事件及这些事件对应的日期，下载这些日期的POLDER卫星二级云产品，包括云量、云光学厚度、云相态和云顶压力数据，组建样本集。用人工目视解译的方法对卫星影像中的云和灰霾区域进行区分，分别对云、灰霾区域的上述云参数结果进行采样统计，确定POLDER云参数识别灰霾的判据T1～T4。此外，下载2013年1月至2013年3月期间华北雾霾天气对应的POLDER二级云产品数据，作为试验数据。步骤S1进一步包括：

S1.1、从PARASOL网站下载2012年一月到2013年一月全球范围内POLDER二级P3L2TRGB数据，分辨率为18km。根据新闻、天气预报筛选出典型霾天气的区域，作为主要试验区域。利用C程序语言读取该数据中的图像高度和宽度，象元经纬度数据以及云量数据。重新生成与输入数据等高和等宽的Geotiff影像，所有象元的初始赋值为1，作为POLDER灰霾范围识别的掩膜。针对待研究区域，筛选该区域范围内灰霾事件及这些事件对应的日期，下载这些日期的卫星云量、云光学厚度、云相态和云顶压力数据，组建样本集。

S1.2、用人工目视解译的方法对样本集中影像的云和灰霾区域进行区分，分别对云、灰霾、晴空区域的上述云参数结果进行采样统计。

S1.3、统计灰霾或云覆盖区域的云量数据、晴空区域的云量数据，确定云量判据T1＝25％。

S1.4、统计灰霾区域的云光学厚度数据、云区域的光学厚度数据，确定云光学厚度判据T2＝0.55。

S1.5、统计灰霾区域各种云相态出现的频率(水云、冰云、冰水混合云、无法确定类型、晴空)，频率超过R＝10％的云相态类型即为T3，实际数据分析，该云相态判据为冰云、冰水混合云、无法确定类型三种类型，即T3＝{冰云、冰水混合云、无法确定类型}。

S1.6、统计灰霾区域的云顶压力数据、云区域的云顶压力数据，确定云顶压力判据T4＝500hPa。

S2、确定云和灰霾共存的范围。使用C语言读取当期POLDER云参数产品中云量数据，并通过和对应的判据T1进行比较，确定云和灰霾共存的区域。步骤S2进一步包括：

S2.1、针对已经确定霾识别判据的区域，下载卫星云参数产品，包括云量、云顶压力、云相态和云光学厚度数据。

S2.2、利用C语言读取云参数产品中的图像高度和宽度，象元经纬度数据以及云量数据。重新生成与输入数据等高和等宽的影像，像素值均为1，作为霾识别的掩膜。这里所说的掩膜又可以称为霾掩膜，是指用来对待识别图像进行霾识别的特定图像，卫星影像和掩膜进行对应像素值相乘处理后，是霾的象元数据被保留，非霾的象元像素值则变为0。

S2.3、读取卫星影像中的云量数据。

S2.3、利用逐个象元的云量结果和对应的判据T1＝25％进行比较，小于或等于T1的象元为云和灰霾共存的象元，该象元在掩膜中对应的值保持不变，而没有通过判据判断的象元的掩膜值乘以0，作为步骤S2的提取结果。

S3、利用云光学厚度数据识别S1结果中灰霾的分布。针对S1选出的区域，利用当期POLDER云光学厚度数据和对应的判据T2进行比较，进一步提取霾分布的区域。S3步骤进一步包括：

S3.1、读取遥感影像中的云光学厚度数据数据。

S3.2、对于步骤S2识别出的象元，将云光学厚度结果和对应的判据T2＝0.55进行比较，大于或等于T2的象元为霾，其对应的掩膜中的象元值保持不变，而没有通过判据判断的象元的掩膜值乘以0，作为步骤S3的提取结果。

S4、利用云相态数据识别S2中的灰霾分布。针对S2获取的区域，利用当期P3L2TRGB云相态数据和对应的判据T3进行比较，进一步提取霾分布的区域。步骤S4进一步包括：

S4.1、利用C语言读取卫星影像中的云相态数据，即每个象元的在以下几类中的归属：水云、冰云、冰水混合云、晴空或者无法确定相态。

S4.2、对于步骤S3识别出的象元，逐个判断其云相态类别。云相态类别与判据T3＝{冰云，冰水混合云，无法确定类型}中的类别一致的象元，其对应的掩膜值保持不变，而没有通过判据判断的象元的掩膜值乘以0，作为步骤S4的提取结果。

S5、利用云顶压力数据识别S3中的灰霾分布。在步骤S3提取的霾区域的基础上，利用当期P3L2TRGB云顶压力数据和对应的判据T4进行比较，进一步提取霾分布的区域。步骤S5进一步包括：

S5.1、读取卫星影像中的云顶压力数据。步骤S5.2对于步骤S4识别出的象元，依次判断其云顶压力结果，大于等于判据T4＝500hPa的象元,其对应的掩膜值保持不变，而没有通过判据判断的象元的掩膜值乘以0，从而得到霾掩膜。

以上步骤中使用的云产品数据都需要经过校正，校正系数可以从云参数产品中读取，具体的做法是将云产品数据乘以校正系数即可。本发明还可应用与其它能够获得云量、云光学厚度、云相态、云顶压力的卫星传感器，例如美国Terra卫星上的MODIS传感器、NPP卫星上的VIIRS传感器等。

以上实施方式仅用于解释本发明的方法，而并非对本发明的限制，有关技术领域的普通技术人员在不脱离本发明的精神和范围的情况下,还可以做出各种变化和变形，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (8)

1.一种基于卫星云参数产品确定灰霾分布的方法，其特征在于，包括如下步骤：

针对待研究区域，下载历史云参数产品中包含卫星云量、云光学厚度、云相态频率和云顶压力的云参数产品组建样本集，对样本集中影像的云和灰霾区域进行区分，分别对云、灰霾、晴空区域进行采样统计，确定对应于所述各个云参数的识别判据，分别为T1、T2、T3、T4；

利用当期云参数产品中云量数据和对应的判据T1进行比较，确定云和灰霾都可能存在的区域；

在其它步骤选出的区域，利用当期云参数产品中云光学厚度数据和对应的判据T2进行比较，进一步提取霾分布的区域；

在其它步骤选出的区域，利用当期云参数产品中云相态数据和对应的判据T3进行比较，进一步提取霾分布的区域；

在其它步骤选出的区域，利用当期云参数产品中云顶压力数据和对应的判据T4进行比较，进一步提取霾分布的区域；

其中，确定所述云参数判据T1、T2、T3、T4的方法是

统计灰霾或云覆盖区域的云量数据，构建“灰霾或云”频率随云量变化的拟合曲线；统计晴空区域的云量数据，构建晴空出现频率随云量变化的拟合曲线,两曲线交点的云量值作为云量判据T1；

统计灰霾区域的云光学厚度数据，构建灰霾频率随云光学厚度变化的拟合曲线；统计云区域的光学厚度数据，构建云出现频率随云光学厚度变化的拟合曲线，两曲线交点的云光学厚度值为作为云光学厚度判据T2；

统计灰霾区域各种云相态出现的频率，超过一定比例R的云相态类型作为识别判据T3；

统计灰霾区域的云顶压力数据，构建灰霾出现频率随云顶压力变化的拟合曲线；统计云区域的云顶压力数据，构建云出现频率随云顶压力变化的拟合曲线，两曲线交点对应的云顶压力值作为云顶压力判据T4。

2.如权利要求1所述基于卫星云参数产品确定灰霾分布的方法，其特征在于，

云量大于或等于T1的象元为云或灰霾；

云光学厚度大于或等于T2的象元为灰霾；

云相态类别与T3一致的象元为灰霾；

云顶压力大于等于判据T4的象元为灰霾。

3.如权利要求1所述基于卫星云参数产品确定灰霾分布的方法，其特征在于，

所述云参数产品需要经过校正，校正系数从云参数产品中读取，将云参数产品的数据乘以校正系数。

4.如权利要求1～3任一项所述基于卫星云参数产品确定灰霾分布的方法，其特征在于，所述云参数产品为美国Terra卫星上的MODIS传感器获得的数据。

5.如权利要求1～3任一项所述基于卫星云参数产品确定灰霾分布的方法，其特征在于，所述云参数产品为NPP卫星上的VIIRS传感器获得的数据。

6.如权利要求1～3任一项所述基于卫星云参数产品确定灰霾分布的方法，其特征在于， 所述云参数产品为POLDER卫星影像，包括成像的几何信息和图像信息，具体为搭载在PARASOL卫星上的POLDER3传感器的CCD数据，该方法使用的是POLDER二级辐射收支产品中的云量、云光学厚度、云相态和云顶压力数据。

7.如权利要求6所述基于卫星云参数产品确定灰霾分布的方法，其特征在于，所述二级辐射收支产品，分辨率为18km，格式为文本格式。

8.如权利要求7所述基于卫星云参数产品确定灰霾分布的方法，其特征在于，所述云参数判据：

T1＝25％；

T2＝0.55；

T3＝{冰云，冰水混合云，无法确定类型}；

T4＝500hPa。