CN106290090A - 评估颗粒物浓度的方法和装置 - Google Patents

Abstract

本发明提供评估颗粒物浓度的方法和装置。该方法包括：a)分别在不同时间对地表区域进行拍摄以晴朗天气影像数据和雾霾天气影像数据；b)分别从晴朗天气影像数据和雾霾天气影像数据提取晴朗天气彩色图像数据和雾霾天气彩色图像数据；c)从晴朗天气彩色图像数据和雾霾天气彩色图像数据分别提取出晴朗天气暗通道数据和雾霾天气暗通道数据；d)计算出暗通道差值数据，并得到引导滤波数据；e)利用引导滤波数据，对标准气溶胶光学厚度值进行重采样，得出每一个像素的重采样的气溶胶光学厚度值；f)利用参数反演方法，确定定量反演关系式，并根据定量反演关系式和每一个像素的重采样的气溶胶光学厚度值，得出每一个像素的颗粒物浓度。

Description

评估颗粒物浓度的方法和装置

技术领域

本发明涉及遥感图像处理技术领域，具体的说是一种使用卫星影像来评估颗粒物浓度的方法和装置，适用于环境分析与处理领域。

背景技术

随着信息技术的发展，遥感技术越来越多的被用来获取和分析各种信息。利用遥感进行颗粒物浓度的监测已经成为研究的重点。

目前颗粒物浓度监测技术主要包括地面监测网络监测和卫星遥感监测。由于卫星遥感图像能够监测更大的连续的空间范围，光谱分辨率也越来越高，能够提取的信息也越来越多，所以使用卫星遥感来监测颗粒物浓度的技术研究有着更为重要的意义。卫星遥感监测主要使用多光谱传感器或中分辨率成像光谱仪获取的可见光及红外波段的数据为基础，处理得到气溶胶光学厚度(AOD)，并进行参数反演得到颗粒物浓度在分析区域的分布。《利用细模态气溶胶光学厚度估计PM2.5》对雾霾污染严重时的气溶胶光学厚度与颗粒物污染直接的相关关系进行了分析，并利用太阳光度计在地面观测得到的气溶胶光学厚度数据进行了模型估算和验证。《Retrieval，Validation，and Application of the 1-km Aerosol OpticalDepth From MODIS Measurements Over Hong Kong》使用MODIS原始影像为初始数据，使用太阳光度计的地面观测数据辅助处理得到1km分辨率的AOD，以弥补标准的AOD产品分辨率太低(10km)，在城市区域监测的不足。

上述从遥感图像对颗粒物进行分析的方法，能够对对颗粒物浓度进行反演而获取，但是对于较小尺度的城市区域，地面观测条件苛刻，难以大量复制使用，而标准AOD产品则分辨率太低，无法呈现颗粒物浓度差异。

发明内容

为了解决基于遥感技术对颗粒物浓度的分析中存在的不易适应较小尺度区域分析的问题，本发明提供了一种评估颗粒物浓度的方法，该方法包括以下步骤：

a)利用卫星传感器，分别在不同时间对地表区域进行拍摄以获得所述地表区域的晴朗天气影像和雾霾天气影像，所述地表区域中的多个像素在所述晴朗天气影像中各自的像素值组成晴朗天气影像数据，所述多个像素在所述雾霾天气影像中各自的像素值组成雾霾天气影像数据；

b)分别从所述晴朗天气影像数据和所述雾霾天气影像数据提取所述多个像素中的每一个像素在三个光通道中的每一个光通道中的像素值，并分别作为晴朗天气彩色图像数据和雾霾天气彩色图像数据；

c)从所述晴朗天气彩色图像数据和所述雾霾天气彩色图像数据分别提取出所述每一个像素在晴朗天气暗通道中的像素值以及所述每一个像素在雾霾天气暗通道中的像素值，并分别作为晴朗天气暗通道数据和雾霾天气暗通道数据；

d)计算所述雾霾天气暗通道数据和所述晴朗天气暗通道数据的差值，作为暗通道差值数据，并通过所述暗通道差值数据和所述雾霾天气暗通道数据得出所述每一个像素在引导滤波后的像素值，作为引导滤波数据；

