分享

ENVI正射校正、统计大气校正、地形起伏度与影像清晰度计算

 米木十一 2011-10-15

ENVI正射校正扩展模块介绍

 

   ENVI正射校正扩展模块随着ENVI4.6.1 在3月份正式发布。

产品概述

  ENVI的正射校正扩展模块采用瑞典Spacemetric公司研制的精确正射校正数学模型,Spacementric公司与卫星和航空的影像数据供应商具有多年合作的经验,能够根据目前最流行的卫星和航天传感器来优化他们的模型。
      ENVI正射校正扩展模块一次可以完成大区域、若干景影像和多传感器的正射校正,并能以镶嵌结果的方式输出,提供接边线、颜色平衡等工具。整个模块是基于流程化的向导式操作方式和工程化管理,经过五个步骤,即输入数据、收集控制点、选择连接点(Ties)、确定接边线和输出结果,即使非专家用户也能轻松快速地得到精确的结果。具有灵活的功能扩展,包括自定义传感器模型和批处理。
但与ArcGIS的整合时,您可以无缝地将结果导入ArcGIS平台中而不需要中断你的工作流程。

ENVI正射校正、统计大气校正、地形起伏度与影像清晰度计算 11.JPG

 


 

基于流程化的向导式操作方式和工程化管理



支持传感器模型

ENVI正射校正扩展模块支持大多数传感器模型,而且还支持自定义传感器模型。

l卫星传感器
ASTER, Landsat, IKONOS, EROS, QuickBird, WorldView-1, SPOT, ERS, RADARSAT-1, IRS/ RESOURCESAT(P6), CARTOSAT-1(P5), FORMOSAT-2, OrbView/GeoEye-1, KOMPSAT, UK-DMC, AlgeriaSat-1, NigeriaSat-1, Beijing-1
l航空传感器
Leica ADS, Intergraph, Z/I, ZI/DMC, Vexcel UltraCam, Trimble/Applanix DSS, Generic Frame Camera, Generic Pushbroom Camera

产品特性

l采用的正射校正方法具有可靠和高精度的特点,并且该方法被行业所认可。
l支持大区域范围内的多幅影像、多传感器的一次正射校正。
l具有镶嵌结果的功能,并提供接边线和颜色平衡辅助工具。
l采用流程化的向导式操作方式和工程化管理。
l自定义传感器模型
l提供接口函数,便于扩展功能。

 

 

 

 

 

 

 

 

 

 

ENVI中基于统计学模型的大气校正方法详解

 

     基于统计学模型的反射率反演的方法主要有平场域法(Flat FieldFF)、对数残差法(Log Residuals)、内部平均法(Internal Average Relative ReflectanceIARR)、经验线性法(Empirical Line)。集中在Basic Tools->Preprocessing-> Calibration Utilities菜单下。

1.       平场域定标(Flat Field Calibration)

Flat Field定标工具通过选择图像中一块具有高反射率、光谱变化平坦的区域,利用这个区域的平均光谱值来模拟飞行时的大气条件下的太阳光谱。将每个像元的DN值除以选择区域的平均光谱值得到相对反射率,以此来消除大气的影响。

在使用这个工具前,需要利用ENVI提供的感兴趣区绘制工具(ROI Tool)在被定标图像上选择感兴趣区作为平场域(Flat Field),感兴趣区可选择沙漠、大块水泥地、沙地等区域。操作过程如下:

(1)           在主菜单中,选择Basic Tools->Preprocessing-> Calibration Utilities-> Flat Field在打开的Calibration Input File对话框中,选择输入文件,单击“OK

(2)           在打开的Flat Field Calibration Parameters面板中(图13.4),在标有“Select ROI for  Calibration”一栏中,选择感兴趣区(只能选择一个),作为平场域定标的平均波谱区。

(3)           选择输出路径及文件名,单击“OK”执行定标处理。

[转载]ENVI中基于统计学模型的大气校正方法详解
13.4 Flat Field Calibration Parameters面板

