青藏高原冬季裸土區(qū) MODIS 溫度反演的完整可運(yùn)行 IDL 代碼

青藏高原冬季裸土區(qū)MODIS溫度反演的完整可運(yùn)行IDL代碼青藏高原冬季裸土區(qū)MODIS溫度反演的完整可運(yùn)行IDL代碼（含MODTRAN大氣校正、海拔修正、冬季參數(shù)適配），以及ENVI中批量處理MODIS溫度反演的實(shí)操流程，覆蓋從數(shù)據(jù)讀取到精度驗(yàn)證的全流程，可直接適配青藏高原冬季場(chǎng)景：一、青藏高原冬季裸土區(qū)MODISLST反演（完整IDL代碼）代碼說(shuō)明適配場(chǎng)景：青藏高原冬季（11-次年2月）、裸土/凍土主導(dǎo)區(qū)域；核心優(yōu)化：海拔氣壓修正、冬季干燥大氣參數(shù)、裸土比輻射率適配、迭代收斂保障；輸入數(shù)據(jù)：MOD021KM（輻射亮度）、MOD03（地理定位）、MOD07（大氣廓線(xiàn)）；輸出數(shù)據(jù)：地表溫度影像（℃）、比輻射率影像、精度驗(yàn)證日志。idl;=====================================環(huán)境初始化=====================================PROQinghai_Tibet_Winter_LST;1.設(shè)置路徑（替換為你的數(shù)據(jù)/輸出路徑）input_path='D:/MODIS_QTP/Winter/';輸入數(shù)據(jù)路徑output_path='D:/MODIS_QTP/Result/';輸出結(jié)果路徑modtran_path='C:/MODTRAN6/';MODTRAN安裝路徑SETENV,'MODTRAN_DIR='+modtran_path;關(guān)聯(lián)MODTRAN路徑FILE_MKDIR,output_path;創(chuàng)建輸出目錄;2.讀取文件列表files_mod02=FILE_SEARCH(input_path,'MOD021KM*.hdf')IFN_ELEMENTS(files_mod02)EQ0THENBEGINPRINT,'未找到MOD021KM數(shù)據(jù)！'RETURNENDIF;3.初始化精度驗(yàn)證日志log_file=FOPEN(output_path+'LST_Validation.log',/WRITE)FPRINT,log_file,'青藏高原冬季裸土區(qū)LST反演日志|日期：'+SYSTIME()FPRINT,log_file,'==========================================';=====================================批量處理=====================================FORf_idx=0,N_ELEMENTS(files_mod02)-1DOBEGINfile_mod02=files_mod02[f_idx]file_name=FILE_BASENAME(file_mod02,'.hdf')PRINT,'處理：',file_name;-------------------------步驟1：讀取MOD021KM數(shù)據(jù)（熱紅外+可見(jiàn)光）-------------------------hdf_id=HDF_OPEN(file_mod02,/READ,ERROR=err)IFerrNE0THENBEGINFPRINT,log_file,'錯(cuò)誤：'+file_name+'-無(wú)法打開(kāi)HDF文件'CONTINUEENDIF;熱紅外Band31（索引10）、Band32（索引11）DN值+定標(biāo)系數(shù)data31=HDFSD_GETDATA(hdf_id,'EV_1KM_Emissive',0,10)data32=HDFSD_GETDATA(hdf_id,'EV_1KM_Emissive',0,11)gain31=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','gain',10)offset31=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','offset',10)gain32=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','gain',11)offset32=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','offset',11);可見(jiàn)光Band1（紅，索引0）、Band2（近紅外，索引1）DN值+定標(biāo)系數(shù)data_b1=HDFSD_GETDATA(hdf_id,'EV_1KM_Reflective',0,0)data_b2=HDFSD_GETDATA(hdf_id,'EV_1KM_Reflective',0,1)gain_b1=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','gain',0)offset_b1=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','offset',0)gain_b2=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','gain',1)offset_b2=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','offset',1)HDF_CLOSE,hdf_id;-------------------------步驟2：輻射定標(biāo)-------------------------;熱紅外：DN→表觀輻射亮度（W·m^-2·sr^-1·μm^-1）L31=data31*gain31+offset31L32=data32*gain32+offset32;可見(jiàn)光：DN→表觀反射率rho_b1=data_b1*gain_b1+offset_b1rho_b2=data_b2*gain_b2+offset_b2;-------------------------步驟3：提取青藏高原大氣/地理參數(shù)-------------------------;3.1從MOD07提取大氣參數(shù)（冬季水汽含量0.5~1.5g/cm2）file_mod07=FILE_SEARCH(input_path,'MOD07*'+STRMID(file_name,9,13)+'*.hdf')IFN_ELEMENTS(file_mod07)GT0THENBEGINhdf07_id=HDF_OPEN(file_mod07[0],/READ)water_vapor=HDFSD_GETDATA(hdf07_id,'Water_Vapor',0,0)water_vapor=MEAN(water_vapor,/NAN);研究區(qū)平均水汽含量HDF_CLOSE,hdf07_idENDIFELSEBEGINwater_vapor=1.0;無(wú)MOD07時(shí)用經(jīng)驗(yàn)值（青藏高原冬季）ENDELSEwater_vapor=MIN([MAX([water_vapor,0.5]),1.5]);限制范圍;3.2從MOD03提取海拔+觀測(cè)天頂角（海拔修正）file_mod03=FILE_SEARCH(input_path,'MOD03*'+STRMID(file_name,9,13)+'*.hdf')IFN_ELEMENTS(file_mod03)GT0THENBEGINhdf03_id=HDF_OPEN(file_mod03[0],/READ)altitude=HDFSD_GETDATA(hdf03_id,'Height',0,0);地表海拔（km）zenith=HDFSD_GETDATA(hdf03_id,'SensorZenith',0,0);傳感器天頂角HDF_CLOSE,hdf03_id;海拔修正天頂角：zenith_corr=zenith-(altitude/6371)*180/!PIzenith_corr=zenith-(altitude/6371)*180/!PIzenith_mean=MEAN(zenith_corr,/NAN);平均天頂角;海拔修正氣壓（標(biāo)準(zhǔn)大氣壓1013.25hPa，海拔4km≈600hPa）pressure=1013.25*EXP(-altitude/8.4);氣壓經(jīng)驗(yàn)公式ENDIFELSEBEGINzenith_mean=30.0;無(wú)MOD03時(shí)默認(rèn)30°pressure=600.0;青藏高原平均氣壓ENDELSE;-------------------------步驟4：MODTRAN大氣校正（冬季大氣模式）-------------------------;4.1配置MODTRAN參數(shù)modtran_pro,/INIT,$ATMOSPHERE='SUBARCTIC_WINTER',$;青藏高原冬季大氣模式ALTITUDE=0,$;地表海拔（后續(xù)修正）ZENITH=zenith_mean,$;修正后天頂角AZIMUTH=0,$WATER_VAPOR=water_vapor,$;水汽含量WAVELENGTH=[11.03],$;Band31中心波長(zhǎng)（μm）/THERMAL,$OUTPUT_FILE=output_path+'modtran_'+file_name+'.dat';4.2運(yùn)行MODTRAN模擬modtran_pro,/RUN,ERROR=mod_errIFmod_errNE0THENBEGINFPRINT,log_file,'錯(cuò)誤：'+file_name+'-MODTRAN運(yùn)行失敗'CONTINUEENDIF;4.3讀取MODTRAN結(jié)果（大氣上行/下行輻射、半球反射率）read_modtran,output_path+'modtran_'+file_name+'.dat',L_up,L_down,r_atm;4.4海拔氣壓修正地表輻射亮度L31_corr=L31*(1013.25/pressure);氣壓修正;計(jì)算地表真實(shí)輻射亮度Lsrho_31=0.02;青藏高原裸土熱紅外反射率（經(jīng)驗(yàn)值）Ls_31=L31_corr-L_up-(L_down*rho_31)/(1-rho_31*r_atm+1e-8);-------------------------步驟5：計(jì)算比輻射率（冬季裸土適配）-------------------------;5.1計(jì)算NDVI（掩膜植被/裸土）ndvi=(rho_b2-rho_b1)/(rho_b2+rho_b1+1e-8)ndvi[ndviLT-1ORndviGT1]=!VALUES.F_NAN;5.2植被覆蓋度（冬季裸土為主，NDVI閾值調(diào)整）ndvi_min=0.0;冬季裸土NDVIndvi_max=0.5;冬季稀疏植被NDVIfv=(ndvi-ndvi_min)/(ndvi_max-ndvi_min+1e-8)fv[fvLT0]=0.0fv[fvGT1]=1.0fv[WHERE(FINITE(fv)EQ0)]=0.0;NaN設(shè)為裸土;5.3裸土比輻射率（青藏高原冬季凍土）eps_v=0.982;高寒草甸eps_s=0.970;凍土裸土d_eps=0.0005;混合修正項(xiàng)eps_31=fv*eps_v+(1-fv)*eps_s+d_eps;水體掩膜（NDVI<0且亮溫>273.15）T31=14388/(11.03*ALOG(3.7418e16/(11.03^5*L31)+1));星上亮溫water_mask=(ndviLT0)AND(T31GT273.15)eps_31[water_mask]=0.995;-------------------------步驟6：普朗克迭代反演LST（冬季收斂?jī)?yōu)化）-------------------------c1=3.7418e-16;W·m2c2=1.4388e-2;m·Klambda31=11.03e-6;mTs=FLTARR(SIZE(Ls_31,/DIMENSIONS))FORi=0,N_ELEMENTS(Ls_31)-1DOBEGINIF~FINITE(Ls_31[i])OR~FINITE(eps_31[i])ORLs_31[i]<=0THENBEGINTs[i]=!VALUES.F_NANCONTINUEENDIF;初始值（冬季限制200~280K）T0=c2/(lambda31*ALOG(c1/(lambda31^5*L31[i])+1))T0=MIN([MAX([T0,200.0]),280.0])T1=T0delta=1.0iter_num=0WHILEABS(delta)GT0.01ANDiter_numLT50DOBEGINLb=c1/(lambda31^5*(EXP(c2/(lambda31*T1))-1))Ls_calc=eps_31[i]*LbdLdT=eps_31[i]*c1*c2/(lambda31^6*(EXP(c2/(lambda31*T1))-1)^2)*EXP(c2/(lambda31*T1))IFABS(dLdT)<1e-10THENBEGINdelta=0.0T1=!VALUES.F_NANBREAKENDIFdelta=(Ls_31[i]-Ls_calc)/dLdTdelta=MIN([MAX([delta,-1.0]),1.0]);單次修正±1KT1=T1+deltaiter_num=iter_num+1ENDWHILETs[i]=(iter_num>=50)?!VALUES.F_NAN:T1ENDFOR;-------------------------步驟7：結(jié)果處理與保存-------------------------Ts_C=Ts-273.15;轉(zhuǎn)換為℃;掩膜異常值（冬季裸土溫度范圍：-30~10℃）Ts_C[Ts_CLT-30ORTs_CGT10]=!VALUES.F_NAN;保存LST影像envi_map=ENVI_GET_MAP_INFO(FILE=file_mod02)ENVI_WRITE_ENVI_FILE,Ts_C,$FILENAME=output_path+file_name+'_LST.dat',$INTERLEAVE=0,$MAP_INFO=envi_map,$DESCRIPTION='青藏高原冬季裸土區(qū)LST（℃）_MODTRAN校正';保存比輻射率ENVI_WRITE_ENVI_FILE,eps_31,$FILENAME=output_path+file_name+'_EPS.dat',$INTERLEAVE=0,$MAP_INFO=envi_map,$DESCRIPTION='Band31比輻射率';-------------------------步驟8：精度驗(yàn)證（實(shí)測(cè)數(shù)據(jù)對(duì)比）-------------------------IFFILE_TEST(input_path+'QTP_meas_data.txt')THENBEGINmeas_data=READ_ASCII(input_path+'QTP_meas_data.txt',DATA_COLUMNS=[3])meas_T=meas_data.DATA[0,*]lat=meas_data.DATA[1,*]lon=meas_data.DATA[2,*];經(jīng)緯度轉(zhuǎn)像元坐標(biāo)ENVI_CONVERT_COORD,envi_map,lon,lat,x,y,/GEO_TO_IMAGEinv_T=FLTARR(N_ELEMENTS(x))FORj=0,N_ELEMENTS(x)-1DOBEGINix=FIX(x[j])iy=FIX(y[j])IFixGE0ANDixLT(SIZE(Ts_C,1))ANDiyGE0ANDiyLT(SIZE(Ts_C,2))THENBEGINinv_T[j]=Ts_C[iy,ix]ENDIFELSEBEGINinv_T[j]=!VALUES.F_NANENDELSEENDFOR;計(jì)算精度指標(biāo)valid_idx=WHERE(FINITE(inv_T)ANDFINITE(meas_T),n_valid)IFn_validGT0THENBEGINrmse=SQRT(MEAN((inv_T[valid_idx]-meas_T[valid_idx])^2))mae=MEAN(ABS(inv_T[valid_idx]-meas_T[valid_idx]))r2=CORRELATE(inv_T[valid_idx],meas_T[valid_idx])^2FPRINT,log_file,file_name+'-驗(yàn)證樣本數(shù)：'+STRTRIM(n_valid,2)FPRINT,log_file,'RMSE：'+FORMAT='(F6.2)',rmse+'℃|MAE：'+FORMAT='(F6.2)',mae+'℃|R2：'+FORMAT='(F6.3)',r2ENDIFELSEBEGINFPRINT,log_file,file_name+'-無(wú)有效驗(yàn)證樣本'ENDELSEENDIF;釋放內(nèi)存HEAP_FREE,data31,data32,L31,L32,Ts,Ts_C,eps_31,Ls_31PRINT,file_name+'處理完成'ENDFOR;=====================================收尾=====================================FPRINT,log_file,'=========================================='FCLOSE,log_filePRINT,'全部處理完成！結(jié)果保存至：',output_pathEND;-------------------------輔助函數(shù)：讀取MODTRAN模擬結(jié)果-------------------------PROread_modtran,modtran_file,L_up,L_down,r_atmL_up=0.0L_down=0.0r_atm=0.0IF~FILE_TEST(modtran_file)THENRETURN;讀取MODTRAN輸出文件（提取大氣上行/下行輻射）lun=GET_LUN()OPENR,lun,modtran_file,/GET_LINEline=''WHILE~EOF(lun)DOBEGINREADF,lun,lineIFSTRPOS(line,'UPWARDRADIANCE')NE-1THENBEGINREADF,lun,lineL_up=FLOAT(STRMID(line,20,15));提取上行輻射值ENDIFIFSTRPOS(line,'DOWNWARDRADIANCE')NE-1THENBEGINREADF,lun,lineL_down=FLOAT(STRMID(line,20,15));提取下行輻射值ENDIFIFSTRPOS(line,'ATMOSPHERICALBEDO')NE-1THENBEGINREADF,lun,liner_atm=FLOAT(STRMID(line,20,15));提取大氣半球反射率ENDIFENDWHILEFREE_LUN,lun;單位轉(zhuǎn)換（MODTRAN輸出單位→W·m^-2·sr^-1·μm^-1）L_up=L_up*1e-6L_down=L_down*1e-6END代碼使用說(shuō)明路徑配置：修改input_path（輸入數(shù)據(jù)路徑）、output_path（輸出路徑）、modtran_path（MODTRAN安裝路徑）；輸入數(shù)據(jù)：MOD021KM、MOD03、MOD07文件需放在同一目錄，命名格式需匹配（如MOD021KM.A2025001.0000.061...和MOD03.A2025001.0000.061...）；實(shí)測(cè)數(shù)據(jù)文件QTP_meas_data.txt格式：經(jīng)度緯度實(shí)測(cè)溫度(℃)，每行一個(gè)測(cè)點(diǎn)；運(yùn)行方式：IDL中編譯代碼后，輸入Qinghai_Tibet_Winter_LST回車(chē)運(yùn)行；輸出文件：*_LST.dat：地表溫度影像（℃）；*_EPS.dat：Band31比輻射率影像；LST_Validation.log：精度驗(yàn)證日志（RMSE/MAE/R2）。二、ENVI中批量處理MODIS溫度反演實(shí)操流程（青藏高原適配）前置準(zhǔn)備安裝插件：ENVI5.3+安裝MODISToolkit+MODTRANLink+BatchProcessingTool；數(shù)據(jù)整理：將MOD021KM/MOD03/MOD07按日期分文件夾，確保同日期數(shù)據(jù)命名匹配；參數(shù)預(yù)設(shè)：新建文本文檔QTP_Winter_Params.txt，寫(xiě)入青藏高原冬季參數(shù)：plaintextATMOSPHERE=SUBARCTIC_WINTERWATER_VAPOR=1.0ZENITH_CORR=TRUEEPS_S=0.970EPS_V=0.982NDVI_MIN=0.0NDVI_MAX=0.5批量處理步驟步驟1：批量格式轉(zhuǎn)換（HDF→ENVI）打開(kāi)ENVI→工具箱→

BatchProcessing→BatchTool；選擇InputFileList：添加所有MOD021KM.hdf文件；選擇Process：RasterManagement→ModisHDFtoENVI；設(shè)置參數(shù)：勾選BandSelection：Band1、Band2、Band31、Band32；OutputDirectory：設(shè)置轉(zhuǎn)換后輸出目錄；勾選AddGeolocation：選擇對(duì)應(yīng)MOD03文件；點(diǎn)擊RunBatch，批量轉(zhuǎn)換為ENVI格式。步驟2：批量云掩膜