# 4.2.2 看图识雾霾

“厚德载雾，自强不吸，霾头苦干，再创灰黄”，雾霾不断肆虐，走上街头越来越多的人开始带上“防毒面具”。雾霾的主要成分是PM2.5，即细颗粒物，它指环境空气中空气动力学当量直径小于等于2.5微米的颗粒物。PM2.5能较长时间悬浮于空气中，其在空气中含量浓度越高，则空气污染越严重。根据中国环境监测总站发布的《2016年11月空气质量状况》，2016年11月，北京市空气质量优良天数比例仅43.3%；其余天数中，重度污染6天，严重污染1天，其首要污染物就是PM2.5。截至2016年11月底，全国共设1436个国家空气质量监测站点，由国家统一运行维护，监测数据直报国家并对外公开。但有检测站是否意味着PM2.5的情况被严格管控了呢？事实上并没有，全国大量的监测站点数据质量其实难以保证。能否通过其他方式来判断雾霾的情况？有人说，每天早上一出门，今天的天气状况也能猜个八九不离十。那有没有可能通过图片来判别雾霾情况呢？这也就是本节要介绍的内容啦！

在上节内容中，我们介绍了图像操作的基本知识，本节我们将结合雾霾图片的实际案例，为大家介绍如何从图片中提取我们想要的特征。在此案例中，我们将主要使用imager包中的函数来完成我们的主要工作。

数据：2016年12月到2017年3月从熊大办公室拍摄的窗外照片共206张照片，每张图片按拍摄时间命名。

    # 清空工作空间
    rm(list = ls())

    # 引入处理图像所需要的包
    # install.packages(knitr)
    require(knitr)
    # install.packages(RgoogleMaps)
    library(RgoogleMaps)
    # install.packages(imager)
    library(imager)
    # install.packages(jpeg)
    library(jpeg)
    # install.packages(ggplot2)
    library(ggplot2)

首先，我们要提取图像特征来作为回归函数中的自变量。如何从图像中提取与雾霾程度有关的特征呢？直观上，人们想到雾霾，就会联想到“伸手不见五指”的灰蒙蒙的景象。因此，我们从衡量图像清晰程度的角度出发，对图像特征进行观察和分析，得到如下解释性变量。为了直观的展示变量是否有效，我们首先将以空气质量状况迥异的两张图片（见下图）为例来展示变量值的差距。其次，我们以206张图片的特征取值为样本绘制了箱线图，使得结果更具有说服力。

左右两边的图片分别是在空气优良和雾霾严重的时候拍摄的，你能分辨出来吗（估计没啥问题）？你的分辨原因是啥？左边天气蓝，层次清晰，右边天气灰，一团浆糊；左边视觉更饱满，右边让人昏昏欲睡……可能这样的理由你可以说一箩筐。那么，从图片中，我们能提取出并量化这样的特征吗？让我们跟着感觉走，找找这些变量吧。

我们会以FigureData中的这两张图片来做一些展示。

    Image1 = load.image('FigureData/201612221230.jpg')  # 路径名＋图片名
    Image2 = load.image('FigureData/201612211230.jpg')

# 提取图像特征

## 1. 灰度差分的方差

第一个变量是灰度差分的方差。通过对图片观察我们可以发现，雾霾影响最明显的地方之一就是远处图物体的边缘，正如图片中的远方的树木的边缘和山的轮廓，天气晴朗时我们能清晰地分辨出这些景物，而到了雾霾时的图像就变得特别模糊。为了提取这种轮廓特征，我们就需要先将灰度差分。在“熊大图像处理”一节中我们已经介绍了灰度图片——灰度是指一张黑白图像中的颜色深度。对于彩色照片，灰度就是色彩（红绿蓝三色）的线性组合，可以认为亮部灰度值高，暗部灰度值低。灰度的变化就可以反映出图像中明暗的变化，雾霾越严重，图像越不清晰，灰度变化也就越小。我们进一步思考，灰度变化最剧烈的地方其实正是这些物体边界的地方（同一个区域的物体（绿树），因为很相似，灰度变化不大；但是绿树和蓝天交界处会有很大的变化）。我们先对灰度差分下一个定义：某一像素点的灰度差分等于下格像素点灰度值减去上格像素点灰度值。这个差分就可以帮助我们提取这种边缘的轮廓信息：

    DIF = function(M){
      # 函数说明：计算图像灰度差分的方差的函数，输入的是图像矩阵，输出图像的灰度差分
      M = grayscale(M, drop = TRUE);  # 转灰度图
      M = as.matrix(M)
      DIF = M[-1, ] - M[-nrow(M), ]  # 下格像素-上格像素
      return(DIF)
    }

    diffvar = function(M){
      # 函数说明：计算图像灰度差分的方差的函数，输入的是图像矩阵，输出图像的灰度差分的方差
      diffvarFeature = var(as.vector(DIF(M)))
      return(diffvarFeature)
    }

