我发论文啦系列之一《融合FY-3C号和FY-4A号卫星数据的积雪面积变化研究—以祁连山区为例》

0.前言

开设这个系列的原因呢，主要是为了学术交流、思想感悟（内心OS：为了记录自己当时有多菜），反正我总有一种当时很牛逼，后来菜到自己都不想看的感觉。这或许就是成长吧。看完以后别取关就行，哈哈哈~

言归正传，这篇文章是我大三（2018年11月）开始做的，但到现在才出刊，可以说是经历坎坷，路途漫长。算是我第一次独立完成的吧，当然也离不开老师的帮助，在此表示感谢！现在想一想每一步都是心酸啊。这篇博文就记录一下所有的过程、数据、代码以及思考和感悟吧！

1.摘要

大家先看看这个内容是不是和自己有关系，再往下看吧。当然对内容不感兴趣的可以去看看我的感悟（不要脸的请求）。
祁连山区积雪类型丰富、判识复杂，是中国积雪研究的典型区域。因此，精确地监测祁连山区积雪面积变化及其时空演变，对祁连山区生态环境和社会经济发展等具有重要意义。FY-3CMULSS 利用多阈值积雪指数模型提供全球日积雪覆盖产品，FY-4A AGRI 传感器每15~60 min提供一景覆盖全球的多光谱影像。基于FY-4A AGRI 高时间分辨率的特征，构建适合于FY-4A号数据的动态多阈值多时相云隙间积雪识别方法，很大程度上减小了云对光学数据识别积雪造成的影响，并结合FY-3C MULSS 积雪覆盖日产品较高空间分辨率的优势，融合得到去除云后的FY3C4 积雪覆盖数据。利用Landsat 8 OLI 卫星数据对融合后的积雪数据进行对比验证，结果表明融合FY-3C 和FY-4A 后的数据能更好地判识祁连山区的积雪覆盖情况。以MODISMOD10A2 积雪产品为真实值，随机检验了2018 年3 月~2019 年3 月融合后数据的积雪判识精度，发现无云情况下方法的总体精度可达到85.25%。进一步研究发现祁连山区积雪面积在海拔、气候和坡向等因素的影响下时空分布极不均匀，总体呈现出冬春季节大于夏秋季节，以及东部积雪面积大于西部积雪面积的特征。

2.研究的目的和意义

随着中国风云系列气象卫星的蓬勃发展及数据量的剧增［16］，有必要结合FY-3C MULSS（多仪器融合数据）积雪覆盖日产品较高空间分辨率的特点和FY-4A AGRI（多通道扫描成像辐射计）数据高时间分辨率的特点，对祁连山区积雪面积变化和时空演变进行研究。一方面有利于发挥我国风云系列气象卫星数据的优势，另一方面弥补了目前大多数研究中仅用单一数据研究积雪面积变化的不足。
啊哈哈哈哈~~最后一句说的我好害羞（正经脸）！

3.研究区

就是它祁连山区，我又想起了小学写的作文，我的家乡坐落在祁连山脚下，这里风景优美~~~不过风景确实挺美的，下面再给几张家乡的风景图。

4.研究数据

主要的风云系列呢，都是国产的数据，资料什么的在官网上都很多我就不详细说明了，怎么说呢，虽然我们国家的风云4卫星已经进入了世界顶尖行列了，但是相对比其他国家的新一代静止轨道卫星依然存在差距，主要表现在时间分辨率、数据预处理流程化等方面，时间分辨率最值得一提的就是其他国家基本都看到稳定的每十分钟一景的数据，而我们国家15-60min不等（宣传说是10分钟，但为什么会这样，原因不明）。

5.数据预处理

那个时候读书和论文少，没文化，仅仅在预处理问题上耗费了我大量的时间经历，还请了很多朋友帮我处理了一部分（自己对号入座一下），在此也表示一一谢过（爱你们哦），想想真的都是好人啊，浪费自己的时间来帮我。现在看来，还不是因为自己菜，哎~~~
FY-4A数据的预处理主要也就是个几何校正和反射率数据的获取，那篇博客在我去年发现新方法后就发在CSDN了，以免大家再浪费时间，经过我的验证那个方法并没有问题，而且不需要编程，要说批量处理的话你可以利用IDL调用那个ENVI的工具就好了。反射率很容易都给出了定标系数计算一下就完了。
（FY-4A数据几何校正处理方法博文链接：静止轨道卫星几何校正——FY-4A几何校正（风云系列静止卫星可以参考））
FY-3C数据预处理也浪费了我的很多时间，当时就什么都不懂啊~各种问各种查，最终找到了一篇讲有关于这种数据的预处理方法才得以解决，简单来说，按照HDF中给出的信息（坐标系，左上角右下角坐标），定义头文件就可以了，当时就很心塞纯手工自己一个又一个弄的（应该有办法批量，呜呜呜）
这里就不得不说一说MODIS和Landsat的好处了，用的人多很容易查，不过GEE确实是现在处理他们的一大利器，当然你不会也无妨有很多“单机版”方法处理。我也写过相关的博客（单机版（不编程）：MODIS数据MRT处理集锦（方法+报错+原因））（GEE：GEE学习记录（三）：Google Earth Engine处理MODIS产品——以MOD10A1积雪产品为例）

