ENVI/IDL （5.3版本）一鍵化實現LST-NDVI的干濕邊方程擬合，并得到TVDI計算結果圖

0 原理介紹

利用IDL將NDVI異常值進行剔除，NDVI取值范圍為0.2~1(植被覆蓋區)，對反演的地表溫度Ts選擇5%的置信區間，構建NDVI-Ts特征空間，進行兩種反演方法(單通道算法、單窗算法)的干濕邊方程的擬合。
其中干濕邊方程擬合公式和TVDI計算公式如下：

其中：其中Tsmin是濕邊方程，Tsmax是干邊方程，TVDI的值域為[0,1]。TVDI越大，土壤濕度越低，TVDI越小，土壤濕度越高。
以NDVI為橫坐標，LST為縱坐標，得到NDVI-LST的散點圖。根據NDVI-LST的散點圖，就能得到每個NDVI值對應的LST的最大最小值，即要求的干邊和濕邊，最后根據干邊和濕邊的散點圖，就能擬合到干邊方程和濕邊方程。

1 打開影像所在文件夾

// 格式為datENVI,/RESTORE_BASE_SAVE_FILES;ENVI批處理模式/RESTORE_BASE_SAVE_FILESENVI_BATCH_INITfile_path=dialog_pickfile( get_path = newPath,$/DIRECTORY,$title='SELECT A DIRECTORY CONTAINS DAT-FILE')file_path_list=file_search(file_path,'*.dat')file_name_list=file_basename(file_path_list); print,file_path_listfilecount=N_ELEMENTS(file_path_list)

2 數據讀取

//************************************** 數據讀取*******************************************ENVI_OPEN_FILE, LSTname, r_fid =fid_LST;print,LSTnameENVI_FILE_QUERY, fid_LST, $nb = nb, $ ;波段數ns = ns, $ ;列數nl = nl ,$;行數dims=dimsmap_info=envi_get_map_info(fid=fid_LST)ENVI_OPEN_FILE, NDVIname, r_fid =fid_NDVILST=ENVI_GET_DATA(fid = fid_LST, dims = dims, pos = 0)+273.15 //℃+273.15---KNDVI=ENVI_GET_DATA(fid = fid_NDVI, dims = dims, pos = 0)

3 直方圖統計

hist=Histogram(LST,min=273.16,Binsize=0.1);統計直方圖，獲取累計像元個數totalCounts = TOTAL(hist, /CUMULATIVE);獲取累積百分比totalPercents = totalCounts/totalCounts[-1]//*************************************** 設定NDVI統計閾值及間隔*******************************NDVI_min=0.2 //NDVI>0.2為植被區域NDVI_max=1INTERVAL=0.01 ;統計間隔NUMS=CEIL((NDVI_max-NDVI_min)/INTERVAL) ;NDVI統計區間數目NDVIs=FINDGEN(NUMS)*INTERVAL+NDVI_MIN ;每個NDVI統計區間的起始值NDVIe=NDVIs+INTERVAL ;每個NDVI統計區間的終止值NDVIm=NDVIs+INTERVAL/2.0 ;;每個NDVI統計區間的均值//********************************* 干濕邊曲線擬合 *******************************************LST_MAX=FLTARR(NUMS)LST_MIN=FLTARR(NUMS) FOR i=0,NUMS-1 DO BEGINW=WHERE(NDVI1 GT NDVIs[i] AND NDVI1 LT NDVIe[i] ,COUNT)if count gt 0 then beginLST_MAX[i]=MAX(LST1[w]) ;當前NDVI區間對應的LST最大值LST_MIN[i]=MIN(LST1[w]) ;當前NDVI區間對應的LST最小值endifENDFOR

4 計算干濕邊方程

A1=REGRESS(NDVIm,LST_MAX,CONST=B1,CORRELATION=R1) ;干邊方程A2=REGRESS(NDVIm,LST_MIN,CONST=B2,CORRELATION=R2) ;TVDI=(LST-A2*NDVI-B2)/(((A1*NDVI+B1)-(A2*NDVI+B2)))

5 擬合干濕邊方程

FLAG=1 if FLAG EQ 1 THEN BEGIN;O_FN=DIALOG_PICKFILE(TITLE='干濕邊擬合圖像保存為：')+'.png'o_fn=outpath+filename+'SW1.png'PLOT_FITTING_FIG,o_fn,$NDVI1,LST1,NDVIm, $LST_MIN,LST_MAX, $A1,B1,R1,A2,B2,R2ENDIF

6 干濕邊方程擬合函數（參考：徐永明老師的課本）

PRO PLOT_FITTING_FIG,OUT_NAME,NDVI,LST,NDVIm, $LST_MIN,LST_MAX,A1,B1,R1,A2,B2,R2xrange=[0.2,1]y_scale=max(lst_max)-min(lst_min)yrange=[min(lst_min)-y_scale*0.5,max(lst_max)+y_scale*0.5]plot_scatter=plot(NDVIm,LST_MAX,LINESTYLE=6,$XRANGE=xrange,YRANGE=yrange,XTITLE='NDVI',YTITLE='LST(K)',$/BUFFER,XMINOR=5,YMINOR=5);MARGIN=[0.2,0.2,0.2,0.2],print,'scatter'PLOT_PT_DRY=PLOT(NDVIm,LST_MAX,SYMBOL=2,SYM_SIZE=0.8, $/SYM_FILLED,COLOR='RED',LINESTYLE=6,XTITLE='NDVI',YTITLE='LST(K)',/OVERPLOT,/BUFFER)PLOT_PT_WET=PLOT(NDVIm,LST_MIN,SYMBOL=2,SYM_SIZE=0.8, $/SYM_FILLED,COLOR='BLUE',LINESTYLE=6,XTITLE='NDVI',YTITLE='LST(K)',/OVERPLOT,/BUFFER)X=[MIN(NDVIm),Max(NDVIM)]PLOT_LINE_DRY=PLOT(X,Y_DRY,COLOR='Black',THICK=1.5,XTITLE='NDVI',YTITLE='LST(K)',/OVERPLOT, $/BUFFER)PLOT_LINE_WET=PLOT(X,Y_WET,COLOR='Black',THICK=1.5,XTITLE='NDVI',YTITLE='LST(K)',/OVERPLOT, $/BUFFER)EQUATION_DRY='Ts = '+STRING(A1,FORMAT='(F6.2)')+'*NDVI+'+STRING(B1,FORMAT='(F6.2)')r1_dry='$R^2$='+STRING(R1^2,FORMAT='(F6.2)')EQUATION_WET='Ts = '+STRING(A2,FORMAT='(F6.2)')+'*NDVI+'+STRING(B2,FORMAT='(F6.2)')r2_wet='$R^2$='+STRING(R2^2,FORMAT='(F6.2)')plot_scatter.save,out_name;,border=0plot_pt_dry.save,out_name,border=0plot_pt_wet.save,out_name,border=0plot_line_dry.save,out_name,border=0plot_line_wet.save,out_name,border=0 END

參考文獻

遙感二次開發語言IDL_徐永明

總結

以上是生活随笔為你收集整理的基于ENVI/IDL 的一键化实现LST-NDVI的干湿边方程拟合，并得到TVDI计算结果图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。