【INSAR形变监测】SENTINEL-1 TOPS INSAR全流程 (ISCE2+MINTPY)

McWuuu

已于 2024-03-14 17:30:56 修改

阅读量1.7k

点赞数 17

文章标签： sentinel python

于 2024-03-14 17:02:40 首次发布

本文链接：https://blog.csdn.net/McWuuu/article/details/136715642

版权

SENTINEL-1 TOPS INSAR全流程 (ISCE2+MINTPY)

博主这段时间接到用InSAR做形变监测的任务，我虽然是遥感专业的，但对SAR、INSAR一窍不通，刚开始使用了ENVI+SARSCAPE，但我就很钻牛角尖，还是想要通过开源的方式实现这些东西，所以发现了ISCE+MintPy的解决方案。国内相关的教程很少，而且都比较碎片化或者不全，我折腾了好久才把这一套走通，踩坑无数，希望这个教程可以对大家有所帮助！

0、INSAR入门

推荐有兴趣的同学查看这个课程，里面有2021、2022、2023年 isce开发团队开的网课资料，包括jupyter notebook和视频，可以学习！

https://github.com/parosen/Geo-SInC/tree/main

一、安装软件

1. 制作镜像/linux安装软件

ISCE2只能运行在linux系统下，所以我们可以使用docker，docker安装在这里不详细赘述了。

官方有ISCE2的镜像，但并没有安装MintPy，所以我们先不用。安装好docker后，可以拉取官方的ubuntu镜像，然后按照张云俊博士（ISCE和MintPy主要贡献者）的步骤进行即可。参考

https://github.com/yunjunz/conda-envs

这里我也提供一个我安装好的docker镜像，可以直接使用，省一堆麻烦，直接去dockerhub pull：

docker pull mcwuuuuuuu/insarenv

镜像创建好（pull好）后，可以这样创建容器，-v是挂载本地文件夹到容器内

这些东西只是一个示范，大家可以自由发挥。

# 创建容器
docker run -it -d --name xxx -v C:\workspace\:/workspace/ mcwuuuuuuu/insarenv /bin/bash
# 进入容器
docker exec -it xxx /bin/bash

因为做insar时序分析非常耗存储空间，所以建议把所有文件都放在挂载的文件夹里，防止docker容器变得异常庞大，占用C盘空间

2. 设置环境变量

我linux比较菜，不知道为什么每次保存完环境变量重启容器后又没了 😦，所以进容器之后输入这行指令

alias load_insar='conda activate insar; source ~/insar_env_setup.rc'

启动这个命令

load_insar

输入load_insar后，会自动进入insar的conda环境，并且将MintPy、ISCE2、PyAps、PySolid、SSARA等软件添加至环境变量，方便直接调用。

本来这个是张云俊博士提供的bash.rc文件，存放在~/tools/conda-envs/insar/config.rc，但如果想要使用stackSentinel.py，需要单独加一行指令，即export PATH=${PATH}:${ISCE_STACK}/topsStack，本来是没有这一行的，因为isce2中几个stack功能是冲突的，如果使用别的stack功能，比如stripmapStack，则需要删除这个环境变量并重新添加，所以本身的config.rc中是没有的，我把它添加进了insar_env_setup.rc，默认大家是想要使用哨兵的stackSentinel.py，如果有朋友需要修改的话可以参照一下网址:

https://github.com/isce-framework/isce2/tree/main/contrib/stack

测试一下是否成功初始化，能运行就是成功了

topsApp.py -h            # test ISCE-2
smallbaselineApp.py -h   # test MintPy
stackSentinel.py -h      # test topsStack

二、下载数据

1. 使用SSARA下载SENTINEL-1数据

首先，在SSARA的配置文件里输入自己ASF的账号密码，如果你使用我的镜像，或者按照上面的教程一模一样安装自己的系统的话，文件在~/tools/utils/SSARA/password_config.py，如果自定义安装，则在自己的SSARA目录下。
在这里插入图片描述
填写中间这个asfuser和asfpass，其实就是NASA EARTHDATA的账号密码。

ASF地址：https://search.asf.alaska.edu/#/，进去也可以注册账号，这里不赘述，自行百度如何创建。