e)利用所述引导滤波数据，对从所述雾霾天气影像数据计算得到的像素的标准气溶胶光学厚度值进行重采样，得出所述每一个像素的重采样的气溶胶光学厚度值；

f)利用参数反演方法，从像素的颗粒物浓度随所述重采样的气溶胶光学厚度值变化的多个函数关系式中，确定一个函数关系式作为定量反演关系式，并根据所述定量反演关系式和所述每一个像素的重采样的气溶胶光学厚度值，得出所述每一个像素的颗粒物浓度。

在步骤b)中，所述三个光通道包括红光通道、蓝光通道和绿光通道。

通过关系式(1)，得到所述晴朗天气暗通道数据和所述雾霾天气暗通道数据，

J^dark(x)＝min_y∈Ω(x)(min_c∈{r,g,b}J^c(y)) 关系式(1)

其中，Ω(x)表示以所述多个像素中的第x个像素为中心的预定的像素区域，y是所述像素区域中的第y个像素，c是所述红光通道r、蓝光通道b和绿光通道g中的一个光通道，J^c(y)表示所述晴朗天气彩色图像数据或所述雾霾天气彩色图像数据中在所述像素区域中的所述第y个像素在光通道c中的像素值，J^dark(x)是所述第x个像素在所述晴朗天气暗通道中的像素值或者在所述雾霾天气暗通道中的像素值，其中，J^dark(x)是所述像素区域中的所有像素在在所述三个光通道中的最小像素值。

如此，晴朗天气暗通道数据和雾霾天气暗通道数据可以作为雾霾分布特征的参考信息。

所述像素区域是根据所述地表区域中的地物的尺寸而预定的。也就是说，根据影像中景物的区域分布特征来设置该像素区域的范围和大小，从而可以使相同类型景物的像素落入同一像素区域中。

在步骤d)中，将所述每一个像素在所述雾霾天气暗通道中的像素值减去所述每一个像素在所述晴朗天气暗通道中的像素值，得到所述每一个像素的暗通道差值，作为所述暗通道差值数据，并且根据以下等式(2)得出所述引导滤波数据，

$q_i=\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}(a_k{I}_i+b_k)$ 等式(2)

其中，q_i表示所述多个像素中的第i个像素在引导滤波后的像素值，ω_k是以所述多个像素中的第k个像素为中心的滤波像素区域，|ω|是所述滤波像素区域中的像素个数，I_i是所述雾霾天气暗通道数据中的所述第i个像素的像素值，a_k和b_k是所述滤波像素区域中的线性因数，并通过以下等式(3)、(4)来分别获得a_k和b_k，

$a_k=\frac{\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}{I}_i{p}_i-{\mathrm{\μ}}_k\stackrel{\‾}{p_k}}{{\mathrm{\σ}}_k^2+\mathrm{\ϵ}}$ 等式(3)

$b_k=\stackrel{\‾}{p_k}-a_k{\mathrm{\μ}}_k$ 等式(4)

其中，p_i是所述暗通道差值数据中所述第i个像素的差值，μ_k和σ_k分别是所述雾霾天气暗通道数据中的所述滤波像素区域中的所有像素值的平均值和方差，是所述暗通道差值数据中的所述滤波像素区域中的所有差值的平均值，ε是调节参数。

由于暗通道的亮度值为一个很小的值，所以雾霾成为暗通道的主要成分、因为雾霾颗粒物会对光线在大气传输时产生较强的漫反射和散射作用，所以晴朗天气的暗通道更符合景物本来的暗通道结果，而雾霾天气暗通道数据则是附加了雾霾的影响。因此，暗通道差值数据可用于表示雾霾的分布差异。

由于在经过相减算法后，暗通道差值数据在景物边缘区域有较多噪声误差，所以，通过上述步骤，可以进行平衡消减，着重在景物边缘区域进行平滑滤波。

标准气溶胶光学厚度值有M×N个，每个标准气溶胶光学厚度值对应于所述多个像素中的m×n个像素,M、N、m、n为大于1的整数，

并且，通过以下等式(5)得到所述每一个像素的重采样的气溶胶光学厚度值，

${\mathrm{AOD}}_i={\mathrm{AOD}}_t*\frac{q_{i\∈w\left(t\right)}}{{\mathrm{\Σ}}_{s=1}^{m\×n}{q}_{s\∈w\left(t\right)}/(m\×n)}$ 等式(5)

