基于Python的Tensorflow卫星数据分类神经网络

  • 2019 年 10 月 4 日
  • 筆記

深度学习已经占据了解决复杂问题的大多数领域,地理空间领域也不例外。文章的标题让您感兴趣,因此希望熟悉卫星数据集 ; 目前,Landsat 5 TM。机器学习(ML)算法如何工作的知识很少,将帮助快速掌握这本动手教程。对于那些不熟悉ML概念的人,简而言之,它是建立一个实体的一些特征(特征或X)与其他属性(值或标签或Y)之间的关系 – 提供了大量的例子(标记数据) )到模型,以便从中学习,然后预测新数据(未标记数据)的值/标签。这对于机器学习来说已经足够理论了!

卫星数据的一般问题:

卫星数据中的两个或更多要素类(例如,建造/贫瘠/采石场)可具有相似的光谱值,这使得该分类在过去的几十年中成为具有挑战性的任务。由于上述问题,传统的监督和无监督方法不能成为完美的分类器,尽管它们可以稳健地执行分类。但是总会有相关的问题。通过以下示例来理解这一点:

在上图中,如果使用垂直线作为分类器并仅沿着x轴移动它,使其将所有图像分类为右侧作为房屋,则答案可能不是直截了当的。这是因为数据的分布是这样的,即不可能只用一条垂直线将它们分开。但是,这并不意味着房屋根本无法归类!

假设使用红线(如上图所示)来分隔这两个要素。在这种情况下,大多数房屋都是由分类器确定的,但房子仍被遗漏,一棵树被误分类为房屋。为了确保不会留下任何一个房子,可以使用蓝线。在这种情况下,分类器将覆盖所有房屋; 这被称为高召回率。然而,并非所有的分类图像都是真正的房屋,这被称为低精度。同样,如果使用绿线,所有分类为房屋的图像都是房屋; 因此,分类器具有高精度。在这种情况下召回的次数会减少,因为还有三所房子被遗漏了。在大多数情况下,这种权衡 在精确度和召回之间保持。

上面展示的房屋和树木问题类似于建筑物,采石场和贫瘠土地的情况。卫星数据的分类优先级可能因目的而异。例如,如果想确保所有的组合单元被归类为组合,没有留下任何东西,并且你更少关心具有类似签名的其他类的像素被归类为组合,那么一个模型与需要高召回率。相反,如果优先级是仅对纯组合像素进行分类而不包括任何其他类像素,并且可以放弃混合的组合像素,则需要高精度分类器。通用模型将使用房屋和树木的红线来保持精确度和召回之间的平衡。

当前范围中使用的数据

在这里,将把Landsat 5 TM的六个波段(波段2 – 波段7)视为特征,并尝试预测二进制构建类。2011年为班加罗尔及其相应的二元建筑层获得的多光谱Landsat 5数据将用于训练和测试。最后,2005年为海德拉巴收购的另一个多光谱Landsat 5数据将用于新的预测。

由于使用标记数据来训练模型,因此这是一种受监督的ML方法。

多光谱训练数据及其相应的二进制构建层

将在Python中使用Google的Tensorflow库来构建神经网络(NN)。需要以下其他库,请确保提前安装它们:

  • pyrsgis – 读写 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')

pyrsgis包的栅格模块分别将GeoTIFF的地理位置信息和数字编号(DN)值作为NumPy数组读取。

打印一下读过的数据的大小。

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

从输出中可以看出,班加罗尔图像中的行数和列数是相同的,并且多光谱图像中的层数是相同的。该模型将基于所有频带上的相应DN值来学习确定像素是否构建,因此,多光谱图像应具有以相同顺序堆叠的相同数量的特征(频带)。

现在将数组的形状更改为二维数组,这是大多数ML算法所期望的,其中每行代表一个像素。pyrsgis包的转换模块将做到这一点。

重组数据的架构

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像素导致的问题,这些像素通常具有极高和极低的值。

现在,将分割数据以进行训练和验证。这样做是为了确保模型没有看到测试数据,并且它对新数据的表现同样出色。否则,模型将过度拟合并且仅在训练数据上表现良好。

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。

包括NN在内的许多ML算法都期望归一化数据。这意味着直方图在一定范围(此处为0到1)之间被拉伸和缩放。将规范化功能以满足此要求。可以通过减去最小值并除以范围来实现归一化。由于Landsat数据是8位数据,因此最小值和最大值分别为0和255(2⁸= 256个值)。

请注意,从标准化数据计算最小值和最大值始终是一个好习惯。为避免复杂性,将在此处坚持使用8位数据的默认范围。

另一个额外的预处理步骤是将特征从二维重塑为三维,使得每行代表单个像素。

# 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构建模型。首先,将使用顺序模型,一个接一个地添加图层。有一个输入层,节点数等于nBands。使用具有14个节点和“ relu ”作为激活功能的一个隐藏层。最后一层包含两个节点,用于二进制构建类,具有' 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)

神经网络架构

如第10行所述,使用' adam '优化器编译模型。(还有其他几个你可以检查。)现在将使用的损失类型是分类 – 稀疏 – 交叉熵。模型性能评估的度量标准是“ 准确性 ”。

最后,使用两个时期(或迭代)在xTrain和yTrain上运行模型。根据数据大小和计算能力,安装模型需要一些时间。模型编译后可以看到以下内容:

预测单独保存的测试数据的值,并执行各种精度检查。

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 nodesn", cMatrix)  print("nP-Score: %.3f, R-Score: %.3f" % (pScore, rScore))

所述SOFTMAX函数产生用于每个类类型单独的列的概率值。只提取第一类(构建),如上面代码片段中的第六行所述。用于地理空间相关分析的模型变得难以评估,因为与其他一般ML问题不同,依赖于广义总结误差是不公平的; 空间位置是获胜模型的关键。因此,混淆矩阵,精度和召回可以更清晰地反映模型的表现。

终端中显示的混淆矩阵,精度和召回

如上面的混淆矩阵所示,有数千个组合像素被分类为非组合,反之亦然,但与总数据大小的比例较小。在测试数据上获得的精度和召回率大于0.8。

总是可以花一些时间并执行一些迭代来找到隐藏层的最佳数量,每个隐藏层中的节点数以及获得准确性的时期数。一些常用的遥感指数,如NDBI或NDWI,也可以在需要时用作特征。达到所需精度后,使用模型预测新数据并导出GeoTIFF。具有微小调整的类似模型可以应用于类似的应用。

predicted = model.predict(feature2005)  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')

请注意,使用预测的概率值导出GeoTIFF,而不是其阈值化的二进制版本。总是可以在以后的GIS环境中对浮点类型图层进行阈值处理,如下图所示。

Hyderabad构建层由模型使用多光谱数据预测

已经精确评估了模型的准确性并进行了调用 – 还可以对新预测的栅格进行传统检查(例如kappa系数)。除了上述卫星数据分类的挑战之外,其他直观的限制包括由于光谱特征的变化,模型无法预测在不同季节和不同区域获得的数据。

在本文中使用的模型是NN的一个非常基本的架构,包括卷积神经网络(CNN)在内的一些复杂模型已经被研究人员证明可以产生更好的结果。这种分类的主要优点是一旦模型被训练就具有可扩展性。

请在此处查找使用的数据和完整脚本。

https://github.com/PratyushTripathy/Landsat-Classification-Using-Neural-Network