AlexeyPechnikov / pygmtsar

PyGMTSAR (Python InSAR): Powerful and Accessible Satellite Interferometry
http://insar.dev/
BSD 3-Clause "New" or "Revised" License
419 stars 91 forks source link

[Help]: PS candidates selection #159

Open jie666-6 opened 1 month ago

jie666-6 commented 1 month ago

In pyGMTSAR, it seems that the PS candidates are selected based on the Amplitude Dispersion Index (ADI), as written in Stack_ps.py: psfunction = (ps.average/(2*ps.deviation))

However, in the paper "Permanent Scatterers in SAR Interferometry," the ADI is calculated as the ratio of the standard deviation to the average, which is: ADI = (ps.deviation/ps.average).

Generally, ps.deviation/ps.average<0.25 can be regarded as PS pointes.

I am wondering if you can provide any references for the calculation ps.average / (2 * ps.deviation), as I would like to check the reference and understand what threshold should be used for PSC selection.

So I am wondering if there is any reference you can provide for ps.average/(2*ps.deviation) calculation since I would like to check the reference and to know which kind of threshold should I select for PSC selection.

Additionally, the paper mentions that before conducting statistical analysis of the amplitude values, images must be radiometrically corrected to ensure comparability. This implies that radiometric correction is necessary before PSC selection, correct? In pyGMTSAR, there is temporal normalization instead of radiometric correction before calculating psfunction. Am I right, or did I miss something?

    # normalize image amplitudes (intensities)
    tqdm_dask(mean := dask.persist(data.mean(dim=['y','x'])), desc='Intensity Normalization')
    # dask.persist returns tuple
    norm = mean[0].mean(dim='date') / mean[0]
    # compute average and std.dev.
    stats = (norm * data).pipe(lambda x: (x.mean(dim='date'), x.std(dim='date')))

Thanks a lot and I am looking forward to your reply.

AlexeyPechnikov commented 1 month ago

A larger persistent scatterers function value indicates a higher probability, which can be directly used to weight pixels for enhanced interferogram calculation without applying any threshold. Anyway, the psfunction and ADI are essentially inverses of each other. Initially, I used 1/ADI, but later I found the definition of the PSFunction and renamed the function accordingly.

You can find more details in the paper here: Fault creep along the southern San Andreas from interferometric synthetic aperture radar, permanent scatterers, and stacking and refer to the corresponding code in the GMTSAR GitHub repository.

Although corrections can be applied, my tests found them unnecessary. Corrections are defined on a low-resolution grid (such as a 40x40 original Sentinel-1 pixel) and do not affect neighboring pixels, making them less significant for this approach. But don't hesitate to experiment!