下面两幅图展示了差分后像素值的分布情况：

    dif1 = DIF(Image1)
    dif2 = DIF(Image2)
    par(mfrow = c(1, 2),family="STKaiti")
    hist(as.cimg(dif1), main = "晴朗天气",
         breaks = c(-1, -0.3, -0.25, -0.2, -0.15, -0.1, 0, 0.1, 0.15, 0.2, 0.25, 0.3, 1), ylim = c(0, 7),
         xlab = "差分后像素值")

    hist(as.cimg(dif2), main = "雾霾天气",
         breaks = c(-1, -0.3, -0.25, -0.2, -0.15, -0.1, 0, 0.1, 0.15, 0.2, 0.25, 0.3, 1), ylim = c(0, 7),
         xlab = "差分后像素值")

非常直观地，好天气图像中由于细节轮廓多，像素点里差分后的值所变动的范围就大（左图框内区域）；而对于雾霾天气，由于其细节会变少，导致大部分像素差分后的值都聚集在0左右。那么如何反应这种区别呢？统计学告诉我们，变量的方差是衡量变动的一个指标，我们只需要求灰度差分的方差即可。遍历图像中所有像素点求其灰度差分，而后计算灰度差分的方差，将其作为我们的一个变量。例如上图中晴朗天气的灰度差分为约3.9×〖10〗(-3)，而雾霾天气的灰度差分仅仅为0.9×〖10〗(-3)。

我们将不同的PM2.5等级对应灰度差分方差用箱线图展现出来。从图中可以看出，随着PM2.5污染程度的增加，灰度差分的方差呈减小趋势。

    ### 描述自变量 ###
    ## 污染水平vs灰度差分的方差 ##
    data = read.csv('Output.csv')   # 读入数据，含有从图片中提取好的自变量
    pm25 = read.csv('PM2.5数据.csv')   # 读入pm2.5数据
    data$pm25 = pm25$PM25
    data$污染水平 = cut(data$pm25,breaks = c(-Inf, 50, 100, 200, Inf), labels = c("优", "良", "轻度污染", "较严重污染"))

    ggplot(data, aes(x = factor(污染水平), y = diffvarB, fill = 污染水平, alpha = 污染水平)) +
      geom_boxplot(aes(fill = 污染水平)) +
      scale_fill_manual(values = c('chartreuse2', 'gold1', 'orange1', 'firebrick3')) + labs(x = "PM2.5", y = "灰度差分的方差") +
      scale_alpha_manual(values = c(0.5, 0.5, 0.5, 0.5)) +
      scale_x_discrete(limits = c("优", "良", "轻度污染", "较严重污染")) +
      theme(axis.title.x = element_text(size = 14), axis.title.y = element_text(size = 14), legend.text = element_text(size = 12), legend.title = element_text(size = 14), axis.text = element_text(size = 12), text = element_text())+
        theme_grey(base_family = "STKaiti")

## 2. 清晰度

第二个解释性变量是清晰度。这一变量来源于计算机视觉图像去雾领域的名作，它表示反射光未衰减的比例。

拍照的原理是物体的反射光被相机所接收，对于像素点x，设反射光为J(x)（不同景物有着不同的反射光），但它在传播过程中会有一定的衰减，我们将最后到达相机的剩余反射光的比例记为t(x)，即为我们定义的清晰度（其本质是透射系数）。最终到达相机的物体的反射光即为J(x)t(x)。除此之外，还有一部分光也会到达相机，这就是大气背景光，我们设其为A，到达的大气背景光占所有接收到的光的比例为(1-t(x))，即到达相机的大气背景光为A(1-t(x))。两者相加为相机收到的总光强，等于J(x)t(x)+A(1-t(x))。简单来说，由于雾霾天气对t(x)的影响较大，雾霾越大，t(x)越小（即景物反射光的比例越低，大气背景光比例越高），因此我们在这里可以估计t(x)，作为其中一个重要解释变量。根据文章中的暗通道先验结果，t(x)可以通过一定的公式估计。

