简介：面向遥感与图像处理初学者的阈值分割水体提取示例包，基于Matlab实现，解决从全色影像中快速识别水体区域的核心问题，适用于GIS分析、环境监测等场景的初步水体检测。压缩包内含2个文件，共2.18MB，其中water.m为Matlab脚本，展示从影像预处理、阈值设定到二值化分割的完整流程；out.bmp为处理后的输出结果，可直观对照算法效果。当前已有1270人学习下载。通过该资源可掌握阈值分割的基本原理、Otsu等动态阈值选取思路，并了解连通域分析、噪声去除等后处理手段，也可作为进一步结合区域生长或机器学习方法提升精度的起点。资料体量小巧，适合课程设计、毕业设计或快速原型验证时参考使用。

1. 全色影像的水体提取：为什么阈值分割是最快的起手式

全色影像的水体提取，落到实操上最先能跑通的思路就是阈值分割——单波段影像没有色彩可用，水体与地物的灰度差异往往比纹理特征更直接。我以前处理过不少灰度遥感影像，水面在晴空下通常偏暗，而道路、裸土和建筑亮度偏高，一个合理的阈值就能把大部分水面圈出来，再补几道后处理就能出成果图。这份专题把从影像读入、直方图分析、阈值计算、二值化到后处理输出的完整链路整理成了一个可复用的 Matlab 流程，包里的 water.m 脚本把整套逻辑串了起来，直接改路径就能跑，出来就是 out.bmp 这样的二值结果图。

它解决的不是“最先进的水体提取算法”问题，而是“在没有额外波段、没有训练样本的情况下，如何快速、可解释地把水面位置标出来”。适合三类人：做 GIS 底图快速预处理的从业者、课程设计里需要跑通几种分割算法的学生、以及刚接触遥感影像想建立完整流程概念的开发者。接下来我从原理讲到代码，再把最常见的翻车点逐条拆开。

2. 阈值分割与全色影像的特性匹配：先理解灰度差异，再动手调阈值

2.1 单波段影像里水体与非水体的灰度差异从哪来

全色影像只有一个波段，记录的是一段较宽波长范围内地表反射太阳辐射的总强度。水体对入射光的吸收能力明显强于大多数陆地地物，特别是当水体比较深、含沙量不高的时候，进入水体的光大部分被吸收掉，能返回传感器的能量很少，所以水体在灰度影像上通常表现为暗色调。而道路、裸土、建筑屋顶、干燥植被反射率高，呈现为亮色调。这个“水暗、地亮”的总体规律，是阈值分割能成立的物理基础。

但事情不会总有这么理想。水面不是完美的朗伯体，太阳角度合适时，水面会出现镜面反射，把阳光直接反射进传感器，这时水面的亮度可能比周围裸地还高。含沙量高的浑浊水体、水华暴发的水面、浅水区露出水底的区域，灰度都会往亮的方向偏移。也就是说，“水体偏暗”只能作为先验知识，不能当作绝对定律。实际操作中我会先用直方图确认当前影像里水体到底落在哪个灰度区间，再决定是用“小于阈值”还是“大于阈值”来切分。

还有一个容易被忽略的现象：山体阴影、建筑物阴影在灰度上和水体非常接近，都属于暗色调目标。阈值分割本身不区分“暗水体”和“暗阴影”，它只按灰度切，所以阴影误提是这个方法最典型的误差来源。这一点后面在避坑章节里会专门展开。

2.2 固定阈值、Otsu 自适应阈值与局部阈值的适用边界

阈值分割的思想是把像素按灰度值划分成两类：低于阈值的归为一类，高于阈值的归为另一类。最朴素的做法是固定阈值，比如肉眼看到直方图谷底在 0.4，就写死level = 0.4。这种方法好在快、可控、结果可复现，坏在换一张影像就得重新调一次，光照条件不同时同一套参数可能完全失效。

