GMTSAR应用实例—基于Sentinel-1在加拿大魁北克省区域的形变监测

已于 2025-03-01 22:37:23 修改
阅读量838 收藏 25
点赞数 18
文章标签： sentinel
于 2025-02-27 16:31:34 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/zhaojinbiao0716/article/details/145853846
版权

  • 前言

学习记录，欢迎交流。后续会对结果进行分析并且调整参数。

  • 设备参数

宿主机：Windows 10 16G运存

虚拟机：VMware Workstation Pro 16

Ubuntu系统：20.04

  • 软件安装

可参考GMTSAR安装教程-CSDN博客

  • 数据准备

1、Sentinel-1A数据：目前下载该数据多使用美国阿拉斯加ASF平台 

   ASF平台：ASF Home | Alaska Satellite Facility

2、精密轨道数据：

   Copernicus Browser：Copernicus Browser

3、DEM数据：GMTSAR提供的网页下载已经停止支持，可以通过以下命令下载

make_dem.csh W E S N [mode]

其中mode 1:SRTM-1s 2:SRTM-3s

使用上述命令下载DEM时需注意：下载范围经纬度最大为4×4（本次实验时将6×4的范围分为两个3×4的区域进行下载）

下载好后的DEM需使用gmt grdpaste拼接，示例命令如下：

gmt grdpaste file_a.grd file_b.grd -Goutfile.grd

注意：DEM覆盖范围一定要足够大，保证可以覆盖处理的条带，不然很容易报错，尤其是有几景影像覆盖区差异偏大的情况下。

  • 操作流程

1、文件夹部署

1.1、按照指导手册依次创建各级文件夹：

1.2、各文件夹放置内容

（1）将解压后的Sentinel-1﹡.SAFE文件放置在data文件夹下

（2）将精密轨道文件﹡.EOF文件放置在orbit文件夹下

（3）将拼接后的DEM数据﹡.grd文件放置在主目录下的topo文件夹下

（4）F1、F2、F3文件夹下应有intf_in、raw、SLC、topo文件夹和batch_tops.config文件

         其中batch_tops.config文件可从/usr/local/GMTSAR/gmtsar/csh中获得

以下脚本可实现自动构建文件夹

#!/bin/bash

# 定义基础目录
base_dir="asc"

# 创建主目录结构
mkdir -p "${base_dir}/data"
mkdir -p "${base_dir}/orbit"
mkdir -p "${base_dir}/SBAS"
mkdir -p "${base_dir}/merge"
mkdir -p "${base_dir}/topo"
mkdir -p "${base_dir}/reframed"

# 创建 F1, F2, F3 目录及其子目录
for folder in F1 F2 F3; do
    mkdir -p "${base_dir}/${folder}/raw"
    mkdir -p "${base_dir}/${folder}/SLC"
    mkdir -p "${base_dir}/${folder}/intf_in"
    mkdir -p "${base_dir}/${folder}/topo"

    # 创建软链接（路径设置最好为系统文件外的文件，不然是只读模式无法修改参数）
    ln -s /home/insar/GMTSAR/gmtsar/csh/batch_tops.config "${base_dir}/${folder}/batch_tops.config"
done

echo "文件夹结构和软链接已成功创建。"

2、文件链接

2.1、dem.grd链接至F1/topo、F2/topo、F3/topo与merge

进入目标文件夹运行命令

ln -s ../../topo/dem.grd

2.2、链接影像

Sentinel-1A数据是由三个IW条带组成，每个条带都会有9个左右的burst，这也是文件夹分成F1、F2、F3的原因，其中数字1、2、3分别对应条带1、2、3，然后我们需要知道自己的研究区落在哪个条带，然后开始链接数据，命令如下