AOD_i是所述第i个像素的重采样的气溶胶光学厚度值，AOD_t是第t个标准气溶胶光学厚度值，w(t)是与所述第t个标准气溶胶光学厚度值对应的所述m×n个像素的区域，q_s是在区域w(t)中的第s个像素在引导滤波后的像素值，1≤i≤M×N×m×n，1≤t≤M×N，s、t为整数。

如此，可以提高空间分辨率。

在步骤f)中，所述参数反演方法如下：将至少三个像素各自的重采样的气溶胶光学厚度值和所述至少三个像素各自的颗粒物浓度的实际测量值，分别带入所述多个函数关系式中进行反演计算，以确定每个函数关系式中的各个系数值，

所述定量反演关系式的确定如下：利用所述至少三个像素各自的颗粒物浓度的实际测量值、根据确定了所述各个系数值的所述每个函数关系式计算出的各个颗粒物浓度值、以及所述至少三个像素各自的颗粒物浓度的实际测量值的平均值，计算每个函数关系式的拟合度，并将拟合度的值最大的一个函数关系式确定作为所述定量反演关系式。

所述多个函数关系式包括线性函数关系式、多项式函数关系式、指数函数关系式、对数函数关系式、幂函数关系式。

本发明还提供一种评估颗粒物浓度的装置，该装置利用如上所述的评估颗粒物浓度的方法来评估所述每一个像素的颗粒物浓度。

通过本发明的评估颗粒物浓度的方法和装置，可以提高颗粒物浓度的分辨率，且该方法和装置能够适用于较小尺寸的城市区域，对于地面观测条件不苛刻，且方便大量复制使用，可以呈现颗粒物浓度的差异。

本发明可以以较低的成本来计算得到能够使用于城市级雾霾分析的具有较高空间分辨率的颗粒物浓度信息。无需经过传统暗像元法的暗区域实地选取及光谱测量，可改善空间分辨率，并且可以显著降低数据获取成本。

附图说明

图1是根据本发明的实施例的评估颗粒物浓度的方法的流程图；

图2(a)显示了根据本发明的实施例的9个像素在红光通道中的像素值；

图2(b)显示了根据本发明的实施例的9个像素在蓝光通道中的像素值；

图2(c)显示了根据本发明的实施例的9个像素在绿光通道中的像素值；

图3显示了根据本发明的实施例的晴朗天气暗通道数据的一部分；

图4显示了根据本发明的实施例的雾霾天气暗通道数据的一部分；

图5显示了根据本发明的实施例的暗通道差值数据的一部分。

具体实施方式

下面结合实施例，对本发明进行详细说明。

图1是根据本发明的实施例的评估颗粒物浓度的方法的流程图。在步骤S11中，利用卫星传感器，分别在不同时间对地表区域进行拍摄以获得所述地表区域的晴朗天气影像和雾霾天气影像，地表区域中的多个像素在晴朗天气影像中各自的像素值组成晴朗天气影像数据，多个像素在雾霾天气影像中各自的像素值组成雾霾天气影像数据。本实施例中，卫星传感器是MODIS，分辨率为1km，晴朗天气影像和雾霾天气影像的大小均为500×500个像素，该地表区域位于东经114.84度至121.20度，北纬32.51度至37.53度。晴朗天气影像的拍摄时间2014年11月23日上午10：45，雾霾天气影像的拍摄时间2013年12月20日上午10：45。影像中该地表区域主要由云层、裸露地表(陆地、海洋)、雾霾覆盖层组成。由于季节时间相差不大，所以地物覆盖情况相差不大。

这里，晴朗天气影像中的每个像素值组成晴朗天气影像数据，雾霾天气影像的每个像素值组成雾霾天气影像数据。也就是说，晴朗天气影像数据和雾霾天气影像数据分别有500×500个像素。

在步骤S12中，分别从晴朗天气影像数据和雾霾天气影像数据提取多个像素中的每一个像素在三个光通道中的每一个光通道中的像素值，并分别作为晴朗天气彩色图像数据和雾霾天气彩色图像数据。本实施例中，三个光通道包括0.65μm红光通道、0.47μm蓝光通道和0.55μm绿光通道。

