基于深度学习的时空间数据融合

hellogis
字数统计: 2.8k   |   阅读时长≈ 11 分钟

基于深度学习的时空间数据融合

以 Landsat-Sentinel NDVI 数据为例，测试利用深度学习进行时空间数据融合的可行性

1、待解决的问题

  1. 融合方法

    深度学习 VS 传统方法

    哪种方法更好

  2. 融合策略

    Blending-then-Index（BI） VS Index-then-Blending（IB）

    先融合再计算 NDVI 还是先计算 NDVI 再融合

  3. 重要性

    融合方法 VS 融合策略

    哪种更重要

2、Pytorch 深度学习环境配置

CUDA + CUDNN + Anaconda + Pytorch 安装

CUDA, CUDNN和pytorch 版本需要一一对应，否则无法运行

这里统一选择 10.1 的版本

2.1、CUDA 安装

CUDA(ComputeUnified Device Architecture)，是显卡厂商NVIDIA推出的运算平台。 CUDA是一种由NVIDIA推出的通用并行计算架构，该架构使GPU能够解决复杂的计算问题。

  1. 下载地址

    https://developer.nvidia.com/cuda-toolkit-archive

    CUDA 版本

  2. 我这里选择 CUDA Toolkit 10.1 update2 Archive 版本

    CUDA 安装选择

    从官网上下载比较困难，需要挂梯子，最好能找到国内镜像。

  3. 安装

    • 双击下载好的 cuda_10.1.243_426.00_win10.exe 执行安装程序
    • 选择自定义安装，也可以选精简安装，我是为了能自定义安装路径（记住自定义的路径），我的驱动程序组件的勾选如下：

    CUDA 自定义安装选项

    • 其它默认

2.2、CUDNN 下载

  1. 下载地址

    https://developer.nvidia.com/rdp/cudnn-archive

  2. 需要注册、登录，我选择了 cuDNN v7.6.5 这版

    cuDNN 下载

  3. 解压，将其中的 bin、include、lib 文件夹复制到 CUDA 安装目录下，直接复制就可以，不会有覆盖的内容

    cuDNN 复制到 CUDA

2.3、安装 Anaconda

  1. 下载

    https://mirrors.tuna.tsinghua.edu.cn/anaconda/archive

    可选择最新的

  2. 安装

    可 for all user 安装，修改安装路径 无法勾选配置系统环境变量的话，后续可自行配置

  3. 验证安装

    1
    conda --version

2.4、配置国内镜像

1
2
3
4
5
6
7
8
9
10
11
# 添加清华anaconda镜像，为 http 的方式，https 方式可能会出现验证问题：
conda config --add channels http://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels http://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --add channels http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --set show_channel_urls yes
# Conda 附加库:
conda config --add channels http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/pytorch/

conda config --show-sources #查看当前使用源
conda config --remove channels 源名称或链接 #删除指定源
conda config --add channels 源名称或链接 #添加指定源

2.5、构建虚拟环境

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 打开 Anaconda Prompt 这一步很重要，不要在命令行操作，在命令行操作可以是 conda init cmd.exe
# 创建虚拟环境
conda create -n dl python=3.8
# 也可以指定创建路径的命令，但是这种方法没有环境名称，只能用路径来代替环境名称
conda create -p C:\\Users\\Admin\\Desktop\\deeplearing-work\\dl python=3.8
# 查看Python环境
conda info --env
# 激活环境
activate dl # 或者
activate C:\\Users\\Admin\\Desktop\\deeplearing-work\\dl
# 命令行中
conda activate dl # 或者
conda activate C:\\Users\\Admin\\Desktop\\deeplearing-work\\dl
# 退出环境 
deactivate
# 删除环境
conda remove -n python35 --all # 或者
conda remove -p C:\\Users\\Admin\\Desktop\\deeplearing-work\\dl --all

# 安装 pytorch CPU 版
pip3 install torch torchvision torchaudio
pip install scipy pandas scikit-image scikit-learn termcolor

2.6、安装 GDAL 环境

  1. 下载安装包

    由于老师 ppt 上的 GDAL 安装包下载路径失效，从别的网站上找到了安装包，但是版本略有差异

    https://www.lfd.uci.edu/~gohlke/pythonlibs/

    我下载的版本为 GDAL-3.4.3-cp38-cp38-win_amd64.whl 最好保存在英文路径下，我存放在 C:\Users\Admin\Desktop\deeplearing-work 路径下

  2. 安装 GDAL 软件包

    1
    pip install C:\\Users\\Admin\\Desktop\\deeplearing-work\\GDAL-3.4.3-cp38-cp38-win_amd64.whl

    GDAL 安装成功