ln -s ../../data/*.SAFE/measurement/*iw1*vv*xml
ln -s ../../data/*.SAFE/annotation/*iw1*vv*tiff
ln -s ../../orbit/*EOF

 逐一链接各文件夹数据较为费时间，所以编写了一个批量链接脚本create_link.sh如下：

#!/bin/bash

# 定义路径
source_dir="/media/insar/Elements/asc/data/*.SAFE"
orbit_dir="/media/insar/Elements/asc/orbit"
target_dir="/media/insar/Elements/asc/F1/raw"
output_file="SAFE"


# 链接dem.grd  ln -s /media/insar/Elements/asc/topo/dem.grd
# 检查目标目录是否存在，如果不存在则创建
if [ ! -d "$target_dir" ]; then
    mkdir -p "$target_dir"
    echo "创建目标目录: $target_dir"
fi

# 清空或创建 SAFE 文件
> "$output_file"

# 遍历每个 .safe 目录
for safe_dir in $source_dir; do
    if [ -d "$safe_dir" ]; then
        # 查找 measurement 目录中的 .tiff 文件
        if [ -d "$safe_dir/measurement" ]; then
            ls "$safe_dir/measurement"/s1a-iw1-slc-vv-*.tiff >> "$output_file" 2>/dev/null
        else
            echo "警告：目录不存在 $safe_dir/measurement"
        fi

        # 查找 annotation 目录中的 .xml 文件
        if [ -d "$safe_dir/annotation" ]; then
            ls "$safe_dir/annotation"/s1a-iw1-slc-vv-*.xml >> "$output_file" 2>/dev/null
        else
            echo "警告：目录不存在 $safe_dir/annotation"
        fi
    else
        echo "警告：目录不存在 $safe_dir"
    fi
done

# 查找 orbit 目录中的 .EOF 文件
if [ -d "$orbit_dir" ]; then
    #ls也可以用 
    #ls "$orbit_dir"/*.EOF >> "$output_file" 2>/dev/null  
    find "$orbit_dir" -type f -name "*.EOF" >> "$output_file"
else
    echo "警告：目录不存在 $orbit_dir"
fi

# 检查 SAFE 文件是否为空
if [ ! -s "$output_file" ]; then
    echo "未找到任何 .tiff 或 .xml 文件！"
    exit 1
fi

# 批量创建符号链接
cat "$output_file" | while read -r file_path; do
    if [ -e "$file_path" ]; then
        ln -sf "$file_path" "$target_dir/"
        echo "已创建链接: $file_path -> $target_dir/"
    else
        echo "文件不存在: $file_path"
    fi
done

echo "符号链接创建完成！所有文件已链接到 $target_dir。"

3、建立影像精轨对应文件data.in

链接完文件后，需要建立一个“data.in”文件，里面放置数据名和精密轨道名对应起来，原指导手册中未提供自动化脚本，自行编写处理脚本如下：

#!/bin/bash

# 定义路径
tiff_dir="/media/insar/Elements/asc/F1/raw/*.tiff"
eof_dir="/media/insar/Elements/asc/orbit/*.EOF"
output_file="/media/insar/Elements/asc/F1/raw/data.in"

# 获取 .tiff 文件列表
tiff_files=($tiff_dir)

# 获取 .EOF 文件列表
eof_files=($eof_dir)

# 检查文件数量是否匹配
if [ ${#tiff_files[@]} -ne ${#eof_files[@]} ]; then
    echo "错误：.tiff 文件和 .EOF 文件数量不匹配！"
    exit 1
fi

# 清空或创建输出文本
> "$output_file"

# 按顺序配对并输出
for ((i=0; i<${#tiff_files[@]}; i++)); do
    # 提取文件名（去掉路径和扩展名）
    tiff_name=$(basename "${tiff_files[$i]}" .tiff)
    eof_name=$(basename "${eof_files[$i]}")

    # 输出配对结果
    echo "$tiff_name:$eof_name" >> "$output_file"
done

echo "配对结果已写入$output_file文本！"

4、选择主影像

4.1、这里需要在raw文件下输入命令行preproc_batch_tops.csh data.in dem.grd 1，生成一个时空极限列表和图像，所有生成文件如下图所示，然后运行mv baseline_table.dat ../将baseline_table.dat文件复制到上一目录F1下。

preproc_batch_tops.csh data.in dem.grd 1
mv baseline_table.dat ../

4.2、时空基线图像如下图所示，我们需要选取一个主影像master，主影像的位置大约在时空基线的中间，此次实验S1_20230704_ALL_F1就比较符合要求。

4.3、然后你需要将你所选的主影像，此次实验为放S1_20210427_ALL_F1在data.in文件的最上方，因为系统会默认为第一个为主影像，如下图。 

4.4、然后输入命令行preproc_batch_tops.csh data.in dem.grd 2，再跑一次模式2，该模式为配置工作，如果有F2或者F3需要，要进行同样的操作。

preproc_batch_tops.csh data.in dem.grd 2

5、建立时空干涉对

5.1、这里需要返回F1文件夹，输入时空参数建立干涉对，本次实验中，设置时间基线150，空间基线100，参数没有固定要求，可按照需求构建，干涉对数量不宜过少，大致为全部数据时间跨度的一半就可以。之后输入命令行如下：

select_pairs.csh baseline_table.dat 150 100

（1）生成时空基线图_baseline.ps

（2）干涉对列表_intf.in

6、干涉

6.1、运行一个干涉图来测试设置

（1）在每个子路径中运行一次测试干涉图，以便确保一切设置正确，并且通过计算一次 topo_ra.grd 文件来优化处理过程，而不是对每个干涉图进行多次计算。

（2）编辑 batch_tops.config 文件，设置首次运行

进入F1文件下，设置batch_tops.config文件，主要参数如下：

master_image = S1_20230704_ALL_F1
Proc_stage = 1
shift_topo = 0
filter_wavelength = 200
range_dec = 8
azimuth_dec = 2
threshold_snaphu = 0
hreshold_geocode = 0
region_cut = 0/10000/10000/13400
#为了减少工作量，这里可以进行研究区裁剪，这里用到的范围是雷达坐标系中的Range和Azimuth，可先进行一个干涉对的全F1条带干涉，然后使用结果图的范围来确定的裁剪范围。

（3）生成测试干涉图和必要的运行文件

我们需要将地形图映射到相位中，供主影像使用，因此编辑 batch_tops.config，设置 proc_stage = 1 并制作单个干涉图。这样就会生成一个大的 trans.dat 文件以及一个 topo_ra.pdf 文件，您可以对其进行检查。这一步大约用时~1小时。

使用命令head -1 intf.in > one.in选取一个干涉对单独放置到one.in文件下，然后使用命令intf_tops.csh one.in batch_tops.config干涉。这一步会在F1文件下的topo文件下生成trans.dat文件，这个后续地理编码的时候会用到，

head -1 intf.in > one.in
intf_tops.csh one.in batch_tops.config

6.2、进行全部干涉对处理

进行完一个干涉对之后，将配置文件batch_tops.config里的proc_stage改为2，然后输入intf_tops.csh intf.in batch_tops.config。如果想要多个进程并行处理可以输入命令intf_tops_parallel.csh intf.in batch_tops.config 3这样就可以三个干涉对同时进行，取决于你的电脑配置，如果提示安装parallel工具，输入命令apt get install parallel就可以。全部跑之后，在intf_all文件下会有生成的干涉对，此次实验有29个干涉对。

intf_tops.csh intf.in batch_tops.config
intf_tops_parallel.csh intf.in batch_tops.config 3
apt get install parallel

由于电脑运行内存较小，本次实验未使用并行处理，逐条处理用时约1~2小时。

7、 跨子路径合并

本次实验未涉及，可参考技术手册sentinel_time_series_5.pdf中的第七步Merging Across Subswaths

8、解缠干涉图

解缠是 InSAR 时间序列处理过程中最耗时、计算量最大的步骤之一，因此我们要想方设法限制对感兴趣区域内不必要区域或像素的处理。在解包之前，需要确定的一个关键问题是所关注区域的相干性如何变化，因为相干性低会导致位移估算精度降低，解缠也会更加困难。

8.1、相干性较差区域解缠

可以采取几个步骤来减少相干性差的像素对解包的影响。第一个步骤是堆叠相干性网格，第二个步骤是褶皱掩膜，掩膜掉相干性低的像素（这种情况经常出现在水体区域）。

（1）堆叠相干性网格

可以通过以下两种方法创建 corr_stack.grd 文件：（1）对输入的部分相干性文件（corr.grd）求平均值，或（2）运行发行版中提供的堆叠脚本。本实验使用堆叠脚本

在上一步处理后的文件夹中运行stack.csh，首先将干涉文件夹中的相干文件corr.grd 路径写入 corr.grd_list

ls 201*/corr.grd > corr.grd_list
stack.csh corr.grd_list 1 corr_stack.grd std.grd
#"1"是缩放因子（本例中不缩放）,"corr_stack.grd"和 "std.grd "是输出文件（我们只使用 "corr_stack.grd"）

