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

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

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

具体实现代码如下：

Compile_opt idl2
    ENVI,/restore_base_save_files
    ENVI_BATCH_INIT
    
    tm_file ='E:\yaoganhomework\ch08+Classification_tm\Classification_tm\tm_data\LT51240362007139_JZ.img'
        
    ENVI_OPEN_FILE, tm_file, r_fid=fid
    ENVI_FILE_QUERY, fid, dims=dims, ns=ns, nl=nl, nb=nb
    
    tm_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 begin
  for j=0,dims[1]-1 do begin
    arr[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 of
      0:begin
         tm_data[i,j,0]=250
         tm_data[i,j,1]=237
         tm_data[i,j,2]=115
          end 
      1:begin 
        tm_data[i,j,0]=33
        tm_data[i,j,1]=138
        tm_data[i,j,2]=33 
        end
      2:begin 
        tm_data[i,j,0]=255
        tm_data[i,j,1]=0
        tm_data[i,j,2]=0 
        end
      3:begin 
        tm_data[i,j,0]=23
        tm_data[i,j,1]=23
        tm_data[i,j,2]=115 
        end
      4:begin 
        tm_data[i,j,0]=189
        tm_data[i,j,1]=189
        tm_data[i,j,2]=189 
        end
    endcase
  endfor
endfor


class_data[*,*,0]=tm_data[*,*,0]
class_data[*,*,1]=tm_data[*,*,1]
class_data[*,*,2]=tm_data[*,*,2]

tv,class_data,true=3,/order

    pos_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 6
    
    out_name ='E:\yaoganhomework\ch08+Classification_tm\Classification_tm\tm_class_jz.img'
    envi_enter_data, class_data, map_info = map_info, r_fid=out_fid
    envi_output_to_external_format, fid=out_fid, out_name = out_name, dims=dims, pos=pos_class, /envi

输出的结果如下：

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