如果只是用sentinel的话，其他两个una和eosso的账号密码可以不填，如果有需要用其他SAR卫星，可以把其他两个也注册了然后填上。

输入账号密码后，切换到需要存储下载SENTINEL数据的文件夹

cd $worksdir/data/sentinel # 自己定义

调用SSARA中的ssara_federated_query.py进行查询，-p是卫星种类，-r是相对轨道，-f是轨道内的第n景，-s开始时间，-e结束时间，--print显示在命令行，下图是ASF手动查询结果，与ssara命令的对照图
在这里插入图片描述

以下是我查询的例子：

$SSARA_HOME/ssara_federated_query.py -p SENTINEL-1A,SENTINEL-1B -r 142 -f 101 -s 2023-01-01 -e 2024-01-15 --print

这里查询的是S1A/S1B，相对轨道142，第101景，从20230101 - 20240115之间的所有影像，结果（部分如下图所示）

ASF,Sentinel-1A,51139,2023-11-09T10:12:05.504081Z,2023-11-09T10:12:32.456535Z,142,647,647,IW,None,ASCENDING,R,VV+VH,https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20231109T101205_20231109T101232_051139_062AE5_F220.zip
ASF,Sentinel-1A,51314,2023-11-21T10:12:05.692476Z,2023-11-21T10:12:32.642875Z,142,647,647,IW,None,ASCENDING,R,VV+VH,https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20231121T101205_20231121T101232_051314_0630F7_A851.zip
ASF,Sentinel-1A,51489,2023-12-03T10:12:05.244565Z,2023-12-03T10:12:32.197019Z,142,647,647,IW,None,ASCENDING,R,VV+VH,https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20231203T101205_20231203T101232_051489_0636F1_3FDF.zip
ASF,Sentinel-1A,51664,2023-12-15T10:12:04.700112Z,2023-12-15T10:12:31.650511Z,142,647,647,IW,None,ASCENDING,R,VV+VH,https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20231215T101204_20231215T101231_051664_063D15_5082.zip
ASF,Sentinel-1A,51839,2023-12-27T10:12:03.661178Z,2023-12-27T10:12:30.611577Z,142,647,647,IW,None,ASCENDING,R,VV+VH,https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20231227T101203_20231227T101230_051839_064321_8E61.zip
ASF,Sentinel-1A,52014,2024-01-08T10:12:03.394552Z,2024-01-08T10:12:30.342895Z,142,647,647,IW,None,ASCENDING,R,VV+VH,https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20240108T101203_20240108T101230_052014_06491A_4C23.zip
...

查看觉得没问题，就可以用相同的命令下载了，把--print改成--download即可，后面可以加上--parallel=n来并行下载

$SSARA_HOME/ssara_federated_query.py -p SENTINEL-1A,SENTINEL-1B -r 142 -f 101 -s 2023-01-01 -e 2024-01-15 --download --parallel=4

在这里插入图片描述

2. 下载轨道数据

还是从ASF的archive下载精密轨道数据，网址：

https://s1qc.asf.alaska.edu/aux_poeorb/

这里的文件命名非常具有迷惑性，第一个日期是轨道数据处理完成时间，第二个日期是轨道数据获取的前一天，第三个日期是获取的后一天，而且网页上显示不全，所以找起来非常麻烦。
在这里插入图片描述

这里提供一个我写好的小脚本，需要修改两个地方，即input_dir和orbit_dir即可调用，主要是我太懒，不想再封装，凑合用。

# ! /usr/bin/env python
# @author: McWuuu
# @time: 2024/3/13

import os
import glob
import time
import requests
from datetime import datetime, timedelta

# archive地址
url = 'https://s1qc.asf.alaska.edu/aux_poeorb/'
# 获取全部信息
txt = requests.get(url).content
# 将文件名转换成列表
lst = [str(i[:77])[2:-1] for i in txt.split(b'href="')]

# TODO 下载好的tops数据存放文件夹
input_dir = '/workspace/insar/data/tops'
# TODO 创建轨道数据存放的文件夹
orbit_dir = os.path.join(os.path.dirname(input_dir), 'orbit')
os.makedirs(orbit_dir, exist_ok=True)
os.chdir(orbit_dir)