图2(a)、图2(b)和图2(c)是晴朗天气彩色图像数据中的一部分，分别显示了500×500个像素中处于相同位置的9个像素分别在三个光通道中的像素值。图2(a)显示了该9个像素在0.65μm红光通道中各自的像素值，图2(b)显示了该9个像素在0.47μm蓝光通道中各自的像素值，图2(c)显示了该9个像素在0.55μm绿光通道中各自的像素值。晴朗天气彩色图像数据包括从晴朗天气影像数据提取的每个像素在每个光通道中的像素值，雾霾天气彩色图像数据包括从雾霾天气影像数据提取的每个像素在每个光通道中的像素值。

接下来，在步骤S13中，从晴朗天气彩色图像数据和雾霾天气彩色图像数据分别提取出每一个像素在晴朗天气暗通道中的像素值以及每一个像素在雾霾天气暗通道中的像素值，并分别作为晴朗天气暗通道数据和雾霾天气暗通道数据。

这里，通过关系式(1)，得到该晴朗天气暗通道数据和雾霾天气暗通道数据，

J^dark(x)＝min_y∈Ω(x)(min_c∈{r,g,b}J^c(y)) 关系式(1)

其中，Ω(x)表示以多个像素中的第x个像素为中心的预定的像素区域，y是该像素区域中的第y个像素，c是红光通道r、蓝光通道b和绿光通道g中的一个光通道，J^c(y)表示晴朗天气彩色图像数据或雾霾天气彩色图像数据中在像素区域中的第y个像素在光通道c中的像素值，J^dark(x)是第x个像素在晴朗天气暗通道中的像素值或者在雾霾天气暗通道中的像素值，其中，J^dark(x)是该像素区域中的所有像素在三个光通道中的最小像素值。

本实施例中，举例说明提取晴朗天气暗通道数据的过程。预定的像素区域Ω(x)例如包含3×3个像素，参见图2(a)、图2(b)和图2(c)，其中显示了在该像素区域Ω(x)中的9个相同位置的像素在红光通道r、蓝光通道b和绿光通道g中各自的像素值。该像素区域Ω(x)是以第520个像素为中心，即图2(a)、图2(b)和图2(c)中第2行第2列的位置处的像素，且x＝520。在像素区域Ω(520)中，一共有9个像素，按从左到右、从上到下的顺序依次为第1、第2、、、第9个像素。从图2(a)、图2(b)和图2(c)可以看出，第1个像素在红光通道r中的像素值为5，在蓝光通道b中的像素值为86，在绿光通道g中的像素值为45；第5个像素(即，多个像素中的第520个像素)在红光通道r中的像素值为26，在蓝光通道b中的像素值为56，在绿光通道g中的像素值为123，等等。如此，可以得到9个像素各自在三个光通道中的像素值。可以看出，图2(a)的红光通道r中的9个像素值中的最小值是5，图2(b)的蓝光通道b中的9个像素值中的最小值是32，图2(c)的绿光通道g中的9个像素值中的最小值是22，那么三个光通道中的最小值为5。

如此，通过关系式(1)可以得出J^dark(520)为5。也就是说，J^dark(520)是像素区域Ω(520)中的9个像素在图2(a)、图2(b)和图2(c)显示的三个光通道中的最小像素值。

按照上述关系式(1)，对其他像素进行提取，以最终获得每个像素在晴朗天气暗通道中的像素值，并作为晴朗天气暗通道数据。同样，按照上述关系式(1)，可以获得每个像素在雾霾天气暗通道中的像素值，并作为雾霾天气暗通道数据。

图3显示了晴朗天气暗通道数据的一部分，显示了与图2(a)、图2(b)和图2(c)中的9个相同位置的像素各自在晴朗天气暗通道中的像素值。图4显示了雾霾天气暗通道数据的一部分，显示了与图2(a)、图2(b)和图2(c)中的9个相同位置的像素各自在雾霾天气暗通道中的像素值。

另外，像素区域Ω(x)是根据地表区域中的地物的尺寸而预定的，且包含(2X+1)×(2X+1)个像素，X是大于等于1的整数。