下面的代码详细阐释了计算的过程：对于输入的图片M，维度为a*b，我们可以计算每个像素点在三个通道上的最小值，得到矩阵minMatrix，再对这个矩阵图像的每个像素点进行遍历操作，即在以每个小像素点(m,n)为中心的31*31小方格范围里找到最小值，并将其保存下来，得到一个新的维数为a*b的矩阵。因此我们称呼这个矩阵为图像的暗通道矩阵。

    clearRate = function(M) {
      # 函数说明：计算图像清晰度的函数，输入的是图像矩阵，输出图像的清晰度

      DIM = dim(M)
      minMatrix = matrix(NA, DIM[1], DIM[2])
      minMatrix = pmin(M[, , 1], M[, , 2], M[, , 3])  # 找到三通道中的最小值
      darkChannel = matrix(NA, DIM[1], DIM[2])
      for(j in 1:DIM[1]) {  # 再在每一个格子的周边31*31的格子中找到最小值
        for(k in 1:DIM[2]) {
          up = j - 15  # 上下左右搜索范围大小是15个像素点
          down = j + 15
          left = k - 15
          right = k + 15
          if(up < 1) {  # 当搜索范围超过图像边界时的处理，取搜索范围和图像范围的交集
            up = 1
         }
          if(down > DIM[1]) {
            down = DIM[1]
          }
          if(left < 1){
            left = 1
          }
          if(right > DIM[2]) {
            right = DIM[2]
          }
          darkChannel[j, k] = min(minMatrix[up:down, left:right])  # 得到（j, k）位置处的暗通道值
        }
      }

      # 利用文章中公式计算，每一点都有清晰度的值，形成分布，找出 70 分位数。
      clearRateFeature = quantile(1 - darkChannel/max(darkChannel), 0.7)  # 将A的估计简化为三通道都为相等常数
                                                                          # 且是暗通道中最大值的一点
      return(clearRateFeature)
    }

当然这个矩阵也可以用灰度图像表示出来，如下图：

暗通道

    clearRate(Image1)
    ##        70%
    ##  0.8676471
    clearRate(Image2)
    ##        70%
    ##  0.7531381

针对我们的展示图片，晴朗天气对应清晰度为0.87，而雾霾天对应清晰度为0.75。

箱线图的结果也证明了这是一个很好的解释变量：

    ## 污染水平vs清晰度 ##
    ggplot(data, aes(x = factor(污染水平), y = clearRate, fill = 污染水平, alpha = 污染水平)) +
      geom_boxplot(aes(fill = 污染水平)) +
      scale_fill_manual(values = c('chartreuse2', 'gold1', 'orange1', 'firebrick3')) + labs(x = "PM2.5", y = "清晰度") +
      scale_alpha_manual(values = c(0.5, 0.5, 0.5, 0.5)) +
      scale_x_discrete(limits = c("优", "良", "轻度污染", "较严重污染")) +
      theme(axis.title.x = element_text(size = 14), axis.title.y = element_text(size = 14), legend.text = element_text(size = 12), legend.title = element_text(size = 14), axis.text = element_text(size = 12), text = element_text())+
        theme_grey(base_family = "STKaiti")

## 3.信息熵

第三个解释性变量是信息熵。众所周知，熵是度量“变化性”的一种指标，当图片颜色变化比较明显时，熵值较大。因此，当天朗气清，惠风和畅时，我们能够看到蓝天白云、红花绿草。而当雾霾严重时，图像中的多数信息被雾霾所掩盖，因此残存的信息较少，我们也就只能看到“一桶浆糊”了。基于此，我们引入了信息熵这一变量，雾霾越严重，图像中的信息越少，则信息熵越小。直观上，空气污染严重时，雾霾掩盖了图像中的部分信息。图像看上去灰蒙蒙一片，难以辨认具体细节。这里我们用信息熵来量化图像的平均信息量。