2.7、Pytorch 安装

1
conda install pytorch torchvision cudatoolkit=10.1

3、实验内容与步骤

3.1、实验内容

  1. 数据准备

    选取农田、森林作为典型研究区，各研究区准备三个时相的 Landsat-Sentinel 影像对

    Landsat-Sentinel 影像对

  2. 数据预处理

    • 部分像元缺失，需要插值
    • Landsat 与 Sentinel 影像存在较大几何偏差，需要人工几何校正

  3. 融合流程

    • BI：先融合NIR-RGB数据，再计算NDVI
    • IB：先计算NDVI，再融合NDVI数据

  4. 精度验证

    • 决定系数：$\large R^2 = 1 - \frac{\sum\limits^n_{i=1}(y_i - p_i)^2}{\sum\limits^n_{i=1}(y_i - \bar{y}_i)^2}$
    • 均方根误差：$\large RMSE = \sqrt{\frac{\sum\limits^n_{i=1}(y_i - p_i)^2}{n}}$
    • 结构相似性：$\large SSIM = \frac{(2\bar{y}*\bar{p}+C_1)(2\sigma_{yp}+C_2)}{(\bar{y}^2+\bar{p}^2+C_1)(\sigma^2_y+\sigma^2_p+C_2)}$

3.2、实验步骤

数据准备：

  • 分别打开cropland 和forest文件夹，可以看到三个时相的Landsat-Sentinel 影像对，数据已经过异常值去除、配准等预处理。其中，第一、第三时相数据用于深度学习模型训练，第一、第二时相数据用于模型预测。
  • 影像共包含四个波段，分别为B-G-R-NIR, 反射率值范围为0~10000，在ArcGIS或Envi中使用波段计算器生成NDVI数据备用

时相数据

  1. 原始数据 NDVI 计算

计算公式：$\large NDVI = \frac{(NIR - R)}{NIR + R}$，NIR 表示近红外波段，R 表示红外波段

ArcGIS 10.8 安装教程：

https://blog.csdn.net/Augenstern_QXL/article/details/123781664

以 ArcGIS 10.8 计算 LB_144.tif 的 NDVI 为例，不保证操作的正确性，因为我已经忘完了我的专业知识了

  • 加载数据

    数据加载

  • 打开影像分析

    影像分析

  • 设置影像分析选项

    NDVI波段设置

  • 进行 NDVI 计算

    开始计算

  • 得到 NDVI 计算的结果

    NDVI计算结果

  • 设置已分类显示影像

    影像分类显示

  • 最终显示效果

    NDVI分类显示效果

    也不知道算没算对

  1. 代码文件
  • SRCNN 文件夹中包含三个文件
    • utils.py：实现了数据处理所需的函数，包含基于 GDAL 库的 geotiff 的读写、日志的输出、统计各 batch 的loss、模型参数的保存和读取、数据增强（旋转、镜像）、NDVA 计算、读取影像、归一化、预处理等
    • srcnn.py：是训练的主要代码
    • assess.py：精度验证的代码

  1. 开始训练

    1
    2
    3
    4
    # 运行深度学习代码
    python srcnn.py --train --test --mode BI --area cropland
    # 其中，mode 参数为融合策略，可选BI和IB，
    # area 参数为研究区，可选 cropland 或 forest 相对路径

    开始训练：

    开始训练

    loss 随着训练时间的增加而迅速减小

    训练完成

    在漫长的等待中，终于训练完成了顶级的电脑性能，是老师电脑训练时间的十几倍。。。

    最终 cropland 在 BI 策略下，训练的精度验证的结果为：R2=0.833378,RMSE=0.08739,SSIM=0.813840；比传统的方法是要好很多的。

    在目录下生成两个 tif 文件，一个是 refer_cropland.tif（应该是第二时相 NDVI 的计算结果），一个是 SRCNN_cropland_BI.tif（模型预测 NDVI 的结果）

    在 ArcMap 中加载如下所示：

    • refer_cropland.tif

      refer_cropland.tif

    • SRCNN_cropland_BI.tif

      SRCNN_cropland_BI.tif

      预测结果和真实结果还是有些差异的

  2. 再训练一个 IB 策略的模型

    1
    python srcnn.py --train --test --mode IB --area cropland

    最终训练的结果如下：

    IB训练完成

    生成 SRCNN_cropland_IB.tif 影像，加载到 ArcMap 中，显示如下：

    SRCNN_cropland_IB.tif

  3. 对比 BI 策略预测生成的影像、 IB 策略预测生成的影像、根据波段计算的 NDVI 影像、第二时相的真实影像，看是 BI 策略的影像效果更好还是 IB 策略的影像效果好

    我觉得是 BI 的好，我比较眼拙，好难分辨

  4. 使用 STARFM 方法，BI 模式

    ENVI 5.6 安装：

    https://zhuanlan.zhihu.com/p/564116451

    • spatiotemporalblending.sav 安装：

      http://www.chen-lab.club/?post_type=products&page_id=18952

    • 打开 ENVI + IDL

      ENVI+IDL

    • BI 模式使用 T1 时刻的 Sentinel-Landsat 影像对，和 T2 时刻的 Landsat 影像进行融合。在 IDL 命令行输入 spatiotemporalblending 打开数据融合工具

      打开数据融合工具

    • 然后，点击Import Coarse Images, 选择 T1 和 T2 时刻的 Landsat 影像；然后选中 T1 时刻的 Landsat 影像，点击 Add Fine Image，选择 T1 时刻的 Sentinel 影像

      数据融合设置

    • 几分钟后，生成ENVI格式的影像文件，打开影像观察融合结果，使用栅格计算器计算NDVI

      • 将融合的影像加载进 ENVI

    加载融合后的图像

    • 计算 NDVI

      打开 NDVI 计算工具

      计算NDVI

      输入波段开始计算

    开始计算