创建 corr_stack.grd 后，可使用命令查看它 

ncview corr_stack.grd

（2）创建相干性掩膜

根据堆叠相干性创建一个水体掩膜。这将屏蔽掉任何不符合我们设置的阈值的数据（这里我们选择了 0.07，但可以随意尝试）。我们选择的值屏蔽了大部分水体，只屏蔽了一小部分陆地。

gmt grdmath corr_stack.grd 0.07 GE 0 NAN = mask_def.grd

（3）对每张干涉图进行解缠

ⅰ.在干涉图目录中，首先列出干涉图列表，将其输入脚本：

ls -d 2023* > intflist

ⅱ.创建以下脚本并将其命名为 "unwrap_intf.csh”

#!/bin/csh -f
# intflist contains a list of all date1_date2 directories.
cd $1
ln -s ../mask_def.grd .
snaphu_interp.csh 0.001 40
cd ..

将该脚本放入GMTSAR环境变量设置的文件夹中，否则运行时会显示未找到该命令，移动脚本时可能会由于权限不足导致无法粘贴，可使用如下命令打开的窗口复制粘贴

sudo nautilus

ⅲ.开始解缠

为了节省时间，还可以选择与 unwrap_parallel.csh 脚本并行运行。这将调用你刚刚创建的脚本（unwrap_intf.csh），因此请确保你正确命名了它。同样，6 是线程数。这意味着多个干涉图将同时解包，从而减少完成所有干涉图所需的时间。例如，逐个运行每个干涉图所需时间约为 1 小时。本次实验只解缠一个条带，运行时间约9小时。