另一方面，我们看到的图片是由R（Red，红色）、G（Green，绿色）以及B（Blue，蓝色）三个色彩通道组成的。虽然一般对于彩色图片，通常的做法是处理成灰度图片处理；但是，对于“看图识雾霾”而言，有没有哪个颜色通道是更加重要的呢？为了保留更多的彩色图像信息，我们首先将预处理图片分到R、G、B三个空间下，并分别进行了测试。结果符合我们的直觉，发现B空间下效果较好（毕竟，没有雾霾时蓝天白云中的蓝色成分还是很重要的），因此这里以B空间下的图像为例进行分析。

以上文中空气质量状况迥异的两张图片为例，第一张图像的空气质量状况为优，PM2.5=7，其信息熵为1.624；第二张图像的空气质量等级达到了较严重污染，PM2.5=436，其信息熵为0.820，明显小于第一张照片。

    entro = function(M) {
      # 函数说明：计算图像信息熵的函数，输入的是图像矩阵，输出图像的信息熵

      y = 30 / 255  # 阈值(可调优参数)

      B1 = M[, , 3]  # 蓝色通道的矩阵
      B2 = M[, , 3]  # 蓝色通道的矩阵
      DIM = dim(B1)  # 蓝色通道的矩阵的维度
      n = DIM[2]  # DIM[2]是蓝色通道的矩阵的列
      m = DIM[1]  # DIM[1]是蓝色通道的矩阵的行

      dif1 = abs(B1[(1:(m - 4)), (1:(n - 4))] - B1[(5:m), (5:n)])
      dif2 = abs(B1[(5:m), (1:(n - 4))] - B1[(1:(m - 4)), (5:n)])
      B2 = B1[(3:(m - 2)), (3:(n - 2))]
      B2[which(dif1 < y & dif2 < y) ] = 0
      temp1 = B1
      temp1[(3:(m - 2)), (3:(n - 2))] = B2
      B2 = temp1  # temp1的作用就只是为了传值

      # 计算信息熵
      t = table(B2)  # 统计每种结果出现多少次
      result = names(t)
      freq = as.numeric(t)
      sum = sum(freq)  # 总次数
      entropy = -sum(log(freq / sum) * (freq / sum))

      return(entropy)
    }

我们展示一下示例图片的信息熵计算：

    entro(Image1)
    ##  [1] 1.623529
    entro(Image2)
    ##  [1] 0.8200629

此外，从箱线图可以看出，随着PM2.5污染程度的增加，信息熵呈减小趋势。

    ## 污染水平vs信息熵 ##
    ggplot(data, aes(x = factor(污染水平), y = entropy, fill = 污染水平, alpha = 污染水平)) +
      geom_boxplot(aes(fill = 污染水平)) +
      scale_fill_manual(values = c('chartreuse2', 'gold1', 'orange1', 'firebrick3')) + labs(x = "PM2.5", y = "信息熵") +
      scale_alpha_manual(values = c(0.5, 0.5, 0.5, 0.5)) +
      scale_x_discrete(limits = c("优", "良", "轻度污染", "较严重污染")) +
      theme(axis.title.x = element_text(size = 14), axis.title.y = element_text(size = 14), legend.text = element_text(size = 12), legend.title = element_text(size = 14), axis.text = element_text(size = 12), text = element_text())+
        theme_grey(base_family = "STKaiti")

## 4. 饱和度

第四个解释性变量是饱和度。饱和度指图像中色彩的鲜艳程度，例如下图这束花的照片，调高它的饱和度，可以看到花朵确实变得娇艳欲滴了。相信一些读者朋友也发现了，其实饱和度与图像中灰色成分所占的比例有关。雾霾会使图像变灰，色彩鲜艳程度减弱；雾霾越严重，图像饱和度越低。

为了获得图片的饱和度，我们可以使用下述代码，这里主要靠的是imager包中的RGBtoHSV()函数，这个函数可以将图像从RGB空间转化为HSV空间。在RGB空间里，每个颜色的度量是看它RGB分量的数值，这是一种对机器很友好的理解方式，但是并不是对人类很适用。想象一下情侣出游，女孩说哇天好蓝啊，男孩回复说是啊，似乎rgb空间分量是(0,191,255)（这个分量指代的颜色是deep sky blue），然后女生表示需要换一个说人话的男朋友。正确的回答方式，应该是“天确实挺蓝，但好像比平时见到的天空颜色更深一些，但还是你最美blablabla”。这其中所蕴含的描述颜色的方式就是HSV空间所做的事情。在这个空间里，每个颜色都可以由色调（H）、饱和度（S）以及明亮程度（V）来表示（类似的空间还有HIS、HSL等等，具体细节还需要大家自己学习）。

