引言
由于笔者从事二分类问题的研究，另外实验室也有此类算法的实现需要，因此从该文开始，笔者将连续讲解基于Python封装函数进行遥感图像特征的分割方法与实现。本文将以NDVI特征为例，进行农田环境下建筑物的分割。本章节代码共分成两部分，第一部分，输出值为图像灰度化后的阈值，第二部分输出为NDVI特征阈值，根据NDVI阈值，在ArcMap中实现结果输出及展示。
主要方法原理
(1)最大类间方差法(OTSU)
简单的说，这种算法假设一副图像由前景色和背景色组成，通过统计学的方法来选取一个阈值，使得这个阈值可以将前景色和背景色尽可能的分开。或者更准确的说是在某种判据下最优，与数理统计领域的 fisher 线性判别算法其实是等价的。
假设取一个最优阈值把原图像分为前景色(A部分)与背景色(B部分)，那么两部分的类间方差越大，说明两部分差别越大，便能有效的分割图像。所以该算法最关键的是找到最优阈值。当数据的灰度直方图呈现双峰分布时，二分类的Otsu算法就能够较好地实现图像分割。原始图像灰度直方图Otsu算法下的分割结果
然而Otsu 方法不是万能的。当目标与背景的大小比例悬殊时，类间方差准则函数可能呈现双峰或多峰，此时效果不好。这时就要考虑其他的办法了，例如改进式多阈值OUST算法，当然这是后话，本节只探讨二分类问题。图像中的多峰情况
数据介绍
原始图像为4波段遥感图像，分别为蓝、绿、红、近红外波段，数据类型为uint16。下图分别显示了本实验所用的真彩色图像和NDVI数据：NDIV图像RGB图像
Python代码实现(一)
import cv2
import numpy as np
import gdal
# Step1. 定义加载图像的函数
def image_open(image):
data = gdal.Open(image)
if data == "None":
print("无法打开数据")
return data
#Step2. 归一化处理
File_Path = r"E:\yync\try\edata\NDVI批处理实验\结果阈值实验遥感图像_NDVI.tif"
image = image_open(File_Path).ReadAsArray()
Normalize_data = np.zeros(image.shape)
for i in range(image.shape[0]):
for j in range(image.shape[1]):
Normalize_data[i][j] = (image[i][j] - np.min(image)) / (np.max(image) - np.min(image))
Normalize_data = Normalize_data*256
data = Normalize_data.astype(np.uint8)
#Step3. 运行封装OTSU函数并输出灰度化后的阈值
ret, th2 = cv2.threshold(data, 0, 256, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
print(ret)
代码关键点说明：
代码中的Step2采用了传统的归一化处理方式，公式如下：公式中，x为需要归一化的元素，而Min和Max分别是图像中的最小值和最大值。
基于上述归一化的结果进行uint8级别的灰度化操作(*256)
Python代码实现(二)
import numpy as np
import gdal
# Step1. 定义加载图像的函数
def image_open(image):
data = gdal.Open(image)
if data == "None":
print("无法打开数据")
return data
#Step2. 读取图像为np数组
File_Path = r"E:\yync\try\edata\NDVI批处理实验\结果阈值实验遥感图像_NDVI.tif"
image = image_open(File_Path).ReadAsArray()
#Step3. 输出NDVI分割阈值(手动输入NDVI_Normalize的分子，分子即为主程序输出的数值)
NDVI_Normalize = 153 / 256
NDVI_Threshold = NDVI_Normalize*(np.max(image) - np.min(image)) + np.min(image)
print(NDVI_Threshold)处理结果
运行完代码块(二)之后，我们会得到一个NDVI的分割阈值，本例中图像的阈值为0.406.
分割结果图如下图所示：黑色为建筑区和道路，白色为田地
根据结果可得，分割效果并不是很好，因此，考虑其他特征或者预处理方法，可增大提取成功的概率。
