[Suggestions]: Some suggestions on the Colab results

[Xiaohua-Eric-Xu]

I was trying to make some suggestions in another thread (a while ago) but maybe this is a better place for such discussions.

1. In your Colab page, things look fine before doing SBAS. Then after SBAS, things start to look strange. e.g., below the time-series has a spatial jump that does not exists in the input data. I am not sure how it happened but it does not look correct as time-series analysis should not produce something that does not exist in the input data.

2. It is mentioned in the Colab that "Remove trend and apply gaussian filter to fix ionospheric effects and solid Earh's tides. " I think I mentioned in prior discussions that this is a dangerous thing to do but it seems that you didn't buy it. Here are some stuffs to consider. a) C-band is not much affected by ionosphere compared to L-band data. b) Deformation contains components from all wavelengths, a simple bandpass may get rid of deformation as well. c) effects of solid earth tides are prominent only when you start to look at large regions and could be easily removed with simple calculation rather than detrending/high-pass filtering. d) atmospheric noise also has components over all wavelengths, and the part below (atm noise produced by gravity waves) remains after high-pass filtering.

3. Personally I love to see advanced users making their own processing chains involving their own flavors, but if you would like to constantly populate your results under the GMTSAR repository, please make sure your results are done in a correct way (that's good for both you and us). Otherwise it may mislead users (especially beginners) and I would recommend that you keep it private. If in any case you don't believe my suggestions, fell free to verify the stuffs above with any other InSAR scientists.

[AlexeyPechnikov]

Thanks for the suggestions! As you remember, I’ve already asked you about the Sentinel-1 SBAS reference example in GMTSAR tickets.

1. Here is 2piN phase jump between the unwrapped phases. Actually, I’ve built GMTSAR Docker image https://hub.docker.com/r/mobigroup/gmtsar to check the case and reproduce the results. See the GMTSAR produced displacement maps ranges below:

As the results, the output map looks as

But you share it as a reference solution, so I provide the exact GMTSAR results in the notebook using GMTSAR sbas tool and add the enhanced ones using detrending and correlation weighted least-squares calculation.

2. For a square area in geographic coordinates ~50km size, the Gaussian filter sigma 12km is 1/4 of the image size. What information lost do you see here? For a real-world data, we use band-pass filtering to remove noise. What is a physical nature of a pixel-scale changes? That’s a noise because the sample rate is not enough (Nyquist’s theorem). What is a physical nature of a raster-size scale changes? That’s a noise because the sample size is not enough (three-sigma rule, but often 4 sigma required for a discrete signals processing). There is no way to measure an unlimited spatial spectrum and Gaussian filtering does the job extracting the valid part of a measured spectrum. There is no way to work here on spatial components like dozens of kilometers (ionospheric effects) and hundreds of kilometers (solid earth tides).

and the part below (atm noise produced by gravity waves) remains 

Yes, because it has the scale in the valid spectrum range. We cannot remove it by this way. Weighted least-squares solution on a per pixel time series can do the job.

a) We still have ionospheric delay and it can be removed using Gaussian filtering with a larger sigma. In case when we already use Gaussian filtering, we don’t need to care about the delay.

b) Gaussian filter does not crop spatial components, but only modulates the amplitudes. While the output displacements can be changed (actually, not) we save all the deformations spatial components. It should be a very hard Gaussian filtering (like to 1km in this case) to affect the SBAS results.

c) Solid earth tides removed using Gaussian filtering even with much larger sigma. It means we son’t need to care about them when we already use Gaussian filter.

[AlexeyPechnikov]

@Xiaohua-Eric-Xu You always or almost always recommend use reference point which approach is unstable by design and produce unpredictable SBAS results. Select different point and the SBAS results are useless. There is no way to fix it in the spatial domain, but that’s easy in (spatial) frequency domain using reference frequency instead of the reference point. We can think about the low-pass filter cut-off frequency as of the reference frequency. This approach is stable and predictable by design. Power spectrum of natural processes is limited and so power decreasing (fast). For an example, topography is fractal and the power vs wavelength decreasing linearly in double logarithmic coordinates. It’s fast enough for all the numerical calculations stability. The power spectrum values mean spatial components dispersion for a matrix calculations. Calculate band-pass Gaussian filter for wavelengths N-1 and N+1 and the raster dispersion is the power spectrum value for the wavelength N. For our case it means that doubling or halving the reference frequency (wavelength) produces the similar SBAS results with the same linear and seasonal trends. That’s following to the theory and easy to check on practice (try different filter wavelength in the example notebook above). The notebook shows the technique and allows to change the wavelength to check the SBAS results stability. Even the more, the reference point cannot be selected automatically, and the reference frequency is defined using easy arithmetics using 2-3-4 sigma rule (1/4 of the image size for 2 sigma and 1/8 for 4 sigma).

