Landsat 5 多光谱数据分类指导手册原作者：Pratyush Tripathy

翻译：荆雪涵

姐妹篇雪涵：【博客翻译】CNN 与中分辨率遥感影像分类​zhuanlan.zhihu.com

深度学习在很多邻域使用，解决了许多复杂的问题，在空间分析领域也不例外。既然你已经点进来了，那估计你对遥感影像以及Landsat 5 TM数据有一定的了解。我们只需要了解一下机器学习的基础知识，便可以快速掌握这个教程的内容。如果你没有接触过机器学习，其实用一句话说就是，建立一系列特征值X与目标值Y的联系。我们通过提供带有标签的训练数据，从而习得X和Y之间的联系，最后我们可以使用训练好的模型预测那些还没有标签的数据。

基于卫星数据学术问题

卫星影像中的一些土地类别是比较难区分的，比如城市用地/裸地/矿区，这些地物的光谱值十分的相似，所以过去几十年，区分相似光谱的地物一直是一个难题。传统的监督分类法或者非监督分类法，虽然在分类领域十分常见，但是也带来一些问题，我们看一下下面的例子。

上图，如果我们用一条垂线来分割这张图片，我们将这条垂线从左往右移动，如果我们规定垂线右边都归为房子类，很显然这个分类方法是有问题的，因为房子和树木的分布并没有一个明确的边界。所以我们是无法用一条直线来进行分类的。但是这并不意味我们无法将上面这张图片正确分类。

在上面这张动图中，我们分别用蓝线，红线，绿线对树木和房子进行分类，我们规定线的右边为房子类别。如果用红线分类的话，绝大多的房子都被正确分类，只有一个房子没有被正确分类，并且还有一颗树木被错误的归为了房子类别。如果我们用蓝线的话，我们就可以正确分类所有的房子，但是会有三棵树被错误的分类。在使用蓝线分类的例子中，蓝线将所有的房子都正确分类，我们称之为高召回(Recall)，但是在房子类别中，有三颗树木被错误分类，我们称为低精度(Precision)。类似的，如果我们用绿线分类，绿线的右边都是房子，我们就称为高精度(Precision)不过召回率就会低一些，因为有三个房子没有被识别到。很多时候，我们都要进行精度率和召回率之间的博弈。

上面树木和房子分类的例子，其实就是我们对遥感影像进行土地分类的一个类别。所以每一幅遥感影像的分类重点是要具体问题具体分析的。举个例子，如果你要保证所有的城市用地都被正确归类，那么你的最终结果中就会出现其他一些和城市用地光谱类似的地物，也就是说你需要一个高召回率的分类器。反过来看，如果你的最终结果不可以出现出城市用地外的其他地物，你就需要一个高精度的分类器。

本文使用的数据

我们将使用Landsat 5 TM数据的波段2-波段7，共计六个波段的数据，我们的目标是从遥感影像中区分出城市用地(Built-up)。2011年印度班加罗尔的Landsat 5 的多光谱数据以及城市用地二值图层，将被用作训练数据和测试数据。我们将使用2011年海得拉巴的数据来进行预测。

因为我们使用了标注数据，所以我们的方法称为监督机器学习。

我们将使用Google的Tensorflow框架，并使用Python语言构建一个神经网络。我们需要以下几个库。pytrsgis 用于读写GeoTIFF格式的数据

scikit-learn 用于数据预处理以及精度检测

numpy 用于基本数组操作

首先将影像放入一个文件夹里，使用如下代码读取GeoTIFF影像图片。

import os

from pyrsgis import raster

os.chdir("E:\\yourDirectoryName")

mxBangalore = 'l5_Bangalore2011_raw.tif'

builtupBangalore = 'l5_Bangalore2011_builtup.tif'

mxHyderabad = 'l5_Hyderabad2011_raw.tif'

# Read the rasters as array

ds1, featuresBangalore = raster.read(mxBangalore, bands='all')