在步骤S14中，计算雾霾天气暗通道数据和晴朗天气暗通道数据的差值，作为暗通道差值数据，并通过暗通道差值数据和雾霾天气暗通道数据得出每一个像素在引导滤波后的像素值，作为引导滤波数据。

其中，将每一个像素在雾霾天气暗通道中的像素值减去每一个像素在晴朗天气暗通道中的像素值，得到每一个像素的暗通道差值，作为暗通道差值数据。

具体的，参见图3和图4，对于本实施例中的这9个像素中的第1个像素，得到该第1个像素的暗通道差值为20-5＝15。对于第9个像素，是1-2＝-1，那么第9个像素的暗通道差值就取0。也就是说，如果某个像素在雾霾天气暗通道中的像素值减去其在晴朗天气暗通道中的像素值的差为负数，那么该像素的暗通道差值为0。

如此，可以得到这9个像素的暗通道差值，如图5所示，且图5显示了暗通道差值数据的一部分。

然后，根据以下等式(2)得出引导滤波数据，

$q_i=\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}(a_k{I}_i+b_k)$ 等式(2)

其中，q_i表示500×500个像素中的第i个像素在引导滤波后的像素值，ω_k是以500×500个像素中的第k个像素为中心的滤波像素区域，|ω|是该滤波像素区域中的像素个数，I_i是雾霾天气暗通道数据中的第i个像素的像素值，a_k和b_k是该滤波像素区域中的线性因数，并通过以下等式(3)、(4)来分别获得a_k和b_k，

$a_k=\frac{\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}{I}_i{p}_i-{\mathrm{\μ}}_k\stackrel{\‾}{p_k}}{{\mathrm{\σ}}_k^2+\mathrm{\ϵ}}$ 等式(3)

$b_k=\stackrel{\‾}{p_k}-a_k{\mathrm{\μ}}_k$ 等式(4)

其中，p_i是暗通道差值数据中第i个像素的差值，μ_k和σ_k分别是雾霾天气暗通道数据中的位于该滤波像素区域ω_k中的所有像素值的平均值和方差，是暗通道差值数据中的位于该滤波像素区域ω_k中的所有差值的平均值，ε是调节参数。

参见图4和图5，本实施例中，例如k取520，即，ω₅₂₀是以第520个像素为中心的滤波像素区域，|ω|＝9。I_i是雾霾天气暗通道数据中的第i个像素的像素值，从图4可以看出，雾霾天气暗通道数据中的第520个像素I₅₂₀的像素值为32；p_i是暗通道差值数据中第i个像素的差值，从图5可以看出，暗通道差值数据中第520个像素的差值p₅₂₀＝27。μ_k和σ_k分别是雾霾天气暗通道数据中的位于该滤波像素区域ω_k中的像素值的平均值和方差。本例中，从图4可以得出，λ₅₂₀＝12.55，和σ₅₂₀＝7.29。在本实施例中，ε例如是0.4。本例中，是暗通道差值数据中的位于滤波像素区域ω₅₂₀中的9个差值的平均值，并从图5的9个差值可以计算得出。如此，根据上述等式(3)和(4)，可以分别计算得出一组a₅₂₀和b₅₂₀。

本例中，第520个像素也被包含在以第519个像素为中心的滤波像素区域ω₅₁₉中，因此，对于I₅₂₀，可以根据图4、图5中的数据，从等式(3)和等式(4)计算得出另一组a₅₁₉和b₅₁₉。本例中，滤波像素区域包含9个像素，那么对于第520个像素(I₅₂₀)，可以计算得出9组a_k和b_k。将这9组a_k和b_k代入上述等式(2)，从而计算得出q₅₂₀的值。按照上述方式，可以计算出500×500个像素中的每个像素在引导滤波后的像素值，从而得到引导滤波数据。

接着，在步骤S15中，利用上述得到的引导滤波数据，对从雾霾天气影像数据计算得到的像素的标准气溶胶光学厚度值进行重采样，得出每一个像素的重采样的气溶胶光学厚度值。

这里，从雾霾天气影像数据计算像素的标准气溶胶光学厚度值的方法和现有技术中的一样，因此这里不再详述。

本例中，标准气溶胶光学厚度值有M×N个，每个标准气溶胶光学厚度(AOD)值对应于500×500个像素中的m×n个像素,M、N、m、n为大于1的整数。