If you don't interested in the theory see the practice. Just for the example see https://www.planetmechanic.net/code

InSAR data provides extremely high resolution velocity measurements, but is subject to significant long-wavelength errors. GPS data provides high accuracy at long wavelengths, but is very low resolution. This method combines the two via a simple interpolation and filtering method, producing one dataset controlled by GPS at long wavelengths, while retaining unaltered InSAR coverage of short-wavelength features.

And again, in GMTSAR tickets you always or almost always recommend use GPS data. It means replacement of the long wavelength power spectrum of the displacement rasters and so it's equal the gaussian filtering in the notebook. Sure, the reference point idea is a far outside of this way. Actually, I automate and explain that you already use.

[Xiaohua-Eric-Xu]

1. 

This is a spatial jump. I did notice a while ago there being an error in README_sbas.txt earlier but did not know whether the fix gets updated to topex download site. i.e., scene.tab should be in a chronological order. None the less, even if I used the prior scene.tab, there should not be any spatial jump. I am not sure what unwrapping parameters you used but there should not be a spatial jump. Below is the same date in the time-series (with the wrong scene.tab) I just got from the example.

Gaussian filtering is a low pass filter that may take out all components below or above a certain wavelength. But it does not mean that anything above/below that wavelength stays intact. They are just less affected with a smaller gain in the filter. Also, result being predictable does not necessarily mean it is the right result. In fact I am not sure why end to end testing does not grant you predictability (we do that often to make sure the software's on track). If you mean instability from unwrapping, then that ambiguity (N2pi) is by nature/definition and referencing to a pixel or solving for the ambiguities should fix that. If you enforce "stability/predictability" with removing parts of data, then real errors, bugs or wrong methods may just get covered up. One good example is the "Baseline Re-estimation". Before everyone seem to believe that the long wavelength ramp comes from orbital errors, and would estimate a baseline associated error with that ramp, update the orbit records and then the ramp's gone. But we all know that if you do so, it actually takes out not just orbit error, but also tropospheric phase, ionospheric phase, tectonic deformation, even oscillator drift (in ENVISAT), etc. Keeping doing this will cover up all these "signals" and prevent from figuring out how to precisely estimate each component of radar phase. I don't think folks do that estimation any more.

You always or almost always recommend use reference point which approach is unstable by design and produce unpredictable SBAS results. Select different point and the SBAS results are useless. There is no way to fix it in the spatial domain, but that’s easy in (spatial) frequency domain using reference frequency instead of the reference point.

If you pick a ref point/area, then deformation is considered relative to that point/area. One main reason that one would not detrend or filter the results is that deformation could be of longer wavelength than your data coverage. I don't agree that you could not recover any signal beyond the measuring bounds, e.g., https://www.nature.com/articles/nature04797. Tectonic motion typically span hundreds to thousands km, and you can definitely recover the right signal of it even with just simply referencing to an area, with a measurement that spans only a 100-200km. If you filter the data, there's no way you could get that in the data.

It means replacement of the long wavelength power spectrum of the displacement rasters and so it's equal the gaussian filtering in the notebook.

Referencing to one point/area is no way equal to gaussian filtering. Indeed there are approaches to remove and restore with gaussian filtering, but you need to have the solution (with many GNSS, or in other words, referencing to a group of points) at long wavelength first.

In fact, folks in the field are already starting to pursue absolute measurement with InSAR, considering we are already directly seeing plate motion both in range and azimuth measures.

[AlexeyPechnikov]

1. I rechecked the processing, and the result was correct. Let's remove unreproducible manual manipulations in the example bash scripts using 3rd party artifacts ("in order to correct for Elevation Antenna Pattern Change") to produce the exact results for the included Sentinel-1 scenes. The linear ramp is expected and GMTSAR SBAS with smooth=1 cannot remove it when PyGMTSAR detrending removes it completely. Actually, I don't understand why do you include the problematic scenes from 2015 while there are fixed in 2016 (when the example created) especially when GMTSAR cannot fix the ramp itself. Anyway, the notebook does all the processing by correct and reproducible way and it works fine for the updated scenes or any other scenes.
2. GMTSAR SBAS atmospheric correction algorithm is an iterative timeseries smothing to minimize the solution dispersion and it's emulate frequency filtering. It means Gaussian filtering with sigma >> timerange does the same job much faster and effective and without unpredictable signal modifications. Did you check the Timofeeva algorithm output signal spectrum?.. It's strange to say about "problem" of Gaussian filtering when you use tension sufraces to damage topography spectrum completely and the iterative timeseries changes to damage the spectrum completely inatead of the well-known frequency domain operations.
3. Using the reference frequency we can separate long-wavelength reference surface and correct it using GPS or other data or use it as is. Removing completely the reference surface we have the correct solution for the task which cannot be solved by the reference point approach. It works great even for absolute measurements.

[AlexeyPechnikov]

If you pick a ref point/area, then deformation is considered relative to that point/area. One main reason that one would not detrend or filter the results is that deformation could be of longer wavelength than your data coverage. I don't agree that you could not recover any signal beyond the measuring bounds, e.g., https://www.nature.com/articles/nature04797. Tectonic motion typically span hundreds to thousands km, and you can definitely recover the right signal of it even with just simply referencing to an area, with a measurement that spans only a 100-200km. If you filter the data, there's no way you could get that in the data.

Let me show you an example. There are two measurements of a harmonic signal: 0 in the 1st point and 1 in the second. Please use any method to detect the wavelength! Obviously, that's impossible. That's all about long wavelength estimation for a small area. We can detect the long-wavelength part of the spatial spectrum using low-pass filtering but it's impossible to detect the exact frequencies. There is no information lost when we use gaussian filtering for it because there is no information about the wavelengths.

The linked paper estimates just a linear trends for a separate points ("slip rate") without any analysis of the spatial frequencies and their amplitudes:

But what shows the picture? Nothing I think. Here would be underground water or magma reservoirs for the areas with higher LOS velocity variations but that's impossible to detect when seasonal and yearly trends missed. Where are the reservoirs, when and how these are filling, how deep buried these are, how the variation is correlated to seismic events... here are no answers but only picture of the surface pixels movements which known a priori. And the 10 years timeseries results are applied to 200-300 years movements... that does not sense of course and 300 years ago the movements could be very different. "It is reasonable to assume" instead of the science analysis, really?

When you select a single-pixel reference point it's unstable because outliers in the rasters stack crash all the results. And when you use an area instead of the single pixel then the area scale is the reference wavelength but there are some border effects which needs to be removed. Using the reference area and gaussian filtering to exclude the border effects we see exactly the same approach as using GPS fitting for long wavelength. Matrix calculation of Gaussian filter in floating window provides better results than just classical Fourier transform of the entire image to compute the power spectrum components but you are free to use any spectrum calculation as you prefer in case you aware about border effects and so on. Anyway, spectral domain operations are robust and noise tolerant and that is known from the theory while spatial (or time) domain operations are not.

And let me repeat again: GMTSAR scripts provide as science tools (like to GPS-related long wavelengths spatial spectrum continuation) as art tools (like to reference point based approach which can be used effectively by experienced geologist but it is intuition based and not reproducible and the results accuracy is not predictable). I'm interested to use and ehnance the 1st and totally exclude the 2nd. I work to provide reproducible pipeline which works for any scenes and any areas without hidden artifacts and scripts and weird ways of parameters selections. All the parameters (like the reference wavelength) should be calculated using exact formulas and the numerical computations are stable and the spectrums are valid.

[AlexeyPechnikov]

With disabled 3rd party artefacts and hidden scripts exactly the same jump produced using GMTSAR processing:

# download and unpack the example
wget -c http://topex.ucsd.edu/gmtsar/tar/S1A_Stack_CPGF_T173.tar.gz
tar xvzf S1A_Stack_CPGF_T173.tar.gz -C .
# clean 3rd party artefacts for pre-processing the Sentinel-1 SLC scenes
find raw_orig -name '2015*_manifest.safe' -exec cp /dev/null '{}' \;
# run the example processing scripts
./README_prep.txt
./README_proc.txt
./README_sbas.txt

I've added the issue description on the GMTSAR DockerHub page https://hub.docker.com/repository/docker/mobigroup/gmtsar