ds2, labelBangalore = raster.read(builtupBangalore, bands=1)

ds3, featuresHyderabad = raster.read(mxHyderabad, bands='all')

代码中的raster模块可以同时读取GeoTIFF影像的地理坐标以及DN值，并且存放于numpy数组中。关于pyrsgis的详细介绍可以参考此页面。

使用如下代码可以将数组的尺寸打印出来。

print("Bangalore multispectral image shape: ", featuresBangalore.shape)

print("Bangalore binary built-up image shape: ", labelBangalore.shape)

print("Hyderabad multispectral image shape: ", featuresHyderabad.shape)

Bangalore multispectral image shape: 6, 2054, 2044

Bangalore binary built-up image shape: 2054, 2044

Hyderabad multispectral image shape: 6, 1318, 1056

从上面的输出结果可以看出班加罗尔的两幅数据的行列数是一致的，都是2054*2044。两幅遥感影像数据波段数也都为6。模型会根据所有波段的DN值来判断像元是否归为城市用地，所以遥感影像的波段数必须保持一致，波段顺序也必须一致。

下面我们要将我们的遥感影像数据转化为一个二维数组，这样才可以输入到机器学习的模型中去。数组的每一行分别代表六个波段的DN值。我们使用pyrsgis中的convert模块可以实现数转化。见下图。

from pyrsgis.convert import changeDimension

featuresBangalore = changeDimension(featuresBangalore)

labelBangalore = changeDimension (labelBangalore)

featuresHyderabad = changeDimension(featuresHyderabad)

nBands = featuresBangalore.shape[1]

labelBangalore = (labelBangalore == 1).astype(int)

print("Bangalore multispectral image shape: ", featuresBangalore.shape)

print("Bangalore binary built-up image shape: ", labelBangalore.shape)

print("Hyderabad multispectral image shape: ", featuresHyderabad.shape)

Bangalore multispectral image shape: 4198376, 6

Bangalore binary built-up image shape: 4198376

Hyderabad multispectral image shape: 1391808, 6

上面代码的第七行，我们提取了数值为1的像元，这样是为了避免数据中出现NoData的像元，NoData通常是极大或者极小的数值，会对训练过程造成影像。

我们会将数据分组，分别用于训练和验证。这样我们才可以说训练出来模型同样适用于预测新的数据。否则，模型就会过适(Overfitting)。

from sklearn.model_selection import train_test_split

xTrain, xTest, yTrain, yTest = train_test_split(featuresBangalore, labelBangalore, test_size=0.4, random_state=42)

print(xTrain.shape)

print(yTrain.shape)

print(xTest.shape)

print(yTest.shape)

(2519025, 6)

(2519025,)

(1679351, 6)

(1679351,)

上面代码中的test_size=0.4代表我们把数据分开，60%为训练数据，40%为验证数据。

许多机器学习模型都需要将数据进行归一化，也就是说我们将DN值的范围变化到0-1之间。

需要注意的是，归一化通常需要计算数据实际的最大最小值来确定数据的范围。不过这里为了简化过程我们使用0-255作为我们数据的范围。

我们还需要将寻来你数据数组从二维变化成三维，新数组的每一行将代表一个像元的信息。

