简介面向遥感与图像处理初学者该专题提取资料聚焦全色影像中的水体识别演示了以阈值分割为核心从单波段影像到二值水体图的完整流程先通过大津法等方法自动选取阈值区分水体与非水体再结合后处理优化提取边缘。压缩包共2个文件包含可运行的Matlab脚本和一张结果图算法脚本负责完整处理流程结果图便于直观对照整体仅2.18MB轻量易用便于直接阅读代码逻辑并对照输出验证效果。已有1270人学习下载适用于环境监测、城市规划等需要快速初步水体识别的场景。脚本覆盖预处理、阈值设定、二值化分割、后处理与结果展示等关键环节既能帮助入门者理解阈值分割原理也为后续改用动态阈值或集成区域生长、机器学习等方法提供了基础。1. 阈值分割提取水体一张全色影像怎么在半小时内给出可用的水体边界全色影像没有近红外波段算不了 NDWI很多人在这一步就卡住了。但真正跑一遍你会发现阈值分割恰恰是应对单波段影像最朴素也最稳定的一招水体在灰度图上稳定偏暗与裸地、植被、屋顶之间的直方图距离足够大一个合理的阈值就能把水面从背景里切出来。这个专题包里的 water.m 就是完成这件事的最小实现脚本读入一张全色影像经过预处理、自动阈值分割、后处理最后输出 out.bmp 这样的二值图白色是水体黑色是背景。它适合两类人一类是手里只有全色影像、需要快速做水体初筛的工程师另一类是想把二值化完整流程跑通、而不是只看教科书例子的新手。2. 全色影像与阈值分割原理先弄清水体为什么能被单波段切出来2.1 全色影像的灰度特征水体在直方图里的位置全色影像是单波段遥感影像记录的是地物在一段宽波段内的综合反射率。水体光谱特性很特殊纯净水对红光到近红外波段的吸收极强反射率往往只有百分之几换算成 8 位灰度就是 0~70 区间的暗像素。相比之下裸土、水泥地、植被的反射率更高灰度值通常在 100 以上。这种灰度分离是阈值分割能做下去的前提。不过“大部分亮、水域暗”只是理想状态。云阴影、建筑物阴影同样是低反射区灰度能压到 10~60和浅水、混浊水高度重叠。还有一个容易忽视的细节是浑浊水体含沙量高或者富营养化的水面反射率会比清水高出一大截灰度值可能抬升到接近裸地水平。这些情况决定了单纯用固定阈值很容易翻车一定要配合直方图观察和后续处理。我一般在一个工程里拿到全色影像第一件事不是直接跑分割而是先让直方图现原形。[counts, grayLevels] imhist(img); figure; stem(grayLevels, counts, Marker, none); xlabel(灰度值); ylabel(像素数);代码逻辑很简单imhist 输出每个灰度级的像素计数stem 把它画出来。看这张图如果水体面积够大直方图左端会出现一个明显的峰右边是背景的峰两峰之间的谷底就是最佳分割点。如果左端峰和右边峰重合在一起、没有谷就说明这个场景不适合全局阈值分割就该考虑用局部自适应阈值或者换特征了。2.2 固定阈值与 Otsu 自动阈值两条路线怎么选阈值分割的数学表达非常简单对每个像素做一个大小比较。设输入灰度 f(i,j)阈值 T水体是暗目标判断式是 f(i,j) ≤ T满足条件的像素标为前景之后再决定是否取反。实现上有两条路线。固定阈值适合光照均匀、灰度分离明显的小幅影像人工从直方图谷底读出 T 即可。缺点是依赖人参与场景一变就要重新读一次批量处理时很不现实。自动阈值里最常用的是 Otsu 方法它遍历所有可能的灰度级找出使背景与前景类间方差最大的那个阈值整个过程不需要人工干预。MATLAB 里用 graythresh 计算返回的是归一化到 0~1 的 level 值实际应用时要么直接传给 imbinarize要么乘以 255 再跟原图比较。Otsu 的思路一句话就能说透让两类像素的灰度均值离得尽量远、各自内部尽量集中能拉开间距的阈值就是好阈值。Otsu 虽然高效但它隐式假设直方图是双峰分布。如果影像里地物类型很多、直方图三个峰以上Otsu 切出来的阈值往往偏向面积更大的类别小水塘可能被漏掉。这种情况可以退回到手动阈值也可以改用局部自适应阈值它在每个像素的邻域里单独计算阈值抗光照不均的能力强但计算代价成倍上涨。做快速初筛Otsu 仍然是第一选择至少它给了一个不需要拍脑袋的起点后面再根据场景微调。2.3 water.m 代码结构与处理流程拆解water.m 是整个专题包的核心处理流程分七步用一段带注释的代码可以描述为% water.m —— 全色影像阈值分割水体提取 clc; clear; close all; % 步骤1读入影像三通道则转灰度 img imread(input.tif); if ndims(img) 3 img rgb2gray(img); end % 步骤2中值滤波去噪3x3 窗口兼顾边缘保留 img_f medfilt2(img, [3 3]); % 步骤3Otsu 自动阈值 level graythresh(img_f); % 步骤4二值化灰度大于等于阈值为前景 bw imbinarize(img_f, level); % 步骤5水体是暗目标亮前景居多就取反 if mean(bw(:)) 0.5 bw ~bw; end % 步骤6形态学后处理填充内部空洞、删除小面积噪点 se strel(disk, 2); bw imclose(bw, se); bw bwareaopen(bw, 100); % 步骤7写结果 imwrite(bw, out.bmp);这段逻辑里最容易被忽略的是第 5 步。全色影像水体是暗目标而 imbinarize 默认把灰度大于阈值的像素置为前景如果不加判断直接使用得到的就是“非水体”图像。这里用前置面积占比做启发式判断如果前景像素超过一半就取反。注意这个启发式在极端场景比如整幅影像几乎全是水下不成立更稳妥的做法是保留判断的同时人工看一眼样本图再确认。第 6 步是工程里真正的坑点。imclose 用半径 2 的圆盘结构元素先做闭运算把水体内部的细小噪点填掉bwareaopen 再把面积小于 100 像素的孤立小斑块删除。这两个参数跟影像的地面分辨率直接相关10 米分辨率影像里 100 像素对应一公顷1 米分辨率里只对应 100 平方米换数据源后必须重新调否则要么留下大量碎斑要么把真正的小水塘删掉。3. 实操从原始全色影像到 out.bmp 的完整跑通步骤3.1 工程目录组织与运行环境把专题包解压后water.m 和 out.bmp 会出现在同一目录下。out.bmp 是样例输出可以当验收参照判断自己机器跑出来的结果在内容上是否一致。我的工程目录一般组织成下面这个样子project/ ├── water.m ├── input.tif ├── out.bmp └── temp/water.m 默认读入当前目录下的 input.tif所以最省事的做法是把待处理的全色影像复制到脚本目录下重命名为 input.tif然后直接运行。如果影像扩展名不是 tif改一下 imread 的参数即可。这里有个建议尽量用 tiff 格式做输入jpg 属于有损压缩灰度值已经偏离原始反射率提取出来的边界会比 tiff 粗糙不少。运行环境方面脚本用到的 medfilt2、graythresh、imbinarize、strel、bwareaopen 都是图像处理工具箱里的基础函数R2018b 及之后的版本基本都兼容。如果本机没有 MATLAB用 Octave 加 image 包也能跑只需要在脚本开头补一行 pkg load image其余代码不用改。我实际跑 Octave 时只遇到过一个差异imbinarize 函数在 Octave 里需要额外确认返回类型一般加一行 logical() 转换就够了。3.2 跑通主流程参数怎么看、怎么调第一次运行如果报错九成是路径问题。直接用 cd 切换到工程目录再运行 water比把脚本路径写死在代码里省事。脚本正常结束后当前目录下会生成 out.bmp打开后白色区域就是提取出来的水体。拿到结果先别急着交付按下面三件事检查。一是看直方图。脚本里没强制画图可以在命令行补一句。[counts, gray] imhist(img_f); figure; plot(gray, counts);如果直方图左峰和右峰没有明显谷底说明 Otsu 的区分度不够结果里大概率混入了阴影或亮噪声。二是看阈值落在哪个灰度范围。graythresh 返回的 level 是归一化值乘以 255 就能看到实际灰度。比如 0.2 对应 51意思是灰度低于 51 的像素被判为水体这个数字符不符合经验一眼就能看出脚本有没有跑偏。三是叠加原图回看这一步最关键。figure; imshow(img, []); hold on; [B, L] bwboundaries(bw); for k 1:length(B) boundary B{k}; plot(boundary(:,2), boundary(:,1), r, LineWidth, 1); endbwboundaries 提取二值图里所有连通域的边界plot 在原图上红色描边水体轮廓直接叠在影像上哪里圈多了、哪里漏了一目了然。我习惯把它固定成流程标配肉眼判断比任何抽象精度指标都直观。参数调整时最常动的是三个位置列成表方便对照。参数在代码中的位置默认值调整方向中值滤波窗口medfilt2(img, [3 3])3×3噪点多时改用 5×5代价是细河道边缘变钝结构元素半径strel(disk, 2)2水体连贯性好可加大到 3~5河道窄则保持 1~2最小面积阈值bwareaopen(bw, 100)100影像分辨率高且碎斑多时调到 500 以上这三个参数是互相牵制的。滤波窗口开大了边缘变钝面积阈值就需要跟着收敛结构元素半径加大能填掉噪声但窄河道被抹平的风险也同步上升。改参数时每次只动一个对比输出后再决定下一步别同时改两三个不然翻车了都不知道是谁干的这算是调这类脚本的一个基本纪律。3.3 输出验证怎么确认 out.bmp 不是一张废图用脚本输出水体像素占比可以快速判断结果合理性。ratio sum(bw(:)) / numel(bw); fprintf(水体像素占比: %.2f%%\n, ratio * 100);对大多数场景水体占整幅影像的比例通常不超过 20%如果输出结果超过 30%就要排查是不是把大片阴影折算进来了。反过来如果占比不足 1%又可能把浑浊水体或小型坑塘漏掉了。注意这只是一个经验范围极端场景比如海岛湿地会被打破具体到某一帧影像时还是要以目视结果为准。更严谨的验证办法是找一个已知水体区域做局部对照在原图上框出一块确定是水面的区域计算这块区域的灰度均值和标准差再看分割结果在相同坐标上是不是白色前景。这个交叉验证做下来基本能确认阈值选取是否可信。如果区域内有明显的道路或船只还要手动确认这些目标不会把水体割裂成不该有的空洞。4. 避坑指南水体阈值分割最常翻车的四类场景与排查方法4.1 阴影被当成水体圈进去现象out.bmp 里除了水面沿高楼边缘、山地北坡出现连成片的暗色区域形状狭长边界僵直肉眼一看就不是自然水面。原因全色影像里建筑物阴影和山地阴影灰度极低和水体落在同一个灰度区间。Otsu 的全局阈值只看直方图统计不理解空间语义灰度重叠的分区就会被一并切进前景。解决先做事后处理用几何特征把形状异常的区域筛掉。水体连通域通常面积大、形状圆润、周长与面积比低阴影则狭长且边缘整齐。可以用 regionprops 统计每个连通域的长宽比把长宽比过大的连通域整体删除。lb bwlabel(bw); stats regionprops(lb, BoundingBox); keep false(numel(stats), 1); for k 1:numel(stats) w stats(k).BoundingBox(3); h stats(k).BoundingBox(4); ar max(w, h) / min(w, h); if ar 5 keep(k) true; end end bw ismember(lb, find(keep));这里重点是 bwlabel 先给每个连通域编号regionprops 返回每个编号的外接框长宽比 ar 大于 5 的判定为阴影特征ismember 只保留符合条件的连通域。删完之后如果还有大片阴影残留就得回到影像层面尝试用上一步的叠加图确认最大几块暗区的分布位置人工判断它们是不是与建筑物边界重合必要时直接手工掩膜。4.2 水体内部空洞密布现象大湖、水库内部出现大量零星黑点水体连成一片但不完整像筛子一样。原因水面自身存在波浪纹理、太阳耀斑或者含沙量不均局部灰度波动越过阈值线导致单个像素或小簇像素被划入背景。解决加大闭运算结构元素半径或者直接对水体区域做孔洞填充。bw imfill(bw, holes);imfill 会填充二值图里所有被前景包围的背景区域对细小孔洞效果极佳。但注意如果水体本身有缺口跟外界连通imfill 填不了必须先做一步闭运算把缺口补上再调用 imfill。我处理这种情况时习惯先闭后填两个操作连着做效果比单纯加大半径稳定。4.3 16 位影像跑出全黑或全白现象输入是 16 位 tiff脚本跑完输出的 out.bmp 全黑或者 imshow 显示时白茫茫一片。原因graythresh 对 uint16 数据处理没有问题但 imbinarize 的结果与显示路径不匹配。16 位影像灰度范围是 0~65535直接用 8 位 BMP 输出时数据被截断暗部细节丢失。还有一种情况是影像里有极少数特别亮的像素拉高了整体显示范围水体区域在画面上看起来接近全黑。解决把原图先线性拉伸到 0~255 再送入分割流程。img double(img); img (img - min(img(:))) / (max(img(:)) - min(img(:))); img uint8(img * 255);这段归一化会把对比度拉开让阈值分割面对的是一个灰度范围完整、分布平稳的影像。之后再跑 graythresh 和 imbinarize结果会稳定得多。注意拉伸之后的灰度直方图形状会改变原图的相对关系还在但绝对灰度值已经不是原始反射率后期如果做定量分析要在报告里注明这一步。4.4 窄河道断裂成虚线现象一条两三米宽的河道在 out.bmp 里断成几段河中间的窄条被分成白色和黑色交替像虚线一样。原因中值滤波的 3×3 窗口在窄河道上相当于一个覆盖多个像素的平滑核灰度值被拉向相邻背景河道中央的像素低于阈值自然就被划成背景。形态学的闭运算半径如果大于河道宽度同样会把河道直接抹掉。解决对细小河道用两种手段组合。一是把中值滤波改为沿河道方向的定向滤波或者干脆跳过中值滤波用原始影像做阈值分割噪声交给后面的面积阈值过滤二是把闭运算半径缩小到 1避免把窄条磨平。如果断裂情况不严重也可以用 bwmorph 对前景做 bridge 运算把断裂处连接起来。bw bwmorph(bw, bridge, inf);bridge 操作会在相邻两个前景像素之间补上中间点专门用于修复细线断裂。不过它有副作用会让图像产生轻微膨胀河岸线条粗一圈用完之后最好再做一次腐蚀抵消否则边界位置会往外偏几个像素。5. 进阶验证与交付把 out.bmp 从能看变成能用一个阈值分割结果真正能交付的前提是对得上一张叠加图。拿到 out.bmp 我会先跑一遍连通域统计按面积排序把最大的几个连通域在原图上标出来。这样做的目的很明确大水体在影像里通常对应最大的暗色连通域如果最大的连通域落在农田或山体阴影上整个结果就不成立不用再花时间做后续精度分析。stats regionprops(bw, Area, BoundingBox); [~, idx] sort([stats.Area], descend); top5 idx(1:min(5, numel(idx))); figure; imshow(img, []); hold on; for k 1:numel(top5) bb stats(top5(k)).BoundingBox; rectangle(Position, bb, EdgeColor, r, LineWidth, 1.5); end这段代码把面积最大的五个连通域外接框画在原图上一眼就能看出最大的水体区域在哪里形状是不是合理。接着把水体像素占比、连通域数量、最大连通域面积这三个统计量打成一串文本跟 out.bmp 一起存进工程记录。换参数重跑时对比这些数字立刻能看出哪次分割更干净比自己肉眼重看快得多。还有一个细节对交付特别关键别把 out.bmp 直接丢给对方就完事。二值图只能说明哪些像素是水没有地理参考和边界光滑度别人没法直接判断可信度。把叠加图、统计文本和 out.bmp 三者放在一起输出基本就能达到可复核的最低标准。矢量化的需求要用到 bwboundaries 提取坐标再导出这一步留给 GIS 端做更合适。第一次把这份脚本接进工作流时我交出去的结果被上游反手退回截图上一栋楼的阴影被圈成了水面。那次翻车让我意识到跑脚本只是开始真正决定交付质量的是跑完之后那十五分钟的检查和统计。从那以后我每次跑完分割都会强制走一遍叠加检查和连通域统计这个习惯帮我拦下了好几次低质量输出。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取