自动化SAR强度分析

在工作的 20 多年中，我处理过许多不同的数据类型和数据模式，但对合成孔径雷达 (SAR) 的了解并不多。直到最近几年，我发现人们对 SAR 数据的兴趣越来越大，这部分归功于免费提供的 Sentinel-1 数据以及更多商业化的高分辨率 SAR 小型卫星。因此，我决定是时候冒险一试，将 SAR 数据添加到我的技能包中。在学习ENVI SARscape 教程中的 Sentinel-1 强度分析时，我想到编写一个自动执行这些步骤的程序将是一个很好的学习练习。

SAR 强度分析用于多种应用，例如检测浮油、研究植被以及使用多个图像随时间监测植被。在提取多张强度图像时，自动提取强度产品特别有用，因为它可以提高处理速度并减少按钮点击次数。

我使用 IDL 8.8.2、ENVI5.6.2 和SARscape5.6.2编写了这个程序。数据可从如下链接获取，下载后解压到本地目录。

链接：https://pan.baidu.com/s/1wsLJO0DQUZTiejuPFvTZUQ?pwd=envi

提取码：envi

SARscape5.6.2软件的试用请浏览：https://www.cnblogs.com/enviidl/p/16266813.html 

运行程序：

• 将下面的代码复制到名为“sar_intensityanalysis.pro”的文件中，在IDL上方工具栏点击编译按钮。
• 在 IDL 命令提示符下执行以下命令：

IDL> sar_intensityanalysis

运行该程序后，系统将提示您选择数据目录。sar_intensityanalysis 有两个输入关键字。有关这些关键字的信息，请参阅pro例程的开头注释。

;------------------------------------------------------------------------------
;+
; This helper routine allows the user to toggle between two ENVIRasterLayers
; in the ENVI display.
;-
pro SAR_ToggleLayers, eLayer, $
  DESCRIPTION=desc, $
  DESTROY=destroy, $
  END_DESCRIPTION=descEnd, $
  TITLE=title, $
  TOP_LEVEL_BASE=tlb
  compile_opt idl2, logical_predicate
  if (n_elements(eLayer) NE 2) then begin
    print, 'Invalid input: ELAYER must contain two ENVIRASTERLAYERs'
    return
  endif
  for i = 0, 1 do eLayer[i]->SetProperty, HIDE=i
  if ((n_elements(tlb) EQ 0) || ~widget_info(tlb, /VALID_ID)) then begin
    e = envi(/CURRENT)
    e->GetProperty, WIDGET_ID=wENVI
    tlb = widget_base(/COLUMN, GROUP_LEADER=wENVI, /FLOATING, $
      TITLE='Intensity Data', TLB_FRAME_ATTR=9)
    wBase = widget_base(tlb, /ALIGN_CENTER, /ROW)
    wLabel = widget_label(wBase, /DYNAMIC_RESIZE, UNAME='label', VALUE=' ')
    wBaseControl = widget_base(tlb, /COLUMN, UNAME='controls', XPAD=0, YPAD=0)
    wBase = widget_base(wBaseControl, /EXCLUSIVE, /ROW)
    wButton = [ $
      widget_button(wBase, VALUE='cross-polarized (VH)', UNAME='button_0', UVALUE=eLayer[0]), $
      widget_button(wBase, VALUE='co-polarized (VV)', UNAME='button_1', UVALUE=eLayer[1])]
    widget_control, wButton[0], /SET_BUTTON
    wBase = widget_base(wBaseControl, /ALIGN_RIGHT, /ROW)
    wButtonClose = widget_button(wBase, UNAME='button_continue', VALUE='Continue')
    widget_control, tlb, /REALIZE
  endif else begin
    widget_control, tlb, /MAP
    wBaseControl = widget_info(tlb, FIND_BY_UNAME='controls')
    wLabel = widget_info(tlb, FIND_BY_UNAME='label')
    wButton = [widget_info(tlb, FIND_BY_UNAME='button_0'), $
      widget_info(tlb, FIND_BY_UNAME='button_1')]
    wButtonClose = widget_info(tlb, FIND_BY_UNAME='button_continue')
  endelse
  widget_control, wBaseControl, /SENSITIVE
  for i = 0, 1 do widget_control, wButton[i], SET_BUTTON=(1-i)
  if (n_elements(title) EQ 1) then widget_control, tlb, BASE_SET_TITLE=title
  widget_control, wLabel, SET_VALUE=((n_elements(desc) GT 0) ? desc[0] : ' ')
  if keyword_set(destroy) then widget_control, wButtonClose, SET_VALUE='Close'
  active = !true
  while (active) do begin
    sEvent = widget_event(tlb, /NOWAIT)
    switch sEvent.id of
      wButton[0]: index = 0
      wButton[1]: begin
        if (sEvent.select EQ 1) then begin
          if (n_elements(index) EQ 0) then index = 1
          eLayer[index]->SetProperty, HIDE=0
          eLayer[1-temporary(index)]->SetProperty, HIDE=1
        endif else [] = temporary(index)
        break
      end

      wButtonClose: begin
        active = !false
        break
      end

      else:
    endswitch
  endwhile
  if keyword_set(destroy) then begin
    widget_control, tlb, /DESTROY
  endif else begin
    widget_control, wBaseControl, SENSITIVE=0
    if (n_elements(descEnd) GT 0) then begin
      widget_control, wLabel, SET_VALUE=descEnd[0]
    endif
  endelse