如上所述，晴朗天气影像和雾霾天气影像的分辨率都为1km。本例中，标准AOD值的分辨率例如为10km，也就是说，一共有50×50个标准AOD值，即，M＝N＝50，且每个标准AOD值对应于10×10个像素，即m＝n＝10。

然后，通过以下等式(5)得到每一个像素的重采样的气溶胶光学厚度值，

${\mathrm{AOD}}_i={\mathrm{AOD}}_t*\frac{q_{i\∈w\left(t\right)}}{{\mathrm{\Σ}}_{s=1}^{m\×n}{q}_{s\∈w\left(t\right)}/(m\×n)}$ 等式(5)

AOD_i是第i个像素的重采样的气溶胶光学厚度值，AOD_t是第t个标准气溶胶光学厚度值，w(t)是与第t个标准气溶胶光学厚度值对应的10×10个像素的区域，q_s是在区域w(t)中的第s个像素在引导滤波后的像素值，1≤i≤M×N×m×n，1≤t≤M×N，s、t为整数。本例中，1≤i≤500×500，1≤t≤50×50。

例如，以第520个像素为例，第2个标准AOD值对应的10×10个像素的区域中包含第520个像素。根据步骤S14中得出的引导滤波数据，并从上述等式(5)可以计算出AOD₅₂₀。如此，可以计算出500×500个像素中的每个像素的重采样的AOD值。

接下来，在步骤S16中，利用参数反演方法，从像素的颗粒物浓度随重采样的气溶胶光学厚度值变化的多个函数关系式中，确定一个函数关系式作为定量反演关系式，并根据定量反演关系式和每一个像素的重采样的气溶胶光学厚度值，得出每一个像素的颗粒物浓度。

多个函数关系式包括线性函数关系式、多项式函数关系式、指数函数关系式、对数函数关系式、幂函数关系式，当然还可以有其他函数关系式。

本实施例中，这些函数关系式例如如下：

线性函数： y＝a₁x+b₁

多项式函数： y＝a₂x²+b₂x+c₂

指数函数： y＝a₃e^x+b₃

对数函数： y＝a₄ln(x)+b₄

幂函数： $y=a_5{x}^{b_5}+c_5$

具体的，该参数反演方法如下：将根据步骤S15得出的至少三个像素各自的重采样的气溶胶光学厚度值和该至少三个像素各自的颗粒物浓度的实际测量值，分别代入上述多个函数关系式中进行反演计算，以确定每个函数关系式中的各个系数值。

这里，当实际测量值大于3组(即，3个像素的重采样的气溶胶光学厚度值和各自的颗粒物浓度的实际测量值)时，可以通过最小二乘法和幂函数回归计算等拟合方法求最优解。实际测量值越多，得到的反演关系越稳定和精确。通过解方程组得到各函数关系式中的系数a、b、c。

本实施例中，计算得到的具体函数关系式如下：

指数函数：y＝74.982×e^0.8422x

线性函数：y＝110.52x+71.588

对数函数：y＝51.599ln(x)+173.62

多项式函数：y＝-26.649x²+148.49x+61.014

幂函数：y＝168.05×x^0.4405

然后，定量反演关系式的确定如下：利用至少三个像素各自的颗粒物浓度的实际测量值、根据确定了各个系数值的每个函数关系式计算出的各个颗粒物浓度值、以及至少三个像素各自的颗粒物浓度的实际测量值的平均值，计算每个函数关系式的拟合度，并将拟合度的值最大的一个函数关系式确定作为所述定量反演关系式。

本实施例中，按照以下等式(6)来计算每个函数关系式中的拟合度。

$R_T=1-\frac{\underset{V=1}{\overset{Z}{\mathrm{\Σ}}}{(y_V-y_{\mathrm{VY}}^{*})}^2}{\underset{V=1}{\overset{Z}{\mathrm{\Σ}}}{(y_V-\stackrel{\‾}{y})}^2}$ 等式(6)

其中，R_T表示上述确定各系数值的5个函数关系式中的第T个函数关系式的拟合度，1≤T≤5，Z表示上述至少三个像素的个数，例如Z＝3，y_V是至少三个像素中的第V个像素的颗粒物浓度的实际测量值，1≤V≤Z，表示将第V像素的重采样的气溶胶光学厚度值作为自变量x代入根据上述确定各系数值的5个关系函数式中的第T个函数关系式计算得出的值，表示至少三个像素各自的颗粒物浓度的实际测量值的平均值。T、V、Z为整数。