Otsu 方法解决的是“阈值能不能自动算”的问题。它的基本思路是遍历所有可能的阈值，计算用该阈值分割后两类像素的类内方差和类间方差，目标是让类间方差最大、类内方差最小。类间方差越大，说明两类像素被分开得越彻底。对应到 Matlab 图像处理工具箱，函数是graythresh，它返回一个归一化到 [0,1] 区间的阈值，可以直接传给imbinarize使用。很多人说起 Otsu 就默认它一定比固定阈值好，实际并非如此。Otsu 假设灰度直方图大致呈现双峰分布，两个峰分别对应目标与背景。如果影像里地物类型复杂，直方图是多峰的，Otsu 往往会把两个相近的地物类切到一起，反而得到一个很奇怪的阈值。

局部阈值适用于光照不均匀的影像。比如一幅大影像左侧是阴影山区、右侧是开阔平原，全局统计的直方图会把不同区域的灰度分布叠在一起，无论是固定阈值还是 Otsu 都很难同时照顾两个区域的差异。常见做法是把影像切成若干小块，对每个小块分别计算阈值，再拼接结果。这个思路我在第 6 章会给出一个可直接改的脚本版本。

2.3 为什么用 Matlab 脚本做第一步验证最顺手

水体提取这个任务，第一步要的不是高并发、高性能工程实现，而是快速验证“当前影像里的水体用灰度能不能分出来”。Matlab 的优势在于图像读写、直方图可视化、形态学操作这些函数开箱即用，Debug 时能看到中间每一个变量的大小和取值，不像编译型工程那样需要在代码里打一堆日志才能定位问题。

water.m 这种脚本的典型形态是顺序执行：读图、转灰度、归一化、看直方图、算阈值、二值化、做后处理、存图。每个步骤的输入输出都是矩阵或者图像，变量名一目了然。对于遥感方向的研究者来说，这比用 OpenCV 写一整套代码再调参要快得多。即使最终要部署到生产系统，先用 Matlab 验证算法流程、确定参数范围，再移植到其他语言，也是最稳妥的路径。

3. 基础流程落地：预处理、直方图分析与分割脚本实战

3.1 去噪与灰度拉伸：让水体和非水体更容易分开

原始影像直接做阈值分割不是不行，但噪点会直接影响阈值计算。全色影像里的噪声来源包括传感器本身的暗电流噪声、压缩编码带来的块效应、以及大气散射引起的均匀亮度抬升。我的习惯是先做两步：第一步用中值滤波或高斯滤波去噪，第二步做灰度拉伸增强对比度。中值滤波对去除孤立噪点很有效，适合保持边缘；高斯滤波对高斯噪声更对症，但边缘会轻微变糊。水体提取场景下我一般优先用中值滤波，因为水体边缘本身就是渐变过渡的，高斯滤波容易让岸线边界更模糊。

% 步骤1：读取影像，统一为double类型的灰度矩阵 I_original = imread('panchromatic_image.tif'); % 读取全色影像 if size(I_original, 3) == 3 % 防止输入是伪彩色的意外情况 I_original = rgb2gray(I_original); end I = im2double(I_original); % 归一化到[0,1] % 步骤2：中值滤波去噪，3x3窗口适合中等分辨率影像 I_denoise = medfilt2(I, [3 3]); % 步骤3：灰度拉伸，把原本较窄的动态范围拉宽 I_stretched = imadjust(I_denoise, stretchlim(I_denoise, [0.02 0.98]), []);

代码背后的逻辑是先用im2double把像素值统一到 0 到 1 之间，这样后续提到的阈值、对比度参数都归一化到一个统一的尺度，换影像时更容易横向比较。medfilt2的第二个参数是滤波窗口，窗口越大去噪越强，但边缘细节丢失越严重，3×3 是起步值，影像分辨率高时可以试 5×5。imadjust配合stretchlim做的是自动对比度拉伸，[0.02 0.98]表示忽略灰度直方图两端各 2% 的极值后再拉伸，这样可以避免个别高光点把整个灰度范围撑开。

3.2 直方图可视化：先肉眼看双峰，再谈自动阈值

程序自动化的前提是人对当前数据有判断。我会在算阈值之前强制自己看一眼灰度直方图，这一步非常重要。直方图横轴是灰度值，纵轴是像素数量。如果水体偏暗、背景偏亮，直方图上会看到两个峰：暗区的峰对应水体与阴影，亮区的峰对应陆地地物，两个峰之间的谷底就是阈值的最佳候选位置。

