4. 单波段阈值法：近红外单波段阈值法、基于直方图分析的自动阈值选取（Otsu算法）
各位同学，今天我们来聊聊水体提取里最经典、也最基础的方法——单波段阈值法。
你可能会问，都什么年代了，还用这么简单的方法？
嗯，我个人的经验是：简单的方法往往最可靠。尤其是在做快速水质监测或者应急响应时，你根本没时间跑复杂的深度学习模型。这时候，单波段阈值法就是你的救命稻草。
4.1 为什么是近红外波段？
先说说原理。水体在近红外波段（NIR，约0.7-2.5μm）有一个非常明显的特征——几乎完全吸收。
你想想看，太阳光打到水面上，近红外光基本被水分子吃掉了，反射回来的信号极弱。而植被、土壤、建筑物这些地物，在近红外波段反射率都很高。
这就形成了一个天然的对比度。说白了，就是水是黑的，其他东西是亮的。
核心公式：
如果 NIR 反射率 < 阈值 T，则判定为水体；否则为非水体。
我在项目中遇到过这样的情况：用Landsat 8的Band 5（近红外）做太湖的水体提取，阈值设在0.1左右，效果出奇的好。当然，不同卫星、不同季节，阈值要微调。
4.2 手动阈值法：简单但需要经验
手动选阈值，说白了就是靠眼睛看、靠经验猜。
具体怎么做？
- 加载近红外波段影像
- 查看影像的灰度分布
- 找一个能区分水体和陆地的值
- 反复试，直到满意为止
听起来很粗糙对吧？但说实话，我做了十年遥感，有时候还是靠这招。尤其是在数据质量不好的时候，自动算法反而容易翻车。
我的小技巧：
先取一个水体样本点，记录它的NIR值。再取一个陆地样本点，记录它的NIR值。阈值就取这两个值的中间数。简单粗暴，但有效。
4.3 Otsu算法：让计算机帮你选阈值
手动选阈值太累了，而且不同影像要调不同参数。这时候，Otsu算法就派上用场了。
Otsu算法，也叫最大类间方差法。它的核心思想是：找到一个阈值，让水体和陆地的像素分布差异最大。
为什么会这样？
你想想看，如果阈值选得好，水体像素都集中在低值区，陆地像素都集中在高值区，两类之间的方差就很大。如果阈值选得不好，两类像素混在一起，方差就小。
Otsu算法就是遍历所有可能的阈值，找到那个让类间方差最大的值。
Otsu算法步骤：
- 计算影像的直方图（每个灰度级的像素个数）
- 遍历所有可能的阈值 t（从0到255）
- 对于每个 t，计算：
- 水体类（灰度 < t）的像素比例 w0 和平均灰度 μ0
- 陆地类（灰度 ≥ t）的像素比例 w1 和平均灰度 μ1
- 类间方差 σ² = w0 × w1 × (μ0 - μ1)²
- 选择使 σ² 最大的 t 作为最优阈值
我记得第一次用Otsu算法时，心里还有点打鼓。结果跑出来一看，阈值选得比我自己调的还准。从那以后，我就成了Otsu的忠实用户。
4.4 代码实现：从理论到实战
光说不练假把式。我们来看看怎么用Python实现Otsu算法。
import numpy as np
import cv2
from matplotlib import pyplot as plt
# 读取近红外波段影像（假设是单波段灰度图）
nir_image = cv2.imread('landsat_band5.tif', cv2.IMREAD_GRAYSCALE)
# 方法一：使用OpenCV内置的Otsu算法
# 注意：cv2.threshold要求输入是8位图像
_, otsu_threshold = cv2.threshold(
nir_image, 0, 255,
cv2.THRESH_BINARY + cv2.THRESH_OTSU
)
print(f'Otsu算法自动选取的阈值: {_}')
# 方法二：手动实现Otsu算法（加深理解）
def otsu_threshold_manual(image):
# 计算直方图
hist, bins = np.histogram(image.flatten(), 256, [0, 256])
# 归一化
hist = hist / hist.sum()
# 初始化
best_threshold = 0
best_variance = 0
# 遍历所有可能的阈值
for t in range(1, 255):
# 水体类（低值区）
w0 = hist[:t].sum()
if w0 == 0:
continue
mu0 = np.sum(np.arange(t) * hist[:t]) / w0
# 陆地类（高值区）
w1 = hist[t:].sum()
if w1 == 0:
continue
mu1 = np.sum(np.arange(t, 256) * hist[t:]) / w1
# 计算类间方差
variance = w0 * w1 * (mu0 - mu1) ** 2
# 更新最优阈值
if variance > best_variance:
best_variance = variance
best_threshold = t
return best_threshold
manual_threshold = otsu_threshold_manual(nir_image)
print(f'手动实现的Otsu阈值: {manual_threshold}')
# 应用阈值
water_mask = nir_image < manual_threshold
# 可视化
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.imshow(nir_image, cmap='gray')
plt.title('近红外波段原始影像')
plt.axis('off')
plt.subplot(1, 3, 2)
plt.hist(nir_image.flatten(), bins=256, range=[0, 256])
plt.axvline(manual_threshold, color='r', linestyle='--', label=f'阈值={manual_threshold}')
plt.title('直方图与Otsu阈值')
plt.legend()
plt.subplot(1, 3, 3)
plt.imshow(water_mask, cmap='Blues')
plt.title('水体提取结果')
plt.axis('off')
plt.tight_layout()
plt.show()
避坑指南：
我曾经在山区项目里直接用Otsu算法提取水体，结果山体阴影全被识别成了水。为什么？因为阴影在近红外波段也很暗。后来我加了一个坡度掩膜，把坡度大于5度的区域排除掉，效果才好了。
所以记住：Otsu算法不是万能的，它只考虑灰度分布，不考虑空间上下文。复杂场景下需要结合其他信息。
4.5 单波段阈值法的局限性
说了这么多优点，也得讲讲缺点。单波段阈值法有几个硬伤：
| 问题 | 原因 | 解决办法 |
|---|---|---|
| 山体阴影误判 | 阴影在NIR波段反射率低，类似水体 | 结合DEM坡度数据 |
| 浑浊水体漏提 | 悬浮物增加NIR反射率 | 改用绿波段或SWIR波段 |
| 薄云干扰 | 云在NIR波段反射率高 | 先做云掩膜 |
| 阈值不通用 | 不同卫星、不同季节阈值不同 | 使用Otsu自适应阈值 |
嗯，这里要注意：单波段阈值法最适合的场景是：大面积、清洁、深水体的提取。比如大型湖泊、水库、海洋。如果是小河沟、浅滩、或者富营养化严重的水体，建议用多波段方法。
4.6 本章知识体系
为了让你更直观地理解本章的知识结构，我画了一张流程图：
总结一下：
单波段阈值法，说白了就是利用水体在近红外波段的强吸收特性，用一个阈值把水陆分开。手动选阈值靠经验，Otsu算法靠数学。我个人建议：日常用Otsu，关键项目手动微调。
记住，没有万能的方法。单波段阈值法虽然简单，但在很多场景下依然是最优解。尤其是当你需要快速出图、或者数据质量一般的时候，它往往比那些花里胡哨的算法更靠谱。
公众号:蓝海资料掘金营，微信deep3321