unwrap_parallel.csh intflist 1
#由于运存限制，逐张解缠，设置并行处理参数为1

9、运行小基线集 (SBAS) 分析

9.1、准备输入文件

GMTSAR 发行版中添加了一个名为 prep_sbas.csh 的脚本，可以帮您完成这项工作。该脚本会创建 scene.tab 文件和 intf.tab 文件，并为用户输出默认的 SBAS 命令结构。

（1）将干涉图列表、基线列表链接入SBAS文件夹

（2）运行脚本

prep_sbas.csh intf.in baseline_table.dat ../merge unwrap.grd corr.grd
# intf.in干涉图列表
# baseline_table.dat 基线列表
# ../merge解缠干涉图目录的路径
# unwrap.grd解缠文件名
# corr.grd相干性文件名

9.2、Run SBAS

在命令行中键入 “sbas”，可阅读有关用法和命令示例。下面是指导手册的命令示例（全部为一行）：

sbas intf.tab scene.tab 194 34 4612 3300 -range 888122 -incidence 40 -wavelength 0.0554658 -smooth 5.0 -rms -dem
Usage: sbas intf.tab scene.tab N S xdim ydim [-atm ni] [-smooth sf] [-wavelength wl] [-incidence inc] [-range -rng] [-rms] [-dem]

# 194是干涉图的数量
# 34是SAR影像数量
# 为了找到xdim (here is 4612) and ydim (3300) (x and y dimension of interferograms)使用以下命令

gmt grdinfo ../merge/date1_date2/unwrap.grd

# x n_columns” is the value for x dim
# y n_rows” is the value for y dim
# wavelength:雷达波的波长（单位：米）。要找到这个值，请在 supermaster.PRM 中查找"radar_wavelength"
# incidence：默认值为 37 度，但这里使用的值基本上无关紧要（大多数情况下使用默认值即可）
# range: 雷达到干涉图中心的距离，可通过以下公式计算

Range = ({[(speed of light) / (rng_samp_rate) / 2] * ((x_min+x_max)/2)} /2 ) + near_range

# Speed of light is ~3 x 10^8 m/s
# “rng_samp_rate” is in supermaster.PRM
# X_min and x_max come from the gmt grdinfo command above
# “near_range” is in supermaster.PRM

9.3、SBAS后处理

（1）运行 SBAS 后，会得到一个 vel.grd 文件（以雷达坐标和毫米/年为单位）。要将其转换为纬度/经度坐标，逐条执行以下命令：

ln -s ../topo/trans.dat .
ln -s ../F1/intf_all/2018127_2018163/gauss_*
proj_ra2ll.csh trans.dat vel.grd vel_ll.grd
gmt grd2cpt vel_ll.grd -T= -Z -Cjet > vel_ll.cpt
grd2kml.csh vel_ll vel_ll.cpt

之后将看到各种形式的位移网格文件 disp_*.grd。这些都是在该时间点测得的位移（以毫米为单位），可以将这些网格文件中的点组成位移时间序列。 

  • 参考博客

基于GMTSAR软件的富士山时序形变探测（一）-CSDN博客

MagpieZ
关注
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值