# 开始检索
for file in glob.glob(os.path.join(input_dir, 'S1A*.zip')):
    date_str = file.split('_')[-5][:8]
    y, m, d = int(date_str[:4]),  int(date_str[4:6]),  int(date_str[6:8])
    dt = datetime(y, m, d)

    prev_dt = dt - timedelta(days=1)
    next_dt = dt + timedelta(days=1)

    prev_dt_str = f'{prev_dt.year}{prev_dt.month:02d}{prev_dt.day:02d}'
    next_dt_str = f'{next_dt.year}{next_dt.month:02d}{next_dt.day:02d}'

    # search!
    for filename in lst:
        # 通过两个日期进行匹配
        if (filename[-35:-27] == prev_dt_str) & (filename[-19:-11] == next_dt_str):
            print(f'found matched orbit file for {file}: {filename}')
            print(f'start downloading {filename}')
            start = time.time()
            download_url = url + filename
            os.system(f'wget {download_url}')
            end = time.time()
            print(f'download finished, time used:{end - start} sec.')

三、ISCE2处理

1. 切换到输出文件夹

创建一个工作文件夹，用来存放所有的结果，并在命令行切换到这个文件夹

mkdir xxxxxx/work
cd xxxxxx/work

2. 准备DEM数据

首先创建一个.netrc文件，可以参考这里https://github.com/isce-framework/isce2的Notes on Digital Elevation Models这一节

cd ~
touch .netrc

然后在这个文件内输入你的NASA EARTHDATA账号密码，和上面需求的ASF的账号密码一致

machine urs.earthdata.nasa.gov
    login your_earthdata_login_name
    password your_earthdata_password

保护隐私

chmod go-rwx .netrc

设置好账号密码后，可以开始下载dem，建议先切换到刚刚创好的工作文件夹中，可以通过自己的研究区域边界下载（这里需要往外扩大至整数），如：

本来我们的研究区是：

32.39 32.90 116.3 118.13 #SNWE

扩大后：

dem.py -a stitch -b 32 33 116 119 -r -s 1 -c

其中，-a 表示操作，可选的下载：download和拼接：stitch，-b是边界经纬度，后面是SNWE，-r表示是否显示下载状态，-s标识DEM源，默认1即可， -c是转换为WGS84，成功下载会像这样，如果选择-a download，则不会直接拼接
在这里插入图片描述

Tips 完成后如果你没有按我说的切换到自己创建的工作文件夹，则需要把拼接好的DEM文件复制到工作文件夹中

3. 运行stackSentinel.py

使用topsStack中的stackSentinel.py进行处理，安装部分可以看上面的步骤。

stackSentinel.py -s ../../data/tops -d ./demLat_N32_N33_Lon_E116_E119.dem.wgs84 -o ../../data/orbit -a ../../data/aux_cal -b '32.39 32.90 116.3 118.13' -c 1 -m 20230101

-d: DEM路径(其实没用，因为如果工作文件夹没有DEM文件的话，程序会报错，应该是有bug，所以必须要把DEM文件复制到工作文件夹中！)

-m 第一天日期

-c 每一景和最近的n景数据组成干涉对，1就是一一对应，比如第1景和第2景，第2景和第3景…，如果是2的话，第1和第2、第1和第3、第2和第3、第2和第4…所以如果设置很高的话会非常占空间，且很慢。

-s: SENTINEL数据存放文件夹路径

-o: 轨道数据文件夹路径

-a: AUX_CAL，辅助文件主要是用来天线方向图校正，IPF(Instrument Processing Facility)仪器处理设备如果是同一个版本的话则不需要辅助文件。主要使用的是2015年3月之前，也就是说当2015年3月之后和之前的图像生成干涉图时会导致相位偏移的变化，因此需要此文件进行校正。但是还是需要文佳夹Aux存在，可以创一个空的文件夹或者直接写别的文件夹即可，这些内容可以参考这篇博客

https://blog.csdn.net/QJBMFCH/article/details/124060861

正常来说运行完上面的指令会在目前所处的文件夹中生成run_files文件夹，并在其中生成16个步骤的指令文件，如下所示：
在这里插入图片描述

可以一个一个运行，试错

bash run_files/run_01_unpack_topo_reference

