像元二分模型计算植被覆盖度

江北20190411

关注

原创最新推荐文章于 2026-05-30 11:44:11 发布 1.2w 阅读 AI 写同款· GEO 优化 ›

一区TOP期刊提出一种植被覆盖度估算新方法参数“β”表示R反射率的比例，而RE1反射率的比例由参数表示(1-'β').当 R2是最高的，对应的权重参数是最好的。基于 NDVI、SAVI和RVI，结合波段相关性确定估算FPV and FNPV的最佳波段，使用RE3和RE2反射率的加权平均值替换NIR波段，而RE1和RE2反射率的加权平均值替换R波段。1）研究植被理化参数对PV和NPV冠层反射率的影响 2）采取新的RE波段前后使用VI估计FPV and FNPV的差异 3）评估4波段红边植被指数对FPV and FNPV的估计能力。阅读详情

前言

关于改进的像元二分模型第一次出现是在李苗苗老师的博士论文中，但她后来发表的一篇期刊文章把这一部分单独摘了出来并做了应用案例，期刊链接在下方。

密云水库上游植被覆盖度的遥感估算

论文要点

计算公式为：Fc =(NDVI -NDVIsoil) (NDVIveg -NDVIsoil)
问题的关键在于NDVIveg 与NDVIsoil如何取值。

其中 ,NDVIsoil 为完全是裸土或无植被覆盖区域的NDVI 值 , NDVIveg 则代表完全被植被所覆盖的像元的NDVI 值, 即纯植被像元的NDVI 值。NDVIveg 代表着全植被覆盖像元的最大值, 由于植被类型的影响, NDVIveg 值也会随着时间和空间而改变。因此 , 计算植被覆盖度时, 既使同一景影像, 对于NDVIsoil 和NDVIveg 值不能取固定值。在实地应用中,植被覆盖类型随土地利用类型而变化, 对于某一土地利用类型, 由于植被类型近似, 其 NDVIveg 值也相近;而对于特定的土壤类型 ,其 NDVIsoil 值也应该是一定的。因此, 土地利用图和土壤图应作为计算 NDVIveg 和NDVIsoil 值的基础

计算过程

辅助数据介绍

全球土壤数据
HWSD全球土壤数据下载处理
全球土地分类数据
MCD12Q1
关于MCDQ1，由官网介绍可知数据中包含5种分类方法得到的结果。
在这里插入图片描述
数据的说明手册目前网址404了。但我们可以查看数据的头文件获取分类代码信息。
这里贴出了type3的分类代码信息，比较符合我们计算植被覆盖度的目的。

Dataset #3: Land_Cover_Type_3  
  Dims: BYTE (2400 x 2400)    
Attribute 3-1: "long_name"  
  "Land_Cover_Type_3"   
Attribute 3-2: "Water Bodies"  
  0   
Attribute 3-3: "Grasslands"  
  1   
Attribute 3-4: "Shrublands"  
  2   
Attribute 3-5: "Broadleaf Croplands"  
  3   
Attribute 3-6: "Savannas"  
  4   
Attribute 3-7: "Evergreen Broadleaf Forests"  
  5   
Attribute 3-8: "Deciduous Broadleaf Forests"  
  6   
Attribute 3-9: "Evergreen Needleleaf Forests"  
  7  
Attribute 3-10: "Deciduous Needleleaf Forests"  
  8   
Attribute 3-11: "Unvegetated"  
  9   
Attribute 3-12: "Urban and Built-up Lands"  
  10   
Attribute 3-13: "Unclassified"  
  255   
Attribute 3-14: "valid_range"  
  0, 10    
Attribute 3-15: "_FillValue"  
  255

使用MCTK对type3进行投影，重采样，至此LULC和HWSD，NDVI数据都准备好了。

IDL实现

pro vfc_cal
  compile_opt idl2
  e = envi()
  workpath = 'F:\FY3D_MERSI_data\Subset'
  cd,workpath
  ndvifile = workpath + '\NDVI'
  lulcfile = workpath + '\LULC'
  hwsdfile = workpath + '\HWSD'
  
  ndvir = e.openraster(ndvifile)
  hwsdr = e.openraster(hwsdfile)
  lulcr = e.openraster(lulcfile)
  
  spatialref = ndvir.spatialRef
  
  ndvi = ndvir.getdata(bands=0)
  hwsd = hwsdr.getdata(bands=0)
  lulc = lulcr.getdata(bands=0)
  
  ndvi[where(finite(ndvi,/nan))]  = mean(ndvi,/nan)
  lulc[where(finite(lulc,/nan))]  = 4
  
  hwsd_class = hwsd[uniq(hwsd,sort(hwsd))]
  lulc_class = lulc[uniq(lulc,sort(lulc))]
  print,hwsd_class,lulc_class
   
  lutmin = make_array(2,n_elements(hwsd_class),type=5)
  for i=0,n_elements(hwsd_class)-1 do begin 
    lutmin[0,i] = hwsd_class[i]
    ndvi_hwsd = ndvi[where(hwsd eq lutmin[0,i])] 
    ndvi_hwsd = ndvi_hwsd[sort(ndvi_hwsd)]
    ndvi_hwsd = ndvi_hwsd[where(hwsd ge -0.1)]
    number = n_elements(ndvi_hwsd)
    lutmin[1,i] =  ndvi_hwsd[round(0.05*number)]  
  endfor

  lutmax = make_array(2,n_elements(lulc_class),type=5)
  for i=0,n_elements(lulc_class)-1 do begin
    lutmax[0,i] = lulc_class[i]
    ndvi_lulc = ndvi[where(lulc eq lutmax[0,i])]
    ndvi_lulc = ndvi_lulc[sort(ndvi_lulc)]
    number = n_elements(ndvi_lulc)
    lutmax[1,i] =  ndvi_lulc[floor(0.95*number)]
  endfor
 
  dim = size(ndvi)
  n = n_elements(ndvi)
  vfc = make_array(dim[1],dim[2],type=5)
  for i=0,n-1 do begin
    for j=0,n_elements(hwsd_class)-1 do begin
      for k=0,n_elements(lulc_class)-1 do begin
        if lulc[i] eq lutmax[0,k] and hwsd[i] eq lutmin[0,j] then begin
          case 1 of 
            ndvi[i] ge lutmax[1,k] : vfc[i] = 1
            ndvi[i] le lutmin[1,j] : vfc[i] = 0
            else: vfc[i] = (ndvi[i]-lutmin[1,j]) / (lutmax[1,k]-lutmin[1,j])
          endcase
        endif
      endfor
    endfor
  endfor
  
  tempfn = e.GetTemporaryFilename()
  vfcr = e.CreateRaster(tempfn, vfc, SpatialRef = spatialref)
  vfcr.save
 
end