回到正题，我们提取其中的S（饱和度，Saturation）空间就能得到图像的饱和度矩阵（对图像的每一个点，都有一个饱和度的度量），我们将这个矩阵转化为向量并取其99.9%分位数和平均值作为我们的特征。

    GetS = function(Img) {
      # 函数说明：计算图像饱和度的函数，输入的是图像矩阵，输出图像饱和度的分位数及均值

      # 把图像由RGB转成HSV空间
      hsvImg = RGBtoHSV(Img);
      # 返回HSV空间的饱和度S的0.001分位数以及均值作为特征
      S = hsvImg[, , 1, 2]
      result = c(quantile(S[, ], 0.999), mean(S[, ]))
      names(result) = c("quantile", "mean")
      return(result)
    }

从箱线图也可以看出，PM2.5的污染程度越高，图像饱和度越小。

    ## 污染水平vs饱和度的均值 ##
    ggplot(data, aes(x = factor(污染水平), y = Smean, fill = 污染水平, alpha = 污染水平)) +
      geom_boxplot(aes(fill = 污染水平)) +
      scale_fill_manual(values = c('chartreuse2', 'gold1', 'orange1', 'firebrick3')) + labs(x = "PM2.5", y = "饱和度的均值") +
      scale_alpha_manual(values = c(0.5, 0.5, 0.5, 0.5)) +
      scale_x_discrete(limits = c("优", "良", "轻度污染", "较严重污染")) +
      theme(axis.title.x = element_text(size = 14), axis.title.y = element_text(size = 14), legend.text = element_text(size = 12), legend.title = element_text(size = 14), axis.text = element_text(size = 12), text = element_text())+
        theme_grey(base_family = "STKaiti")

## 5. 高频含量

第五个解释性变量是高频含量。如果我们把图片看成一个二维信号，将这个信号在频域上进行分解，可以得到图像频率分布。

这种频率分布其实反映了图像的像素灰度在空间中变化的情况。例如，一面墙壁的图像，由于灰度值分布平坦，其低频成分就较强，而高频成分较弱。换句话说，高频成分决定了图像细节部分。因此，可以很直观地想到，天气晴朗时图像细节非常多，因此高频成分较多；而污染严重时，大多数图像细节被雾霾掩盖，高频成分会大幅减少。因此我们此处定义高频含量为频率分布的99.8分位数减去频率分布的0.2分位数。

在上文所述的例图中，当PM2.5=7时，天气晴朗，细节丰富，可以分辨远山、建筑、天空，此时高频含量为821.20；而PM2.5=436时，雾霾很严重，我们只能看清图像中物体的轮廓，难以分辨更多的细节，高频含量为536.97。对于上述所展示的两张图片，他们的频率分布展示如下：

    grayI1 = grayscale(Image1, drop = TRUE);
    f1 = fft(grayI1);
    F01 = (abs(f1))
    df = data.frame(F0 = as.vector(F01))
    ggplot(df, aes(x = F0)) +
      scale_x_log10(limits = c(0.1, 1000), breaks = c(as.vector(quantile(F01, 0.002)), as.vector(quantile(F01, 0.998))), labels = c("0.98\n0.1%分位数", "822.2\n99.9%分位数")) +
      geom_histogram(binwidth = 0.1, color = "steelblue", fill = "steelblue") +
      labs(x = "晴朗天气频率分布", y = "频数") +
      theme(text = element_text())+
      theme_grey(base_family = "STKaiti")

    grayI2 = grayscale(Image2, drop = TRUE);
    f2 = fft(grayI2);
    F02 = (abs(f2))
    df = data.frame(F0 = as.vector(F02))
    ggplot(df,aes(x = F0)) +
      scale_x_log10(limits = c(0.1, 1000), breaks = c(as.vector(quantile(F02, 0.002)), as.vector(quantile(F02, 0.998))), labels = c("0.51\n0.2%分位数", "537.5\n99.8%分位数")) +
      geom_histogram(binwidth = 0.1, color = "steelblue", fill = "steelblue") +
      labs(x = "雾霾天气频率分布", y = "频数") +
      theme(text = element_text())+
      theme_grey(base_family = "STKaiti")