6.方法、代码和数据

1.方法：总的来说就是论文里面的这样。

简单来说：识别FY-4A积雪，融合两种结果，获得最终结果。

2.FY-4A积雪识别
因为选择了2KM的数据，波段数目就比较少，因此看了一些文献就打算用一些简单的方法识别一下。这个应该是我在全网找到最最最最简单的一种积雪识别的方法了，毕竟受到波段限制和研究区域的限制。

选择样本确定阈值，由于之前才疏学浅没有看到半经验模型的角度效应校正方法，因此用同一种阈值总是出现奇奇怪怪的现象，所以我就只能从成像时间段下手，分阈值。后续大家要是想弄，我觉得先做一下角度效应校正应该也就不需要阈值分段了（那时候我还觉得自己很牛批，搞了一种新方法，还起了个响亮的名字。。。）又是被自己打败系列。

代码1：代码是我现在才写的，之前太low就在ENVI里面用决策树一个一个搞的，正好那时候实习好像空闲时间比较多，保研面试也都结束了，就一个一个慢慢弄了。。。。。。。。这时候又开始想到了“工欲善其事必先利其器”，哎。。。代码里面的数据改一改路径和if处的文件名匹配一直就成了（按照时间一一对应），看不懂的话先去学python就好了。考虑到后面的融合只需要用到FY-4A的积雪像元直接重分类为两类雪（1）和非雪（2），哦对了也记得微调阈值

#批量影像分类和重分类
#2021/1/29
#Jonathan

import gdal
import os
import numpy as np
import glob

list_tif = glob.glob(r"f:\JGT\*.dat")
list_tif_NDSI = glob.glob(r"f:\JGT\*_NDSI.tif")
out_path = 'f:/JGT/'

'''for tif in list_tif_NDSI:
    print(tif[7:19])'''

for tif in list_tif:
    for tif_NDSI in list_tif_NDSI:
        if tif[51:63] == tif_NDSI[7:19]:
            print(tif[51:63])
            print(tif_NDSI[7:19])
            open_tif = gdal.Open(tif)
            open_tif_NDSI = gdal.Open(tif_NDSI)
            (filepath, fullname) = os.path.split(tif)
            (prename, suffix) = os.path.splitext(fullname)
            obtain_tif_NDSI = open_tif_NDSI.GetRasterBand(1).ReadAsArray().astype(float)
            obtain_tif4 = open_tif.GetRasterBand(4).ReadAsArray().astype(float)
            obtain_tif5 = open_tif.GetRasterBand(5).ReadAsArray().astype(float)
            result_tif = np.zeros(obtain_tif4.shape, dtype=np.float)
            Height = obtain_tif4.shape[0]
            Width = obtain_tif4.shape[1]
            #考虑到后面的融合只需要用到FY-4A的积雪像元直接重分类为两类雪（1）和非（2）
            for i in range(Height):
                for j in range(Width):
                    if obtain_tif4[i][j] > 0.078:
                        result_tif[i][j] = 0
                    elif obtain_tif_NDSI[i][j] > 0.18 and obtain_tif5[i][j] > 0.1:
                        result_tif[i][j] = 1
                    else:
                        result_tif[i][j] = 0
            driver = gdal.GetDriverByName("GTiff")
            datatype = gdal.GDT_Float32
            dataset = driver.Create(out_path + tif[51:63] + '_result.tif', Width, Height, 1, datatype)
            if dataset != None:
                dataset.SetGeoTransform(open_tif.GetGeoTransform())  # 写入仿射变换参数
                dataset.SetProjection(open_tif.GetProjectionRef())
            dataset.GetRasterBand(1).WriteArray(result_tif)
            dataset.GetRasterBand(1).SetNoDataValue(0)
            del dataset
            print(tif[51:63] + '完成')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

