简介:基于Qgis二次开发的遥感图像云检测及地物恢复系统,是一份聚焦遥感图像处理的Python毕业设计工程。系统综合运用地理信息系统、数字图像处理与深度学习技术,面向高校计算机、人工智能、通信工程等专业学生及开发者,既可用于毕设答辩,也能作为课设或项目立项演示。资源包共74个文件,包含21个Python脚本及对应pyc字节码文件、5个UI界面文件、XML/QRC配置、多张图标与示例图像、一个TIF测试影像和requirements依赖清单,压缩包整体约51.21MB;目录划分了主程序、算法模块、界面层及Qgis工具类,便于按功能定位代码。源码均实测运行成功,作者评审平均分达96分,附有README说明,可快速打开工程;现有大气校正、波段合成、图层属性等模块,读者可在此基础上修改或扩展云检测与地物恢复算法。目前已吸引254人下载学习,适合需要完整参考案例的中级及以上Python/GIS学习者。
1. 基于Qgis二次开发的遥感图像云检测及地物恢复系统:毕业设计到底在做什么
假设你手上有一批Landsat 8影像,想算植被指数或做土地利用分类,结果每一期都有云块遮挡。云区算出的指数异常离谱,分类模型甚至把云当成了地类。这个标题讲的是在Qgis这个开源GIS里,用Python写一套自动云检测、再把云遮挡下的地物信息恢复出来的流程。它不是空泛的论文框架,而是能落到代码和插件层面的毕业设计实现方案。适合遥感、GIS方向的毕设,也适合想做批量影像去云预处理的人。下面我会把选型理由、最小跑通方法、核心代码和踩坑记录一次说清楚。
2. 为什么选Qgis做二次开发:选型理由与云检测/恢复的底层逻辑
2.1 选Qgis而不是ArcGIS或纯GDAL:二次开发的三个理由
先回答一个逃不开的问题:为什么不用ArcGIS?ArcGIS 的 ArcPy 在栅格处理和分析上很强,但它有两个硬伤。第一是授权和部署,正版费用高,学校机房装一套要折腾半天,答辩现场经常因为许可弹窗导致演示失败;第二是闭环,ArcPy 脚本确实能处理数据,但要让评委直观看到“原图—云掩膜—恢复图”的对比,还得再去别的软件排版出图。Qgis 是开源桌面 GIS,安装包小,启动快,内置 Python 控制台,打开就能写脚本,这对毕业设计来说省下的时间足够你去调好几轮阈值。
第二个理由是 PyQGIS 的底层就是 GDAL/OGR。云检测绕不开栅格读写、地理变换和投影处理,而 PyQGIS 把 GDAL 封装在 Qgis 的 API 里。你在插件中既可以用QgsRasterLayer来显示影像,也可以继续用osgeo.gdal做底层像元读取。两种方式混用基本没问题。这比纯 GDAL 脚本更友好,因为处理完一个中间结果,你立刻把结果添加为 Qgis 图层,用“识别要素”点开像元值看看到底对不对,比在命令行里反复 print 数字直观太多。很多新手拿纯 GDAL 写了一大段脚本,跑完发现掩膜偏移了一个像元,找错找半天,而这种可视化迭代能很快定位。
第三个理由是 Qgis 的布局(Print Layout)出图能力。答辩需要一张图同时展示假彩色影像、云掩膜、恢复结果和精度指标。Qgis 的布局里可以加地图框、图例、比例尺、指北针,直接导出 300 dpi 的 PNG 放进论文。如果只用 GDAL 脚本,最后还得把中间 tif 拉到别的软件里排版,多走一步就多一次格式转换的麻烦。所以这套系统的架构定位就是:Qgis 做容器与展示,Python 写算法,GDAL 做底层数据读写。这也正好覆盖了“数据处理+算法实现+平台展示”三个毕业设计考核点。
2.2 云检测的两条主流技术路线:光谱阈值与语义分割
云检测在遥感图像处理里并不是新问题,主线分为两类。第一类是光谱阈值法。云在可见光波段的反射率显著高于普通地表,尤其在蓝波段;而在热红外波段,云顶亮温明显低于陆地表面。利用这些物理特征,可以构造多通道判断规则:蓝波段反射率大于某值、归一化差异云指数 NDTI 大于某值、亮温小于某值,三个条件求交集得到云掩膜。这个方法最大的优点是零训练成本,不需要标注样本,跑一景 8000×8000 的 Landsat8 影像只要几十秒。每个阈值都有物理含义,论文里好解释,答辩时评委问“为什么这个值设成 0.15”,你可以从云的光谱曲线讲起。
第二类是语义分割路线,也就是把云检测当成遥感图像语义分割任务做逐像元分类。典型模型是 U-Net,输入真彩色或多波段影像切片,输出每个像元“云/非云”的概率。语义分割的上限高,能学习云的纹理和上下文信息,对薄云、碎云更鲁棒;缺点是数据成本太高。你需要先手工标注一批云掩膜,而遥感图像标注本身就是一件苦差事,云和雪、云和白色建筑物、云和干燥盐碱地在图像上人眼都容易混淆,标注质量直接决定模型表现。此外还要安装 PyTorch、配置 CUDA、做切片增强,这些对毕业设计进度是很大的挑战。
我一般建议的路线是:系统主线用光谱阈值加形态学后处理,把语义分割作为扩展对比实验。这样你有稳定可交付的流程,也能在论文里写“未来可用深度学习方法替换云检测模块”,反而体现你对技术选型有判断。如果基本功扎实,想冲击高分,可以在阈值法跑通之后,再用同一套标注数据训练一个 U-Net,算出两套模型的 IoU 对比。但不要把深度学习作为唯一实现路径,万一训练不收敛,连系统都交不出来。
2.3 地物恢复的核心思路:时序填充、邻域插值和多角度替补
检测到云之后,难点转移到“云底下到底是什么”。地物恢复有三种工程策略,我按优先级排个序。
第一种是多时相替换。卫星对同一区域有重访周期,Landsat8 是 16 天,Sentinel-2 是 5 天。你今天拍的影像上某块有云,另一天拍的同一位置可能晴空无云。把有云影像记为基准,找同区域的其他时相影像,把云像元替换成参考时相同一位置的像元值。这个策略要求两期影像的辐射差异不大,且地物本身没有剧烈变化。对水体、建筑、裸地这类稳定地物,效果最好;对农田,最好选择物候相近的时相,否则替换上去的是不同生长季的植被,分类结果还是错。
第二种是邻域插值修复。如果只有单景影像,找不到多时相参考,那就只能利用云周围的无云像元去做空间插值。最常用的实现是 OpenCV 的 inpainting 算法,或者自己写半径加权平均。这个方法只适用于小块云和孤立云。云区直径超过几十个像元时,插值出来的区域是平滑估计,真实地表纹理全丢了,后续做分类会把这些位置误判成均质土地。所以邻域插值在这套系统里定位成兜底方案,只在无参考时相时启用。
第三种是多角度替补,利用卫星侧摆能力从另一视角拍摄同一地点,看到云遮挡背后的地物。这条对未来商业星座来说可行,但 Landsat 和 Sentinel 这类常规光学卫星很难在短时间内获得同一日期不同角度的配对影像,毕业设计数据上不好支撑,一般只在论文综述里提一句即可。因此整套恢复流程采用“多时相替换为主,邻域插值为辅”的组合策略,下面的实现部分会详细给代码。
3. 搭好PyQGIS开发环境:从Qgis安装到第一个云检测脚本跑通
3.1 先搞清楚Qgis内置Python和独立Python的关系
很多同学照网络教程先在系统里装了一个 Python 3.10,又装了 Qgis,然后在系统 Python 里执行pip install qgis,结果报错一屏。这是最典型的起步坑。Qgis 是自带 Python 解释器的,这个解释器位于 Qgis 安装目录下(比如C:\Program Files\QGIS 3.28\apps\Python\python.exe),Qgis 本体就是用这个解释器启动的。系统里另外装的 Python 和 Qgis 没有关系,在系统 Python 里装 PyQGIS 需要处理 Qt 和 Qgis 库的依赖版本,复杂度很高,没有必要去碰。
正确做法是直接使用 Qgis 内置环境:启动 Qgis,打开菜单栏的“插件→Python 控制台”或“处理→Python 控制台”,在里面输入from qgis.core import QgsRasterLayer,正常情况下不会报错。第三方库 nginx、numpy、scipy、opencv 默认不一定齐全,需要在 Qgis 的内置 Python 环境里安装。打开命令行,切换到apps\Python目录,执行python.exe -m pip install numpy scipy opencv-python,注意这里用的是 Qgis 自己的 python.exe,而不是桌面上的 Python 快捷方式。如果安装后控制台里 import 还是失败,检查一下是否装到了用户目录下的 site-packages,Qgis 环境可能不自动扫描那里。
判断环境是否正常有一个小技巧:在 Qgis 的 Python 控制台里运行import sys; print(sys.executable),它会打印当前解释器的路径。如果路径显示的是 Qgis 安装目录下的 apps/Python,那说明环境正确。如果不是,说明控制台用了错误的解释器,需要检查 Qgis 的 Python 插件设置。
3.2 在Qgis的Python控制台里跑通第一个云检测
环境确认后,写一个最小可运行的云检测脚本,目的是验证数据链路没有问题。这里使用 Landsat8 蓝波段表观反射率产品,注意必须是 TOA 反射率,不能用原始 DN 值。
from osgeo import gdal import numpy as np # 打开Landsat8蓝波段TOA反射率文件 raster_path = r"D:\landsat\LC08_B2_TOA.TIF" ds = gdal.Open(raster_path) if ds is None: raise IOError("无法打开文件,请检查路径是否包含中文或空格") band = ds.GetRasterBand(1) blue = band.ReadAsArray().astype(np.float32) # 最简单的云检测:蓝波段反射率大于0.15认为可能是云 cloud_mask = blue > 0.15 # 统计云量 cloud_ratio = float(cloud_mask.sum()) / cloud_mask.size print("影像大小: {} x {}".format(blue.shape[1], blue.shape[0])) print("云像元占比: {:.2f}%".format(cloud_ratio * 100)) # 保存掩膜为GTiff,并拷贝原影像的地理参考 driver = gdal.GetDriverByName("GTiff") out_ds = driver.Create(r"D:\landsat\cloud_mask.tif", blue.shape[1], blue.shape[0], 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(cloud_mask.astype(np.uint8)) out_ds.FlushCache() out_ds = None逻辑说明:先用 GDAL 打开蓝波段反射率,转成 float32 numpy 数组;blue > 0.15得到布尔掩膜,True 表示云像元;然后统计占比;最后创建一个新 tif 保存掩膜。在创建输出文件时,SetGeoTransform和SetProjection是把原始影像的地理仿射变换和投影信息拷贝给输出,否则 Qgis 里无法将掩膜和原图叠在一个坐标系上。
参数说明:0.15 这个阈值只是最小验证门限,它对应 Landsat8 蓝波段云反射率的下界。实际使用时不同季节、不同太阳高度角下,薄云反射率可能只有 0.1,地面雪或云可能到 0.3 以上,所以阈值必须根据影像动态调整。后续第 4 章会换成多特征组合,并把所有阈值做成可配置参数。
3.3 把检测脚本做成Qgis插件入口
控制台里跑通只能算实验脚本,毕业设计要交的是一个“系统”。最常见的方式是把算法封装成 Qgis 插件。插件是一个目录,至少包含三个文件:
cloud_detect_plugin/ ├── __init__.py ├── metadata.txt └── cloud_detect_plugin.pymetadata.txt写插件名称、版本、Qgis 最低版本;__init__.py定义classFactory(iface)返回插件实例;核心算法放在cloud_detect_plugin.py中。如果你不想从零写,可以在 Qgis 的插件管理器里安装 Plugin Builder,它会生成完整的模板,然后你把自己的代码填进去。
插件里有两类错误要避免。一是import osgeo必须放在算法函数内部,而不是模块顶层。因为 Qgis 加载插件时会在主线程执行模块顶层代码,此时 GDAL 的环境可能还没初始化完全,很容易出现加载失败。二是不要给插件文件起qgis.py、gdal.py这类名字,会遮蔽同名模块导致 import 异常。插件写好之后,放到用户目录下的AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins,重启 Qgis 后就能在“插件”菜单里看到。插件的好处是可以在界面上加按钮、加进度条,也可以注册成 Processing 算法,让用户从处理工具箱里调用。
4. 核心实现:把云检测和地物恢复写成一个可交付的流程
4.1 多特征云检测:蓝波段+NDTI+表面温度阈值的组合拳
单一蓝波段阈值会把雪、水体边缘、亮色屋顶都误判成云。更稳妥的做法是用蓝波段反射率、归一化差异云指数 NDTI 和亮温三个特征组合判断。
from osgeo import gdal import numpy as np from scipy import ndimage def read_band(path): """读取单波段影像,返回数组和数据集对象""" ds = gdal.Open(path) if ds is None: raise IOError("打开文件失败: " + path) return ds.GetRasterBand(1).ReadAsArray().astype(np.float32), ds def cloud_detect(blue_path, swir_path, temp_path, out_path, blue_th=0.15, ndti_th=0.05, temp_max=300.0): """ 多特征云检测 blue_path: 蓝波段TOA反射率 swir_path: 短波红外1波段TOA反射率 temp_path: 热红外亮温(单位K) """ blue, ds_blue = read_band(blue_path) swir, _ = read_band(swir_path) temp, _ = read_band(temp_path) # 归一化差异云指数:云在蓝光反射高、短波红外反射低 ndti = (blue - swir) / (blue + swir + 1e-6) # 云像元需要同时满足三个条件 cloud = (blue > blue_th) & (ndti > ndti_th) & (temp < temp_max) # 形态学后处理:开运算去掉孤立噪声,闭运算填充云内部空洞 cloud = ndimage.binary_opening(cloud, iterations=2) cloud = ndimage.binary_closing(cloud, iterations=2) # 保存掩膜 driver = gdal.GetDriverByName("GTiff") out_ds = driver.Create(out_path, ds_blue.RasterXSize, ds_blue.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds_blue.GetGeoTransform()) out_ds.SetProjection(ds_blue.GetProjection()) out_ds.GetRasterBand(1).WriteArray(cloud.astype(np.uint8)) out_ds.FlushCache() out_ds = None return cloud # 调用示例 cloud_detect( r"D:\landsat\LC08_B2_TOA.TIF", r"D:\landsat\LC08_B6_TOA.TIF", r"D:\landsat\LC08_B10_Temp.tif", r"D:\landsat\cloud_mask.tif" )逻辑说明:先读入三个波段。ndti是归一化差异云指数,云的蓝光反射率高、短波红外反射率低,所以 NDTI 会明显大于 0;水体的蓝光反射率虽高,但短波红外也很高,两者相减后 NDTI 会被压小。亮温条件temp < 300用来排除高温地表。最后用 scipy 的ndimage.binary_opening和binary_closing做形态学操作,能滤掉零星误检的孤立像元,也能把云块内部由于薄云漏检产生的空洞填上。
参数说明:常用参数区间如下表。
| 参数 | 典型值 | 调参方向 |
|---|---|---|
| blue_th | 0.15 | 误检雪/白地多时升到0.20;薄云漏检时降到0.12 |
| ndti_th | 0.05 | 误检水体时升到0.10;漏检卷云时降到0.02 |
| temp_max | 300 K | 去除高温地表误检时降到295;保留暖云时提高到305 |
| iterations | 2 | 噪声大时开运算提高;大云被削边时降低 |
这些阈值一定不要死记,每期影像都要做一次分位数诊断。我在开发时习惯先打印蓝波段反射率的 2% 和 98% 分位数,再判断阈值范围,比如 98% 分位数才 0.12,说明这期影像整体偏暗,0.15 的阈值就会把大量云漏掉。
4.2 地物恢复第一招:基于多时相影像的最小替换
云掩膜生成后,如果研究区有多期影像,优先用多时相替换。假设你已经准备了 n 期同一区域的影像,并且用上面cloud_detect对每一期生成了云掩膜。现在需要把每一期的对应波段数组放在一个列表里,做逐像元替换。
def reconstruct_from_series(mask_list, band_arrays): """ 用多时相无云像元恢复基准影像 mask_list: 每期影像的云掩膜列表,True表示云 band_arrays: 每期影像的同一波段的反射率数组列表 返回: result: 恢复后的波段数组 filled: 布尔数组,标记哪些像元已被填上 """ # 以第一期作为基准,初始化结果 result = band_arrays[0].copy() filled = ~mask_list[0] # 第一期非云的位置算已填 for i in range(1, len(band_arrays)): arr = band_arrays[i] # 这一期非云,且当前结果还没填好的位置 valid = ~mask_list[i] & ~filled result[valid] = arr[valid] filled[valid] = True return result, filled逻辑说明:初始化结果等于第一期影像,filled记录第一期已经无云的位置。然后从第二期开始遍历,凡是当期不是云、且之前没填过的像元,就用当期值覆盖。这样做保留了最早可用的有效像元;如果你想优先使用时间上最接近的参考影像,可以调整 mask_list 的排序。返回的filled数组里False的位置表示所有时相都是云,是下一步要处理的空洞。
参数说明:多时相替换有一个容易被忽略的前置条件,就是影像必须先做裁剪和对齐。如果不同期影像的边界不完全重合,存在偏移,直接按像元位置替换会串位。我在处理前先运行 GDAL 的gdal.Warp把所有期影像重采样到基准影像的网格,并统一投影、统一空间分辨率。另外,替换后往往需要做直方图匹配,否则两期影像的辐射差异会在云边界形成明显的拼接缝。直方图匹配的简化做法是统计两期影像在无云区的均值和标准差,对参考期影像做线性变换:ref_new = (ref - ref_mean) / ref_std * base_std + base_mean,做完这一步,拼缝基本能压下去。
4.3 地物恢复第二招:没有多时相数据时的邻域修复
当只剩单景影像,或者多时相替换后还有filled=False的残留空洞,就轮到邻域插值兜底。OpenCV 的 inpainting 算法在遥感小空洞修复里很好用,因为它不是简单拉平,而是沿边缘的梯度方向向洞里递推颜色信息,对几到几十个像元的孔洞效果不错。
import cv2 def neighbor_inpaint(image, cloud_mask, radius=5): """ image: 多波段数组,形状为 HxWxC,值域可为0-1反射率 cloud_mask: 云掩膜,True表示云/空洞 radius: 修复半径,决定参考邻域范围 """ # OpenCV的inpaint要求mask是0和255的uint8格式 mask = (cloud_mask * 255).astype(np.uint8) # 对浮点波段先缩放到0-255,修复后再还原,这里用简单归一化 images_u8 = np.zeros(image.shape, dtype=np.uint8) reduced = np.zeros(image.shape, dtype=np.float32) restored = np.zeros(image.shape, dtype=np.float32) for b in range(image.shape[2]): ch = image[:, :, b] ch_min, ch_max = ch.min(), ch.max() ch_u8 = (((ch - ch_min) / (ch_max - ch_min + 1e-6)) * 255).astype(np.uint8) # inpaint是单通道处理,逐波段来 repaired = cv2.inpaint(ch_u8, mask, radius, cv2.INPAINT_TELEA) restored[:, :, b] = repaired.astype(np.float32) / 255 * (ch_max - ch_min) + ch_min return restored逻辑说明:OpenCV 的 inpaint 只接受 uint8 单通道。因此对每个波段先做最小值最大值线性拉伸到 0-255,修复后再拉伸回原值域。mask 中 255 表示需要修复的区域。cv2.inpaint会从空洞边缘开始,用 TELEA 算法逐步向内填充。radius表示每个像元搜索邻域半径,半径越大能修越大的洞,但平滑程度也随之增加。
参数说明:radius=5是一个比较折中的值。如果云块比较碎,可以用 3,保留更多纹理细节;如果云块直径较大,可以试 10 或 15。超过 15 之后,修补区域的纹理基本消失,建议改用多时相方案。此外,逐波段独立修复会破坏波段间的光谱关系,所以只把它作为最后手段。实际流程中,我在执行neighbor_inpaint前,会对filled=False的位置生成一个单独的空洞掩膜,只修这些位置,不把原本的无云区再修一遍,避免无谓的细节损失。
5. 避坑与排查:云检测和地物恢复最容易翻车的5个地方
5.1 云却检不出来:阈值设置和影像校准不一致
现象:同一套阈值在一个研究区效果很好,换到另一个区域的影像,云漏掉一大半,或者在整景影像上检测出大量条带状噪声。
原因:不同影像可能没有做辐射定标,有些是原始 DN 值,有些是表面反射率产品;另外不同季节太阳高度角不同,表观反射率整体偏移。
解决:第一步先把所有输入统一到 TOA 反射率。Landsat8 的 Level-1 产品可以用 MTL 文件里的RADIANCE_MULT_BAND_x和RADIANCE_ADD_BAND_x将 DN 转成辐亮度,再除以余弦太阳天顶角得到反射率。第二步在跑算法之前打印蓝波段反射率的 2% 和 98% 分位数,确认影像动态范围正常。如果 98% 分位数只有 0.1,说明这期影像整体很暗,需要把 blue_th 下调到 0.08 再试。
5.2 雪地和亮色建筑被误判成云,导致地物恢复把好像素覆盖掉
现象:云掩膜在山区大片标白,城市里的白色屋顶也被圈进来。恢复阶段把这些非云的好像素当成云,用其他时相替换,结果把原本正确的建筑颜色改没了。
原因:雪和云在可见光反射率很接近,白屋顶和薄云在 NDTI 上也没有明显区分。只看可见光波段,无法排除这些表面类型。
解决:在检测条件中加入亮温或近红外反射率的区分。云顶亮温普遍低于 280K,而雪和建筑地表的亮温通常高于 285K,把temp_max调到 285 而不是 300,能过滤掉大部分低温地物误检。如果研究区确实有积雪,则再加一条“近红外反射率大于 0.3 则不算云”的规则,因为云的近红外反射率比雪低。最后,对小面积连通域做过滤,面积小于几十个像元的斑点直接置为非云,可以进一步压住零散误检。
5.3 多时相替换后出现拼接缝和跳变
现象:恢复图像像拼图,云边界一圈有明显的补丁色调,甚至能看到参考影像的辐射差异。
原因:两期影像的太阳高度角、气溶胶、传感器响应状态不同,直接硬替换必然存在辐射不一致。简单拷贝像素值没有做匀色。
解决:替换之前先做直方图匹配。在两期影像的无云重叠区分别计算均值和标准差,然后对参考影像做线性辐射归一化,使参考影像的均值和方差对齐基准影像。做完这一步,再执行替换。另外在云边界做一个羽化过渡,比如对云掩膜做 5 到 10 个像元的缓冲,在缓冲区内距离加权混合基准值和参考值,拼缝几乎看不出来。
5.4 PyQGIS脚本在插件里import一直失败
现象:在 Qgis 控制台里运行import osgeo正常,但做成插件后一加载就报ModuleNotFoundError,插件菜单里看不到自己的工具。
原因:插件模块顶层引用了osgeo或cv2,Qgis 加载插件时主线程没有预先设置好相关环境路径,或者插件文件名与标准库冲突,导致导入异常。
解决:把from osgeo import gdal、import numpy、import cv2全部移入到算法函数内部,模块顶层只保留from qgis.core import ...等 Qgis API。这牺牲了一点加载速度,但换来了稳定。另外,检查插件目录名和文件名是否包含qgis、gdal、numpy等关键词,如果有,改名。这个问题在插件打包给别人用时尤其常见,提前规避能省很多答疑环节。
5.5 大影像处理内存溢出,程序静默崩溃
现象:处理一景完整 Landsat8 影像,脚本执行到一半 Qgis 突然失去响应,几秒后整个程序退出,控制台没有任何错误输出。
原因:将三到四个全波段一次性ReadAsArray成 float32,再加上中间 NDTI、阈值判断、形态学膨胀变量,内存占用轻松超过 8 GB。Qgis 进程内存不足时,操作系统直接将其杀掉,而不是抛出 Python 异常。
解决:改成分块处理。用band.ReadAsArray(xoff, yoff, xsize, ysize, bufxsize, bufysize)每次读 512×512 像元,处理完写回输出文件,再继续下一块。块之间要留出 morph 算子所需的重叠区域,比如形态学 iterations=2 时需要额外扩展 4 到 6 行。另一个更省事的途径是先把影像下采样或者用gdal_translate -outsize 10% 10%做粗略预览,在缩略图上调好阈值,再处理全分辨率,好过直接让进程崩溃。
6. 进阶验证:把毕设做出能答辩的精度报告和对比图
6.1 用人工标注样本算IoU与像素精度
云检测不能只拿眼睛看“差不多”。答辩评委都会问“准确率怎么算”。我习惯的做法是在研究区随机选 3 到 5 个 1024×1024 的子区,用 Qgis 手动绘制云多边形,栅格化成真实云掩膜,然后和算法输出做逐像元对比。
from sklearn.metrics import confusion_matrix def cloud_metrics(truth_mask, pred_mask): """ truth_mask: 人工标注的云掩膜,1表示云 pred_mask: 算法输出的云掩膜,1表示云 """ t = truth_mask.ravel() p = pred_mask.ravel() tn, fp, fn, tp = confusion_matrix(t, p, labels=[0, 1]).ravel() iou = tp / (tp + fp + fn + 1e-6) pa = (tp + tn) / (tp + tn + fp + fn) return { "IoU": round(iou, 4), "PixelAccuracy": round(pa, 4), "误检像元": int(fp), "漏检像元": int(fn), }代码逻辑很直接:把两期掩膜拉平为一维数组,用混淆矩阵算出四类统计。IoU是云区预测正确像素数与“预测云或真实云的总数”之比,IoU 大于 0.7 算合格,低于 0.5 就需要重新调阈值。PixelAccuracy是整体正确率,它容易被大面积非云像素拉高,所以主要看 IoU。把多组子区的结果汇总成表,既是论文里的定量分析,也是答辩时最有力的证据。
6.2 用Qgis布局出图做对比展示
数值指标之外,可视化对比很关键。在 Qgis 的布局里加三个地图框,第一个框显示原始假彩色影像,第二个框叠加云掩膜(用红色半透明);第三个框显示恢复后的影像。再放一个表格显示 IoU 和像素精度。导出设 200 dpi,分辨率足够放进论文。这张图的作用比十段文字都大,它直接告诉评委:你这个系统的输入、输出、误差是可视化的闭环。
最后讲一个我自己的教训:刚开始做云检测时,我只盯着 IoU 数字,不看图上结果,结果发现算法把一处大面积蓝绿色屋顶当成云标成 1,IoU 依然高,因为那个屋顶只占子区面积的 1%。从那以后我给自己定了个习惯——每一组数据跑完,必须把结果叠加进 Qgis 目视一圈,再算指标。数字会骗人,图不会。希望这个习惯能帮你少熬几个改参数的夜,也希望这篇笔记能让你把系统做扎实,顺利走进答辩教室。希望帮到你。
本文还有配套的精品资源,点击获取