整体上看，基本满足污染程度越高，则高频含量越小。

    GetFourier = function(Img) {
      # 函数说明：计算图像高频变量的函数，输入的是图像矩阵，输出图像频率分布的分位数之差

      # 傅里叶变换
      grayI = grayscale(Img, drop = TRUE);
      f = fft(grayI);
      F0 = (abs(f));
      # 返回频率分布的的99.8分位数和0.2分位数之差
      result = (quantile(F0[, ], 0.998) - quantile(F0[, ], 0.002))
      return(result)
    }

从箱线图的结果来看，基本满足污染程度越高，则高频含量越小的特点，证明我们这个变量具有解释效力。

    ## 污染水平vs高频含量 ##
    ggplot(data, aes(x = factor(污染水平), y = Fourier, fill = 污染水平, alpha = 污染水平)) +
      geom_boxplot(aes(fill = 污染水平)) +
      scale_fill_manual(values = c('chartreuse2', 'gold1', 'orange1', 'firebrick3')) + labs(x = "PM2.5", y = "高频含量") +
      scale_alpha_manual(values = c(0.5, 0.5, 0.5, 0.5)) +
      scale_x_discrete(limits = c("优", "良", "轻度污染", "较严重污染")) +
      theme(axis.title.x = element_text(size = 14), axis.title.y = element_text(size = 14), legend.text = element_text(size = 12), legend.title = element_text(size = 14), axis.text = element_text(size = 12), text = element_text())+
        theme_grey(base_family = "STKaiti")

# 批量处理图像