或者参考我这样写一个脚本把16个依次跑一遍

# ! /usr/bin/env python
# @author: McWuuu
# @time: 2024/3/13

import os
import time
import glob

cwd = os.getcwd()

scripts = glob.glob(os.path.join(cwd, 'run_files', 'run_*'))

if len(scripts) == 16:
    for i in range(1, 17):
        for script in scripts:
            if f'run_{i:02d}' in script:
                print(f'start run {script}')
                start = time.time()
                os.system(f'bash {script}')
                end = time.time()
                print(f'{script} done, time:{end - start} sec')

开始跑，如果很多景的话会很慢的，建议挂机做别的事
在这里插入图片描述

跑完之后的结果会像这个结构，忽略这个mintpy文件夹，这是下一步使用mintpy时，我们创建的文件夹。

在这里插入图片描述

建议按照mintpy的要求查看是否有这些文件：https://github.com/insarlab/MintPy/blob/main/docs/dir_structure.md，当然日期肯定不用按照这个来，主要看merged里的东西是否齐全，否则无法使用mintpy进行分析。

$DATA_DIR/GalapagosSenDT128
├── baselines
│   ├── 20141213_20141225
│   │   └── 20141213_20141225.txt
│   └── ...
├── reference
│   ├── data.rsc    #generated by prep_isce.py
│   ├── IW1
│   ├── IW1.xml
│   ├── IW2
│   └── IW2.xml
├── merged
│   ├── geom_reference
│   │   ├── hgt.rdr
│   │   ├── hgt.rdr.full.vrt
│   │   ├── hgt.rdr.full.xml
│   │   ├── hgt.rdr.vrt
│   │   ├── hgt.rdr.xml
│   │   ├── lat.rdr
│   │   ├── lat.rdr.full.vrt
│   │   ├── lat.rdr.full.xml
│   │   ├── lat.rdr.vrt
│   │   ├── lat.rdr.xml
│   │   ├── lon.rdr
│   │   ├── lon.rdr.full.vrt
│   │   ├── lon.rdr.full.xml
│   │   ├── lon.rdr.vrt
│   │   ├── lon.rdr.xml
│   │   ├── los.rdr
│   │   ├── los.rdr.full.vrt
│   │   ├── los.rdr.full.xml
│   │   ├── los.rdr.vrt
│   │   ├── los.rdr.xml
│   │   ├── shadowMask.rdr
│   │   ├── shadowMask.rdr.full.vrt
│   │   ├── shadowMask.rdr.full.xml
│   │   ├── shadowMask.rdr.vrt
│   │   └── shadowMask.rdr.xml
│   └── interferograms
│       ├── 20141213_20141225
│       │   ├── filt_fine.cor
│       │   ├── filt_fine.cor.vrt
│       │   ├── filt_fine.cor.xml
│       │   ├── filt_fine.unw
│       │   ├── filt_fine.unw.conncomp
│       │   ├── filt_fine.unw.conncomp.vrt
│       │   ├── filt_fine.unw.conncomp.xml
│       │   ├── filt_fine.unw.vrt
│       │   ├── filt_fine.unw.xml
│       │   ├── ...
│       ├── 20141213_20150307
│       └── ...
├── secondarys
│   ├── 20141225
│   │   ├── IW1
│   │   ├── IW1.xml
│   │   ├── IW2
│   │   └── IW2.xml
│   ├── 20150307
│   └── ...
└── mintpy
    └── GalapagosSenDT128.txt

四、MintPy做时序分析

INSAR时序分析软件MintPy：https://github.com/insarlab/MintPy

1. 编写inputConfig文件

默认情况下，MintPy会读取~/tools/MintPy/src/mintpy/defaults/smallbaselineApp.cfg里的设置运行，这个文件里设置非常繁琐，所以我们可以选择一些我们需要改动的东西，存到另一个设置文件中，如下文所示，mintpy会优先读取我们的新设置，并将新设置的参数覆盖原设置参数。

首先，在我们存放ISCE2结果的文件夹中创建mintpy的文件夹，用来存放mintpy的结果，然后切换到mintpy文件夹。

mkdir mintpy
cd mintpy

在这里面创建一个cfg或者txt文件，将以下内容复制进去：