% 步骤4：绘制灰度直方图，观察分布形态 figure; histogram(I_stretched, 256); xlabel('灰度值'); ylabel('像素数量'); title('预处理后的灰度直方图'); grid on; % 顺带输出灰度分位数，帮助判断阈值大致范围 pcts = prctile(I_stretched(:), [10 25 50 75 90]); fprintf('灰度分位数 10%%: %.3f, 25%%: %.3f, 50%%: %.3f, 75%%: %.3f, 90%%: %.3f\n', pcts);

这一步不是走形式。直方图如果呈现出明显的双峰形态，说明阈值分割在这个场景下大概率可行；如果直方图只有一个大峰加一个缓坡，说明水体面积占比太小或者水体与背景灰度严重重叠，这时候强行做全局阈值分割，结果基本靠猜。看直方图的另一层用意是判断水体到底在亮的还是暗的一侧。大部分水体是暗的，但镜面反射严重的水面会形成高亮峰，不看图直接写IBW = I < level就会把整个水面漏掉。

3.3 用 graythresh 自动算阈值，并生成初始分割图

如果直方图双峰清晰，可以直接用graythresh计算 Otsu 阈值。返回的阈值是归一化数值，比如 0.42 代表原图中灰度值 107（255 × 0.42）作为分界点。得到阈值后要做的不是急着存图，而是先确认分割方向：水体暗就保留小于阈值的像素，水体亮就保留大于阈值的像素。这个方向搞反了，输出的结果图就是陆地而不是水体。

% 步骤5：Otsu自动计算阈值并生成初始二值图 level = graythresh(I_stretched); % 返回[0,1]之间的阈值 fprintf('Otsu自动阈值: %.3f\n', level); % 默认假设水体比背景暗，所以保留低于阈值的区域 BW_initial = I_stretched < level; % 如果直方图显示水体在高亮侧，则改用大于号 % BW_initial = I_stretched > level;

代码里level是唯一的关键参数。graythresh的输入可以是二维灰度矩阵，输出是归一化数值，不需要人为指定灰度分级数，这是它的方便之处。BW_initial是一个逻辑矩阵，值为 1 的像素代表被判定为水体的区域。这里要特别提醒：Otsu 算出来的阈值不是金标准。当水体面积相对于整幅影像非常小时，背景像素占绝对多数，Otsu 会把阈值推向背景那一侧，导致水体区域被严重缩小。遇到这种情况，我会在 Otsu 阈值基础上人工偏移，比如把阈值改成level * 0.9来保留更多暗部区域，具体偏移量要根据直方图谷底位置修正。

4. 后处理与结果输出：从二值图到一份能用的 out.bmp

4.1 小噪斑清理与孔洞填充的关键参数

初始二值图上通常有两种噪声：一种是零散分布在陆地区域的孤立白点，可能是亮色地物的小反射区被误判为水体；另一种是水面内部的小黑点，可能是波纹、悬浮物造成的灰度异常。前者用面积过滤，后者用孔洞填充。

% 步骤6：删除面积小于50像素的孤立噪斑 minArea = 50; BW_area = bwareaopen(BW_initial, minArea); % 步骤7：填充水体内部的细小孔洞 BW_filled = imfill(BW_area, 'holes');

bwareaopen的作用是把连通域面积小于minArea的白色区域全部置为 0。这个值不是固定的——影像分辨率越高，一个真实水体斑块的面积越大，minArea可以设置得越大而不担心误删；如果影像里存在细长的沟渠这类真实水体，minArea设置太大反而会把这些线性水体整段删掉。我的经验是先设 50 看效果，统计一下图里连通域面积的分布再调。imfill的'holes'模式只填充白色区域内部的黑色孔洞，不会扩展到外部背景，所以不用担心把河岸两侧的陆地填进来。

4.2 形态学操作把水体边缘修圆滑

分割边缘往往带有锯齿或碎刺，这在后期出图时非常难看，也会让面积统计偏离真实值。形态学闭运算是水体提取后处理里的标准动作：先膨胀再腐蚀，它的特点是能弥合狭窄的断裂、填平边缘小凹陷，同时大体保持原区域的面积大小。开运算则是先腐蚀再膨胀，适合去掉零散的小突起。针对水面区域，我一般优先用闭运算。