本实施例中，根据上述等式(6)，得到5个函数关系式的拟合度分别为：指数函数的拟合度为0.7501，线性函数的拟合度为0.7972，对数函数的拟合度为0.7909，多项式函数的拟合度为0.9203，幂函数的拟合度为0.7426。

其中，多项式函数的拟合度的值最大，所以多项式函数是最优反演关系式，并将该多项式函数确定作为定量反演关系式。

然后，根据该多项式函数：y＝-26.649x²+148.49x+61.014，将每个像素的重采样的气溶胶光学厚度值作为自变量x代入，从而得出每个像素的颗粒物浓度。

另外，本发明还提供一种使用上述评估颗粒物浓度的方法来评估每个像素的颗粒物浓度的装置。

通过本发明的评估颗粒物浓度的方法和装置，可以提高颗粒物浓度的分辨率，且该方法和装置能够适用于较小尺寸的城市区域，对于地面观测条件不苛刻，且方便大量复制使用，可以呈现颗粒物浓度的差异。

虽然本发明的特定实施例已被描述，但这些实施例只通过实例的方式进行表述，并不意欲限制本发明的范围。实际上，本文描述的创新方法可以通过各种其他形式实施；此外，也可以进行对本文描述的方法和***的各种省略、替代和改变而不背离本发明的精神。附后的权利要求及其等同内容的目的是涵盖落入本发明的范围和精神内的这样的各种形式或修改。

Claims (10)

1.一种评估颗粒物浓度的方法，其特征在于，所述方法包括以下步骤：

a)利用卫星传感器，分别在不同时间对地表区域进行拍摄以获得所述地表区域的晴朗天气影像和雾霾天气影像，所述地表区域中的多个像素在所述晴朗天气影像中各自的像素值组成晴朗天气影像数据，所述多个像素在所述雾霾天气影像中各自的像素值组成雾霾天气影像数据；

b)分别从所述晴朗天气影像数据和所述雾霾天气影像数据提取所述多个像素中的每一个像素在三个光通道中的每一个光通道中的像素值，并分别作为晴朗天气彩色图像数据和雾霾天气彩色图像数据；

c)从所述晴朗天气彩色图像数据和所述雾霾天气彩色图像数据分别提取出所述每一个像素在晴朗天气暗通道中的像素值以及所述每一个像素在雾霾天气暗通道中的像素值，并分别作为晴朗天气暗通道数据和雾霾天气暗通道数据；

d)计算所述雾霾天气暗通道数据和所述晴朗天气暗通道数据的差值，作为暗通道差值数据，并通过所述暗通道差值数据和所述雾霾天气暗通道数据得出所述每一个像素在引导滤波后的像素值，作为引导滤波数据；

e)利用所述引导滤波数据，对从所述雾霾天气影像数据计算得到的像素的标准气溶胶光学厚度值进行重采样，得出所述每一个像素的重采样的气溶胶光学厚度值；

f)利用参数反演方法，从像素的颗粒物浓度随所述重采样的气溶胶光学厚度值变化的多个函数关系式中，确定一个函数关系式作为定量反演关系式，并根据所述定量反演关系式和所述每一个像素的重采样的气溶胶光学厚度值，得出所述每一个像素的颗粒物浓度。

2.如权利要求1所述的方法，其特征在于，在步骤b)中，所述三个光通道包括红光通道、蓝光通道和绿光通道。

3.如权利要求2所述的方法，其特征在于，在步骤c)中，通过关系式(1)，得到所述晴朗天气暗通道数据和所述雾霾天气暗通道数据，

J^dark(x)＝min_y∈Ω(x)(min_c∈{r,g,b}J^c(y)) 关系式(1)