end

;------------------------------------------------------------------------------
;+
; This program performs all of the steps done in the Sentinel-1 Intensity
; Analysis in ENVI SARscape Tutorial
;
; :Keywords:
;   DIRECTORY: in, optional, type="string"
;     Set this to the full path of the directory containing the tutorial data.
;     If this is not set, a dialog will be presented for the user to select
;     this folder.
;   SHOW_TOGGLE_DIALOG: in, optional, type="boolean"
;     Set this keyword to have a dialog displayed that will cause the halt
;     execution at the end of each step.  The dialog will allow the user to
;     toggle between the VH and VV data in the ENVI display.
;
;-
pro SAR_IntensityAnalysis, $
  DIRECTORY=dirRoot, $
  SHOW_TOGGLE_DIALOG=showDialog
  compile_opt idl2, logical_predicate
  err = 0
  catch, err
  if (err NE 0) then begin
    catch, /CANCEL
    help, /LAST_MESSAGE
    return
  endif
  if (n_elements(dirRoot) EQ 0) then begin
    dirRoot = dialog_pickfile(/DIRECTORY, /MUST_EXIST, $
      TITLE='Select the Tutorial Data Directory')
  endif
  ;
  ; Open and display the raw intensity data.
  ;
  file = filepath('sentinel1_83_20171212_095750732_IW_SIW1_D_'+['VH','VV']+'_slc_list_pwr', $
    ROOT_DIR=dirRoot, SUBDIR='2017-12-12-SLC-Power')
  if ~(total(file_test(file, /REGULAR)) EQ 2) then begin
    print, 'Input directory is missing data files'
    return
  endif
  e = envi(/CURRENT)
  if ~obj_valid(e) then e = envi()
  eView = e->GetView()
  eRaster = [e->OpenRaster(file[0]), e->OpenRaster(file[1])]
  eLayer = [eView->CreateLayer(eRaster[0]), eView->CreateLayer(eRaster[1])]
  eView->Zoom, LAYER=eLayer[0], /FULL_EXTENT
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayer, END_DESCRIPTION='Creating subset...', $
      TITLE='Raw Intensity Data', TOP_LEVEL_BASE=tlb
  endif
  ;
  ; Create spatial subsets of the VH and VV intensity images.
  ;
  fileAOI = filepath('StudySite.shp', ROOT_DIR=dirRoot)
  fileCut = filepath(file_basename(file)+'_cut', ROOT_DIR=dirRoot)
  eTaskCut = envitask('SARsToolsManualSelection')
  eTaskCut.input_sarscapedata = file
  eTaskCut.cut_shape_file = fileAOI
  eTaskCut.output_sarscapedata = fileCut
  eTaskCut.cut_is_georeferenced = !false
  eTaskCut.USE_only_shape_corners = !true
  eTaskCut.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  eTaskCut->Execute
  eRasterCut = [e->OpenRaster(fileCut[0]), e->OpenRaster(fileCut[1])]
  eLayerCut = [eView->CreateLayer(eRasterCut[0]), $
    eView->CreateLayer(eRasterCut[1])]
  eView->Zoom, LAYER=eLayerCut[0], /FULL_EXTENT
  for i = 0, 1 do eLayer[i]->Close
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayerCut, END_DESCRIPTION='Filtering data...', $
      TITLE='Cut Intensity Data', TOP_LEVEL_BASE=tlb
  endif
  ;
  ; Apply a filter to reduce speckling noise in the VH and VV images.
  ;
  fileFilter = filepath(file_basename(fileCut)+'_fil', ROOT_DIR=dirRoot)
  eTaskFilter = envitask('SARsWF_ToolsGenericFilterSingleImage')
  eTaskFilter.filt_type = 'Gamma APM'
  eTaskFilter.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  for i = 0, 1 do begin
    eTaskFilter.input_sarscapedata = fileCut[i]
    eTaskFilter.output_sarscapedata = fileFilter[i]
    eTaskFilter->Execute
  endfor
  eRasterFilter = [e->OpenRaster(fileFilter[0]), e->OpenRaster(fileFilter[1])]
  eLayerFilter = [eView->CreateLayer(eRasterFilter[0]), $
    eView->CreateLayer(eRasterFilter[1])]
  for i = 0, 1 do eLayerCut[i]->SetProperty, HIDE=1
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayerFilter, $
      END_DESCRIPTION='Geocoding and radiometric correction...', $
      TITLE='Filtered Cut Data', TOP_LEVEL_BASE=tlb
  endif
  ;
  ; Create a DEM file for georeferencing.
  ;
  fileDEM = filepath(file_basename(fileFilter[0])+'_srtm3_dem', ROOT_DIR=dirRoot)
  eTaskDEM = envitask('SARsToolsDEMExtractionSRTM4')
  eTaskDEM.output_cartographic_system = ['GEO-GLOBAL','','GEO','','WGS84','','0.00']
  eTaskDEM.reference_sarscapedata = fileFilter
  eTaskDEM.output_sarscapedata = fileDEM
  eTaskDEM.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  eTaskDEM->Execute
  ;
  ; Georeference the VH and VV images to a standard map projection, also called
  ; geocoding. As part of this step, radiometric calibration and normalization
  ; also will be done.
  ;
  fileGeo = filepath(file_basename(fileFilter)+'_geo', ROOT_DIR=dirRoot)
  eTaskGeo = envitask('SARsBasicGeocoding')
  eTaskGeo.input_sarscapedata = fileFilter
  eTaskGeo.output_sarscapedata = fileGeo
  eTaskGeo.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  eTaskGeo.dem_sarscapedata = fileDEM
  eTaskGeo->Execute
  eRasterGeo = [e->OpenRaster(fileGeo[0]), e->OpenRaster(fileGeo[1])]
  eLayerGeo = [eView->CreateLayer(eRasterGeo[0]), $
    eView->CreateLayer(eRasterGeo[1])]
  eView->Zoom, LAYER=eLayerGeo[0], /FULL_EXTENT
  foreach eLayerClose, [eLayerCut,eLayerFilter] do eLayerClose->Close
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayerGeo, /DESTROY, $
      TITLE='Georeferenced Images', TOP_LEVEL_BASE=tlb
  endif
  ;
  ; Create a color composite image from the VH and VV images.
  ;
  fileColor = filepath(file_basename(fileGeo[0])+'_color', ROOT_DIR=dirRoot)
  vh = eRasterGeo[0]->GetData()
  vv = eRasterGeo[1]->GetData()
  dim = size(vh, /DIMENSIONS)
  rgb = fltarr(dim[0],dim[1],3)
  rgb[*,*,0] = vv
  rgb[*,*,1] = vh
  rgb[*,*,2] = vv/vh
  eRasterColor = enviraster(rgb, URI=fileColor)
  eRasterColor->Save
  eLayerColor = eView->CreateLayer(eRasterColor)
  eView->Zoom, LAYER=eLayerColor, /FULL_EXTENT
  for i = 0, 1 do eLayerGeo[i]->SetProperty, HIDE=1
end

此示例演示了以编程方式调用ENVI SARscape任务以自动运行强度分析。如果您正在寻找 SAR 数据的自动化工作流程，ENVI SARscape Analytics将引导用户完成使用 SAR 数据执行的一些最常见的处理任务。而且，如果您认为可以使用 SAR 数据解决问题，请联系我们——我们可以与您合作，根据您的特定项目要求创建自定义工作流。

作者：Daryl Atencio

翻译、代码校验：马宇龙