给出一个示例的结果2019年1月16日数据：

3.FY-4A结果融合
注意：可以参加试验的数据数量根据实际情况调整（可以比10-12景少一点），一定要保证云比较少的数据参加运算，数量可以不多但是质量要高。

简单来说，先把数据分类，一种就是有必要融合的（晴间云数据），另一种就是没有必要融合的（云数据），多说无益，第一幅是晴间云第二幅是云数据2018.03.12（我懒，没做2019.1.16，意思就这个意思大家康康就行）。之后呢，对于晴间云数据进行融合，有5景（这是一个经验值根据实际情况可以进行调整，不局限在2-3景，2-3景只是个平均数，比如在冬季云雪误判的比例比较大就需要增多景数来控制误差例如此时需要5景以上来确定）以上的FY-4A结果判为积雪的像元就认为当日是积雪，最后得到结果。这里没有代码的原因就是2-3景是经验值需要人为再去验证例如这里还是用2-3景的话就会出现大量误差，这也是一个很大的弊端。

4.FY-3C和FY-4A融合

融合方法如下图所示，简而言之呢，就是把FY-3C识别的和FY-4A识别的积雪一致的区域保留下来，其他的积雪像元都用FY-4A的，其实相当于没怎么用FY-3C的积雪识别结果。

废话不多说上代码：
在此之前请注意原始FY-4A的积雪是1，FY-3C有自己的代表值意义，去看数据说明就有了，这里在用的时候并没有都写到，因为有些类型的地物在祁连山区确实也不存在。融合前先给FY-4A重采样一下，不然不可能会失败。

#批量融合影像
#2019/09/09
#Jonathan

import gdal
import os
import numpy as np
import glob

list_tif_FY_4 = glob.glob(r"f:\JGT\*._FY-4A")
list_tif_FY_3 = glob.glob(r"f:\JGT\*_FY-3A.tif")

out_path = 'f:/JGT/'

'''for tif in list_tif_NDSI:
    print(tif[7:19])'''

for tif4 in list_tif_FY_4:
    for tif3 in list_tif_FY_3:
        if tif4[51:63] == tif3[7:19]:
            print(tif4[51:63])
            print(tif3[7:19])
            open_tif4 = gdal.Open(tif4)
            open_tif3 = gdal.Open(tif3)
            (filepath, fullname) = os.path.split(tif4)
            (prename, suffix) = os.path.splitext(fullname)
            obtain_tif4 = open_tif4.GetRasterBand(1).ReadAsArray().astype(float)
            obtain_tif3 = open_tif3.GetRasterBand(1).ReadAsArray().astype(float)
            result_tif = np.zeros(obtain_tif4.shape, dtype=np.float)
            Height = obtain_tif4.shape[0]
            Width = obtain_tif4.shape[1]
            for i in range(Height):
                for j in range(Width):
                    if obtain_tif4[i][j] == 1:
                        result_tif[i][j] = obtain_tif4[i][j]
                    elif obtain_tif3[i][j] == 0 or obtain_tif3[i][j] == 200:
                        result_tif[i][j] = obtain_tif4[i][j]
                    else:
                        result_tif[i][j] = obtain_tif3[i][j]
            driver = gdal.GetDriverByName("GTiff")
            datatype = gdal.GDT_Float32
            dataset = driver.Create(out_path + tif4[51:63] + '_result.tif', Width, Height, 1, datatype)
            if dataset != None:
                dataset.SetGeoTransform(open_tif4.GetGeoTransform())  # 写入仿射变换参数
                dataset.SetProjection(open_tif4.GetProjectionRef())
            dataset.GetRasterBand(1).WriteArray(result_tif)
            dataset.GetRasterBand(1).SetNoDataValue(0)
            del dataset
            print(tif4[51:63] + '完成')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

融合结果：提取一下像元数量，面积就很容易得到了。白色是雪(1)，红色是云(50)，黄色是无雪，浅蓝色是水。

提取面积：属性表看积雪像元数量40544个，但是像元单位为0.01°，经过换算面积为40381平方公里。

7.验证

MOD10A2验证，由于MOD10A2分辨率是8天算上可以用的FY3C4数据一共能用于验证的数据就更少了，所以做一下验证大致总体精度OA能在85%，当然由于无雪地表的影响精度会偏高。和Landsat8做了对比发现结果一些结果有所改善。

8.关于这个研究的思考（有兴趣继续的可以看看）