% 步骤8：用半径为3的圆形结构元素做闭运算 se = strel('disk', 3); BW_smooth = imclose(BW_filled, se); % 再次清理闭运算可能产生的小连接（用一次开运算兜底） BW_final = bwareaopen(BW_smooth, minArea);

strel('disk', 3)生成一个半径为 3 像素的圆形结构元素，半径决定边缘平滑的程度。半径太小，起不到弥合缝隙的作用；半径太大，会把两条距离较近的岸线之间的陆地误填成水体，导致两条河汊被并成一条。3 是一个相对保守的起步值，对大多数中等分辨率影像都适用。闭运算之后我没有直接结束，而是再做一次面积过滤，这个操作防止闭运算把原本不相连的噪声斑块桥接成大面积区域。

4.3 输出结果与效果评估：怎么判断提取是靠谱的

out.bmp是最终保存的二值结果图。Matlab 里逻辑矩阵保存成 BMP 时，1 和 0 会被写成 255 和 0，显示效果就是白色水体、黑色背景，这正好符合遥感制图里目标高亮显示的习惯。

% 步骤9：保存结果，同时输出统计信息 imwrite(BW_final, 'out.bmp'); % 统计水体区域面积与斑块数量 stats = regionprops(BW_final, 'Area', 'Perimeter', 'BoundingBox'); areas = [stats.Area]; fprintf('检测到水体斑块 %d 个\n', numel(stats)); fprintf('水体总面积: %d 像素，最大斑块面积: %d 像素\n', sum(areas), max(areas));

保存之后不能直接收工。我会把BW_final以半透明方式叠加到原图上：用imshowpair对比原图和二值图，或者用hold on把水体边界勾勒在原图上，这样可以直观看到每个水体斑块与真实岸线的重合程度。一个比较实用的经验是观察提取面积占影像总面积的比例，如果比例明显高于正常水域占比，比如一幅以山地为主的影像提取出 30% 的面积是“水体”，那基本可以断定是把阴影整片划进来了。

5. 避坑与排查：水体提取最常见的四个翻车现场

5.1 阴影被当成水体：灰度相近不是充分的判据

现象：分割结果里出现大片不规则暗色区域，明显是山体阴影或者建筑物阴影，而非水体。在丘陵城市影像里尤其严重，水体面积统计被放大数倍。

原因：水体与阴影在灰度上高度重叠，简单的阈值分割只能看到灰度数值，无法区分“低反射的平坦水面”和“低反射的遮挡区域”。阴影区域没有纹理，但水体表面其实也存在微弱的灰度波动，只是阈值分割没用到这个信息。

解决：最直接的办法是在后处理里增加位置约束。阴影通常与高大地物相邻，区域形状呈条带状，而水体边缘相对平滑、形状更自然。用regionprops输出每个连通域的Eccentricity和Solidity，把细长且形状不规则的暗色区域筛掉。另外，如果影像覆盖范围内有数字高程数据，直接按坡向和太阳方位角生成阴影掩膜来排除，这个方案最稳健。

5.2 亮水面被漏提：镜面反射让水体亮度反超陆地

现象：分割结果里水体区域出现大面积空洞，只在边缘提取出一圈轮廓，中心区域被判定为陆地。特别容易发生在太阳高度角低、水面平静的影像上。

原因：平静水面在特定角度下产生镜面反射，反射光的强度远大于漫反射的陆地地物。水体区域直方图出现高亮峰，而全局阈值又落在暗部谷底，亮水面自然被划进“非水体”一类。

解决：先把影像做一次方向性滤波或者局部对比度增强，把水面的高光区域与周边水体在灰度上拉近。更直接的做法是检查直方图：如果水体区域内同时存在暗部和亮部两个峰，全局阈值无法一次覆盖，可以考虑对分割结果做“暗水体 ∪ 亮水体”的并集操作，即分别提取I < level_dark和I > level_bright的区域，再做并集和孔洞填充。

5.3 自动阈值与人工判断冲突：Otsu 的阈值是不是一定靠谱

现象：graythresh算出的阈值明明让分割结果走向失控，例如把大范围暗色裸土、阴影全部切进来，肉眼可见地不合理。