其中，Ω(x)表示以所述多个像素中的第x个像素为中心的预定的像素区域，y是所述像素区域中的第y个像素，c是所述红光通道r、蓝光通道b和绿光通道g中的一个光通道，J^c(y)表示所述晴朗天气彩色图像数据或所述雾霾天气彩色图像数据中在所述像素区域中的所述第y个像素在光通道c中的像素值，J^dark(x)是所述第x个像素在所述晴朗天气暗通道中的像素值或者在所述雾霾天气暗通道中的像素值，其中，J^dark(x)是所述像素区域中的所有像素在在所述三个光通道中的最小像素值。

4.如权利要求3所述的方法，其特征在于，所述像素区域是根据所述地表区域中的地物的尺寸而预定的。

5.如权利要求4所述的方法，其特征在于，在步骤d)中，将所述每一个像素在所述雾霾天气暗通道中的像素值减去所述每一个像素在所述晴朗天气暗通道中的像素值，得到所述每一个像素的暗通道差值，作为所述暗通道差值数据，并且根据以下等式(2)得出所述引导滤波数据，

$q_i=\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}(a_k{I}_i+b_k)$ 等式(2)

其中，q_i表示所述多个像素中的第i个像素在引导滤波后的像素值，ω_k是以所述多个像素中的第k个像素为中心的滤波像素区域，|ω|是所述滤波像素区域中的像素个数，I_i是所述雾霾天气暗通道数据中的所述第i个像素的像素值，a_k和b_k是所述滤波像素区域中的线性因数，并通过以下等式(3)、(4)来分别获得a_k和b_k，

$a_k=\frac{\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}{I}_i{p}_i-{\mathrm{\μ}}_k\stackrel{\‾}{p_k}}{{\mathrm{\σ}}_k^2+\mathrm{\ϵ}}$ 等式(3)

$b_k=\stackrel{\‾}{p_k}-a_k{\mathrm{\μ}}_k$ 等式(4)

其中，p_i是所述暗通道差值数据中所述第i个像素的差值，μ_k和σ_k分别是所述雾霾天气暗通道数据中的所述滤波像素区域中的所有像素值的平均值和方差，是所述暗通道差值数据中的所述滤波像素区域中的所有差值的平均值，ε是调节参数。

6.如权利要求5所述的方法，其特征在于，在步骤e)中，所述标准气溶胶光学厚度值有M×N个，每个标准气溶胶光学厚度值对应于所述多个像素中的m×n个像素,M、N、m、n为大于1的整数，

并且，通过以下等式(5)得到所述每一个像素的重采样的气溶胶光学厚度值，

${\mathrm{AOD}}_i={\mathrm{AOD}}_t*\frac{q_{i\∈w\left(t\right)}}{{\mathrm{\Σ}}_{s=1}^{m\×n}{q}_{s\∈w\left(t\right)}/(m\×n)}$ 等式(5)

AOD_i是所述第i个像素的重采样的气溶胶光学厚度值，AOD_t是第t个标准气溶胶光学厚度值，w(t)是与所述第t个标准气溶胶光学厚度值对应的所述m×n个像素的区域，q_s是在区域w(t)中的第s个像素在引导滤波后的像素值，1≤i≤M×N×m×n，1≤t≤M×N，s、t为整数。

7.如权利要求6所述的方法，其特征在于，在步骤f)中，所述参数反演方法如下：将至少三个像素各自的重采样的气溶胶光学厚度值和所述至少三个像素各自的颗粒物浓度的实际测量值，分别带入所述多个函数关系式中进行反演计算，以确定每个函数关系式中的各个系数值，

所述定量反演关系式的确定如下：利用所述至少三个像素各自的颗粒物浓度的实际测量值、根据确定了所述各个系数值的所述每个函数关系式计算出的各个颗粒物浓度值、以及所述至少三个像素各自的颗粒物浓度的实际测量值的平均值，计算每个函数关系式的拟合度，并将拟合度的值最大的一个函数关系式确定作为所述定量反演关系式。

8.如权利要求7所述的方法，其特征在于，所述滤波像素区域是根据所述地表区域中的地物的尺寸而预定的。

9.如权利要求8所述的方法，其特征在于，所述多个函数关系式包括线性函数关系式、多项式函数关系式、指数函数关系式、对数函数关系式、幂函数关系式。

10.一种评估颗粒物浓度的装置，其特征在于，所述装置利用如权利要求1-9中任一项所述的评估颗粒物浓度的方法来评估所述每一个像素的颗粒物浓度。