基于深度学习的多光谱图像语义分割 —— 从原理到实践（第四节）

经过前三节的数据准备与分块处理，我们现在手握预训练的 U‑Net 模型（或自己训练的模型），终于可以对测试图像进行真正的语义分割了。
本节将完整演示测试图像推理、后处理去噪、结果可视化以及植被覆盖率计算，带你走完实际应用的全链路。

---

4. 测试图像推理、后处理与植被覆盖率计算

4.1 测试图像同样需要分块推理

测试图像 test_data 的尺寸为 12446×7654×7（宽×高×通道），整张图像有近 1 亿像素，加上 6 个通道，数据量超过 1.3 GB。直接加载会耗尽内存，更别提送入 GPU。因此推理时也必须采用分块策略。

不过推理时我们可以使用更大的补片（例如 1024×1024），因为推理不涉及反向传播，显存占用比训练小得多，大补片能减少块数，加快速度。

我们只取前 6 个光谱通道作为输入（第 7 通道是掩膜，不用于预测），并创建分块图像对象：

patchSize = [1024 1024];
bimTest = blockedImage(test_data(:,:,1:6), BlockSize=patchSize);

matlab
bimSeg = apply(bimTest, @(bs) semanticseg(bs.Data, net, Outputtype="uint8"), ...
    PadPartialBlocks=true, PadMethod=0);

4.2 使用 semanticseg 对每个块进行分割

semanticseg 是 Computer Vision Toolbox 提供的深度学习分割函数，它接受图像块和训练好的网络，输出每个像素的类别索引（1 或 2）。我们利用 apply 函数将这个操作应用到分块图像的每一个块上：

参数解释：

@(bs)semanticseg(bs.Data, net, Outputtype=“uint8”)：匿名函数，对每个块数据 bs.Data 调用 semanticseg，输出类型设为 uint8（节省内存）。

PadPartialBlocks=true：当图像边缘不足一个完整块时，自动填充（填充值由 PadMethod 指定）。

PadMethod=0：填充值为 0（对类别索引无影响，因为我们的类别是 1 和 2）。

apply 返回一个新的 blockedImage 对象 bimSeg，它记录了每个块的分割结果，但数据仍分散在磁盘或内存中。

4.3 将所有块拼接成完整图像

使用 gather 函数从分块图像中读取所有块并拼接成一整张图像：

matlab
segmentedImage = gather(bimSeg);

现在 segmentedImage 是一个 12446×7654 的二维矩阵（单通道标签图），每个像素值为 1（非植被）或 2（植被）。

4.4 应用掩膜——只保留有效区域

测试图像的第 7 通道是掩膜，标记了哪些像素属于有效分析区域（值为 1），哪些是图像边界或无效区域（值为 0）。我们需要将无效区域的预测结果屏蔽掉，以免干扰后续统计。

matlab
maskTest = test_data(:,:,7);   % 提取掩膜
segmentedImage = segmentedImage .* uint8(maskTest ~= 0);

这样，掩膜为 0 的位置在 segmentedImage 中变为 0，后续统计时这些像素将被忽略（nnz 不会计数 0）。

4.5 观察原始分割结果（含噪声）

先显示一下未去噪的分割图，看看效果：

matlab
imshow(segmentedImage, []);
title("Segmented Image (Before Denoising)");

你会看到分割图中散布着许多孤立的小点、细碎区域，这就是典型的椒盐噪声。噪声产生的原因通常是模型在局部上下文不一致时做出了随机性误判。

4.6 中值滤波去除椒盐噪声

中值滤波是一种非线性滤波，它将每个像素的值替换为其邻域窗口内所有像素的中值。这种操作对孤立噪声点非常有效，同时能较好地保留边缘。

我们采用 7×7 的窗口（窗口大小应大于噪声颗粒尺寸，但不宜过大以免模糊真实边界）：

matlab
segmentedImage = medfilt2(segmentedImage, [7 7]);
imshow(segmentedImage, []);
title("Segmented Image with Noise Removed");

滤波后，分割图变得平滑干净，区域边界更加清晰。

4.7 将分割结果叠加到 RGB 图像上可视化

为了更直观地评估分割质量，我们将去噪后的标签图叠加在测试图像的彩色（RGB）版本上。rgbTest 我们在第二节已经生成（经过直方图均衡化）。