在上文提到的5种特征的基础上，我们将对图像逐个分析，提取变量，这里会用到我们之前介绍过的对文件夹中图像做批量处理的知识。

    # 读取图像数据名称
    fileName = dir(path = "FigureData")
    dataList = vector("list", 1)

    # 每个图像数据依次读入并且计算出图片对应的自变量
    for (i in 1:length(fileName)) {
      # 请确保您的图像都是Jpg格式，否则请换用相应的图像读取函数
      a = Sys.time()
      Image = load.image(paste("FigureData", fileName[i], sep = '/'))  # 路径名＋图片名
      v1 = diffvar(Image)
      v2 = clearRate(Image)
      v3 = entro(Image)
      v4 = GetS(Image)
      v5 = GetFourier(Image)
      dataList[[i]] = list(fileName = fileName[i], diffvarB = v1, clearRate = v2,
                                   entropy = v3, Squan = v4[1], Smean = v4[2], Fourier = v5)
      rm(Image)
      b = Sys.time()
      print(b - a)  # 由于一次程序运行时间过长，我们输出每次循环的运行时间，既是对程序的检查，也可作为程序保持运行的标志
    }
    ##  Time difference of 15.92419 secs
    ##  Time difference of 16.54898 secs
    ##  Time difference of 16.42315 secs
    ##  Time difference of 16.46133 secs
    ##  Time difference of 15.73884 secs
    ##  Time difference of 15.48412 secs
    ##  Time difference of 15.02853 secs
    ##  Time difference of 15.26224 secs
    ##  Time difference of 15.19045 secs
    ##  Time difference of 15.05841 secs
    ##  Time difference of 14.95099 secs
    ##  Time difference of 16.27709 secs
    ##  Time difference of 17.70153 secs
    ##  Time difference of 16.69625 secs
    ##  Time difference of 16.44896 secs
    ##  Time difference of 17.25257 secs
    ##  Time difference of 16.84369 secs
    ##  Time difference of 16.85229 secs
    ##  Time difference of 16.84731 secs
    ##  Time difference of 15.28556 secs
    ##  Time difference of 15.13625 secs
    ##  Time difference of 14.78328 secs
    ##  Time difference of 14.806 secs
    ##  Time difference of 14.61024 secs
    ##  Time difference of 18.95375 secs
    ##  Time difference of 17.74332 secs
    ##  Time difference of 16.98731 secs
    ##  Time difference of 17.33427 secs
    ##  Time difference of 16.49257 secs
    ##  Time difference of 20.47617 secs
    ##  Time difference of 15.9348 secs
    ##  Time difference of 17.75488 secs
    ##  Time difference of 22.76869 secs
    ##  Time difference of 18.78608 secs
    ##  Time difference of 16.70432 secs
    ##  Time difference of 15.96799 secs
    ##  Time difference of 16.49764 secs
    ##  Time difference of 16.43725 secs
    ##  Time difference of 15.97202 secs
    ##  Time difference of 15.69032 secs
    ##  Time difference of 15.94555 secs
    ##  Time difference of 15.83344 secs
    ##  Time difference of 17.76137 secs
    ##  Time difference of 17.37338 secs
    ##  Time difference of 16.99412 secs
    ##  Time difference of 16.1108 secs
    ##  Time difference of 16.08341 secs
    ##  Time difference of 16.64453 secs
    ##  Time difference of 15.68177 secs
    ##  Time difference of 17.12513 secs
    ##  Time difference of 15.72473 secs
    ##  Time difference of 15.70522 secs
    ##  Time difference of 16.34324 secs
    ##  Time difference of 15.65355 secs
    ##  Time difference of 17.54957 secs
    ##  Time difference of 18.71052 secs
    ##  Time difference of 16.97534 secs
    ##  Time difference of 16.42148 secs
    ##  Time difference of 16.43705 secs
    ##  Time difference of 16.53729 secs
    ##  Time difference of 15.8265 secs
    ##  Time difference of 16.58908 secs
    ##  Time difference of 15.44073 secs
    ##  Time difference of 15.76447 secs
    ##  Time difference of 15.69024 secs
    ##  Time difference of 15.14942 secs
    ##  Time difference of 15.598 secs
    ##  Time difference of 18.12931 secs
    ##  Time difference of 15.74765 secs
    ##  Time difference of 15.32713 secs
    ##  Time difference of 15.53763 secs
    ##  Time difference of 15.43392 secs
    ##  Time difference of 15.88 secs
    ##  Time difference of 15.50641 secs
    ##  Time difference of 15.01433 secs
    ##  Time difference of 15.3125 secs
    ##  Time difference of 15.4079 secs
    ##  Time difference of 15.4017 secs
    ##  Time difference of 15.34148 secs
    ##  Time difference of 14.94748 secs
    ##  Time difference of 15.3656 secs
    ##  Time difference of 15.36706 secs
    ##  Time difference of 15.32753 secs
    ##  Time difference of 15.80034 secs
    ##  Time difference of 15.99865 secs
    ##  Time difference of 15.17265 secs
    ##  Time difference of 15.04236 secs
    ##  Time difference of 17.07283 secs
    ##  Time difference of 15.61664 secs
    ##  Time difference of 16.85171 secs
    ##  Time difference of 16.21294 secs
    ##  Time difference of 16.77883 secs
    ##  Time difference of 17.76379 secs
    ##  Time difference of 15.55624 secs
    ##  Time difference of 15.1698 secs
    ##  Time difference of 15.11495 secs
    ##  Time difference of 14.87736 secs
    ##  Time difference of 15.04062 secs
    ##  Time difference of 15.0941 secs
    ##  Time difference of 15.08058 secs
    ##  Time difference of 17.19763 secs
    ##  Time difference of 14.96451 secs
    ##  Time difference of 15.42124 secs
    ##  Time difference of 15.31433 secs
    ##  Time difference of 15.17624 secs
    ##  Time difference of 15.2079 secs
    ##  Time difference of 15.17041 secs
    ##  Time difference of 15.04048 secs
    ##  Time difference of 15.17527 secs
    ##  Time difference of 15.30013 secs
    ##  Time difference of 15.42181 secs
    ##  Time difference of 14.97502 secs
    ##  Time difference of 15.04043 secs
    ##  Time difference of 15.17839 secs
    ##  Time difference of 15.12678 secs
    ##  Time difference of 15.12052 secs
    ##  Time difference of 15.41905 secs
    ##  Time difference of 15.49695 secs
    ##  Time difference of 16.11805 secs
    ##  Time difference of 19.10498 secs
    ##  Time difference of 18.42919 secs
    ##  Time difference of 17.56066 secs
    ##  Time difference of 17.27953 secs
    ##  Time difference of 17.06164 secs
    ##  Time difference of 17.23434 secs
    ##  Time difference of 16.9896 secs
    ##  Time difference of 16.58675 secs
    ##  Time difference of 17.58042 secs
    ##  Time difference of 16.71057 secs
    ##  Time difference of 16.99359 secs
    ##  Time difference of 16.45428 secs
    ##  Time difference of 16.66697 secs
    ##  Time difference of 16.46655 secs
    ##  Time difference of 16.61738 secs
    ##  Time difference of 15.94547 secs
    ##  Time difference of 15.6587 secs
    ##  Time difference of 16.16573 secs
    ##  Time difference of 15.91045 secs
    ##  Time difference of 15.5864 secs
    ##  Time difference of 14.92693 secs
    ##  Time difference of 16.48011 secs
    ##  Time difference of 16.22959 secs
    ##  Time difference of 15.7154 secs
    ##  Time difference of 16.38313 secs
    ##  Time difference of 15.6085 secs
    ##  Time difference of 17.11937 secs
    ##  Time difference of 15.12234 secs
    ##  Time difference of 15.30148 secs
    ##  Time difference of 16.78499 secs
    ##  Time difference of 15.30657 secs
    ##  Time difference of 15.80185 secs
    ##  Time difference of 16.16764 secs
    ##  Time difference of 17.49454 secs
    ##  Time difference of 15.61177 secs
    ##  Time difference of 16.37329 secs
    ##  Time difference of 14.7184 secs
    ##  Time difference of 14.94318 secs
    ##  Time difference of 15.37248 secs
    ##  Time difference of 15.30247 secs
    ##  Time difference of 15.26182 secs
    ##  Time difference of 16.10957 secs
    ##  Time difference of 16.25879 secs
    ##  Time difference of 15.95095 secs
    ##  Time difference of 14.86068 secs
    ##  Time difference of 15.03337 secs
    ##  Time difference of 16.30319 secs
    ##  Time difference of 15.75683 secs
    ##  Time difference of 15.95851 secs
    ##  Time difference of 14.93456 secs
    ##  Time difference of 14.77218 secs
    ##  Time difference of 14.93489 secs
    ##  Time difference of 14.99591 secs
    ##  Time difference of 14.96275 secs
    ##  Time difference of 16.18051 secs
    ##  Time difference of 15.71754 secs
    ##  Time difference of 15.09511 secs
    ##  Time difference of 15.00476 secs
    ##  Time difference of 14.99111 secs
    ##  Time difference of 15.00643 secs
    ##  Time difference of 16.99831 secs
    ##  Time difference of 16.65316 secs
    ##  Time difference of 15.97225 secs
    ##  Time difference of 15.82448 secs
    ##  Time difference of 16.76548 secs
    ##  Time difference of 17.06682 secs
    ##  Time difference of 17.60652 secs
    ##  Time difference of 16.10572 secs
    ##  Time difference of 15.09421 secs
    ##  Time difference of 15.03156 secs
    ##  Time difference of 15.47989 secs
    ##  Time difference of 16.26912 secs
    ##  Time difference of 15.23154 secs
    ##  Time difference of 15.06071 secs
    ##  Time difference of 15.21822 secs
    ##  Time difference of 15.45617 secs
    ##  Time difference of 15.39119 secs
    ##  Time difference of 15.20212 secs
    ##  Time difference of 15.10998 secs
    ##  Time difference of 15.49305 secs
    ##  Time difference of 15.41319 secs
    ##  Time difference of 15.56492 secs
    ##  Time difference of 15.16238 secs
    ##  Time difference of 15.17531 secs
    ##  Time difference of 15.18835 secs
    ##  Time difference of 15.42153 secs
    ##  Time difference of 15.41463 secs

    # 保存
    save(dataList, file="....RData")
    dataOutput = data.frame(matrix(unlist(dataList), nrow = length(fileName), byrow = T))
    colnames(dataOutput) = names(dataList[[1]])

    # 输出csv文件
    # 请将对应时间的PM2.5值补充到表中，并将该列命名为PM25，再用第二部分的代码来做进一步的回归分析
    write.csv(dataOutput, file = 'Output.csv')