跳转至

Phase unwrapping

Notebook

Converted from nbs/Tutorials/CLI/CampiFlegrei/05_unwrap.ipynb (text and code cells; the outputs are not stored in the repository). Run it in Jupyter or VS Code for the interactive maps.

Data: Campi Flegrei. Sentinel-1A IW, descending track 22, VV, 92 dates from 2019-01-02 to 2021-12-29 (every 12 days, reference 2020-07-07), cropped to the caldera (Pozzuoli and the Gulf of Pozzuoli) at full resolution. A bell-shaped uplift of tens of centimetres in the centre; the town gives many PS. The data and how they were made are described in data/CampiFlegrei/README.md.

In this tutorial, we unwrap the interferograms on the refined measurement points (PS+DS) with the extended minimum cost flow (EMCF): a redundant network of image pairs (every image with the next three) is unwrapped together, and the loops of image pairs that still do not add up to zero are corrected by phase closure, which gives the unwrapped phase of every image. More details and the expected results: docs/workflows/05_unwrap.md.

import moraine.cli as mc
import moraine as mr
import zarr
import numpy as np
import toml
logger = mc.get_logger()
!mkdir phase_unwrapping -p

Inputs: the phase history and the hilbert index of the refined points, the raster shape and the pixel spacings (in meters, they convert the grid index to distances).

pc_ph = './pixel_refinement/hix/pc/pc_ph.zarr'
pc_hix = './pixel_refinement/hix/pc/pc_hix.zarr'
pc_e = './pixel_refinement/hix/pc/pc_e.zarr/'
pc_n = './pixel_refinement/hix/pc/pc_n.zarr/'

pc_gix = './phase_unwrapping/pc_gix.zarr'
pc_unw = './phase_unwrapping/pc_unw.zarr'
pc_ts = './phase_unwrapping/pc_ts.zarr'
pc_misclosure = './phase_unwrapping/pc_misclosure_fraction.zarr'
pc_change = './phase_unwrapping/pc_change_fraction.zarr'
pc_region = './phase_unwrapping/pc_region.zarr'

meta = toml.load('./load_data/meta.toml')
spacing = dict(range_pixel_spacing=meta['range_pixel_spacing'], azimuth_pixel_spacing=meta['azimuth_pixel_spacing'])
shape = zarr.open('./load_data/rslc.zarr', mode='r').shape[:2]
nimages = zarr.open(pc_ph, mode='r').shape[1]

Grid index (azimuth, range) of the points from their hilbert index:

mc.pc_gix(pc_hix, pc_gix, shape=shape)

The network of image pairs: every image with the next three (it closes loops, e.g. (a, b), (b, c), (a, c)):

image_pairs = mr.TempNet.from_bandwidth(nimages, bandwidth=3).image_pairs
image_pairs.shape

EMCF unwrapping of all interferograms of the network:

mc.emcf_pc(pc_gix, pc_ph, image_pairs, pc_unw, **spacing)

Phase closure correction: where the unwrapped interferograms of a loop of image pairs do not add up to zero, every region of points (connected without long edges) is corrected by whole cycles in the interferograms that disagree with the majority; where that does not close every loop of a point, the interferograms fitting most of its loops are kept. The output ts is the unwrapped phase of every image relative to the first image (ref); the interferogram (a, b) is ts[:, a] - ts[:, b], so every loop closes. misclosure_fraction is, per point, the fraction of interferograms that did not fit the loops before the correction, change_fraction is the share of its interferograms that the correction changed (a point quality for masking); region is -1 for points of small regions, corrected one by one.

mc.unwrap_correct_closure_pc(pc_gix, pc_ph, pc_unw, image_pairs, pc_ts, **spacing, ref=0,
                             misclosure_fraction=pc_misclosure, change_fraction=pc_change, region=pc_region)
mc.pc_pyramid(
    pc_ts,
    './phase_unwrapping/pc_ts_geo_pyramid',
    x = pc_e,
    y = pc_n,
    ras_resolution=20,
)
mc.pc_pyramid(
    pc_misclosure,
    './phase_unwrapping/pc_misclosure_geo_pyramid',
    x = pc_e,
    y = pc_n,
    ras_resolution=20,
)
mc.pc_pyramid(
    pc_change,
    './phase_unwrapping/pc_change_geo_pyramid',
    x = pc_e,
    y = pc_n,
    ras_resolution=20,
)

The wrapped sequential interferograms next to the unwrapped ones, the differences of the image phases ts[:, image] - ts[:, image + 1]; the slider image is shared.

dates = './load_data/meta.toml'
wrapped = mc.view('./pixel_refinement/hix/pc/pc_ph_pyramid', show='intf_seq', dates=dates, label='Wrapped phase')
unwrapped = mc.view('./phase_unwrapping/pc_ts_geo_pyramid', show=lambda v, image: v[..., image] - v[..., image + 1],
                    clim=(-10, 10), label='Unwrapped phase')
wrapped + unwrapped

Where the interferograms disagreed before the correction, and how many interferograms of every point the correction changed:

mc.view('./phase_unwrapping/pc_misclosure_geo_pyramid', clim=(0, 0.5), label='Misclosure fraction') + \
    mc.view('./phase_unwrapping/pc_change_geo_pyramid', clim=(0, 0.2), label='Change fraction')