多种机器学习算法实现遥感影像监督分类(LR,KNN,NB,SVM,RF)

引言

遥感影像的监督分类是用被确认类别的样本像元去识别其他未知类别像元的过程。它是在分类之前通过目视判读和野外调查，对每一种类别选取一定数量的训练样本，通过数学模型对其进行训练，根据训练结果将待分数据划分到与其最相似的样本类，以此完成对整个图像的分类。

本文基于逻辑回归、最近邻、朴素贝叶斯、支持向量机、随机森林5种分类算法，选择一幅Landsat8 OLI多光谱影像，为每类地物选择一定数量的样本（shp格式）进行监督分类，将不同算法的分类结果导出为带有空间参考的tif影像。

结果展示

准备遥感影像

首先准备一幅用于分类的遥感影像，本文用到的影像为Landsat8 OLI多光谱影像，可从地理空间数据云下载，对下载后的影像进行辐射定标、大气校正等操作，并适当裁剪，本文影像裁剪后大小为5000×5000，水平分辨率为30m。由于需要对影像进行目视解译以选取训练样本，应选择各种地物区分度较好的假彩色图像作为分类影像。根据Landsat8 OLI不同波段组合的应用场景遥感影像解译，如何更好的识别地物类型？ ，本文选择了SWIR1、NIR、Blue波段的假彩色影像，将其另存为tif格式数据。

预处理后的遥感影像（真彩色）

制作训练样本

打开tif格式的遥感影像，通过目视判别选取不同地物类型的样本，本文所使用的影像地理位置为昆明市周边，该区域为亚热带气候，地表植被类型丰富，且生长周期不尽相同，因此地物判别有一定的难度。我们结合地物光谱曲线，将这幅影像的地物分为6类：城镇用地、天然林地、水体、农业用地1、农业用地2、农业用地3，农业用地被划分成了3类，具体植被类型不是很清楚，但在假彩色影像中确实有明显区别。使用GIS软件制作训练样本，在影像中为每类地物新建点要素，不同地物类别的样本总数不应相差太远，把选好的地物样本保存为shp文件。注意：遥感影像和训练样本数据的坐标参考系应保持一致（最好是投影坐标系）。

选取地物样本（假彩色）

提取训练样本

我们需要将shp文件中点要素对应的栅格值提取出来，并对每类地物的值进行标记，作为监督分类的训练数据和标签。对遥感影像某点数值的提取通过geopandas和rasterio库实现，geopandas官方文档中有一个采样栅格的点数据的示例教程Using GeoPandas with Rasterio to sample point data，我们直接模仿示例代码对tif影像进行操作。

首先通过rasterio库读取tif遥感影像，得到一个rasterio.io.DatasetReader对象，该对象的很多属性记录了遥感影像的元数据信息rasterio 1.4b1.dev0 documentation，接下来在导出tif格式分类结果的时候会用到它们。

src = rasterio.open(rasterpath)

使用geopandas打开所有地物类别的shp文件，它以GeoDataFrame的结构来表示空间矢量数据，包括index、data、geometry三部分，geometry记录了每条数据的矢量结构（可以是点、线、面要素）。通过属性geometry.x和geometry.y可以访问到每个矢量点的xy坐标。

gdf = gpd.read_file(filepath)

geometry存储了点要素信息

将shp文件中所有点要素的xy坐标导出到列表中。

coord_list = [(x, y) for x, y in zip(gdf.geometry.x, gdf.geometry.y)]

使用DatasetReader的sample()方法将栅格影像对应xy坐标的值提取出来，由于本文使用的影像有三个波段，因此每个坐标都对应着3个像元值。

X_list += [x for x in raster.sample(coord_list)]

最后将每个样本的标签也使用列表保存起来，地物类型使用0、1、2．．．数字表示。

Y_list += [idx] * len(coord_list)

本文遥感影像分类的输入特征数为3，即R、G、B通道3个波段，可以用三维散点图展示训练样本在每个特征中的分布情况。

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.scatter(X_train[:, 0], X_train[:, 1], X_train[:, 2], c=Y_train, s=50)
plt.show()

可以看出6种地物类型的输入特征还是有较好的区分度的，接下来读取遥感影像三个波段的完整数据，用作测试集，使用DatasetReader的read()方法读取每个波段的栅格数组，属性indexes记录了每个波段的索引序号，一般序号默认为1、2、3．．．。

bands = [np.ravel(raster.read(band)) for band in range(1, 4)]

将3个波段的数组合并为X_test，再将数组转置。

X_test = np.vstack((bands[0], bands[1], bands[2]))
X_test = X_test.T

遥感影像分类

本文准备了LR、KNN、NB、SVM、RF5种算法对影像影像进行分类，由于最近邻、支持向量机、线性回归等算法基于“距离”构建解算器，需要对训练和预测数据进行标准化处理，从而避免方差较大的特征主导目标函数，导致解算器难以正确的从其它特征中进行学习。对逻辑回归也应用标准化能够提高模型的收敛速度，朴素贝叶斯和随机森林则不需要数据的标准化操作不能不用也不可乱用的标准化和归一化处理。

提前将模型参数优化好后，利用pipeline工具构建预处理和分类估算器。

model = {'LR': make_pipeline(StandardScaler(), LogisticRegression(solver='lbfgs', multi_class='multinomial')),
         'KNN': make_pipeline(StandardScaler(), KNeighborsClassifier(n_neighbors=3)),
         'NB': GaussianNB(),
         'SVM': make_pipeline(StandardScaler(), NuSVC(nu=0.21, kernel='rbf')),
         'RF': RandomForestClassifier(n_estimators=57, random_state=0)}