2.       对数残差(Log  Residuals

对数残差定标工具将数据除以波段几何均值,后再除以像元几何均值,可以消除光照、大气传输、仪器系统误差、地形影响和星体反照率对数据辐射的影响。定标结果的值在1附近。操作过程如下:

(1)           在主菜单中,选择Basic Tools->Preprocessing->Calibration Utilities->Log Residuals,在Log Residuals Calibration Input File对话框中,选择输入文件,单击“OK”。

(2)           在打开的Log Residuals Calibration Parameters面板中,选择输出路径及文件名,单击OK执行定标处理

3.       内部平均(IAR)反射率定标

IAR (Internal Average Relative) Reflectance定标工具假定整幅图像的平均光谱基本代表了大气影响下的太阳光谱信息。把图像DN值与整幅图像的平均辐射光谱值相除,得到的结果为相对反射率。该工具特别适用于没有植被的干旱区域。操作过程如下:

(1)           在主菜单中,选择Basic Tools->Preprocessing->Calibration Utilities-> IAR Reflectance在打开的Calibration Input File对话框中,选择输入文件,单击“OK”。

(2)           在打开的IARR Calibration Parameters面板中,选择输出到“File”或“Memory”。点击“OK”执行定标处理。

4.       经验线性定标(Empirical Line Calibration)

Empirical Line 定标方法是假设图像DN值与反射率之间存在线性关系:

反射率= 增益 * DN+ 偏移

利用两个已知点的地面反射光谱值,再计算图像上对应像元点的平均DN值,然后利用线性回归求出增益和偏移值,建立DN值与反射率之间的相互关系式,进行反射率的定标。消除了太阳辐亮度和大气程辐射。

    ENVIEmpirical Line定标工具要求至少需要一个已知区域的地面反射光谱值(Field Spectra)作为参照波谱,以及图像上对应像元点的波谱曲线(Data Spectra)。它们可以来自波谱剖面或波谱曲线、波谱库、感兴趣区、统计文件和ASCII文件。输入的波谱将自动被重采样,以与选择的数据波长相匹配。也可以用已经存在的系数对数据集进行定标。

1)      计算系数并定标(Compute Factors and Calibrate

Empirical Line 定标工具时,一般可以在图像上选择一个暗区和一个亮区作为已知区域(假定这些区域中的参照波谱是可以获得的)。使用越多的已知波谱也可以提高定标精度,至少需要一组已知区域的波谱。操作过程如下:

(1)           选择Basic Tools->Preprocessing-> Calibration Utilities-> Empirical Line-> Compute Factors and Calibrate。在打开的Empirical Line Input File对话框中,选择输入文件。单击“OK”,

(2)           在打开的Empirical Line Spectra面板中(图13.5),需要选择地面反射光谱值(Field Spectra)以及图像上对应像元点的波谱值(Data Spectra)。

 

[转载]ENVI中基于统计学模型的大气校正方法详解13.5 Empirical Line Spectra面板

l  选择数据(图像)波谱

Empirical Line Spectra对话框中,单击“Import Spectra”按钮,打开Data Spectral Collection对话框(图13.6)。

Data Spectral Collection对话框中,选择Import->from ROI/EVF from input file,选择定义好的感兴趣区文件,点击“Apply”,波谱名被输入到Empirical Line Spectra面板中。点击“Cancel”,关闭Data Spectra Collection对话框。

 

[转载]ENVI中基于统计学模型的大气校正方法详解

13.6 Data Spectral Collection对话框

l  选择参照波谱

Empirical Line Spectra对话框中,单击“Import Spectra”按钮,打开Feild Spectra Collection对话框,这个对话框与Data Spectra Collection对话框类似。

Feild Spectra Collection对话框中,选择Import->from ASD binary file,选择对应图像波谱区域用ASD波谱仪测量波谱文件。点击“Apply”,波谱名被输入到Empirical Line Spectra面板中。点击“Cancel”,关闭Feild Spectra Collection对话框。

(3)           回到Empirical Line Spectra面板中,在顶部的列表内点击波谱名选择数据波谱。在底部的列表中,点击相应的参照波谱名,单击击“Enter Pair”按钮使两个波谱相关联,相关联的波谱将被列在“Selected Pairs”文本框中。(说明:如果错误的选择了关联波谱,在“Selected Pairs”文本框中单击关联波谱可移除)

(4)           重复(3)(4)步骤可选择关联波谱。

(5)           单击击“OK”。打开Empirical Line Calibration Parameters对话框。

(6)           Empirical Line Calibration Parameters对话框中,选择输出定标结果文件路径及文件名,将定标系数保存在ASCII文件中,在“Output Calibration Filename”文本框中键入第二个文件名,定标系数文件的默认扩展名是 .cff

(7)           单击OK执行定标过程。

2)      使用现有系数定标(Calibrating Using Existing Factors)

