IDL实战：从零实现遥感影像的极大似然监督分类

1. 遥感影像分类基础入门第一次接触遥感影像分类时我盯着屏幕上五彩斑斓的卫星图完全无从下手。后来才明白这就像教计算机玩找不同游戏——只不过我们要找的是地表上的农田、森林、水域等地物类型。极大似然监督分类就是其中最经典的游戏规则之一。传统分类方法主要分两大类监督和非监督分类。监督分类就像老师手把手教学需要我们先圈出样本区域比如这片红色是城市那片绿色是森林计算机通过这些样本学习判别规则。而非监督分类则是让计算机自己找规律适合我们对地物类型不太了解的场景。为什么选择极大似然法实测下来它的分类精度通常比最小距离、马氏距离等方法更稳定。原理其实很直观假设每类地物的光谱特征都服从正态分布当一个像元的光谱值在森林类中出现的概率最大时就把它划为森林。就好比看到一只动物有羽毛和喙我们判断它是鸟而不是鱼因为鸟的特征组合出现概率更高。在ENVI等商业软件中这个算法已经被封装成工具按钮。但就像只会用计算器做算术无法成为数学家一样要真正掌握分类精髓就得从底层实现它。接下来我将用IDL语言带大家从零搭建这个经典算法。2. 数据准备与预处理先准备好我们的实验材料——Landsat TM影像数据。这里我用的是LT51240362007139_JZ.img这个文件包含6个波段的多光谱数据。在IDL中读取数据的代码如下ENVI_OPEN_FILE, tm_file, r_fidfid ENVI_FILE_QUERY, fid, dimsdims, nsns, nlnl, nbnb tm_data intarr(ns,nl,nb) tm_data[*,*,0] envi_get_data(fidfid, dimsdims, pos0) ;...依次读取其他波段读取后建议做个简单检查用print,tm_data[3,3,*]查看某个像元各波段值是否正常。遇到过数据异常的情况比如全是0或65535溢出值这种像元需要后续过滤处理。关键的一步是获取训练样本统计值。我们需要为每类地物准备两个核心参数均值向量各类样本在各波段的平均值协方差矩阵反映各波段间的相关性比如农田类的参数可能是mi_crop [87.8,41.9,44.2,81.8,75.2,33.9] ;6个波段的均值 cov_crop [[3.84,1.41,...], [...], ...] ;6×6协方差矩阵获取这些参数有两种途径在ENVI中圈选样本区域自动计算手动测量典型地物样本适合小区域验证建议新手先用方法1记得检查样本可分离性——用ENVI的ROI Separability工具Jeffries-Matusita指数大于1.9才合格。曾经因为样本区分度不够导致分类结果一片混乱这个坑大家一定要避开。3. 极大似然算法核心实现算法核心在于计算每个像元属于各类别的概率。公式看起来复杂但拆解后很容易理解P(Class_i|x) -\frac{1}{2}[\ln|\Sigma_i| (x-\mu_i)^T\Sigma_i^{-1}(x-\mu_i)]用IDL实现这个公式时重点注意三个操作矩阵求逆INVERT(cov_matrix)矩阵乘法#运算符行列式计算DETERM()具体实现如下For i0,ns-1 do begin For j0,nl-1 do begin x Transpose(tm_data[i,j,*]) ;当前像元光谱向量 ;计算属于城市类的概率 diff x - mi_city p_city -(alog(determ(cov_city)) $ transpose(diff)#invert(cov_city)#diff)/2 ;...同理计算其他类别概率 ;找出最大概率对应的类别 pro_array [p_city, p_forest, p_water, p_crop, p_nudation] max_p max(pro_array) ;给像元赋类别值 case 1 of (p_city eq max_p): class_data[i,j,*] [255,0,0] ;红色 (p_forest eq max_p): class_data[i,j,*] [33,138,33] ;绿色 ;...其他类别 endcase endfor endfor调试时遇到过两个典型问题协方差矩阵奇异导致求逆失败——解决方法是在对角线加微小扰动如1e-6概率计算溢出——对数值太大时改用对数运算4. 分类结果输出与优化得到分类结果后可以用ENVI的导出功能保存为图像envi_enter_data, class_data, map_infomap_info, r_fidout_fid envi_output_to_external_format, fidout_fid, out_nameresult.img但直接分类的结果往往会有椒盐噪声——零星散落的错分像元。这时需要后处理众数滤波用3×3窗口内的多数类别替换中心像元去除小图斑删除面积小于阈值的孤立区域IDL中实现众数滤波的代码片段kernel [1,1,1,1,1,1,1,1,1] ;3×3矩形核 filtered CONVOL(class_data, kernel, /NAN) class_data (filtered ge 5) ? 1 : 0 ;阈值设为5我曾用北京郊区影像测试未滤波的总体精度只有78%滤波后提升到85%。但要注意滤波过度会导致地块边界模糊需要在保持边界和去除噪声间权衡。5. 精度验证与误区规避分类结果不能只看视觉效果必须定量评估。常用方法是混淆矩阵类别实测农田实测城市实测水域用户精度分类农田855292.4%分类城市372194.7%分类水域429093.8%制图精度92.4%91.1%96.8%总体精度89.3%实现代码;假设有验证样本truth和分类结果result confusion INTARR(nclass,nclass) FOR i0, nclass-1 DO BEGIN FOR j0, nclass-1 DO BEGIN confusion[i,j] TOTAL((truth EQ i) AND (result EQ j)) ENDFOR ENDFOR新手容易踩的坑样本偏差训练样本仅来自影像某区域导致其他区域分类差过拟合样本太少时协方差矩阵估计不准未考虑时相变化夏季植被茂盛时训练的模型不适合冬季影像有个项目曾因训练样本全来自旱季导致雨季影像中水体被错分为阴影。后来通过增加多时相样本解决了这个问题。6. 进阶优化方向当基础版本跑通后可以尝试以下优化波段选择用J-M距离评估各波段区分度剔除冗余波段jm 2*(1-exp(-bhattacharyya_distance))引入纹理特征结合GLCM矩阵的对比度、同质性等指标并行加速将影像分块用MPI_LAUNCH并行计算测试发现加入NDVI指数作为额外波段能使农田/草地区分精度提升7%。但要注意各特征量纲不同需要标准化处理。对于高光谱数据建议先做MNF变换降维。曾经处理过200波段的AVIRIS数据直接分类需要3小时降维到30个主成分后只需15分钟精度反而提高了2%。7. 与其他算法的对比实验在相同数据上对比几种算法效果算法总体精度Kappa系数运行时间最大似然89.3%0.862.1s支持向量机91.2%0.8915.7s随机森林90.8%0.888.3s1D-CNN92.1%0.9023.5s虽然深度学习更精准但极大似然的优势在于物理意义明确无需大量训练样本计算效率高特别是在应急监测场景下当需要快速获取初步分类结果时这个经典算法依然是首选。去年参与某次洪灾评估时用极大似然法10分钟就完成了受灾范围估算为救援决策争取了宝贵时间。