原因：Otsu 是纯统计方法，它假设直方图双峰且两类像素占比相对均衡。当水体面积极小时，暗部峰太小，Otsu 计算的类间方差最大值点会偏向亮部峰一侧，导致阈值偏低，把大量中间灰度地物并入暗部。这时不是代码问题，而是方法前提不满足。

解决：在graythresh输出基础上增加一个人工修正参数。定义一个偏移量offset，最终阈值取threshold = level + offset。offset为正时提高阈值、提取更多暗部区域，为负时降低阈值、收紧水体范围。这个参数不强求数学推导，根据直方图谷底的位置微调就行。另外可以用graythresh的第二个输出或者直接看level是否落在直方图谷底附近来辅助判断。

5.4 全局阈值在大幅影像上失效：光照不均导致阈值漂移

现象：同一幅影像里，一部分水域被正确提取，另一部分水域在同一阈值下完全断裂或者根本没提取出来。影像拼接痕迹明显时尤其常见。

原因：大幅影像往往由多景数据拼接而成，各景的成像时间、太阳角度、大气条件不同，灰度基准不一致。全影像统计出一个全局阈值，对亮的那一景合适，对暗的那一景就偏了。

解决：分块处理是工程上的常见做法。把影像切成 256×256 或 512×512 的块，每个块独立计算阈值，再拼接结果。分块后每块内的光照条件近似均匀，Otsu 的假设在每个块内更可能成立。我在下一章给出这个方案的可运行版本。

6. 进阶：分块自适应阈值与边缘约束的快速验证

分块自适应阈值是全局阈值失效时的首选升级方案。思路很直白：把影像切块，每块计算自己的阈值，再拼接回完整二值图。与全局阈值相比，它能让每一块都用更适合局部统计特征的灰度分界点。

% 分块自适应阈值处理 blockSize = 256; [rows, cols] = size(I_stretched); nRows = ceil(rows / blockSize); nCols = ceil(cols / blockSize); BW_block = false(rows, cols); % 初始化输出矩阵 for i = 1:nRows for j = 1:nCols rStart = (i - 1) * blockSize + 1; rEnd = min(i * blockSize, rows); cStart = (j - 1) * blockSize + 1; cEnd = min(j * blockSize, cols); block = I_stretched(rStart:rEnd, cStart:cEnd); % 块内面积太小或灰度方差太小时，沿用全局阈值 if numel(block) < 100 || std2(block) < 0.05 blockLevel = level; else blockLevel = graythresh(block); end BW_block(rStart:rEnd, cStart:cEnd) = block < blockLevel; end end

分块大小blockSize是这里的关键参数。块太大时，局部光照优势不明显；块太小（比如 64×64）时，每块内水体占比可能过低，Otsu 在单块内失效。256 是一个比较稳的起步值。代码里增加了一个保护条件：块内灰度标准差小于 0.05 时说明这块区域几乎是一种地物，没必要算阈值，直接用全局阈值兜底。

分块处理有一个明显的副作用——块与块之间的阈值不连续，结果图上会出现块状边界。为了消除斑块痕迹，我会在拼接后的二值图上再跑一次 3×3 中值滤波，或者用形态学闭运算做平滑。

除了分块，边缘约束对阈值分割的改进也很明显。水体边界往往是平滑的闭合曲线，而阴影边界通常与地形轮廓伴生。可以先对原图用edge(I, 'canny')提取边缘信息，然后检查每个水体斑块与边缘图像的重合程度，把那些内部边缘过于密集的暗色区域判定为阴影。这个思路本质上是给阈值分割加上一个纹理过滤条件。

验证方法上，如果没有参考水体标签图，我的习惯是把分割结果叠加到原图上逐块目视检查。一个像素数占总面积比例、斑块数量、最大斑块面积这些统计量，也能快速判断结果是否可信。

从那次处理大幅拼接影像被阈值漂移坑过一次之后，我每次拿到全色影像都强制先走一遍直方图检查，再决定用全局阈值还是分块阈值。这个习惯把阈值分割从玄学变成了一个可解释、可复现的工程流程。希望帮到你。

本文还有配套的精品资源，点击获取