IDL监督分类——极大似然估计

代码为HPU遥感专业作业，仅供参考。

 需要解释一下，X为单个像素六个波段的灰度值数组，mi是采集的样本六个波段中每一个波段的平均灰度值，是协方差阵，是协方差阵的逆，是协方差阵的模。

具体实现代码如下：

Compile_opt idl2ENVI,/restore_base_save_filesENVI_BATCH_INITtm_file ='E:\yaoganhomework\ch08+Classification_tm\Classification_tm\tm_data\LT51240362007139_JZ.img'ENVI_OPEN_FILE, tm_file, r_fid=fidENVI_FILE_QUERY, fid, dims=dims, ns=ns, nl=nl, nb=nbtm_data = intarr(ns,nl,nb)tm_data[*,*,0] = envi_get_data(fid=fid, dims=dims, pos=0)tm_data[*,*,1] = envi_get_data(fid=fid, dims=dims, pos=1)tm_data[*,*,2] = envi_get_data(fid=fid, dims=dims, pos=2)tm_data[*,*,3] = envi_get_data(fid=fid, dims=dims, pos=3)tm_data[*,*,4] = envi_get_data(fid=fid, dims=dims, pos=4)tm_data[*,*,5] = envi_get_data(fid=fid, dims=dims, pos=5)map_info = envi_get_map_info(fid=fid)dims=size(tm_data[*,*,0], /DIMENSIONS)class_data=intarr(512,512,3)cov_crop=[[3.846530612,1.415510204,1.928163265,3.993469388,4.413877551,1.16],$[1.415510204,1.469795918,1.680816327,2.423673469,3.172653061,1.283265306],$[1.928163265,1.680816327,2.45877551,3.095510204,4.59755102,1.711020408],$[3.993469388,2.423673469,3.095510204,12.42612245,9.560816327,3.580408163],$[4.413877551,3.172653061,4.59755102,9.560816327,11.88938776,4.392653061],$[1.16,1.283265306,1.711020408,3.580408163,4.392653061,2.667755102]]cov_forest=[[8.704081633,3.946938776,5.746938776,0.842857143,8.967346939,6.395918367],$[3.946938776,2.434285714,3.353469388,0.293061224,5.103673469,3.760816327],$[5.746938776,3.353469388,5.348571429,-1.186938776,7.173877551,5.537142857],$[0.842857143,0.293061224,-1.186938776,34.62897959,10.22367347,-0.196734694],$[8.967346939,5.103673469,7.173877551,10.22367347,19.85469388,9.634285714],$[6.395918367,3.760816327,5.537142857,-0.196734694,9.634285714,8.026122449]]cov_urban=[[32.39020408,15.44,21.72734694,12.6122449,26.56,25.11102041],$[15.44,19.41387755,28.21061224,19.12244898,38.62693878,31.30040816],$[21.72734694,28.21061224,42.77714286,27.89795918,56.42204082,46.19755102],$[12.6122449,19.12244898,27.89795918,36.97959184,62.08163265,41.06122449],$[26.56,38.62693878,56.42204082,62.08163265,141.6995918,94.65877551],$[25.11102041,31.30040816,46.19755102,41.06122449,94.65877551,70.36285714]]cov_water=[[4.865714286,2.68,2.080816327,-0.034285714,0.219591837,0.758367347],$[2.68,2.667755102,2.105306122,-0.148571429,0.672653061,0.864489796],$[2.080816327,2.105306122,2.662857143,1.064489796,2.038367347,1.226122449],$[-0.034285714,-0.148571429,1.064489796,5.198367347,6.121632653,2.742040816],$[0.219591837,0.672653061,2.038367347,6.121632653,8.760816327,3.601632653],$[0.758367347,0.864489796,1.226122449,2.742040816,3.601632653,2.687346939]]cov_barren=[[26.08530612,12.47346939,12.31959184,-6.809795918,-12.78979592,-2.391428571],$[12.47346939,7.387755102,8.497959184,-2.085714286,-2.469387755,0.506122449],$[12.31959184,8.497959184,16.58,2.486530612,9.757142857,5.435510204],$[-6.809795918,-2.085714286,2.486530612,24.85061224,24.55510204,6.666938776],$[-12.78979592,-2.469387755,9.757142857,24.55510204,62.62244898,27.74897959],$[-2.391428571,0.506122449,5.435510204,6.666938776,27.74897959,16.82]]mo_crop=determ(cov_crop)
mo_forest=determ(cov_forest)
mo_urban=determ(cov_urban)
mo_water=determ(cov_water)
mo_barren=DETERM(cov_barren)ni_crop=MATRIX_POWER(cov_crop,-1)
ni_forest=MATRIX_POWER(cov_forest,-1)
ni_urban=MATRIX_POWER(cov_urban,-1)
ni_water=MATRIX_POWER(cov_water,-1)
ni_barren=MATRIX_POWER(cov_barren,-1)mi=[[83.52,37.14,34.52,82.32,54.22,22.16],$[83.9,38.88,35.72,84.06,79.68,32.88],$[117.24,57.12,66.28,65,94.88,62.38],$[110.54,54.84,50.52,28.84,24.12,13.92],$[107.42,55.8,69.54,84.08,122.7,70.58]]arr=intarr(6,1)out=intarr(5,1)for i=0,dims[0]-1 do beginfor j=0,dims[1]-1 do beginarr[0]=tm_data[i,j,0]arr[1]=tm_data[i,j,1]arr[2]=tm_data[i,j,2]arr[3]=tm_data[i,j,3]arr[4]=tm_data[i,j,4]arr[5]=tm_data[i,j,5]out[0]=-(arr-mi[*,0])##ni_crop##TRANSPOSE(arr-mi[*,0])-alog(mo_crop)out[1]=-(arr-mi[*,1])##ni_forest##TRANSPOSE(arr-mi[*,1])-alog(mo_forest)out[2]=-(arr-mi[*,2])##ni_urban##TRANSPOSE(arr-mi[*,2])-alog(mo_urban)out[3]=-(arr-mi[*,3])##ni_water##TRANSPOSE(arr-mi[*,3])-alog(mo_water)out[4]=-(arr-mi[*,4])##ni_barren##TRANSPOSE(arr-mi[*,4])-alog(mo_barren)a=max(out,index)case index of0:begintm_data[i,j,0]=250tm_data[i,j,1]=237tm_data[i,j,2]=115end 1:begin tm_data[i,j,0]=33tm_data[i,j,1]=138tm_data[i,j,2]=33 end2:begin tm_data[i,j,0]=255tm_data[i,j,1]=0tm_data[i,j,2]=0 end3:begin tm_data[i,j,0]=23tm_data[i,j,1]=23tm_data[i,j,2]=115 end4:begin tm_data[i,j,0]=189tm_data[i,j,1]=189tm_data[i,j,2]=189 endendcaseendfor
endforclass_data[*,*,0]=tm_data[*,*,0]
class_data[*,*,1]=tm_data[*,*,1]
class_data[*,*,2]=tm_data[*,*,2]tv,class_data,true=3,/orderpos_class=indgen(3)   ; 3 denote band numbers  if 1, means band number of outdata is 1;                        if 3, means band number of outdata is 3;                        if 6, means band number of outdata is 6out_name ='E:\yaoganhomework\ch08+Classification_tm\Classification_tm\tm_class_jz.img'envi_enter_data, class_data, map_info = map_info, r_fid=out_fidenvi_output_to_external_format, fid=out_fid, out_name = out_name, dims=dims, pos=pos_class, /envi

输出的结果如下：

 主要存在的问题是样本点选的过于集中导致分类效果差强人意。