# vim: set filetype=cfg:
# mintpy.load.processor      = aria  #[isce, aria, hyp3, gmtsar, snap, gamma, roipac], auto for isce
#---------interferogram datasets:
# mintpy.load.unwFile        = ../stack/unwrapStack.vrt
# mintpy.load.corFile        = ../stack/cohStack.vrt
# mintpy.load.connCompFile   = ../stack/connCompStack.vrt

#---------geometry datasets:
mintpy.load.processor        = isce
##---------for ISCE only:
mintpy.load.metaFile         = ../reference/IW*.xml
mintpy.load.baselineDir      = ../baselines
##---------interferogram datasets:
mintpy.load.unwFile          = ../merged/interferograms/*/filt_*.unw
mintpy.load.corFile          = ../merged/interferograms/*/filt_*.cor
mintpy.load.connCompFile     = ../merged/interferograms/*/filt_*.unw.conncomp
##---------geometry datasets:
mintpy.load.demFile          = ../merged/geom_reference/hgt.rdr
mintpy.load.lookupYFile      = ../merged/geom_reference/lat.rdr
mintpy.load.lookupXFile      = ../merged/geom_reference/lon.rdr
mintpy.load.incAngleFile     = ../merged/geom_reference/los.rdr
mintpy.load.azAngleFile      = ../merged/geom_reference/los.rdr
mintpy.load.shadowMaskFile   = ../merged/geom_reference/shadowMask.rdr

# 地面控制点经纬度，这个可以不指定，Mintpy会随机选择一个范围内的点当作控制点（直接删除这一行或者注释掉）
mintpy.reference.lalo        = 32.52, 116.92

mintpy.unwrapError.method = bridging
########## 6. correct_troposphere (optional but recommended)
# 如果是auto,会用PyAPS下载ERA5的在分析数据进行矫正，这个下载逻辑有大问题，所以先不搞了。。。
mintpy.troposphericDelay.method = no

如果你和我一样在ISCE2的结果文件夹中创建了mintpy并进入此文件夹，可以直接照抄这个设置，其实配置文件非常复杂，如果有需求，可以详细查看mintpy本身自带的超详细的配置文件：

~/tools/MintPy/src/mintpy/defaults/smallbaselineApp.cfg

或者直接

smallBaselineApp.py -h # 查看帮助

如果你像我一样只想做基础的，那直接复制即可。

2. 命令行跑MintPy

确保自己在mintpy文件夹，直接运行，耐心等待结果即可。

smallBaselineApp.py inputConfig.cfg

如果想要一步一步运行，参考这个notebook https://github.com/parosen/Geo-SInC/blob/main/EarthScope2023/5.2_Intro_to_MintPy/smallbaselineApp_aria.ipynb

在这里插入图片描述

在这里插入图片描述
跑完后长这样：

/pic

在这里插入图片描述

/geo

在这里插入图片描述

3. 转TIFF

使用save_gdal.py，不想打字了，直接上图

在这里插入图片描述

save_gdal.py geo_velocity.h5 -d velocity -o velocity.tif --of GTiff

在这里插入图片描述

结果展示

在这里插入图片描述

McWuuu

关注

17
点赞
踩
40

收藏

觉得还不错? 一键收藏
19
评论
【INSAR形变监测】SENTINEL-1 TOPS INSAR全流程 (ISCE2+MINTPY)

【INSAR形变监测】SENTINEL-1 TOPS INSAR全流程 (ISCE2+MINTPY)
复制链接

扫一扫

【INSAR形变监测】SENTINEL-1 TOPS INSAR全流程 (ISCE2+MINTPY)

SENTINEL-1 TOPS INSAR全流程 (ISCE2+MINTPY)

0、INSAR入门

一、安装软件

1. 制作镜像/linux安装软件

2. 设置环境变量

二、下载数据

1. 使用SSARA下载SENTINEL-1数据

2. 下载轨道数据

三、ISCE2处理

1. 切换到输出文件夹

2. 准备DEM数据

3. 运行stackSentinel.py

四、MintPy做时序分析

1. 编写inputConfig文件

2. 命令行跑MintPy

3. 转TIFF

“相关推荐”对你有帮助么？