Calibrate Using Existing Factors工具可以使用另一个定标过程中存储的纠正系数来运行经验行定标功能。

(1)           在主菜单中,选择 Basic Tools- >Preprocessing-> Calibration Utilities-> Empirical Line->Calibrate Using Existing Factors。选择输入文件,点击“OK”。

(2)           在打开的Enter Calibration Factors Filename对话框中,选择一个之前定标过程中创建的定标系数文件(.cff)。点击“OK”。

在打开的Empirical Line Calibration Parameters对话框中,选择输出路径及文件名。单击“OK”执行定标过程。

 

摘自《ENVI遥感图像处理方法》科学出版社

 

 

 

 

 

利用IDL实现ENVI下计算地形起伏度

 

地形起伏度:是指在一个特定的区域内,最高点海拔高度与最低点海拔高度的差值。它是描述一个区域地形特征的一个宏观性的指标。

    在很多具有栅格分析的软件中都没有提供计算地形起伏度的功能,虽然可以根据现有的工具进行组合计算,但是还是比较麻烦。ENVI/IDL的栅格运算功能强大,效率高。本文利用IDL制作了计算地形起伏度的程序,并且能集成在ENVI中使用。

注意事项:

    将dixingqifudu.pro文件放到ENVI的安装目录\ProgramFiles\ITT\IDLxx\products\enviXX\save_add下,启动ENVI+IDL,出现在ENVI的topographic菜单下:(或者编译成sav文件,只打开ENVI即可。)

利用IDL实现ENVI下计算地形起伏度        利用IDL实现ENVI下计算地形起伏度

源码:

;;===界面居中的函数======================

PRO CENTERTLB, tlb, x, y, NoCenter=nocenter

  COMPILE_OPT StrictArr
 
  geom = WIDGET_INFO(tlb, /Geometry)
 
  IF N_ELEMENTS(x) EQ 0 THEN xc = 0.5 ELSE xc = FLOAT(x[0])
  IF N_ELEMENTS(y) EQ 0 THEN yc = 0.5 ELSE yc = 1.0 - FLOAT(y[0])
  center = 1 - KEYWORD_SET(nocenter)
  ;
  oMonInfo = OBJ_NEW('IDLsysMonitorInfo')
  rects = oMonInfo -> GetRectangles(Exclude_Taskbar=exclude_Taskbar)
  pmi = oMonInfo -> GetPrimaryMonitorIndex()
  OBJ_DESTROY, oMonInfo
 
  screenSize =rects[[2, 3], pmi]
 
  ; Get_Screen_Size()
  IF screenSize[0] GT 2000 THEN screenSize[0] = screenSize[0]/2 ; Dual monitors.
  xCenter = screenSize[0] * xc
  yCenter = screenSize[1] * yc
 
  xHalfSize = geom.Scr_XSize / 2 * center
  yHalfSize = geom.Scr_YSize / 2 * center
 
  XOffset = 0 > (xCenter - xHalfSize) < (screenSize[0] - geom.Scr_Xsize)
  YOffset = 0 > (yCenter - yHalfSize) < (screenSize[1] - geom.Scr_Ysize)
 
  WIDGET_CONTROL, tlb, XOffset=XOffset, YOffset=YOffset
END
;;;================================================

;自动添加envi菜单
PRO dixingqifudu_define_buttons, buttonInfo
  ENVI_DEFINE_MENU_BUTTON, buttonInfo, VALUE = '地形起伏度计算', $
    uValue = '', $
    event_pro ='dixingqifudu', $
    REF_VALUE = 'Topographic',position=2
end

pro qifudu_result, pp, flag, out_name=out_name ;pp是窗口大小,flag是存贮方式
  compile_opt idl2
   envi_batch_init
  pp=pp;窗口大小
  p1=(pp-1)/2;用于建立窗口大小
  envi_select,fid=fid,pos=pos,dims=dims

  if fid eq -1 then begin             

      envi_batch_exit                  

    return                            

  endif                                

  ENVI_FILE_QUERY, fid, bnames=bnames                                       

  nl=dims[4]-dims[3]+1
  ns=dims[2]-dims[1]+1
  nb=size(pos)
  nb=nb[1]