matlab
cmap = [1 0 1; 0 1 0];   % 第一行：品红（非植被），第二行：绿色（植被）
B = labeloverlay(rgbTest, segmentedImage, Transparency=0.8, Colormap=cmap);
imshow(B, cmap);
title("Labeled Segmented Image (Overlay on RGB)");
numClasses = 2;
ticks = 1/(numClasses*2):1/numClasses:1;
colorbar(TickLabels=cellstr(classNames), Ticks=ticks, TickLength=0, TickLabelInterpreter="none");

现在你可以看到，植被区域被绿色覆盖，非植被区域被品红覆盖，叠加在原图上，便于肉眼验证模型的识别能力。

4.8 计算植被覆盖率——解决生态学实际问题

语义分割不只是为了好看，其最终目的是提取有价值的信息。例如，环保机构常需要知道某区域的植被覆盖百分比，用于生态评估或变化监测。

我们统计分割结果中所有属于植被的像素数，再除以有效区域的总像素数（即 ROI 大小）：

matlab
vegetationPixels = ismember(segmentedImage(:), vegetationClassIDs);
numVegetationPixels = sum(vegetationPixels(:));
numROIPixels = nnz(segmentedImage);   % 有效区域内被预测为 1 或 2 的像素总数
percentVegetationCover = (numVegetationPixels / numROIPixels) * 100;
disp("The percentage of vegetation cover is " + percentVegetationCover + "%");

输出示例：

text
The percentage of vegetation cover is 60.7028%

这说明该测试区域约有 60.7% 的面积被植被覆盖，与 Hamlin Beach 州立公园的真实情况吻合。

4.9 本节完整代码汇总

以下是本节所有代码的完整顺序（从推理到计算覆盖率

matlab
% 1. 分块推理
patchSize = [1024 1024];
bimTest = blockedImage(test_data(:,:,1:6), BlockSize=patchSize);
bimSeg = apply(bimTest, @(bs) semanticseg(bs.Data, net, Outputtype="uint8"), ...
    PadPartialBlocks=true, PadMethod=0);
segmentedImage = gather(bimSeg);

% 2. 应用掩膜
maskTest = test_data(:,:,7);
segmentedImage = segmentedImage .* uint8(maskTest ~= 0);

% 3. 显示原始分割
figure;
imshow(segmentedImage, []);
title("Segmented Image (Before Denoising)");

% 4. 中值滤波去噪
segmentedImage = medfilt2(segmentedImage, [7 7]);
figure;
imshow(segmentedImage, []);
title("Segmented Image with Noise Removed");

% 5. 叠加到RGB可视化
rgbTest = histeq(test_data(:,:,[3 2 1]));  % 若尚未生成
cmap = [1 0 1; 0 1 0];
B = labeloverlay(rgbTest, segmentedImage, Transparency=0.8, Colormap=cmap);
figure;
imshow(B, cmap);
title("Labeled Segmented Image (Overlay on RGB)");
numClasses = 2;
ticks = 1/(numClasses*2):1/numClasses:1;
colorbar(TickLabels=cellstr(classNames), Ticks=ticks, TickLength=0, TickLabelInterpreter="none");

% 6. 计算植被覆盖率
vegetationPixels = ismember(segmentedImage(:), vegetationClassIDs);
numVegetationPixels = sum(vegetationPixels(:));
numROIPixels = nnz(segmentedImage);
percentVegetationCover = (numVegetationPixels / numROIPixels) * 100;
disp("The percentage of vegetation cover is " + percentVegetationCover + "%");

4.10 本节小结

在本节中，我们完成了以下关键任务：

对测试图像进行分块推理，利用预训练的 U‑Net 生成标签图；

拼接所有块，得到完整的预测结果；

使用掩膜屏蔽无效区域；

应用中值滤波去除椒盐噪声，提升分割质量；

将分割结果叠加在 RGB 图像上，实现可视化验证；

基于分割图计算植被覆盖率，展示深度学习在生态监测中的实际价值。

4.11 下一步：训练自己的 U‑Net

到目前为止，我们一直使用预训练模型。但实际项目中，你可能需要针对自己的数据集重新训练模型，或者调整类别数量。从下一节开始，我们将进入训练环节，详细讲解如何利用我们构建的数据存储 dsTrain，结合自定义损失函数和训练选项，从零开始训练 U‑Net。

敬请期待第五节《U‑Net 构建、训练与损失函数设计》。