多波段遥感影像数据增广怎么做？教你用PaddleSeg处理多波段遥感任务
多波段遥感影像数据增广怎么做？教你用PaddleSeg处理多波段遥感任务
【飞桨开发者说】陈奕州，飞桨开发者技术专家（PPDE），PP SIG Model-CV 成员，重庆交通大学测绘科学与技术在读研究生。
项目背景
随着近几年AI的快速发展，除了其理论、模型、结构不断推陈出新以外，AI在各个垂直行业领域也开始遍地开花。在卫星遥感领域，遥感卫星每天产生大量的数据，依靠人工处理效率极低，而搭上AI列车一起发展，也是遥感数据处理未来的发展方向之一。目前，在语义分割中，大多数针对图像的数据增强，是对自然图像的数据增强，例如随机翻转、随机旋转、随机裁剪等。针对遥感图像成像、传输、使用等各个阶段的特点，我们可以有针对性地将其特点用于数据预处理、数据增强等。
图1 遥感工作原理示意图
如图，在卫星获取图像阶段受到太空环境和地球大气的影响，与地面拍摄图像相比，会被噪声、薄云等影响。同时，一些扫描仪出现问题，也可能出现条带状的噪声。而现在的卫星传感器种类多样，且普遍搭载多光谱传感器，相比普通拍摄设备获取的RGB图像，遥感图像可能会有上百个波段（通道）。不同成像方式、不同波段、不同分辨率、不同尺度及数据类型，使得这些异构、多源、海量的遥感数据，注定与普通的自然图像处理存在一定区别。那么在语义分割中，怎么使常规的训练适应于遥感图像的分割呢？
技术方案
针对遥感图像定制网络、数据处理方法等，是未来深度学习在遥感方面应用的重大研究方向。受限于有限的知识，其中最基础的就是针对遥感图像的数据增强或预处理方法。本项目旨在结合遥感数据的特点，修改和定制transforms方法，并将此transforms方法用于遥感图像分割项目，并对比传统transforms方法。项目方案如下：
- 根据特点，思考可能用到的增强方法rs_transforms；
- 基于PaddleSeg的遥感垂类任务，重新打造了rs_transforms方法；
- 使用PaddelSeg进行训练，对比传统的transforms方法。
经过简单的思考，初步想到了以下17种数据增强的方案，如下表：
表1 数据增强的方案
在PaddleSeg中，得益于程序员小哥哥小姐姐的工作，transforms中的代码非常易懂。所以，基于PaddleSeg提供的transforms，我们可以对其进行修改和定制以达到所需的要求。下面可以看一下官方提供的随机垂直翻转的代码，感受一下是不是很容易看懂和修改。
class RandomVerticalFlip:
"""以一定的概率对图像进行垂直翻转。当存在标注图像时，则同步进行翻转。
Args:
prob (float): 随机垂直翻转的概率。默认值为0.1。
"""
def __init__(self, prob=0.1):
self.prob = prob
def __call__(self, im, im_info=None, label=None):
"""
Args:
im (np.ndarray): 图像np.ndarray数据。
im_info (dict): 存储与图像相关的信息。
label (np.ndarray): 标注图像np.ndarray数据。
Returns:
tuple: 当label为空时，返回的tuple为(im, im_info)，分别对应图像np.ndarray数据、存储与图像相关信息的字典；
当label不为空时，返回的tuple为(im, im_info, label)，分别对应图像np.ndarray数据、
存储与图像相关信息的字典和标注图像np.ndarray数据。
"""
if random.random() < self.prob:
im = vertical_flip(im)
if label is not None:
label = vertical_flip(label)
if label is None:
return (im, im_info)
else:
return (im, im_info, label)
基于上面的代码，完成自定义的transforms后，就可以将其Clone到AI Studio的项目中，与PaddleSeg一起使用。具体可移步Gitee:
https://gitee.com/Geoyee/rs_aug
模型训练
数据准备
数据使用L8SPARCS数据集，此数据集由USGS（美国地质勘探局）发布，在PaddleSeg的遥感垂类中也是采用这个数据集做演示项目。下面是USGS关于此数据集的介绍：
This collection of cloud validation data contains 80, 1000 x 1000-pixel, subsets of Landsat 8 OLI/TIRS scenes in .tif format, each with both a manual cloud truth mask and a color composite preview image, in .png format. Each product is provided as a thematic raster, and includes the classes of “cloud”, “cloud shadow”, “snow/ice”, ”water”, and ”flooded”. While these validation images were subjectively digitized by a single analyst, they provide useful information for quantifying the accuracy of clouds flagged by various cloud masking algorithms. The mask and accompanying color composite preview image are provided in PNG image format, and includes all bands from the original Landsat Level-1 data prouduct (GeoTIFF) band-interleaved by pixel (BIP), Landsat Level-1 Quality Assessment (QA)band (GeoTIFF), and its associated Level-1 metadata (MTL.txt file).
简单来说就是这是80张1000x1000的Landsat8 OLI/TIRS数据，共有10个波段，分为阴影、水面阴影、水、冰雪、陆地、云和淹没区7大类（英语不好，如果搞错了请别介意）。考虑到由于淹没区和水面阴影2个类别占比仅为1.8%和0.24%，PaddleSeg将其进行合并，淹没区归为陆地，水面阴影归为阴影，合并后标注包含5个类别。具体请看查看表2：
表2 标签像素类别及颜色
查看一下这个数据集：
# 数据路径
npy_path = 'L8SPARCS/data_npy'
vis_path = 'L8SPARCS/data_vis'
lab_path = 'L8SPARCS/mask'
# 查看
vis_names = os.listdir(vis_path)
random.shuffle(vis_names) # 随机打乱一下
plt.figure(figsize=(20, 4))
for i in range(10):
vis_img = cv2.cvtColor(cv2.imread(os.path.join(vis_path, vis_names[i])), cv2.COLOR_BGR2RGB)
lab_img = cv2.cvtColor(cv2.imread(os.path.join(lab_path, vis_names[i].replace('photo', 'mask'))), cv2.COLOR_BGR2RGB)
plt.subplot(2, 10, i+1);plt.imshow(vis_img);plt.xticks([]);plt.yticks([])
if i == 0: plt.ylabel('Image')
plt.subplot(2, 10, i+11);plt.imshow(lab_img);plt.xticks([]);plt.yticks([])
if i == 0: plt.ylabel('Label')
plt.show()
图2 数据集展示
这里，还需要介绍一下L8SPARCS的波段组成。Landsat8 OLI/TIRS数据共有11个波段，但由于第11波段有点问题，所以数据集中仅有10个波段，每个波段的信息介绍如下表3：
表3 Landsat8 OLI/TIRS波段信息
查看一下图像、增广图像及标签。这里需要注意的有以下几点。首先非几何变化的增加，如对比度、饱和度等变化，只能作用于图像原始的波段上，不能作用于计算出NDVI等，否则就失去了计算的正确性。同时，我们根据官方提供的标签的颜色重新定义了配色表。
# 部分数据增强
transforms = rsg.Compose([
rsg.NDVI(r_band=3, nir_band=4), # 计算NDVI指数并叠在通道11
rsg.NDWI(g_band=2, nir_band=4), # 计算NDWI指数并叠在通道12
rsg.NDBI(nir_band=4, mir_band=5), # 计算NDBI指数并叠在通道13
# 以下的非旋转、缩放等操作仅针对前10个通道
rsg.RandomColor(prob=0.1, alpha_range=[0.8, 1.2], beta_range=[0, 100], band_num=10), # 以0.1的概率随机改变图像的亮度、对比度
rsg.RandomFog(prob=0.1, fog_range=[0.03, 0.56], band_num=10), # 以0.1的概率随机添加薄雾的效果
rsg.RandomSharpening(prob=0.1, laplacian_mode='8-1', band_num=10), # 以0.1的概率随机使用拉普拉斯算子进行图像锐化
rsg.RandomBlur(prob=0.1, ksize=3, band_num=10), # 以0.1的概率随机使用3x3大小的核进行高斯模糊
rsg.RandomSplicing(prob=0.5, direction='Vertical', band_num=10), # 随机以0.5的概率添加竖直方向的拼接色差
rsg.RandomStrip(prob=0.5, strip_rate=0.05, direction='Horizontal', band_num=10), # 随机以0.5的概率添加水平方向的条带噪声
rsg.RandomFlip(prob=0.5, direction='Both'), # 随意以0.5的概率进行两轴翻转
rsg.RandomRotate(prob=0.4, ig_pix=0), # 随机以0.4的概率进行0-90度的旋转
rsg.RandomEnlarge(prob=0.5, min_clip_rate=[0.7, 0.7]), # 随机以0.5的概率放大图像局部到原大小
rsg.RandomNarrow(prob=0.2, min_size_rate=[0.5, 0.5], ig_pix=0), # 随机以0.2的概率缩小全图并填充回原大小
rsg.Normalize(mean=([0] * 10), std=([1] * 10), bit_num=16, band_num=10), # 标准化
rsg.Resize(target_size=1000, interp='NEAREST') # 定义尺寸
])
# 自定义colormap配色表
def colormap(): # 与官方png相同的配色表
return mpl.colors.LinearSegmentedColormap.from_list('cmap', ['#FFFFFF', '#000000', '#00FEFE', '#0C00FB','#969696'], 256)
# 显示
plt.figure(figsize=(20, 8))
for i in range(10):
vis_img = cv2.cvtColor(cv2.imread(os.path.join(vis_path, vis_names[i])), cv2.COLOR_BGR2RGB)
npy_img_path = os.path.join(npy_path, vis_names[i].replace('png', 'npy').replace('photo', 'data'))
lab_img_path = os.path.join(lab_path, vis_names[i].replace('photo', 'mask'))
npy, _, lab = transforms(img=npy_img_path, label=lab_img_path)
npy_RGB = cv2.merge([npy[:, :, 3], npy[:, :, 2], npy[:, :, 1]]) # 官方彩色图
plt.subplot(4, 10, i+1);plt.imshow(vis_img);plt.xticks([]);plt.yticks([]) # 真彩色原图
if i == 0: plt.ylabel('Image')
plt.subplot(4, 10, i+11);plt.imshow(show_img(npy_RGB));plt.xticks([]);plt.yticks([]) # 变化图
if i == 0: plt.ylabel('AUG_Image')
plt.subplot(4, 10, i+21);plt.imshow(show_img(npy[:, :, -3])[:, :, 0], cmap=plt.cm.jet);plt.xticks([]);plt.yticks([]) # NDVI
if i == 0: plt.ylabel('NDVI')
plt.subplot(4, 10, i+31);plt.imshow(show_img(lab)[:, :, 0], cmap=colormap());plt.xticks([]);plt.yticks([]) # 标签图
if i == 0: plt.ylabel('AUG_Label')
plt.show()
图3 数据增强展示
网络训练
网络选用的是PaddleSeg实现的UNet++，加载模型。这里的in_channels=12表示波段数+2，有NDVI与NDWI。
from paddleseg.models import UNetPlusPlus
unet_pp = UNetPlusPlus(
in_channels=12,
num_classes=5,
)
构建训练集，使用自定义的数据增强，验证集同理。仅数据增强，只需要NDVI、NDWI、Resize和Normalize，保证波段数和大小一致。这里from datasets import Dataset拷贝于官方的的Dataset，并有所修改。主要是参考原版，使用了Compose。
from datasets import Dataset
import rs_aug.rs_aug as T
train_transforms = [
T.NDVI(r_band=3, nir_band=4),
T.NDWI(g_band=2, nir_band=4),
T.RandomColor(prob=0.1, alpha_range=[0.8, 1.2], beta_range=[0, 100], band_num=10),
T.RandomFlip(prob=0.5, direction='Both'),
T.RandomRotate(prob=0.4, ig_pix=255),
T.RandomEnlarge(prob=0.5, min_clip_rate=[0.7, 0.7]),
T.RandomNarrow(prob=0.2, min_size_rate=[0.5, 0.5], ig_pix=255),
T.Resize(target_size=(512, 512), interp='NEAREST'),
T.Normalize(mean=([0] * 10), std=([1] * 10), bit_num=16, band_num=10)
# 这里不用计算NDBI，因为没有建筑，数据质量较好，也不用加入条带、拼接色差什么的
]
## 构建训练集
train_dataset = Dataset(
transforms=train_transforms,
dataset_root='L8SPARCS',
num_classes=5,
train_path='L8SPARCS/train.txt',
separator=' ',
mode='train'
)
然后，设置好参数，开始训练。这里使用了Momentum优化器以及学习率多项式衰减策略。损失函数采用的是CrossEntropy损失函数。
import paddle
from paddleseg.models.losses import CrossEntropyLoss
from paddleseg.core import train
# 设置学习率
base_lr = 0.01
lr = paddle.optimizer.lr.PolynomialDecay(base_lr, power=0.9, decay_steps=1000, end_lr=0)
# 设置优化器
optimizer = paddle.optimizer.Momentum(lr, parameters=unet_pp.parameters(), momentum=0.9, weight_decay=4.0e-5)
# 损失
losses = {}
losses['types'] = [CrossEntropyLoss()]
losses['coef'] = [1]
训练过程
相比使用传统数据增强方法的同参数UNet++，两个任务的训练Loss和预测的mIoU、Acc如图4，其中蓝色是transforms，绿色是rs_transforms。可以看到在Loss下降上差距不大，绿色折线稍微低点；而在mIoU和Acc上，绿色折线先低于蓝色折线，在中间部分超过了蓝色折线，最后取得的效果也更胜一筹。表明使用更符合遥感数据的数据增强能够有效的提高模型的精度，而精度先低后高的原因可能是开始时数据变化和通道过多，网络还没有学习到更加有效的参数，而在后期网络学习渐入佳境后，更多增强的效果才显现出来。但是其中还有个问题是（图上看不出来），由于增加了通道和更多的数据增强方法，绿色折线的训练时间显著地大于蓝色折线的训练时间。
图4 训练过程VDL可视化
模型预测
利用训练好的模型，可以预测一张图像，并与原图、标签进行对比。首先，我们将用于展示的图片、用于训练的图片和标签图片的路径准备好。
vis_path = 'L8SPARCS/data_vis/LC81480352013195LGN00_32_photo.png'
img_path = 'L8SPARCS/data_npy/LC81480352013195LGN00_32_data.npy'
lab_path = 'L8SPARCS/mask/LC81480352013195LGN00_32_mask.png'
对图像和标签进行变换。这里主要是需要为图像增加上两个通道的NDVI和NDWI，以保证和模型输入的12通道相同。
img, _, lable = transforms(img=img_path, label=lab_path)
接下来就可以进行预测了：
## 预测
img = paddle.to_tensor(img.transpose((2, 0, 1))[np.newaxis, :]) # 增加一个N通道
pre = infer.inference(model, img)
pred = paddle.argmax(pre, axis=1, dtype='int32').numpy().reshape((512, 512))
这里根据官方提供的标签的颜色重新定义了配色表：
## 自定义colormap配色表
def colormap():
return mpl.colors.LinearSegmentedColormap.from_list('cmap', ['#FFFFFF', '#000000', '#00FEFE', '#0C00FB','#969696'], 256)
最后对结果和原图像及标签进行显示：
## 显示
vis_img = cv2.cvtColor(cv2.imread(vis_path), cv2.COLOR_BGR2RGB)
lab_img = cv2.cvtColor(cv2.imread(lab_path), cv2.COLOR_BGR2RGB)
plt.figure(figsize=(15, 45))
plt.subplot(131);plt.imshow(vis_img);plt.xticks([]);plt.yticks([]);plt.title('Image')
plt.subplot(132);plt.imshow(lab_img);plt.xticks([]);plt.yticks([]);plt.title('Label')
plt.subplot(133);plt.imshow(pred, cmap=colormap());plt.xticks([]);plt.yticks([]);plt.title('Prediction') # 设置配色表
plt.show()
图5 预测结果
后记
在与一般的transforms训练比较中，可以看到使用rs_transforms：
- 训练的轮数有些少，继续训练还能取得更好的效果；
- 效果差距不是很大，因为数据质量较好，没有拼接、条带、使用不同传感器的图像（分辨率不同，模糊的增益应该会更大）等等。大家可以在更杂的数据上进行训练、预测和分析；
- 使用RS_AUG的训练时间，是自带增强训练时间的2倍还多。因为增强的操作较多，同时添加了NDVI和NDWI，所以网络也相应输入更多的波段和参数，大家可以通过PCA进行降维（这里会涉及Normalize会有点问题，因为最小最大值已经改变）。
AI要在每个领域内有所发展，不只需要计算机，也需要该领域的人员一起参与。AI在遥感领域，应该根据遥感图像获取、传输和应用等特点，设计专门的数据增强和网络结构。奈何我还不知道如何修改新的结构？也不知道如何使用更针对性的数据增强方法？以上只是我浅显的尝试。希望未来的时间里，PaddleSeg可以在大家的贡献下，拥有更完善的遥感垂类，不只是数据增强，新的Loss、新的网络能不断涌现，让大家更方便地训练、预测和部署。欢迎大家Fork这个项目！
AI Studio链接：
https://aistudio.baidu.com/aistudio/projectdetail/1485445
如在使用过程中有问题，可加入官方QQ群进行交流：778260830。
如果您想详细了解更多飞桨的相关内容，请参阅以下文档。
·飞桨官网地址·
https://www.paddlepaddle.org.cn/
·飞桨开源框架项目地址·
GitHub: https://github.com/PaddlePaddle/Paddle
Gitee: https://gitee.com/paddlepaddle/Paddle
飞桨(PaddlePaddle)以百度多年的深度学习技术研究和业务应用为基础，是中国首个开源开放、技术领先、功能完备的产业级深度学习平台，包括飞桨开源平台和飞桨企业版。飞桨开源平台包含核心框架、基础模型库、端到端开发套件与工具组件，持续开源核心能力，为产业、学术、科研创新提供基础底座。飞桨企业版基于飞桨开源平台，针对企业级需求增强了相应特性，包含零门槛AI开发平台EasyDL和全功能AI开发平台BML。EasyDL主要面向中小企业，提供零门槛、预置丰富网络和模型、便捷高效的开发平台；BML是为大型企业提供的功能全面、可灵活定制和被深度集成的开发平台。
END