我先来谈谈这项研究的优势吧：思路比较新颖：市面上大多数研究90%以上吧，都是先把静止轨道卫星数据进行融合，然后做的积雪识别工作，而这个研究正好反过来做，去融合结果；
再来谈谈以后可以改进的地方（ps:完全是亲情奉献，因为我以后不做这个方面的研究了，我觉得还有很多可以挖掘，我一并给出来大家，希望有兴趣的朋友能继续做做）：
（1）数据的选择：这项研究用的是2KM的数据，波段太少。但是没办法研究区域小，再用4KM的数据意义就不大了。因此，以后可以扩增研究区选择4KM的数据，4KM的数据波段很多，识别更加准确。我现在完成的FY-4A全圆盘数据的积雪识别工作中（如果顺利现在投出去，今年应该可以再与大家见面），发现4KM的数据真的就很好用！！！
（2）研究区：现在看来，当初走入了一个歧途，限制在祁连山区做这项工作，实际上大可以对全圆盘先做识别再去裁剪出研究区域，就会减少很多问题。
（3）角度效应校正：前面说过分段阈值的事情，但是通过角度效应校正完全可以把这个问题给避免了，Roujean的模型效果就很好，我已经在新的研究中使用了，大家也可以尝试，文献链接：Roujean J L , Leroy M , Deschamps P Y . A bidirectional reflectance model of the Earth"s surface for the correction of remote sensing data[J]. Journal of Geophysical Research, 1992, 97(D18):20455.。
（4）研究的自动化程度不高，人为干扰比较大：由于当时真的是个编程能力很Low的本科生所以很多都是靠大量劳动完成的，也就是俺的同学们。其次呢，算法里面那个对于云数据和晴间云数据的区别，全靠我自己的意识判别的，以及融合数据时所需要用到的判断条件都需要人为查看再决定，在这里应该可以有所改进。
（5）验证数据选择：这里按理说应该用MOD10A1数据验证会好一点，但是那时候想到MOD10A2也是个合成产品，到最后发现算上那个可以用的晴间云数据，也没有多少能于验证的数据，所以以后还是用MOD10A1或者站点数据吧。
大概就是这些内容吧，希望后来的小伙伴能有进一步好的结果。

9.通过这段时间的感悟

其实这个工作开始到结束一共11个月吧，但是出刊用了一年半（等到都没有耐心了），保研结束后我就进入现在的课题组了，没有再花费什么经历了，进入新的课题组后我也做了一个相关的研究，目前基本结束要投稿了，然后就不再涉及相关内容的研究了。有这么几点心得和建议吧，与各位小伙伴分享一下：
（1）学好编程：之前感觉自己还行，但是实战能力真的不足，到现在不能说编程很强，最起码的实现还是没啥大的问题了~工欲善其事必先利其器！！！
（2）意外收获就是读文献但不要被文献禁锢思想：这次很意外的就是由于没有读到很多文献，结果还有了新思路“先识别，再融合（融合结果）”，如果看到文献都是先融合再识别，我一定不会想到新思路。
（3）多与人交流：无论是与导师或者是做相同方向的同学，一定要多交流，多讲会发现很多问题，得到很多启发。
（4）要接受不怎么好的实验结果，新思路和好结果可能不容易兼得，纵观各种研究都是从一个很弱鸡的东西，慢慢变来的，不可能一弄就很好很优秀。就像我这个，我承认精度确实不如那种很成熟的方法，但是难保以后不会变得成熟，对吧。
（5）多想办法解决问题：编程不行有编程不行的办法，只要做就能做完（无非你可能会耗费很多时间），但是通过这次学会了，下次不就更容易了吗？

10.后记

我打算一直把这个系列做下去给出一些思考，我主要是希望大家能一起讨论，我也不怕大家看到什么新的思路抢在我前面弄出个啥（毕竟我说的你可能早就想到了，哈哈哈~），但是请讲武德哦，尊重版权。

论文原文链接：（知网或者期刊网站均可）
融合FY-3C号和FY-4A号卫星数据的积雪面积变化研究——以祁连山区为例
融合FY-3C号和FY-4A号卫星数据的积雪面积变化研究——以祁连山区为例

附录

数据的话我给出了博客里面用的2019年1月16日的数据，包括原始的和识别的结果。
百度网盘链接：（如果失效的话再私信问我要就好了）
链接：https://pan.baidu.com/s/1Z-aNaCQ4PHRXREgD_lqIrQ
提取码：x5kx