青藏高原冬季裸土区 MODIS 温度反演的完整可运行 IDL 代码

青藏高原冬季裸土区MODIS温度反演的完整可运行IDL代码青藏高原冬季裸土区MODIS温度反演的完整可运行IDL代码（含MODTRAN大气校正、海拔修正、冬季参数适配），以及ENVI中批量处理MODIS温度反演的实操流程，覆盖从数据读取到精度验证的全流程，可直接适配青藏高原冬季场景：一、青藏高原冬季裸土区MODISLST反演（完整IDL代码）代码说明适配场景：青藏高原冬季（11-次年2月）、裸土/冻土主导区域；核心优化：海拔气压修正、冬季干燥大气参数、裸土比辐射率适配、迭代收敛保障；输入数据：MOD021KM（辐射亮度）、MOD03（地理定位）、MOD07（大气廓线）；输出数据：地表温度影像（℃）、比辐射率影像、精度验证日志。idl;=====================================环境初始化=====================================PROQinghai_Tibet_Winter_LST;1.设置路径（替换为你的数据/输出路径）input_path='D:/MODIS_QTP/Winter/';输入数据路径output_path='D:/MODIS_QTP/Result/';输出结果路径modtran_path='C:/MODTRAN6/';MODTRAN安装路径SETENV,'MODTRAN_DIR='+modtran_path;关联MODTRAN路径FILE_MKDIR,output_path;创建输出目录;2.读取文件列表files_mod02=FILE_SEARCH(input_path,'MOD021KM*.hdf')IFN_ELEMENTS(files_mod02)EQ0THENBEGINPRINT,'未找到MOD021KM数据！'RETURNENDIF;3.初始化精度验证日志log_file=FOPEN(output_path+'LST_Validation.log',/WRITE)FPRINT,log_file,'青藏高原冬季裸土区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数据（热红外+可见光）-------------------------hdf_id=HDF_OPEN(file_mod02,/READ,ERROR=err)IFerrNE0THENBEGINFPRINT,log_file,'错误：'+file_name+'-无法打开HDF文件'CONTINUEENDIF;热红外Band31（索引10）、Band32（索引11）DN值+定标系数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);可见光Band1（红，索引0）、Band2（近红外，索引1）DN值+定标系数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：辐射定标-------------------------;热红外：DN→表观辐射亮度（W·m^-2·sr^-1·μm^-1）L31=data31*gain31+offset31L32=data32*gain32+offset32;可见光：DN→表观反射率rho_b1=data_b1*gain_b1+offset_b1rho_b2=data_b2*gain_b2+offset_b2;-------------------------步骤3：提取青藏高原大气/地理参数-------------------------;3.1从MOD07提取大气参数（冬季水汽含量0.5~1.5g/cm²）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);研究区平均水汽含量HDF_CLOSE,hdf07_idENDIFELSEBEGINwater_vapor=1.0;无MOD07时用经验值（青藏高原冬季）ENDELSEwater_vapor=MIN([MAX([water_vapor,0.5]),1.5]);限制范围;3.2从MOD03提取海拔+观测天顶角（海拔修正）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);平均天顶角;海拔修正气压（标准大气压1013.25hPa，海拔4km≈600hPa）pressure=1013.25*EXP(-altitude/8.4);气压经验公式ENDIFELSEBEGINzenith_mean=30.0;无MOD03时默认30°pressure=600.0;青藏高原平均气压ENDELSE;-------------------------步骤4：MODTRAN大气校正（冬季大气模式）-------------------------;4.1配置MODTRAN参数modtran_pro,/INIT,$ATMOSPHERE='SUBARCTIC_WINTER',$;青藏高原冬季大气模式ALTITUDE=0,$;地表海拔（后续修正）ZENITH=zenith_mean,$;修正后天顶角AZIMUTH=0,$WATER_VAPOR=water_vapor,$;水汽含量WAVELENGTH=[11.03],$;Band31中心波长（μm）/THERMAL,$OUTPUT_FILE=output_path+'modtran_'+file_name+'.dat';4.2运行MODTRAN模拟modtran_pro,/RUN,ERROR=mod_errIFmod_errNE0THENBEGINFPRINT,log_file,'错误：'+file_name+'-MODTRAN运行失败'CONTINUEENDIF;4.3读取MODTRAN结果（大气上行/下行辐射、半球反射率）read_modtran,output_path+'modtran_'+file_name+'.dat',L_up,L_down,r_atm;4.4海拔气压修正地表辐射亮度L31_corr=L31*(1013.25/pressure);气压修正;计算地表真实辐射亮度Lsrho_31=0.02;青藏高原裸土热红外反射率（经验值）Ls_31=L31_corr-L_up-(L_down*rho_31)/(1-rho_31*r_atm+1e-8);-------------------------步骤5：计算比辐射率（冬季裸土适配）-------------------------;5.1计算NDVI（掩膜植被/裸土）ndvi=(rho_b2-rho_b1)/(rho_b2+rho_b1+1e-8)ndvi[ndviLT-1ORndviGT1]=!VALUES.F_NAN;5.2植被覆盖度（冬季裸土为主，NDVI阈值调整）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设为裸土;5.3裸土比辐射率（青藏高原冬季冻土）eps_v=0.982;高寒草甸eps_s=0.970;冻土裸土d_eps=0.0005;混合修正项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（冬季收敛优化）-------------------------c1=3.7418e-16;W·m²c2=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：结果处理与保存-------------------------Ts_C=Ts-273.15;转换为℃;掩膜异常值（冬季裸土温度范围：-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='青藏高原冬季裸土区LST（℃）_MODTRAN校正';保存比辐射率ENVI_WRITE_ENVI_FILE,eps_31,$FILENAME=output_path+file_name+'_EPS.dat',$INTERLEAVE=0,$MAP_INFO=envi_map,$DESCRIPTION='Band31比辐射率';-------------------------步骤8：精度验证（实测数据对比）-------------------------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,*];经纬度转像元坐标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;计算精度指标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+'-验证样本数：'+STRTRIM(n_valid,2)FPRINT,log_file,'RMSE：'+FORMAT='(F6.2)',rmse+'℃|MAE：'+FORMAT='(F6.2)',mae+'℃|R²：'+FORMAT='(F6.3)',r2ENDIFELSEBEGINFPRINT,log_file,file_name+'-无有效验证样本'ENDELSEENDIF;释放内存HEAP_FREE,data31,data32,L31,L32,Ts,Ts_C,eps_31,Ls_31PRINT,file_name+'处理完成'ENDFOR;=====================================收尾=====================================FPRINT,log_file,'=========================================='FCLOSE,log_filePRINT,'全部处理完成！结果保存至：',output_pathEND;-------------------------辅助函数：读取MODTRAN模拟结果-------------------------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;单位转换（MODTRAN输出单位→W·m^-2·sr^-1·μm^-1）L_up=L_up*1e-6L_down=L_down*1e-6END代码使用说明路径配置：修改input_path（输入数据路径）、output_path（输出路径）、modtran_path（MODTRAN安装路径）；输入数据：MOD021KM、MOD03、MOD07文件需放在同一目录，命名格式需匹配（如MOD021KM.A2025001.0000.061...和MOD03.A2025001.0000.061...）；实测数据文件QTP_meas_data.txt格式：经度纬度实测温度(℃)，每行一个测点；运行方式：IDL中编译代码后，输入Qinghai_Tibet_Winter_LST回车运行；输出文件：*_LST.dat：地表温度影像（℃）；*_EPS.dat：Band31比辐射率影像；LST_Validation.log：精度验证日志（RMSE/MAE/R²）。二、ENVI中批量处理MODIS温度反演实操流程（青藏高原适配）前置准备安装插件：ENVI5.3+安装MODISToolkit+MODTRANLink+BatchProcessingTool；数据整理：将MOD021KM/MOD03/MOD07按日期分文件夹，确保同日期数据命名匹配；参数预设：新建文本文档QTP_Winter_Params.txt，写入青藏高原冬季参数：plaintextATMOSPHERE=SUBARCTIC_WINTERWATER_VAPOR=1.0ZENITH_CORR=TRUEEPS_S=0.970EPS_V=0.982NDVI_MIN=0.0NDVI_MAX=0.5批量处理步骤步骤1：批量格式转换（HDF→ENVI）打开ENVI→工具箱→

BatchProcessing→BatchTool；选择InputFileList：添加所有MOD021KM.hdf文件；选择Process：RasterManagement→ModisHDFtoENVI；设置参数：勾选BandSelection：Band1、Band2、Band31、Band32；OutputDirectory：设置转换后输出目录；勾选AddGeolocation：选择对应MOD03文件；点击RunBatch，批量转换为ENVI格式。步骤2：批量云掩膜