From 3df6eb53fcdcad40e2d6f0f35923b82da7991b6e Mon Sep 17 00:00:00 2001 From: Harmon Date: Tue, 14 Apr 2026 01:39:51 +0800 Subject: [PATCH] Initial commit --- DinsarWorkflow_ALL.pro | 373 ++++++++++++++++++++++++++ README.md | 316 +++++++++++++++++++++++ Run_Batch_Task_Import.pro | 160 ++++++++++++ batch_dinsarworkflow_all.pro | 488 +++++++++++++++++++++++++++++++++++ 4 files changed, 1337 insertions(+) create mode 100644 DinsarWorkflow_ALL.pro create mode 100644 README.md create mode 100644 Run_Batch_Task_Import.pro create mode 100644 batch_dinsarworkflow_all.pro diff --git a/DinsarWorkflow_ALL.pro b/DinsarWorkflow_ALL.pro new file mode 100644 index 0000000..3a6ddd7 --- /dev/null +++ b/DinsarWorkflow_ALL.pro @@ -0,0 +1,373 @@ +;+ +; PURPOSE: +; This script automates a full D-InSAR workflow, from interferogram +; generation to geocoding (Steps 1-6). +; +; This is the final, complete, and integrated version resulting from +; our collaborative development and debugging. +;- + +;================================================================================ +; ### HELPER FUNCTION: READ SML PARAMETER ### +;================================================================================ +FUNCTION Read_SML_Parameter, sml_file, param_name_to_find + IF N_ELEMENTS(sml_file) EQ 0 OR sml_file EQ '' THEN BEGIN + PRINT, 'Error: An empty or undefined filename was passed to Read_SML_Parameter.' + RETURN, '' + ENDIF + IF NOT FILE_TEST(sml_file) THEN BEGIN + PRINT, 'Error: SML file not found -> ' + sml_file + RETURN, '' + ENDIF + lun = -1 + ON_ERROR, 2 + OPENR, lun, sml_file, /GET_LUN, ERROR=err + IF (err NE 0) THEN BEGIN + PRINT, 'Error: Could not open SML file for reading -> ' + sml_file + IF lun NE -1 THEN FREE_LUN, lun + RETURN, '' + ENDIF + found_value = '' + line = '' + start_tag = '<' + STRUPCASE(param_name_to_find) + '>' + end_tag = '' + WHILE NOT EOF(lun) DO BEGIN + READF, lun, line + clean_line = STRTRIM(line, 2) + line_upper = STRUPCASE(clean_line) + start_pos = STRPOS(line_upper, start_tag) + IF start_pos NE -1 THEN BEGIN + end_pos = STRPOS(line_upper, end_tag, start_pos) + IF end_pos NE -1 THEN BEGIN + start_extract = start_pos + STRLEN(start_tag) + value_length = end_pos - start_extract + found_value = STRMID(clean_line, start_extract, value_length) + BREAK + ENDIF + ENDIF + ENDWHILE + IF lun NE -1 THEN FREE_LUN, lun + RETURN, found_value +END + +;================================================================================ +; ### HELPER FUNCTION: CREATE GCPs FROM COHERENCE ### +;================================================================================ +FUNCTION CREATE_GCPS_FROM_COHERENCE, coherence_file, output_shp_file, COH_THRESHOLD=coh_threshold, NUM_POINTS=num_points + + ON_ERROR, 2 + IF FILE_TEST(output_shp_file) THEN FILE_DELETE, output_shp_file, /QUIET + + IF N_ELEMENTS(coh_threshold) EQ 0 THEN coh_threshold = 0.7 + IF N_ELEMENTS(num_points) EQ 0 THEN num_points = 100 + + e = ENVI(/HEADLESS) + + PRINT, 'Opening coherence file: ' + coherence_file + oRaster = e.OpenRaster(coherence_file) + + ns = oRaster.NCOLUMNS + nl = oRaster.NROWS + data = oRaster.GetData(BANDS=[0]) + + PRINT, 'Searching for GCPs with coherence > ' + STRTRIM(coh_threshold, 2) + + grid_dim = CEIL(SQRT(num_points)) + x_step = FLOOR(ns / grid_dim) + y_step = FLOOR(nl / grid_dim) + + pixel_points_x = [] & pixel_points_y = [] + + FOR j=0L, grid_dim-1 DO BEGIN + FOR i=0L, grid_dim-1 DO BEGIN + x_start = i * x_step & y_start = j * y_step + x_end = (i+1) * x_step -1 < (ns-1) + y_end = (j+1) * y_step -1 < (nl-1) + IF (x_start GE ns) OR (y_start GE nl) THEN CONTINUE + cell_data = data[x_start:x_end, y_start:y_end] + max_val = MAX(cell_data, max_idx) + IF max_val GE coh_threshold THEN BEGIN + max_y_cell = max_idx / (x_end - x_start + 1) + max_x_cell = max_idx MOD (x_end - x_start + 1) + pixel_points_x = [pixel_points_x, x_start + max_x_cell] + pixel_points_y = [pixel_points_y, y_start + max_y_cell] + ENDIF + ENDFOR + ENDFOR + + point_count = N_ELEMENTS(pixel_points_x) + IF point_count EQ 0 THEN BEGIN + PRINT, 'Error: No GCPs found.' & OBJ_DESTROY, oRaster & RETURN, 0 + ENDIF + + PRINT, 'Found ' + STRTRIM(point_count, 2) + ' suitable GCPs.' + PRINT, 'Creating shapefile using final authoritative methodology...' + + oShape = OBJ_NEW('IDLffShape', output_shp_file, ENTITY_TYPE=1, /UPDATE) + + oShape->AddAttribute, 'SHP_ID', 3, 10 + oShape->AddAttribute, 'GCP_LABEL', 7, 20 + oShape->AddAttribute, 'GCP_TYPE', 7, 20 + oShape->AddAttribute, 'GCP_COLUMN', 5, 24, PRECISION=6 + oShape->AddAttribute, 'GCP_ROW', 5, 24, PRECISION=6 + oShape->AddAttribute, 'GCP_OTHER_', 7, 10 ; User-verified field name + + attr_template = oShape->GetAttributes(/ATTRIBUTE_STRUCTURE) + + entNew = {IDL_SHAPE_ENTITY} + entNew.SHAPE_TYPE = 1 + + FOR i=0, point_count-1 DO BEGIN + x_float = FLOAT(pixel_points_x[i]) + y_float = FLOAT(pixel_points_y[i]) + + entNew.ISHAPE = i + entNew.BOUNDS = [x_float, y_float, 0.0, 0.0, x_float, y_float, 0.0, 0.0] + oShape->PutEntity, entNew + + attr = attr_template + attr.ATTRIBUTE_0 = i + attr.ATTRIBUTE_1 = 'GCP_' + STRTRIM(i+1, 2) + attr.ATTRIBUTE_2 = 'undefined' + attr.ATTRIBUTE_3 = x_float + attr.ATTRIBUTE_4 = y_float + attr.ATTRIBUTE_5 = '' + oShape->SetAttributes, i, attr + ENDFOR + + oShape->Close + OBJ_DESTROY, [oShape, oRaster] + + PRINT, 'Successfully created GCP shapefile: ' + output_shp_file + RETURN, 1 +END + +;================================================================================ +; ### MAIN WORKFLOW PROCEDURE ### +;================================================================================ +PRO DinsarWorkflow_ALL + + ; --- Initialize SARscape Environment --- + temporary_directory = 'D:\Sarscape_IDL_Area' + SARscape_Batch_Init,Temp_Directory=temporary_directory + e = ENVI(/HEADLESS) + PRINT, 'ENVI 和 SARscape 环境初始化成功！' + PRINT, '' + + ; =================================================================== + ; ### GLOBAL CONFIGURATION ### + ; =================================================================== + + master_base_file = 'G:/Code/InSAR/IDLTools/Task_20250504_20250601_1/LuTan-1_LT1A_39_20250504_091105615_STRIP1_A_HH_slc' + slave_base_file = 'G:/Code/InSAR/IDLTools/Task_20250504_20250601_1/LuTan-1_LT1B_47_20250601_091101714_STRIP1_A_HH_slc' + dem_base_file = 'D:/SRTM30m/SRTMDEM_RSP_SARscape' + workflow_root_name = 'G:/Code/InSAR/IDLTools/test/test' ; Use a consistent root name + + target_ground_resolution_m = 10.0 + filter_method = 'GOLDSTEIN' + unwrapping_coh_threshold = 0.05 + gcp_coh_threshold = 0.7 + gcp_number = 100 + geocoding_coh_threshold = 0.0 + geocoding_pixel_size_m = 10.0 + + ; =================================================================== + ; ### STEP 1: INTERFEROGRAM GENERATION ### + ; =================================================================== + + print, '======================================================' + print, '### STEP 1: INTERFEROGRAM GENERATION ###' + print, '======================================================' + + master_sml = master_base_file + '.sml' & slave_sml = slave_base_file + '.sml' + m_rg_sp_str = Read_SML_Parameter(master_sml, 'PixelSpacingRg') + m_az_sp_str = Read_SML_Parameter(master_sml, 'PixelSpacingAz') + m_inc_ang_str = Read_SML_Parameter(master_sml, 'IncidenceAngle') + s_rg_sp_str = Read_SML_Parameter(slave_sml, 'PixelSpacingRg') + s_az_sp_str = Read_SML_Parameter(slave_sml, 'PixelSpacingAz') + s_inc_ang_str = Read_SML_Parameter(slave_sml, 'IncidenceAngle') + + IF m_rg_sp_str EQ '' OR s_rg_sp_str EQ '' THEN GOTO, WORKFLOW_EXIT_FAILURE + + m_rg_sp = FLOAT(m_rg_sp_str) & m_az_sp = FLOAT(m_az_sp_str) & m_inc_ang = FLOAT(m_inc_ang_str) + s_rg_sp = FLOAT(s_rg_sp_str) & s_az_sp = FLOAT(s_az_sp_str) & s_inc_ang = FLOAT(s_inc_ang_str) + + avg_az_spacing = (m_az_sp + s_az_sp) / 2.0 & azimuth_looks = (long(target_ground_resolution_m / avg_az_spacing)) > 1 + m_ground_rg_sp = m_rg_sp / SIN(m_inc_ang * !DPI / 180.0) & s_ground_rg_sp = s_rg_sp / SIN(s_inc_ang * !DPI / 180.0) + avg_ground_rg_sp = (m_ground_rg_sp + s_ground_rg_sp) / 2.0 & range_looks = (long(target_ground_resolution_m / avg_ground_rg_sp)) > 1 + + task1_obj = obj_new('SARscapeBatch', Module='InSARInterferogramGeneration') + prefix_1 = 'MAIN_INSAR_INTERFEROGRAM_GENERATION_CMD.' + task1_obj->SetParam, prefix_1 + 'INPUT_REFERENCE_FILE_NAME' , master_base_file + task1_obj->SetParam, prefix_1 + 'INPUT_SECONDARY_FILE_NAME' , slave_base_file + task1_obj->SetParam, prefix_1 + 'DEM_FILE_NAME' , dem_base_file + task1_obj->SetParam, prefix_1 + 'OUTPUT_ROOT_FILE_NAME' , workflow_root_name + task1_obj->SetParam, prefix_1 + 'RG_LOOKS_NBR' , string(range_looks) + task1_obj->SetParam, prefix_1 + 'AZ_LOOKS_NBR' , string(azimuth_looks) + task1_obj->SetParam, 'COREGISTRATION_CMD.COREGISTRATION_WITH_DEM_FLAG' , 'OK' + task1_obj->SetParam, 'FLAT_CMD.MAKE_FLATTENING_FLAG', 'OK' + ok_step1 = task1_obj->Execute() + obj_destroy, task1_obj + IF NOT ok_step1 THEN GOTO, WORKFLOW_EXIT_FAILURE + print, '### STEP 1: SUCCESS!' + + ; =================================================================== + ; ### STEP 2: FILTERING AND COHERENCE GENERATION ### + ; =================================================================== + + print, '======================================================' + print, '### STEP 2: FILTERING AND COHERENCE GENERATION ###' + print, '======================================================' + task2_obj = obj_new('SARscapeBatch', Module='InSARFilterAndCoherence') + prefix_main2 = 'MAIN_INSAR_FILTER_COHERENCE_CMD.' + prefix_filt2 = 'FILTERING_CMD.' + task2_obj->SetParam, prefix_main2 + 'INPUT_INTERF_FILE_NAME', workflow_root_name + '_dint' + task2_obj->SetParam, prefix_main2 + 'INPUT_REFERENCE_FILE_NAME', workflow_root_name + '_reference_pwr' + task2_obj->SetParam, prefix_main2 + 'INPUT_SECONDARY_FILE_NAME', workflow_root_name + '_secondary_pwr' + task2_obj->SetParam, prefix_main2 + 'OUT_ROOT_FILE', workflow_root_name + task2_obj->SetParam, prefix_filt2 + 'FILTERING_METHOD', filter_method + task2_obj->SetParam, prefix_main2 + 'COHERENCE_FLAG', 'OK' + task2_obj->SetParam, prefix_main2 + 'INTERF_FILT_FLAG', 'OK' + ok_step2 = task2_obj->Execute() + obj_destroy, task2_obj + IF NOT ok_step2 THEN GOTO, WORKFLOW_EXIT_FAILURE + print, '### STEP 2: SUCCESS!' + + ; =================================================================== + ; ### STEP 3: ORBITAL TREND REMOVAL ### + ; =================================================================== + print, '======================================================' + print, '### STEP 3: ORBITAL TREND REMOVAL ###' + print, '======================================================' + + reflat_output_root = workflow_root_name + '_rrpf' + task3_obj = obj_new('SARscapeBatch', Module='InSARRemoveResidualPhaseFrequency') + prefix_3 = 'MAIN_INSAR_REMOVE_RESIDUAL_PHASE_FREQUENCY.' + task3_obj->SetParam, prefix_3 + 'INTERF_FILE_NAME', workflow_root_name + '_fint' + task3_obj->SetParam, prefix_3 + 'COHERENCE_FILE_NAME', workflow_root_name + '_cc' + task3_obj->SetParam, prefix_3 + 'OUT_ROOT_FILE_NAME', reflat_output_root + ok_step3 = task3_obj->Execute() + obj_destroy, task3_obj + IF NOT ok_step3 THEN GOTO, WORKFLOW_EXIT_FAILURE + print, '### STEP 3: SUCCESS!' + + ; =================================================================== + ; ### STEP 4: PHASE UNWRAPPING ### + ; =================================================================== + print, '======================================================' + print, '### STEP 4: PHASE UNWRAPPING ###' + print, '======================================================' + + reflat_upha_output_file = workflow_root_name + '_upha' + task4_obj = obj_new('SARscapeBatch', Module='InSARPhaseUnwrapping') + prefix_main4 = 'MAIN_INSAR_PHASE_UNWRAPPING_CMD.' + prefix_upha4 = 'UPHA_CMD.' + task4_obj->SetParam, prefix_main4 + 'INFILE_NAME', reflat_output_root + '_fint' + task4_obj->SetParam, prefix_main4 + 'COHERENCEFILE_NAME', workflow_root_name + '_cc' + task4_obj->SetParam, prefix_upha4 + 'UPHA_COH_THRESHOLD', string(unwrapping_coh_threshold) + task4_obj->SetParam, prefix_main4 + 'OUTFILE_NAME', reflat_upha_output_file + ok_step4 = task4_obj->Execute() + obj_destroy, task4_obj + IF NOT ok_step4 THEN GOTO, WORKFLOW_EXIT_FAILURE + print, '### STEP 4: SUCCESS!' + + ; =================================================================== + ; ### STEP 5: REFINEMENT AND REFLATTENING ### + ; =================================================================== + + ; --- STEP 5A: Automatic GCP Generation --- + print, '======================================================' + print, '### STEP 5A: 正在自动生成GCP... ###' + print, '======================================================' + auto_gcp_shp = workflow_root_name + '_auto_gcp.shp' + gcp_success = CREATE_GCPS_FROM_COHERENCE(workflow_root_name + '_cc', auto_gcp_shp, COH_THRESHOLD=gcp_coh_threshold, NUM_POINTS=gcp_number) + + IF gcp_success EQ 0 THEN BEGIN + PRINT, '错误: 自动生成GCP shapefile失败。流程终止。' + GOTO, WORKFLOW_EXIT_FAILURE + ENDIF + PRINT, '### STEP 5A: SUCCESS!' + + ; --- STEP 5B: Refinement and Reflattening Execution --- + print, '======================================================' + print, '### STEP 5B: 正在执行 REFINEMENT AND REFLATTENING ###' + print, '======================================================' + ref_output_root = workflow_root_name + '_reflat' + task_refine = obj_new('SARscapeBatch', Module='InSARRefinementAndReflattening') + prefix_refine = 'MAIN_INSAR_REFINEMENT_AND_REFLATTENING_CMD.' + + task_refine->SetParam, prefix_refine + 'INPUT_UPHA_FILE_NAME', reflat_upha_output_file + task_refine->SetParam, prefix_refine + 'INPUT_REFERENCE_FILE_NAME', workflow_root_name + '_reference_pwr' + task_refine->SetParam, prefix_refine + 'INPUT_SECONDARY_FILE_NAME', workflow_root_name + '_secondary_pwr' + task_refine->SetParam, prefix_refine + 'SLANT_RANGE_DEM_FILE_NAME', workflow_root_name + '_srdem' + task_refine->SetParam, prefix_refine + 'SYNTHETIC_FILE_NAME', workflow_root_name + '_sint' + task_refine->SetParam, prefix_refine + 'COHERENCE_FILE_NAME', workflow_root_name + '_cc' + task_refine->SetParam, prefix_refine + 'OUTPUT_ROOT_NAME', ref_output_root + task_refine->SetParam, prefix_refine + 'DEM_FILE_NAME', dem_base_file + task_refine->SetParam, 'REFINEMENT_CMD.REFINEMENT_GCP_FILE_NAME', auto_gcp_shp + + ok_step5 = task_refine->Execute() + obj_destroy, task_refine + IF NOT ok_step5 THEN GOTO, WORKFLOW_EXIT_FAILURE + print, '### STEP 5B: SUCCESS!' + + ; =================================================================== + ; ### STEP 6: PHASE TO DISPLACEMENT AND GEOCODING ### + ; =================================================================== + print, '======================================================' + print, '### STEP 6: PHASE TO DISPLACEMENT AND GEOCODING ###' + print, '======================================================' + + final_upha_file = ref_output_root + '_upha' + geocoded_output_root = workflow_root_name + '_geo' + + task6_obj = obj_new('SARscapeBatch', Module='InSARPhaseToDisplacement') + + IF (~OBJ_VALID(task6_obj)) THEN BEGIN + PRINT, '错误: 创建 SARscapeBatch 对象失败: InSARPhaseToDisplacement' + GOTO, WORKFLOW_EXIT_FAILURE + ENDIF + + ; --- Set Main Parameters --- + prefix_main6 = 'MAIN_INSAR_PHASE_TO_DISPLACEMENT_CMD.' + task6_obj->SetParam, prefix_main6 + 'INPUT_FILE_NAME', final_upha_file + task6_obj->SetParam, prefix_main6 + 'OUTPUT_FILE', geocoded_output_root + task6_obj->SetParam, prefix_main6 + 'COHERENCE_FILE', workflow_root_name + '_cc' + task6_obj->SetParam, prefix_main6 + 'DEM_FILE_NAME', dem_base_file + task6_obj->SetParam, prefix_main6 + 'COHERENCE_THRESHOLD', STRING(geocoding_coh_threshold) + + ; --- Set Other Required Parameters --- + task6_obj->SetParam, prefix_main6 + 'INTERPOL_TYPE' , '4th_order_cc' + task6_obj->SetParam, prefix_main6 + 'ALLOW_SKIP_REFINEMENT' , 'NotOK' + task6_obj->SetParam, prefix_main6 + 'SARSCAPEENVIRONMENT' , 'IDL_ENVI_ENV' + + ; --- Set Geocoding Parameters --- + prefix_geocode6 = 'GEOCODE_CMD.' + task6_obj->SetParam, prefix_geocode6 + 'GEOCODE_RG_GRID_SIZE', STRING(geocoding_pixel_size_m) + task6_obj->SetParam, prefix_geocode6 + 'GEOCODE_AZ_GRID_SIZE', STRING(geocoding_pixel_size_m) + task6_obj->SetParam, prefix_geocode6 + 'GEOCODE_INTERPOL_BOX_SIZE', '7' + + ; --- Set Displacement Product Generation Parameters --- + prefix_disp6 = 'DISPLACEMENT_PROJECTION_CMD.' + task6_obj->SetParam, prefix_disp6 + 'GENERATE_LOS_FLAG', 'OK' + task6_obj->SetParam, prefix_disp6 + 'GENERATE_VERTICAL_FLAG', 'NotOK' + task6_obj->SetParam, prefix_disp6 + 'GENERATE_MAX_SLOPE_FLAG', 'NotOK' + + print, '正在执行地理编码任务...' + ok_step6 = task6_obj->Execute() + obj_destroy, task6_obj + IF NOT ok_step6 THEN GOTO, WORKFLOW_EXIT_FAILURE + print, '### STEP 6: SUCCESS!' + + ; --- Final Success Message --- + print, '=======================================================' + print, ' FULL D-InSAR WORKFLOW (STEPS 1-6) COMPLETED SUCCESSFULLY!' + print, '=======================================================' + GOTO, WORKFLOW_EXIT_NORMAL + +WORKFLOW_EXIT_FAILURE: + print, 'Workflow terminated due to a failure in one of the steps.' +WORKFLOW_EXIT_NORMAL: + SARscape_Batch_Exit + +end \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..da60f67 --- /dev/null +++ b/README.md @@ -0,0 +1,316 @@ +# ENVI + IDL D-InSAR 批处理流水线 + +本仓库是一套基于 ENVI / IDL / SARscape 的 D-InSAR 自动化处理脚本，面向陆探一号等 SAR 数据的批量导入与差分干涉处理。代码包含两个主要阶段： + +1. 批量导入原始数据，将 `*.meta.xml` 导入为 SARscape 可继续处理的产品。 +2. 对已经完成导入的主辅影像执行完整 D-InSAR 流程，直到位移反演与地理编码输出。 + +仓库当前核心脚本只有 3 个： + +- `Run_Batch_Task_Import.pro`：批量导入 `Task_*` 目录中的原始数据。 +- `batch_dinsarworkflow_all.pro`：批量执行 D-InSAR 流程。 +- `DinsarWorkflow_ALL.pro`：单景测试版 / 调试版完整流程。 + +## 1. 依赖环境 + +运行本项目需要以下软件环境已经正确安装并可在 IDL 中调用： + +- ENVI +- IDL +- SARscape +- 对应 SARscape Batch 模块与 ENVI Task 环境 + +脚本中直接使用了以下对象 / 任务： + +- `ENVI(/HEADLESS)` +- `SARscape_Batch_Init` +- `SARscapeBatch` +- `ENVITask('SARsImportLuTan1')` +- `IDLffShape` + +因此该仓库不是纯 IDL 通用脚本，必须在装有 ENVI + SARscape 的环境中运行。 + +## 2. 仓库结构 + +```text +ENVI_IDL_DInSAR/ +├─ Run_Batch_Task_Import.pro +├─ batch_dinsarworkflow_all.pro +├─ DinsarWorkflow_ALL.pro +└─ README.md +``` + +各脚本职责如下。 + +### 2.1 `Run_Batch_Task_Import.pro` + +用途：批量扫描根目录下的 `Task_*` 文件夹，并分别处理其中的 `master` 与 `slave` 子目录。 + +处理逻辑： + +- 在 `master` / `slave` 目录中查找 `*.meta.xml` +- 使用 `ENVITask('SARsImportLuTan1')` 执行导入 +- 输出直接写回当前 `master` 或 `slave` 目录 +- 如果目录中已经存在 `*.sml`，则认为该目录已导入，自动跳过 + +适合放在整个流程的第一步，用来把原始数据准备成后续 D-InSAR 可识别的输入。 + +### 2.2 `batch_dinsarworkflow_all.pro` + +用途：批量执行完整 D-InSAR 流程。 + +处理逻辑： + +- 扫描根目录下所有 `Task_*` 文件夹 +- 检查每个任务目录下是否存在 `master` / `slave` +- 在 `master`、`slave` 中各取第一个 `*.sml` +- 去掉 `.sml` 后缀，得到 SARscape 输入基名 +- 在任务目录内创建 `dinsar_results` 输出目录 +- 调用 `Execute_Single_Dinsar_Workflow` 执行完整处理 + +这是生产使用时最重要的脚本。 + +### 2.3 `DinsarWorkflow_ALL.pro` + +用途：单组主辅影像的完整 D-InSAR 处理脚本。 + +特点： + +- 直接在脚本内手工指定 `master_base_file`、`slave_base_file`、`dem_base_file` +- 更适合单景调试、流程验证、参数试验 +- 逻辑上与批处理版本一致，但不是批量入口 + +## 3. 建议的数据组织方式 + +根据代码，批处理脚本假定数据目录组织如下： + +```text +/ +├─ Task_20250504_20250601_1/ +│ ├─ master/ +│ │ ├─ *.meta.xml +│ │ ├─ ... 导入后生成的 *.sml、影像及配套文件 +│ ├─ slave/ +│ │ ├─ *.meta.xml +│ │ ├─ ... 导入后生成的 *.sml、影像及配套文件 +│ └─ dinsar_results/ +│ └─ workflow* +├─ Task_xxx/ +└─ ... +``` + +处理顺序通常是： + +1. 先运行 `Run_Batch_Task_Import.pro` +2. 确认每个 `master` / `slave` 目录中已经生成 `*.sml` +3. 再运行 `batch_dinsarworkflow_all.pro` + +## 4. D-InSAR 主流程说明 + +`batch_dinsarworkflow_all.pro` 中的 `Execute_Single_Dinsar_Workflow` 实现了完整 6 步流程： + +### Step 1. 干涉图生成 + +- 从主影像和辅影像的 `.sml` 中读取： + - `PixelSpacingRg` + - `PixelSpacingAz` + - `IncidenceAngle` +- 根据目标地面分辨率自动计算多视数： + - `RG_LOOKS_NBR` + - `AZ_LOOKS_NBR` +- 调用 `InSARInterferogramGeneration` + +### Step 2. 滤波与相干图生成 + +- 调用 `InSARFilterAndCoherence` +- 使用设定的滤波方法，例如 `GOLDSTEIN` +- 生成滤波干涉图与相干图 + +### Step 3. 轨道残余相位频率去除 + +- 调用 `InSARRemoveResidualPhaseFrequency` +- 输入为上一步生成的滤波干涉图和相干图 + +### Step 4. 相位解缠 + +- 调用 `InSARPhaseUnwrapping` +- 使用相干阈值 `UPHA_COH_THRESHOLD` + +### Step 5. 精化与再平地 + +分为两步： + +1. 自动 GCP 生成 +2. Refinement and Reflattening + +自动 GCP 生成的逻辑为： + +- 打开相干图 +- 按网格均匀划分影像 +- 在每个网格中选取相干系数最大的像元 +- 当该像元相干值高于阈值时，将其写入 shapefile +- 最终生成 `*_auto_gcp.shp` + +随后调用 `InSARRefinementAndReflattening`，并将自动生成的 GCP 文件作为输入。 + +### Step 6. 相位转位移与地理编码 + +- 调用 `InSARPhaseToDisplacement` +- 使用最终解缠结果生成 LOS 位移产品 +- 执行地理编码 +- 当前脚本仅开启 LOS 输出： + - `GENERATE_LOS_FLAG = OK` + - `GENERATE_VERTICAL_FLAG = NotOK` + - `GENERATE_MAX_SLOPE_FLAG = NotOK` + +## 5. 关键输入参数 + +以下参数需要在脚本中按实际工程修改。 + +### 5.1 导入脚本参数 + +文件：`Run_Batch_Task_Import.pro` + +- `temporary_directory`：SARscape 临时目录 +- `rootDir`：待扫描的任务根目录 +- `numToProcess`：最多处理多少个任务，`0` 表示全部处理 + +### 5.2 批处理 D-InSAR 参数 + +文件：`batch_dinsarworkflow_all.pro` + +- `temporary_directory`：SARscape 临时目录 +- `rootDir`：包含多个 `Task_*` 目录的根目录 +- `dem_base_file`：DEM 基础文件路径 +- `target_ground_resolution_m`：目标地面分辨率，控制自动多视 +- `filter_method`：滤波方法，默认 `GOLDSTEIN` +- `unwrapping_coh_threshold`：解缠相干阈值 +- `gcp_coh_threshold`：自动 GCP 相干阈值 +- `gcp_number`：目标 GCP 数量 +- `geocoding_coh_threshold`：地理编码时使用的相干阈值 +- `geocoding_pixel_size_m`：地理编码输出像元大小 +- `numToProcess`：最多处理多少个任务，`0` 表示全部处理 + +### 5.3 单景调试参数 + +文件：`DinsarWorkflow_ALL.pro` + +- `master_base_file` +- `slave_base_file` +- `dem_base_file` +- `workflow_root_name` +- 其余处理参数与批处理版本基本一致 + +注意：`master_base_file` 和 `slave_base_file` 填写的是“基名”，不是 `.sml` 文件全名。脚本内部会自动拼接 `.sml` 去读取参数。 + +## 6. 输出结果说明 + +批处理脚本默认在每个任务目录下创建： + +```text +Task_xxx/ +└─ dinsar_results/ + ├─ workflow_dint + ├─ workflow_reference_pwr + ├─ workflow_secondary_pwr + ├─ workflow_fint + ├─ workflow_cc + ├─ workflow_rrpf_fint + ├─ workflow_upha + ├─ workflow_auto_gcp.shp + ├─ workflow_reflat_upha + ├─ workflow_geo* + └─ ... +``` + +其中常见关键中间结果包括： + +- `*_dint`：初始干涉图 +- `*_fint`：滤波后干涉图 +- `*_cc`：相干图 +- `*_upha`：解缠结果 +- `*_auto_gcp.shp`：自动生成的 GCP 矢量文件 +- `*_geo*`：地理编码后的位移产品 + +实际生成的附属文件数量和后缀会受 SARscape 版本及模块设置影响。 + +## 7. 典型使用方式 + +### 7.1 第一步：批量导入 + +在 `Run_Batch_Task_Import.pro` 中修改： + +- `temporary_directory` +- `rootDir` +- `numToProcess` + +然后在 IDL 命令行或脚本环境中运行： + +```idl +.compile Run_Batch_Task_Import.pro +Run_Batch_Task_Import +``` + +### 7.2 第二步：批量执行 D-InSAR + +在 `batch_dinsarworkflow_all.pro` 中修改： + +- `temporary_directory` +- `rootDir` +- `dem_base_file` +- 各处理参数 + +然后运行： + +```idl +.compile batch_dinsarworkflow_all.pro +Batch_DinsarWorkflow_ALL +``` + +### 7.3 单景测试 + +在 `DinsarWorkflow_ALL.pro` 中直接修改主影像、辅影像、DEM 和输出路径，然后运行： + +```idl +.compile DinsarWorkflow_ALL.pro +DinsarWorkflow_ALL +``` + +## 8. 代码实现要点 + +### 8.1 自动读取 SML 参数 + +`Read_SML_Parameter` 会直接从 `.sml` 文件中读取指定 XML 标签内容，用于自动计算多视参数。这意味着： + +- 输入数据必须已经正确导入并生成 `.sml` +- `.sml` 内部需要包含 `PixelSpacingRg`、`PixelSpacingAz`、`IncidenceAngle` 等字段 + +### 8.2 自动 GCP 生成策略 + +`CREATE_GCPS_FROM_COHERENCE` 不是随机取点，而是“规则网格 + 每格取最大相干值”的策略。这样做的优点是： + +- GCP 空间分布更均匀 +- 能避免点位过度集中 +- 能通过相干阈值控制稳定区域 + +### 8.3 批处理脚本的选择策略 + +当前批处理版本在每个 `master` / `slave` 目录中只取“第一个找到的 `.sml`”作为输入，因此建议： + +- 每个 `master` 目录内只保留一组待处理主影像 +- 每个 `slave` 目录内只保留一组待处理辅影像 +- 若目录中存在多份 `.sml`，需要先确认 `FILE_SEARCH` 返回顺序是否符合预期 + +## 9. 注意事项 + +1. 路径需要按本机环境手动修改，仓库中的路径示例均为作者本地路径。 +2. 脚本默认使用 `ENVI(/HEADLESS)`，适合批处理环境。 +3. `temporary_directory` 需要保证有足够磁盘空间。 +4. `gcp_coh_threshold` 过高时，可能导致无法生成足够 GCP，流程会在 Step 5A 失败。 +5. `unwrapping_coh_threshold`、`geocoding_coh_threshold` 需要结合研究区相干性调整。 +6. 本仓库目前没有统一的日志输出、配置文件和异常恢复机制，失败后通常需要根据控制台信息排查。 + + + + + diff --git a/Run_Batch_Task_Import.pro b/Run_Batch_Task_Import.pro new file mode 100644 index 0000000..be7256b --- /dev/null +++ b/Run_Batch_Task_Import.pro @@ -0,0 +1,160 @@ +PRO Run_Batch_Task_Import + + ;=========================================================================== + ; --- 用户配置区域 --- + ; 1. **数据根目录 (rootDir)** + ; 指定包含所有Task_开头文件夹的根目录路径。 + ; !!! 重要: 路径末尾不要带斜杠，脚本会自动添加。 + ; - 正确示例: 'C:\MyData\LuTan1' 或 '\\SERVER\Share\Data' + + ; 2. **处理数量 (numToProcess)** + ; 指定要处理的Task文件夹的最大数量。 + ; - 如果设置为 0 或一个负数, 程序将处理所有找到的文件夹。 + + ;=========================================================================== + temporary_directory = 'D:/Sarscape_IDL_Area' + SARscape_Batch_Init,Temp_Directory=temporary_directory + rootDir = 'K:/T2' + numToProcess = 0 + ;=========================================================================== + ; --- 脚本执行区域 --- + ;=========================================================================== + pathSep = '' + IF (!VERSION.OS_FAMILY EQ 'Windows') THEN pathSep = '\' ELSE pathSep = '/' + e = ENVI(/HEADLESS) + + IF rootDir EQ '' OR FILE_TEST(rootDir, /DIRECTORY) EQ 0 THEN BEGIN + PRINT, '错误: "rootDir" 参数未设置或指定的目录不存在: ' + rootDir + RETURN + ENDIF + + PRINT, '==============================================' + PRINT, '开始批量导入陆探一号数据...' + PRINT, '扫描目录: ' + rootDir + PRINT, '==============================================' + + ; 查找所有Task_开头的文件夹 - 使用兼容的方法 + PRINT, '正在扫描Task_开头的文件夹...' + + ; 兼容旧版本IDL的方法：先搜索所有匹配项，然后筛选出文件夹 + allItems = FILE_SEARCH(rootDir + pathSep + 'Task_*', COUNT=nTotal) + taskFolders = [] + nTasks = 0 + + IF nTotal GT 0 THEN BEGIN + FOR i = 0, nTotal - 1 DO BEGIN + IF FILE_TEST(allItems[i], /DIRECTORY) THEN BEGIN + taskFolders = [taskFolders, allItems[i]] + nTasks = nTasks + 1 + ENDIF + ENDFOR + ENDIF + + IF nTasks EQ 0 THEN BEGIN + PRINT, '在指定目录下没有找到任何Task_开头的文件夹。' + RETURN + ENDIF + + PRINT, '共找到 ' + STRING(nTasks) + ' 个Task文件夹，准备开始处理...' + PRINT, '' + + processed_count = 0 + + ; 使用标准的FOR循环替代FOREACH以确保兼容性 + FOR taskIndex = 0, nTasks - 1 DO BEGIN + taskFolder = taskFolders[taskIndex] + + IF (numToProcess GT 0) AND (processed_count GE numToProcess) THEN BEGIN + PRINT, '' + PRINT, '已达到指定的处理数量 (' + STRING(numToProcess) + ')，程序终止。' + BREAK + ENDIF + + taskName = FILE_BASENAME(taskFolder) + PRINT, '=== 正在处理任务文件夹: ' + taskName + ' ===' + + ; 检查master和slave子文件夹 + masterDir = taskFolder + pathSep + 'master' + slaveDir = taskFolder + pathSep + 'slave' + + ; 处理master文件夹 + IF FILE_TEST(masterDir, /DIRECTORY) THEN BEGIN + PRINT, '--- 处理master文件夹 ---' + processed_count = processed_count + Process_LuTan_Folder(masterDir, pathSep) + ENDIF ELSE BEGIN + PRINT, '警告: ' + taskName + ' 中未找到master文件夹' + ENDELSE + + ; 处理slave文件夹 + IF FILE_TEST(slaveDir, /DIRECTORY) THEN BEGIN + PRINT, '--- 处理slave文件夹 ---' + processed_count = processed_count + Process_LuTan_Folder(slaveDir, pathSep) + ENDIF ELSE BEGIN + PRINT, '警告: ' + taskName + ' 中未找到slave文件夹' + ENDELSE + + PRINT, '完成处理任务文件夹: ' + taskName + PRINT, '' + ENDFOR + + PRINT, '==============================================' + PRINT, '所有任务已完成！' + PRINT, '总共成功处理了 ' + STRING(processed_count) + ' 个数据文件夹。' + PRINT, '==============================================' + +END + +; 处理单个陆探一号数据文件夹的子函数 +FUNCTION Process_LuTan_Folder, folderPath, pathSep + COMPILE_OPT STRICTARR + + baseName = FILE_BASENAME(folderPath) + PRINT, '-> 正在检查文件夹: ' + baseName + + ; 直接使用master/slave文件夹作为输出目录，不再创建envi_import子目录 + outputDir = folderPath + + ; 检查是否已经处理过（通过检查是否存在.sml文件） + check_sml_file = FILE_SEARCH(outputDir + pathSep + '*.sml', COUNT=sml_count) + IF sml_count GT 0 THEN BEGIN + PRINT, '-> 跳过: ' + baseName + ' 文件夹中已存在处理结果。' + RETURN, 0 + ENDIF + + ; 查找元数据文件 - 查找文件夹中所有的.meta.xml文件 + metaFiles = FILE_SEARCH(folderPath + pathSep + '*.meta.xml', COUNT=nMetaFiles) + + IF nMetaFiles EQ 0 THEN BEGIN + PRINT, '-> 警告: 在文件夹 ' + baseName + ' 中未找到元数据文件(*.meta.xml)，跳过。' + RETURN, 0 + ENDIF + + ; 处理找到的元数据文件 + success_count = 0 + FOR i = 0, nMetaFiles - 1 DO BEGIN + inputFile = metaFiles[i] + metaFileName = FILE_BASENAME(inputFile) + PRINT, '-> 准备导入: ' + metaFileName + + ; 执行导入任务 + oSB = obj_new('SARscapeBatch') + task = ENVITask('SARsImportLuTan1') + task.INPUT_FILE_LIST = [inputFile] + task.ROOT_URI_FOR_OUTPUT = outputDir ; 直接输出到master/slave文件夹 + + CATCH, error_status + IF error_status NE 0 THEN BEGIN + PRINT, '-> 错误: 在处理 ' + metaFileName + ' 时发生严重错误。' + PRINT, !ERROR_STATE.MSG + CATCH, /CANCEL + CONTINUE + ENDIF + + task.Execute + CATCH, /CANCEL + PRINT, '-> 成功: 数据已导入到 ' + baseName + ' 文件夹中。' + success_count = success_count + 1 + ENDFOR + + RETURN, success_count +END \ No newline at end of file diff --git a/batch_dinsarworkflow_all.pro b/batch_dinsarworkflow_all.pro new file mode 100644 index 0000000..99b2731 --- /dev/null +++ b/batch_dinsarworkflow_all.pro @@ -0,0 +1,488 @@ +;================================================================================ +; ### HELPER FUNCTION: READ SML PARAMETER ### +;================================================================================ +FUNCTION Read_SML_Parameter, sml_file, param_name_to_find + IF N_ELEMENTS(sml_file) EQ 0 OR sml_file EQ '' THEN BEGIN + PRINT, 'Error: An empty or undefined filename was passed to Read_SML_Parameter.' + RETURN, '' + ENDIF + IF NOT FILE_TEST(sml_file) THEN BEGIN + PRINT, 'Error: SML file not found -> ' + sml_file + RETURN, '' + ENDIF + lun = -1 + ON_ERROR, 2 + OPENR, lun, sml_file, /GET_LUN, ERROR=err + IF (err NE 0) THEN BEGIN + PRINT, 'Error: Could not open SML file for reading -> ' + sml_file + IF lun NE -1 THEN FREE_LUN, lun + RETURN, '' + ENDIF + found_value = '' + line = '' + start_tag = '<' + STRUPCASE(param_name_to_find) + '>' + end_tag = '' + WHILE NOT EOF(lun) DO BEGIN + READF, lun, line + clean_line = STRTRIM(line, 2) + line_upper = STRUPCASE(clean_line) + start_pos = STRPOS(line_upper, start_tag) + IF start_pos NE -1 THEN BEGIN + end_pos = STRPOS(line_upper, end_tag, start_pos) + IF end_pos NE -1 THEN BEGIN + start_extract = start_pos + STRLEN(start_tag) + value_length = end_pos - start_extract + found_value = STRMID(clean_line, start_extract, value_length) + BREAK + ENDIF + ENDIF + ENDWHILE + IF lun NE -1 THEN FREE_LUN, lun + RETURN, found_value +END + +;================================================================================ +; ### HELPER FUNCTION: CREATE GCPs FROM COHERENCE ### +;================================================================================ +FUNCTION CREATE_GCPS_FROM_COHERENCE, coherence_file, output_shp_file, COH_THRESHOLD=coh_threshold, NUM_POINTS=num_points + + ON_ERROR, 2 + IF FILE_TEST(output_shp_file) THEN FILE_DELETE, output_shp_file, /QUIET + + IF N_ELEMENTS(coh_threshold) EQ 0 THEN coh_threshold = 0.7 + IF N_ELEMENTS(num_points) EQ 0 THEN num_points = 100 + + e = ENVI(/HEADLESS) + + PRINT, 'Opening coherence file: ' + coherence_file + oRaster = e.OpenRaster(coherence_file) + + ns = oRaster.NCOLUMNS + nl = oRaster.NROWS + data = oRaster.GetData(BANDS=[0]) + + PRINT, 'Searching for GCPs with coherence > ' + STRTRIM(coh_threshold, 2) + + grid_dim = CEIL(SQRT(num_points)) + x_step = FLOOR(ns / grid_dim) + y_step = FLOOR(nl / grid_dim) + + pixel_points_x = [] & pixel_points_y = [] + + FOR j=0L, grid_dim-1 DO BEGIN + FOR i=0L, grid_dim-1 DO BEGIN + x_start = i * x_step & y_start = j * y_step + x_end = (i+1) * x_step -1 < (ns-1) + y_end = (j+1) * y_step -1 < (nl-1) + IF (x_start GE ns) OR (y_start GE nl) THEN CONTINUE + cell_data = data[x_start:x_end, y_start:y_end] + max_val = MAX(cell_data, max_idx) + IF max_val GE coh_threshold THEN BEGIN + max_y_cell = max_idx / (x_end - x_start + 1) + max_x_cell = max_idx MOD (x_end - x_start + 1) + pixel_points_x = [pixel_points_x, x_start + max_x_cell] + pixel_points_y = [pixel_points_y, y_start + max_y_cell] + ENDIF + ENDFOR + ENDFOR + + point_count = N_ELEMENTS(pixel_points_x) + IF point_count EQ 0 THEN BEGIN + PRINT, 'Error: No GCPs found.' & OBJ_DESTROY, oRaster & RETURN, 0 + ENDIF + + PRINT, 'Found ' + STRTRIM(point_count, 2) + ' suitable GCPs.' + PRINT, 'Creating shapefile using final authoritative methodology...' + + oShape = OBJ_NEW('IDLffShape', output_shp_file, ENTITY_TYPE=1, /UPDATE) + + oShape->AddAttribute, 'SHP_ID', 3, 10 + oShape->AddAttribute, 'GCP_LABEL', 7, 20 + oShape->AddAttribute, 'GCP_TYPE', 7, 20 + oShape->AddAttribute, 'GCP_COLUMN', 5, 24, PRECISION=6 + oShape->AddAttribute, 'GCP_ROW', 5, 24, PRECISION=6 + oShape->AddAttribute, 'GCP_OTHER_', 7, 10 ; User-verified field name + + attr_template = oShape->GetAttributes(/ATTRIBUTE_STRUCTURE) + + entNew = {IDL_SHAPE_ENTITY} + entNew.SHAPE_TYPE = 1 + + FOR i=0, point_count-1 DO BEGIN + x_float = FLOAT(pixel_points_x[i]) + y_float = FLOAT(pixel_points_y[i]) + + entNew.ISHAPE = i + entNew.BOUNDS = [x_float, y_float, 0.0, 0.0, x_float, y_float, 0.0, 0.0] + oShape->PutEntity, entNew + + attr = attr_template + attr.ATTRIBUTE_0 = i + attr.ATTRIBUTE_1 = 'GCP_' + STRTRIM(i+1, 2) + attr.ATTRIBUTE_2 = 'undefined' + attr.ATTRIBUTE_3 = x_float + attr.ATTRIBUTE_4 = y_float + attr.ATTRIBUTE_5 = '' + oShape->SetAttributes, i, attr + ENDFOR + + oShape->Close + OBJ_DESTROY, [oShape, oRaster] + + PRINT, 'Successfully created GCP shapefile: ' + output_shp_file + RETURN, 1 +END + +;================================================================================ +; ### 执行单个D-InSAR工作流的函数 ### +;================================================================================ +FUNCTION Execute_Single_Dinsar_Workflow, master_base_file, slave_base_file, dem_base_file, workflow_root_name, $ + target_ground_resolution_m, filter_method, unwrapping_coh_threshold, $ + gcp_coh_threshold, gcp_number, geocoding_coh_threshold, geocoding_pixel_size_m + + ; =================================================================== + ; ### STEP 1: INTERFEROGRAM GENERATION ### + ; =================================================================== + + PRINT, '======================================================' + PRINT, '### STEP 1: INTERFEROGRAM GENERATION ###' + PRINT, '======================================================' + + master_sml = master_base_file + '.sml' + slave_sml = slave_base_file + '.sml' + + m_rg_sp_str = Read_SML_Parameter(master_sml, 'PixelSpacingRg') + m_az_sp_str = Read_SML_Parameter(master_sml, 'PixelSpacingAz') + m_inc_ang_str = Read_SML_Parameter(master_sml, 'IncidenceAngle') + s_rg_sp_str = Read_SML_Parameter(slave_sml, 'PixelSpacingRg') + s_az_sp_str = Read_SML_Parameter(slave_sml, 'PixelSpacingAz') + s_inc_ang_str = Read_SML_Parameter(slave_sml, 'IncidenceAngle') + + IF m_rg_sp_str EQ '' OR s_rg_sp_str EQ '' THEN RETURN, 0 + + m_rg_sp = FLOAT(m_rg_sp_str) & m_az_sp = FLOAT(m_az_sp_str) & m_inc_ang = FLOAT(m_inc_ang_str) + s_rg_sp = FLOAT(s_rg_sp_str) & s_az_sp = FLOAT(s_az_sp_str) & s_inc_ang = FLOAT(s_inc_ang_str) + + avg_az_spacing = (m_az_sp + s_az_sp) / 2.0 + azimuth_looks = (LONG(target_ground_resolution_m / avg_az_spacing)) > 1 + m_ground_rg_sp = m_rg_sp / SIN(m_inc_ang * !DPI / 180.0) + s_ground_rg_sp = s_rg_sp / SIN(s_inc_ang * !DPI / 180.0) + avg_ground_rg_sp = (m_ground_rg_sp + s_ground_rg_sp) / 2.0 + range_looks = (LONG(target_ground_resolution_m / avg_ground_rg_sp)) > 1 + + task1_obj = OBJ_NEW('SARscapeBatch', Module='InSARInterferogramGeneration') + prefix_1 = 'MAIN_INSAR_INTERFEROGRAM_GENERATION_CMD.' + task1_obj->SetParam, prefix_1 + 'INPUT_REFERENCE_FILE_NAME', master_base_file + task1_obj->SetParam, prefix_1 + 'INPUT_SECONDARY_FILE_NAME', slave_base_file + task1_obj->SetParam, prefix_1 + 'DEM_FILE_NAME', dem_base_file + task1_obj->SetParam, prefix_1 + 'OUTPUT_ROOT_FILE_NAME', workflow_root_name + task1_obj->SetParam, prefix_1 + 'RG_LOOKS_NBR', STRING(range_looks) + task1_obj->SetParam, prefix_1 + 'AZ_LOOKS_NBR', STRING(azimuth_looks) + task1_obj->SetParam, 'COREGISTRATION_CMD.COREGISTRATION_WITH_DEM_FLAG', 'OK' + task1_obj->SetParam, 'FLAT_CMD.MAKE_FLATTENING_FLAG', 'OK' + ok_step1 = task1_obj->Execute() + OBJ_DESTROY, task1_obj + IF NOT ok_step1 THEN RETURN, 0 + PRINT, '### STEP 1: SUCCESS!' + + ; =================================================================== + ; ### STEP 2: FILTERING AND COHERENCE GENERATION ### + ; =================================================================== + + PRINT, '======================================================' + PRINT, '### STEP 2: FILTERING AND COHERENCE GENERATION ###' + PRINT, '======================================================' + + task2_obj = OBJ_NEW('SARscapeBatch', Module='InSARFilterAndCoherence') + prefix_main2 = 'MAIN_INSAR_FILTER_COHERENCE_CMD.' + prefix_filt2 = 'FILTERING_CMD.' + task2_obj->SetParam, prefix_main2 + 'INPUT_INTERF_FILE_NAME', workflow_root_name + '_dint' + task2_obj->SetParam, prefix_main2 + 'INPUT_REFERENCE_FILE_NAME', workflow_root_name + '_reference_pwr' + task2_obj->SetParam, prefix_main2 + 'INPUT_SECONDARY_FILE_NAME', workflow_root_name + '_secondary_pwr' + task2_obj->SetParam, prefix_main2 + 'OUT_ROOT_FILE', workflow_root_name + task2_obj->SetParam, prefix_filt2 + 'FILTERING_METHOD', filter_method + task2_obj->SetParam, prefix_main2 + 'COHERENCE_FLAG', 'OK' + task2_obj->SetParam, prefix_main2 + 'INTERF_FILT_FLAG', 'OK' + ok_step2 = task2_obj->Execute() + OBJ_DESTROY, task2_obj + IF NOT ok_step2 THEN RETURN, 0 + PRINT, '### STEP 2: SUCCESS!' + + ; =================================================================== + ; ### STEP 3: ORBITAL TREND REMOVAL ### + ; =================================================================== + + PRINT, '======================================================' + PRINT, '### STEP 3: ORBITAL TREND REMOVAL ###' + PRINT, '======================================================' + + reflat_output_root = workflow_root_name + '_rrpf' + task3_obj = OBJ_NEW('SARscapeBatch', Module='InSARRemoveResidualPhaseFrequency') + prefix_3 = 'MAIN_INSAR_REMOVE_RESIDUAL_PHASE_FREQUENCY.' + task3_obj->SetParam, prefix_3 + 'INTERF_FILE_NAME', workflow_root_name + '_fint' + task3_obj->SetParam, prefix_3 + 'COHERENCE_FILE_NAME', workflow_root_name + '_cc' + task3_obj->SetParam, prefix_3 + 'OUT_ROOT_FILE_NAME', reflat_output_root + ok_step3 = task3_obj->Execute() + OBJ_DESTROY, task3_obj + IF NOT ok_step3 THEN RETURN, 0 + PRINT, '### STEP 3: SUCCESS!' + + ; =================================================================== + ; ### STEP 4: PHASE UNWRAPPING ### + ; =================================================================== + + PRINT, '======================================================' + PRINT, '### STEP 4: PHASE UNWRAPPING ###' + PRINT, '======================================================' + + reflat_upha_output_file = workflow_root_name + '_upha' + task4_obj = OBJ_NEW('SARscapeBatch', Module='InSARPhaseUnwrapping') + prefix_main4 = 'MAIN_INSAR_PHASE_UNWRAPPING_CMD.' + prefix_upha4 = 'UPHA_CMD.' + task4_obj->SetParam, prefix_main4 + 'INFILE_NAME', reflat_output_root + '_fint' + task4_obj->SetParam, prefix_main4 + 'COHERENCEFILE_NAME', workflow_root_name + '_cc' + task4_obj->SetParam, prefix_upha4 + 'UPHA_COH_THRESHOLD', STRING(unwrapping_coh_threshold) + task4_obj->SetParam, prefix_main4 + 'OUTFILE_NAME', reflat_upha_output_file + ok_step4 = task4_obj->Execute() + OBJ_DESTROY, task4_obj + IF NOT ok_step4 THEN RETURN, 0 + PRINT, '### STEP 4: SUCCESS!' + + ; =================================================================== + ; ### STEP 5: REFINEMENT AND REFLATTENING ### + ; =================================================================== + + ; --- STEP 5A: Automatic GCP Generation --- + PRINT, '======================================================' + PRINT, '### STEP 5A: 正在自动生成GCP... ###' + PRINT, '======================================================' + + auto_gcp_shp = workflow_root_name + '_auto_gcp.shp' + gcp_success = CREATE_GCPS_FROM_COHERENCE(workflow_root_name + '_cc', auto_gcp_shp, COH_THRESHOLD=gcp_coh_threshold, NUM_POINTS=gcp_number) + + IF gcp_success EQ 0 THEN BEGIN + PRINT, '错误: 自动生成GCP shapefile失败。流程终止。' + RETURN, 0 + ENDIF + PRINT, '### STEP 5A: SUCCESS!' + + ; --- STEP 5B: Refinement and Reflattening Execution --- + PRINT, '======================================================' + PRINT, '### STEP 5B: 正在执行 REFINEMENT AND REFLATTENING ###' + PRINT, '======================================================' + + ref_output_root = workflow_root_name + '_reflat' + task_refine = OBJ_NEW('SARscapeBatch', Module='InSARRefinementAndReflattening') + prefix_refine = 'MAIN_INSAR_REFINEMENT_AND_REFLATTENING_CMD.' + + task_refine->SetParam, prefix_refine + 'INPUT_UPHA_FILE_NAME', reflat_upha_output_file + task_refine->SetParam, prefix_refine + 'INPUT_REFERENCE_FILE_NAME', workflow_root_name + '_reference_pwr' + task_refine->SetParam, prefix_refine + 'INPUT_SECONDARY_FILE_NAME', workflow_root_name + '_secondary_pwr' + task_refine->SetParam, prefix_refine + 'SLANT_RANGE_DEM_FILE_NAME', workflow_root_name + '_srdem' + task_refine->SetParam, prefix_refine + 'SYNTHETIC_FILE_NAME', workflow_root_name + '_sint' + task_refine->SetParam, prefix_refine + 'COHERENCE_FILE_NAME', workflow_root_name + '_cc' + task_refine->SetParam, prefix_refine + 'OUTPUT_ROOT_NAME', ref_output_root + task_refine->SetParam, prefix_refine + 'DEM_FILE_NAME', dem_base_file + task_refine->SetParam, 'REFINEMENT_CMD.REFINEMENT_GCP_FILE_NAME', auto_gcp_shp + + ok_step5 = task_refine->Execute() + OBJ_DESTROY, task_refine + IF NOT ok_step5 THEN RETURN, 0 + PRINT, '### STEP 5B: SUCCESS!' + + ; =================================================================== + ; ### STEP 6: PHASE TO DISPLACEMENT AND GEOCODING ### + ; =================================================================== + + PRINT, '======================================================' + PRINT, '### STEP 6: PHASE TO DISPLACEMENT AND GEOCODING ###' + PRINT, '======================================================' + + final_upha_file = ref_output_root + '_upha' + geocoded_output_root = workflow_root_name + '_geo' + + task6_obj = OBJ_NEW('SARscapeBatch', Module='InSARPhaseToDisplacement') + + IF (~OBJ_VALID(task6_obj)) THEN BEGIN + PRINT, '错误: 创建 SARscapeBatch 对象失败: InSARPhaseToDisplacement' + RETURN, 0 + ENDIF + + ; --- Set Main Parameters --- + prefix_main6 = 'MAIN_INSAR_PHASE_TO_DISPLACEMENT_CMD.' + task6_obj->SetParam, prefix_main6 + 'INPUT_FILE_NAME', final_upha_file + task6_obj->SetParam, prefix_main6 + 'OUTPUT_FILE', geocoded_output_root + task6_obj->SetParam, prefix_main6 + 'COHERENCE_FILE', workflow_root_name + '_cc' + task6_obj->SetParam, prefix_main6 + 'DEM_FILE_NAME', dem_base_file + task6_obj->SetParam, prefix_main6 + 'COHERENCE_THRESHOLD', STRING(geocoding_coh_threshold) + + ; --- Set Other Required Parameters --- + task6_obj->SetParam, prefix_main6 + 'INTERPOL_TYPE' , '4th_order_cc' + task6_obj->SetParam, prefix_main6 + 'ALLOW_SKIP_REFINEMENT' , 'NotOK' + task6_obj->SetParam, prefix_main6 + 'SARSCAPEENVIRONMENT' , 'IDL_ENVI_ENV' + + ; --- Set Geocoding Parameters --- + prefix_geocode6 = 'GEOCODE_CMD.' + task6_obj->SetParam, prefix_geocode6 + 'GEOCODE_RG_GRID_SIZE', STRING(geocoding_pixel_size_m) + task6_obj->SetParam, prefix_geocode6 + 'GEOCODE_AZ_GRID_SIZE', STRING(geocoding_pixel_size_m) + task6_obj->SetParam, prefix_geocode6 + 'GEOCODE_INTERPOL_BOX_SIZE', '7' + + ; --- Set Displacement Product Generation Parameters --- + prefix_disp6 = 'DISPLACEMENT_PROJECTION_CMD.' + task6_obj->SetParam, prefix_disp6 + 'GENERATE_LOS_FLAG', 'OK' + task6_obj->SetParam, prefix_disp6 + 'GENERATE_VERTICAL_FLAG', 'NotOK' + task6_obj->SetParam, prefix_disp6 + 'GENERATE_MAX_SLOPE_FLAG', 'NotOK' + + PRINT, '正在执行地理编码任务...' + ok_step6 = task6_obj->Execute() + OBJ_DESTROY, task6_obj + IF NOT ok_step6 THEN RETURN, 0 + PRINT, '### STEP 6: SUCCESS!' + + ; --- Final Success Message --- + PRINT, '=======================================================' + PRINT, ' FULL D-InSAR WORKFLOW (STEPS 1-6) COMPLETED SUCCESSFULLY!' + PRINT, '=======================================================' + + RETURN, 1 + +END + +;================================================================================ +; ### 批量D-InSAR工作流主程序 ### +;================================================================================ +PRO Batch_DinsarWorkflow_ALL + + ; --- 初始化SARscape环境 --- + temporary_directory = 'D:\Sarscape_IDL_Area' + SARscape_Batch_Init, Temp_Directory=temporary_directory + e = ENVI(/HEADLESS) + PRINT, 'ENVI 和 SARscape 环境初始化成功！' + PRINT, '' + + ; =================================================================== + ; ### 全局配置 ### + ; =================================================================== + rootDir = 'G:/Province_Result/new' ; 修改为您的根目录 + dem_base_file = 'F:/Data/SRTM30m/SRTMDEM_RSP_SARscape' ; DEM文件路径 + + ; 处理参数配置 + target_ground_resolution_m = 10.0 + filter_method = 'GOLDSTEIN' + unwrapping_coh_threshold = 0.05 + gcp_coh_threshold = 0.7 + gcp_number = 100 + geocoding_coh_threshold = 0.0 + geocoding_pixel_size_m = 10.0 + + ; 最大处理数量 (0表示处理所有) + numToProcess = 0 + + ; =================================================================== + ; ### 查找并处理所有Task文件夹 ### + ; =================================================================== + + pathSep = '' + IF (!VERSION.OS_FAMILY EQ 'Windows') THEN pathSep = '\' ELSE pathSep = '/' + + IF rootDir EQ '' OR FILE_TEST(rootDir, /DIRECTORY) EQ 0 THEN BEGIN + PRINT, '错误: "rootDir" 参数未设置或指定的目录不存在: ' + rootDir + RETURN + ENDIF + + PRINT, '==============================================' + PRINT, '开始批量D-InSAR处理...' + PRINT, '扫描目录: ' + rootDir + PRINT, '==============================================' + + ; 查找所有Task_开头的文件夹 + PRINT, '正在扫描Task_开头的文件夹...' + taskFolders = FILE_SEARCH(rootDir + pathSep + 'Task_*', /FOLDERS, COUNT=nTasks) + + IF nTasks EQ 0 THEN BEGIN + PRINT, '在指定目录下没有找到任何Task_开头的文件夹。' + RETURN + ENDIF + + PRINT, '共找到 ' + STRING(nTasks) + ' 个Task文件夹，准备开始处理...' + PRINT, '' + + processed_count = 0 + + ; 循环处理每个Task文件夹 + FOREACH taskFolder, taskFolders DO BEGIN + IF (numToProcess GT 0) AND (processed_count GE numToProcess) THEN BEGIN + PRINT, '' + PRINT, '已达到指定的处理数量 (' + STRING(numToProcess) + ')，程序终止。' + BREAK + ENDIF + + taskName = FILE_BASENAME(taskFolder) + PRINT, '=== 正在处理任务文件夹: ' + taskName + ' ===' + + ; 检查master和slave子文件夹 + masterDir = taskFolder + pathSep + 'master' + slaveDir = taskFolder + pathSep + 'slave' + + IF (NOT FILE_TEST(masterDir, /DIRECTORY)) OR (NOT FILE_TEST(slaveDir, /DIRECTORY)) THEN BEGIN + PRINT, '警告: ' + taskName + ' 中缺少master或slave文件夹，跳过。' + PRINT, '' + CONTINUE + ENDIF + + ; 查找master和slave文件夹中的SML文件 + masterSmlFiles = FILE_SEARCH(masterDir + pathSep + '*.sml', COUNT=nMaster) + slaveSmlFiles = FILE_SEARCH(slaveDir + pathSep + '*.sml', COUNT=nSlave) + + IF (nMaster EQ 0) OR (nSlave EQ 0) THEN BEGIN + PRINT, '警告: ' + taskName + ' 中未找到足够的SML文件 (master: ' + STRING(nMaster) + ', slave: ' + STRING(nSlave) + ')，跳过。' + PRINT, '' + CONTINUE + ENDIF + + ; 取第一个找到的SML文件作为主从影像 + master_base_file = masterSmlFiles[0] + slave_base_file = slaveSmlFiles[0] + + ; 移除.sml扩展名，得到基础文件名 + master_base_file = STRMID(master_base_file, 0, STRLEN(master_base_file) - 4) + slave_base_file = STRMID(slave_base_file, 0, STRLEN(slave_base_file) - 4) + + ; 设置工作流输出根名称（在Task文件夹内创建结果） + workflow_root_name = taskFolder + pathSep + 'dinsar_results' + pathSep + 'workflow' + + ; 创建输出目录 + outputDir = FILE_DIRNAME(workflow_root_name) + IF NOT FILE_TEST(outputDir, /DIRECTORY) THEN BEGIN + FILE_MKDIR, outputDir + PRINT, '-> 创建输出目录: ' + outputDir + ENDIF + + PRINT, '-> Master文件: ' + FILE_BASENAME(master_base_file) + PRINT, '-> Slave文件: ' + FILE_BASENAME(slave_base_file) + PRINT, '-> 输出路径: ' + workflow_root_name + PRINT, '' + + ; 执行单个D-InSAR工作流 + success = Execute_Single_Dinsar_Workflow(master_base_file, slave_base_file, dem_base_file, workflow_root_name, $ + target_ground_resolution_m, filter_method, unwrapping_coh_threshold, $ + gcp_coh_threshold, gcp_number, geocoding_coh_threshold, geocoding_pixel_size_m) + + IF success THEN BEGIN + PRINT, '### ' + taskName + ': D-InSAR工作流完成! ###' + processed_count = processed_count + 1 + ENDIF ELSE BEGIN + PRINT, '### ' + taskName + ': D-InSAR工作流失败! ###' + ENDELSE + + PRINT, '' + ENDFOREACH + + PRINT, '==============================================' + PRINT, '批量D-InSAR处理完成！' + PRINT, '总共成功处理了 ' + STRING(processed_count) + ' 个Task文件夹。' + PRINT, '==============================================' + + SARscape_Batch_Exit + +END \ No newline at end of file