```python

# Normalise the data

xTrain = xTrain / 255.0

xTest = xTest / 255.0

featuresHyderabad = featuresHyderabad / 255.0

# Reshape the data

xTrain = xTrain.reshape((xTrain.shape[0], 1, xTrain.shape[1]))

xTest = xTest.reshape((xTest.shape[0], 1, xTest.shape[1]))

featuresHyderabad = featuresHyderabad.reshape((featuresHyderabad.shape[0], 1, featuresHyderabad.shape[1]))

# Print the shape of reshaped data

print(xTrain.shape, xTest.shape, featuresHyderabad.shape)

(2519025, 1, 6) (1679351, 1, 6) (1391808, 1, 6)

上述准备工作完成后，我们可以使用Keras库建立模型。我们会使用一个顺序模型(Sequential Model)，顾名思义就是将神经网络层依次按顺序连接。我们有一个输入数据层，节点数等于波段数。然后是一个节点为14的reLu激活函数层。最后是节点数为2的softmax激活函数层，因为我们最终的预测结果为两类：城市用地，非城市用地。激活函数的更多内容可以参考此处。

from tensorflow import keras

# Define the parameters of the model

model = keras.Sequential([

keras.layers.Flatten(input_shape=(1, nBands)),

keras.layers.Dense(14, activation='relu'),

keras.layers.Dense(2, activation='softmax')])

# Define the accuracy metrics and parameters

model.compile(optimizer="adam", loss="sparse_categorical_crossentropy", metrics=["accuracy"])

# Run the model

model.fit(xTrain, yTrain, epochs=2)

从代码中的第十行看出，我们的模型使用了adam优化器(其他的优化器可以参考这里)，损失函数类型选择categorical-sparse-crossentropy(类别稀疏交叉熵)，具体含义可以参考这里。最后我们使用准确度(accuracy)来衡量模型的表现。

最后我们使用训练数据，训练两个循环。模型的训练时间取决于数据量以及电脑性能。训练中可以看到如下消息。

训练完成后，我们将使用测试数据来验证模型的预测精度。代码如下，

from sklearn.metrics import confusion_matrix, precision_score, recall_score

# Predict for test data

yTestPredicted = model.predict(xTest)

yTestPredicted = yTestPredicted[:,1]

# Calculate and display the error metrics

yTestPredicted = (yTestPredicted>0.5).astype(int)

cMatrix = confusion_matrix(yTest, yTestPredicted)

pScore = precision_score(yTest, yTestPredicted)

rScore = recall_score(yTest, yTestPredicted)

print("Confusion matrix: for 14 nodes\n", cMatrix)

print("\nP-Score:%.3f, R-Score:%.3f" % (pScore, rScore))

softmax函数会计算每一类别的概率。之前代码的第六行我们只提取了城市用地类别。空间分析模型的优劣不像其他机器学习模型那么的好评判，我们通常要将位置信息也考虑进去，而不是简单的把误差叠加。因此，我们将使用混淆矩阵，召回率，精确率等一系列参数来评估模型表现。

从混淆矩阵中可以看出，有上千个城市用地像元被归类为非城市用地类别，但是这一部分的比例相对数据的总量来说是微不足道的。测试数据的精确度和召回率超过了0.8.

你可以尝试花时间将模型多训练几次，从而找到合适的隐藏层个数，合适的隐藏层节点个数合适的训练循环次数。你还可以使用NDBI或者NDVI数据作为输入训练数据。当我们训练好一个合适精度的模型后，我们可以使用这个模型来预测一张新的GeoTIFF遥感图。我们后续只需将模型稍作修改，还可以用于其他的应用。

predicted = model.predict(featuresHyderabad)

predicted = predicted[:,1]

#Export raster

prediction = np.reshape(predicted, (ds.RasterYSize, ds.RasterXSize))

outFile = 'Hyderabad_2011_BuiltupNN_predicted.tif'

raster.export(prediction, ds3, filename=outFile, dtype='float')

需要注意的是，预测后的结果是连续的概率值，并不是一张0或1的二值图。不过使用GIS软件，设置一个阈值，可以很快的将概率数值二值化。得到下面右图的结果。

模型的精度已经在上文评估过了，不过你可以对新数据的预测结果在进行一次评估。遥感影像分类除了之前提到的一些难点，季节，地域，光谱特征多样性也同样会对模型带来一些不确定的影像。

本文介绍了很简单的神经网络，较为复杂的卷积神经网同样被用来做遥感分类，你可以阅读本文的姊妹篇。

文中的代码数据可以打开此链接下载。