print,nb
  tfid=indgen(nb)
  posss=intarr(nb);因为后面生成的都是单波段文件及r_fid,所以都是0表示第一个波段。
  FOR g=0,nb-1 DO BEGIN
    data1 = ENVI_GET_DATA(DIMS=dims, FID=fid , POS=pos[g])
    inherit = envi_set_inheritance(fid, dims, pos[g], /full)
    data2=data1
    FOR i =p1,ns-p1-1 DO BEGIN

      FOR j =p1, nl-p1-1  DO BEGIN

        ;获取一个pp*pp的窗口
        t= data1[[i-p1]:[i+p1],[j-p1]:[j+p1]]
        tmax=max(t,min=tmin)
        data2[i,j]=tmax-tmin
      ENDfor
    ENDFOR
    
      ENVI_WRITE_ENVI_FILE, data2,inherit=inherit ,r_fid=t,out_name=out_name+bnames[pos[g]] ;,/NO_OPEN
     ENVI_FILE_MNG, id=t, /REMOVE
      tfid[g]=t  
  ENDfor
    if flag eq 1 then begin
       envi_doit, 'cf_doit', fid=tfid, pos=posss, dims=dims, $
       remove=0, out_name=out_name, r_fid=r_fid
    endif else begin
       envi_doit, 'cf_doit', fid=tfid, pos=posss, dims=dims, $
       remove=0, /IN_MEMORY, r_fid=r_fid
    endelse
 
end

;事件响应程序

pro dixingqifudu_event,event
  widget_control,event.top,get_uvalue=pstate

   if TAG_NAMES(event,/STRUCTURE_NAME) eq 'WIDGET_KILL_REQUEST' then begin     
      WIDGET_CONTROL,event.top,/destroy
      ptr_free,pstate
      RETURN
   endif
  
  widget_control,event.id,get_uvalue=uval
  case uval of
    'list':  begin
               (*pstate).chuangkou=fix(event.str)  
               print, fix(event.str)        
             if event.index eq -1 then a=fix(event.str)
             help,fix(event.str)
             end
    'file_m':begin
              if event.value eq '内存' then begin
               (*pstate).flag=0
               widget_control,(*pstate).i_base52,map=0
              endif else begin
               (*pstate).flag=1
               widget_control,(*pstate).i_base52,map=1
              endelse
             end 
    'xuanze':begin
              out_name=Dialog_Pickfile(Path=current, /NoConfirm, $
                            Get_Path=path, Filter=['*.*'],title='选择输出文件名')
                  if out_name eq '' then return
                  widget_control,(*pstate).dir_text,set_value=out_name
                  if (*pstate).flag eq 1 then (*pstate).out_name=out_name
             end
    'dir'   :begin             
              widget_control,event.id,get_value=y
              (*pstate).out_name=y
             end
    'queding':begin
   
              if ((*pstate).flag eq 1) and ((*pstate).out_name eq '') then begin
                  tem=DIALOG_MESSAGE('指定文件存储路径')
                  return
              endif
              qifudu_result,(*pstate).chuangkou, (*pstate).flag, out_name=(*pstate).out_name
              envi_batch_exit
             end
    'quxiao':begin
              widget_control,event.top,/DESTROY
              ptr_free,pstate
             end       
    else:
  endcase
end

pro dixingqifudu,event
compile_opt idl2
envi,/restore_base_save_files
envi_select,fid=fid,pos=pos,dims=dims,/band_only
print,dims

 

  tlb=widget_base(title='地形起伏度计算',/COLUMN)
  tbase0=widget_base(tlb,/COLUMN,/frame)
  tbase1=widget_base(tbase0,/row)
  lab1=widget_label(tbase1,value='窗口大小(奇数):    ')
  valu=['3','5','7','9','11']
  list=WIDGET_COMBOBOX(tbase1,value=valu,/EDITABLE,xsize=100,uvalue='list')
  lab2=widget_label(tbase0,value='说明:可键入参数,回车。')