分别对每个模型进行训练和预测。

# 使用逻辑回归算法进行分类学习
model['LR'].fit(X_train, Y_train)
Y_pred = model['LR'].predict(X_test)

导出分类结果

得到分类结果后，将不同地物类型的标签转化为R、G、B值，分别用3个数组存储，进而导出成带有空间参考的tif格式数据，以便在遥感影像处理软件中进行比较和分析。

原始影像与支持向量机分类结果

使用rasterio库以写入形式打开一个新的tif文件，指定mode='w'，由于是对栅格文件进行写入，因此影像的头文件信息需要设置，否则会报错。遥感影像的大小、分辨率等属性在分类过程中没有发生变化，直接将输入影像的头文件信息赋给结果影像即可，包括driver（影像存储格式”GTiff”或”JPEG”）、width和height（影像的长和宽）、count（波段数）、crs（坐标参考系）、transform（设置像元空间和地理空间的转换关系）等。

with rasterio.open(path, mode='w', driver='GTiff', width=raster.width, height=raster.height, count=3, crs=raster.crs, transform=raster.transform, dtype='uint8') as dst:
    dst.write(band1.reshape(raster.width, raster.height), 1)
    dst.write(band2.reshape(raster.width, raster.height), 2)
    dst.write(band3.reshape(raster.width, raster.height), 3)

关于分类结果的验证

通常一幅资源环境遥感影像（30m分辨率）有数千万个像元，我们不可能人工标注每一个像元所对应的地物类别，可以通过目视判别或实地观测在遥感影像中重新选择一部分地物样本（验证样本中应不含训练样本），并标记所属的类别，提取分类结果影像中与之对应的标签值，计算分类结果在验证样本中的错分率作为用户精度。也可以将分类效果较好的地物影像作为近似真值，用于评价其他模型或参数所预测的分类结果。

完整代码

import os
import re
from datetime import datetime as dt
import numpy as np
from numpy.random import shuffle
import matplotlib.pyplot as plt
import rasterio
import geopandas as gpd
from sklearn.pipeline import make_pipeline
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC, NuSVC
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler


def rastersample(raster, shpdir):
    idx = 0
    X_list = []
    Y_list = []
    filelist = os.listdir(shpdir)
    for filename in filelist:
        if re.split('\.', filename)[-1] == 'shp':
            filepath = os.path.join(shpdir, filename)
            gdf = gpd.read_file(filepath)
            coord_list = [(x, y) for x, y in zip(gdf.geometry.x, gdf.geometry.y)]
            X_list += [x for x in raster.sample(coord_list)]
            Y_list += [idx] * len(coord_list)
            idx += 1

    X_train = np.array(X_list)
    Y_train = np.array(Y_list)

    idx = np.arange(len(X_list))
    shuffle(idx)
    X_train = X_train[idx]
    Y_train = Y_train[idx]

    return X_train, Y_train


def rasterCLF(raster, X_train, Y_train):
    bands = [np.ravel(raster.read(band)) for band in range(1, 4)]

    X_test = np.vstack((bands[0], bands[1], bands[2]))
    X_test = X_test.T

    model = {'LR': make_pipeline(StandardScaler(), LogisticRegression(solver='lbfgs', multi_class='multinomial')),
             'KNN': make_pipeline(StandardScaler(), KNeighborsClassifier(n_neighbors=3)),
             'NB': GaussianNB(),
             'SVM': make_pipeline(StandardScaler(), NuSVC(nu=0.21, kernel='rbf')),
             'RF': RandomForestClassifier(n_estimators=57, random_state=0)}

    model['LR'].fit(X_train, Y_train)

    Y_pred = model['LR'].predict(X_test)

    return Y_pred


def write2raster(raster, arr, path):
    gridcount = raster.width * raster.height
    band1 = np.empty(gridcount)
    band2 = np.empty(gridcount)
    band3 = np.empty(gridcount)

    color_list = [[235, 51, 36], [119, 67, 66], [117, 249, 77], [50, 130, 246], [240, 134, 80], [115, 43, 245]]
    for item in range(len(color_list)):
        band1[arr == item] = color_list[item][0]
        band2[arr == item] = color_list[item][1]
        band3[arr == item] = color_list[item][2]

    with rasterio.open(path, mode='w', driver='GTiff', width=raster.width, height=raster.height, count=3, crs=raster.crs, transform=raster.transform, dtype='uint8') as dst:
        dst.write(band1.reshape(raster.width, raster.height), 1)
        dst.write(band2.reshape(raster.width, raster.height), 2)
        dst.write(band3.reshape(raster.width, raster.height), 3)


def main():
    rasterpath = r"C:\Users\XXX\PycharmProjects\MLtest\Raster\Landsat8OLI.tif"
    shpdir = r"C:\Users\XXX\PycharmProjects\MLtest\SHP"

    src = rasterio.open(rasterpath)

    X_train, Y_train = rastersample(src, shpdir)
    pred = rasterCLF(src, X_train, Y_train)

    # 设置分类结果影像的文件名
    tifdir = os.path.dirname(rasterpath)
    newtif = 'LR.tif'
    newpath = os.path.join(tifdir, newtif)

    write2raster(src, pred, newpath)

    src.close()


if __name__ == '__main__':
    dt1 = dt.now()
    main()
    print(f'Duration: {dt.now() - dt1}')