然后在 ENVI 或者 ArcGIS 中计算 NDVI，融合后的图像貌似波段混乱了，依据 band4 为 NIR band3 为 R 算出来不对劲，然后选了个能算的波段，发现效果挺差的。

  1. 使用 STARFM 方法，IB 模式

    IB 模式使用T1时刻的Sentinel-Landsat NDVI pairs，和T2时刻的Landsat NDVI进行融合。

    • 首先计算 T1、T2 时刻的 Landsat 图像的 NDVI 和 T1 时刻的 Sentinel 图像的 NDVI

    • 参考步骤 6 对计算出来的 NDVI 图像进行融合，点击Import Coarse Images, 选择T1和T2时刻的Landsat NDVI；然后选中T1时刻的Landsat NDVI ，点击Add Fine Image，选择T1 时刻的Sentinel影像；注意设置合适的Min Value和 Max Value。

      IB模式融合设置

    • 融合的图像结果如下：

      融合结果

    • 在 ArcGIS 中显示的效果如下：

      分类显示

      如果我的实验步骤没有错误的话，我认为在使用 STARFM 方法下，IB 模式的效果要比 BI 模式的效果好很多

  2. 剩下的 ESTARFM 方法和 FSDAF 方法应该都差不多，也就不再进行了

4、实验结果与总结

4.1、问题

  1. 深度学习 VS 传统方法的优劣？各算法生成的图像有什么特征？主要错误/误差发生在哪些区域？可能的原因？

    • 总的来说深度学习的方法是要优于传统方法的
    • 传统方法生成的图像圆润一点吧，连续一点，深度学习的话有许多散点的情况
    • 说不上来。。。

  2. BI VS IB 哪种策略更适合NDVI融合？

    • 深度学习方法中，可能是 BI 策略好
    • 传统方法中，可能是 IB 策略好

  3. 融合方法 VS 融合策略哪个对融合精度的提升哪个更大？

    融合方法吧

4.2、思考

  1. 上述结论是否适用于其它传感器/其它植被指数的时空间融合？有哪些潜在的应用方向？
  2. 尝试调整模型结构、训练参数、优化器、LOSS函数，对比融合精度。
  3. 你认为SRCNN方法存在哪些不足？如何改进这些不足？

5、总结

  1. 很久没有用过专业软件了，感觉有个两年多了，这个做到后面都没耐心了，这图都看不懂……退化严重，工作时候也就是用QGIS看看矢量文件，连一下PG库，GIS的专业软件都太重了，电脑不想装，这一次作业可把这三大件都凑齐了，QGIS最新版、ArcGIS 10.8、ENVI 5.6，啧啧啧，还好是研习室的电脑，想开了，然后又补了一个VS 2022，2022版好用很多阿，巨硬进步很大。
  2. 第一次跑遥感方面的深度学习代码，代码不能说完全看不懂，略懂，感叹很神奇，一行代码定义的网络和卷积能跑出这样的效果，牛。
  3. 这代码训练一次，在老师的电脑两百多秒，研习室电脑是他的接近二十倍，纯性能差距，这就是配的8k的电脑的性能:happy:
  4. 总算也是干了点专业稍相关的事情了，几大软件都用上了，轮番上阵，hhh。。。
  • Copyrights © 2022-2023 hqz

请我喝杯咖啡吧~

支付宝
微信