lab2=widget_label(tbase0,value='',ysize=10)
  lab2=widget_label(tlb,value='',ysize=10)
 
  i_base5=WIDGET_BASE(tlb,/COLUMN,xsize=230,xpad=4,/FRAME)
  k_lable2=widget_label(i_base5,value='选择保存方式:',/ALIGN_LEFT)
  file_m= CW_BGROUP(i_base5, /row, /EXCLUSIVE, /NO_REL, /RETURN_name,UVALUE='file_m',$
    ['文件','内存'], SET_VALUE=0)
  i_base52=WIDGET_BASE(i_base5,xsize=250,/row)
  label3=widget_label(i_base52,value='输出文件',/align_left)
  xuanze=widget_button(i_base52,value='选择',uvalue='xuanze')
  dir_text=WIDGET_TEXT(i_base52,ysize=1,/EDITABLE,UVALUE='dir',/ALL_EVENTS)
  ibase6=widget_base(tlb,/row,/ALIGN_CENTER)
  queding=WIDGET_BUTTON(ibase6,value='确定',uvalue='queding')
  queding=WIDGET_BUTTON(ibase6,value='取消',uvalue='quxiao')
  CENTERTLB, tlb
  widget_control,tlb,/realize
  pstate=ptr_new({dir_text:dir_text,$
    chuangkou:3,$
    flag:1,out_name:'',$
    i_base52:i_base52},/no_copy)
  widget_control,tlb,set_uvalue=pstate
  xmanager,'dixingqifudu',tlb,/NO_BLOCK
 
end

    以上的代码,包含如何搭建界面、事件响应、代码定义ENVI菜单、ENVI功能函数调用等,遥感业务化系统可以按照这种模式打造。

sav文件下载见:http://bbs./ESRI/viewthread.php?tid=87724&extra=

 

 

 

 

利用IDL+ENVI实现遥感影像清晰度计算

 

清晰度是图像细节边缘变化的敏锐程度。在图像细节的边缘处,光学密度或亮度随位置的变化越敏锐(变化快)、越剧烈(反差大),则细节的边缘就越清晰,可辨程度越高。其计算公式可采用改进后的点锐度法【8】来表示:  式中:M、N分别为图像的行数和列数,df为灰度变化幅值,dx为像元间的距离增量,a表示像素i周围的像素数。

利用IDL+ENVI实现遥感影像清晰度计算

PRO Articulation

  COMPILE_OPT idl2
    ENVI, /restore_base_save_files
    envi_batch_init
  infile = DIALOG_PICKFILE(title = '请选择要打开的文件:')
  ENVI_OPEN_FILE, infile, r_fid=fid
  IF (fid EQ -1) THEN BEGIN
    ENVI_BATCH_EXIT
    RETURN
  ENDIF
  ;获取数据的基本信息
  ENVI_FILE_QUERY, fid, dims=dims, nb=nb ,nl=nl,ns=ns
 
  proj = ENVI_GET_PROJECTION(FID = fid,pixel_size=Dx)
  print,dx
  dx=dx[0]
data1 = MAKE_ARRAY(nl,ns)
  FOR g=0,nb-1 DO BEGIN
    data1 = ENVI_GET_DATA(DIMS=dims, FID=fid , POS=g)

    ;
    ;;清晰度的计算
    ;

    Den = nl*ns

    FOR i =1,nl-2 DO BEGIN
      FOR j = 1,ns-2 DO BEGIN
        ;获取一个3*3的窗口
        t= data1[[i-1]:[i+1],[j-1]:[j+1]]
            help,t
        ;对每个窗口做计算
        d1=ABS((t[0,0]-t[1,1])/(Dx*1.414))
        d2=ABS(t[0,1]-t[1,1])/Dx*1.0
        d3=ABS(t[0,2]-t[1,1])/(Dx*1.414)
        d4=ABS(t[1,0]-t[1,1])/Dx*1.0
        d5=ABS(t[1,2]-t[1,1])/Dx*1.0
        d6=ABS(t[2,0]-t[1,1])/(Dx*1.414)
        d7=ABS(t[2,1]-t[1,1])/Dx*1.0
        d8=ABS(t[2,2]-t[1,1])/(Dx*1.414)
        ;计算每个窗口的df/dx
        dd=d1+d2+d3+d4+d5+d6+d7+d8
        
      ENDFOR
    ENDFOR
    outdata2 = dd/Den

   ;输出每个波段图像的清晰度
    PRINT,outdata2
  ENDFOR
END

个人觉得这个代码的优点在于:

       1、三窗口矩阵的获取,(----多窗口的获取)

       2、获取后指定一变量t,关于矩阵内部的运算直接针对t就行了,这样避免了操作窗口时数组下标混乱。

 

    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多