Thermal Denoising of Cross-Polarized Sentinel-1 Data in Interferometric and Extra Wide Swath Modes

An updated algorithm for removal of thermal noise in Sentinel-1 synthetic aperture radar (SAR) Level-1 Ground Range Detected (GRD) data in cross-polarization is presented. The algorithm is comprised of two steps: correction of the annotated thermal noise magnitude (previously proposed in park2018) and a novel correction of the annotated thermal noise range dependence. The magnitude of the annotated thermal noise is corrected by applying scale and offset coefficients tuned on a few hundred of Sentinel-1 data acquired over surfaces with low backscatter in Interferometric Wide and Extra Wide swath modes in HV and VH polarizations. The values of coefficients for all modes and polarizations for data processed with Instrument Processing Facility (IPF) version 3.1 -3.3 are provided. The range dependence is corrected by minimizing a cost function between the annotated range profiles of thermal noise and antenna pattern gain (APG). An objective validation metric based on comparison of averaged backscatter at inter-swath boundaries is proposed. Validation is performed on hundreds of Sentinel-1 scenes acquired over open ocean, doldrums, deserts and sea ice. It shows that the new algorithm outperforms the standard thermal noise removal algorithm proposed by European Space Agency in almost all cases. Analysis shows that the new algorithm worsens noise correction in cases when the range dependence of the annotated APG does not match with the observed signal, indicating either problems with signal processing on IPF or imprecise annotation of APG.

The European Space Agency (ESA) Sentinel-1 mission is 41 equipped with a C-band SAR instrument which can operate in 42 Terrain Observation with Progressive Scans SAR (TOPSAR) 43 mode and provide Extra Wide (EW) or Interferometric Wide 44 (IW) observations useful for sea ice, ocean and land moni-45 toring. The data in Level-1 Ground Range Detected (GRD) 46 products consist of focused SAR data that has been detected, 47 multi-looked and projected to ground range using the Earth 48 ellipsoid model, and include both the SAR signal intensity 49 and a look-up table (LUT) with values of noise power varying 50 in the range (here and throughout the manuscript under "the 51 range" we mean ground range unless slant range is specified) 52 and azimuth direction. The Sentinel-1 documentation [2] rec-53 ommends applying thermal noise correction by subtraction of 54 noise equivalent sigma zero (NESZ) from radar backscatter: where σ 0 SN is the denoised radar backscatter, and σ 0 is the 56 calibrated signal power provided in the S1 GRD file (DN ) 57 [3], [4]: where A i is is a calibration constant for each sub-swath i. 59 The choice of calibrated sigma zero (instead of e.g., beta zero) 60 is dictated by scientific applications of Sentinel-1 SAR data 61 requiring calibrated ground range data. 62 NESZ (σ 0 N ) is a calibrated product of the range (N A Ri ) and 63 azimuth (N A Ai ) components of the noise power provided in the 64 look-up table: In several studies [1], [5], [6], [7] it was shown that the 66 magnitude of the provided noise values is not sufficiently 67 precise because the noise dynamically depends on factors that 68 are not included in the noise vectors, and even after noise 69 subtraction the artefacts in the signal are still observed. A 70 method for correction of the NESZ values by applying scaling 71 and offset coefficient was proposed in [1]: where K ns and K pb are subswath dependent noise scaling 73 and power balancing coefficients. Initially this method was 74 developed for Level-1 GRD data from Sentinel-1 A and B 75 developed criterion on their homogeneity for the thermal noise 85 as well as residual noise (i.e., multiplicative noise) removal [8], 86 [9]. These adaptive algorithms, however, utilize radar signal as spreading loss (RSL) gains, but in many cases it not so.

96
According to [10], the annotated slant range noise power is 97 computed as: 98 N A R (n s ; R, η) = σ 2 n (n s , η) · G 2 tot (n s ; R, η) where G 2 tot (n s ; R, η) is total gain computed as: 99 G tot (n s ; R, η) = G ds (η)G pg (n s ; η) G eap (n s ; R)G rsl (R) k noise k proc (6) where R is slant range, η is slow time, n s is sub-swath 100 number, G ds (η) is de-scalloping gain, G pg (n s ; η) is PG-101 correction factor, G eap (n s ; R) is two-way elevation antenna 102 pattern gain, G rsl (R) is the range spreading loss, k noise is 103 the noise calibration factor, and k proc is the processor scaling 104 factor.

105
The only slant range-varying variables are G 2 eap (n s ; R) and 106 G 2 rsl (R) and, therefore, the range component of the thermal 107 noise should be proportional to a function that can be called 108 "antenna range gain", G ar : However, it was observed that quite frequently there is an 110 obvious shift between N A R (n s ; R) and G ar in the range direc-111 tion as shown on Fig.1(a). Moreover, the same shift is observed thermal noise reduction. Instead of adaptive correction using 124 radar signal, the new algorithm for thermal noise removal 125 utilizes extra information from the annotated metadata (i.e. 126 G ar ). Thus it combines both the range shift and the magnitude 127 correction of thermal noise vectors.

128
In addition the algorithm was extended for application on 129 Sentinel-1 A/B data acquired in EW and IW mode, in HV and 130 VH polarization and IPF version 3.1, 3.2 and 3.3. We also sug-131 gest a new quantitative metric for comprehensive evaluation of 132 the new algorithm performance. The improvements in thermal 133 noise removal in several surface types including sea ice, ocean 134 and deserts are illustrated.

135
The paper is organized as follows: first we characterize 136 the data used for training and validation, then we present 137 the algorithms for the thermal range noise shift / magnitude 138 correction and the new validation metric, and finally we 139 provide the coefficients that can be used for implementation 140 of our algorithm of NESZ modification and show and discuss 141 the validation results.  Table I. independent Sentinel-1 data were used as specified in Table II.  The overall scheme of the developed algorithm is presented   Next, a vector of range spread loss gain is computed:

150
where Rref and T sr are vectors of reference range and slant 202 range time interpolated from the GGADS on the destination 203 coordinates, and c is the speed of light (see Eq. 9-20 in [10]). 204 Next, a one-dimensional linear spline interpolator using 205 'interp1d' class from the SciPy package (for details see 206 http://docs.scipy.org) for the annotated range noise is trained 207 on vectors of annotated range noise (N A R ) and destination 208 column coordinates (x) from the NADS LUT: The shift of the annotated range noise vector (∆X) is found 210 by minimizing the cost function: where (G ar ) is the vector of "antenna range gain" computed 212 from vectors G eap and G rsl using Eq. 7, x denotes column 213 coordinates, and × is the Pearson's correlation coefficient 214 between G ar and range noise values interpolated on the 215 column coordinates offset by ∆X.

216
Finally, the shifted noise is computed by applying the 217 trained interpolator on the pixel coordinates with the offset: As illustrated on Fig.1(a) the shifted noise matches better 219 with G ar . The signal corrected with shifted NESZ appears to 220 be less varying with range, Fig.1 Interpolate θ elevation and θ roll from GGADS on the destination coordinates x; Compute vector of θ boresight on the destination coordinates; Interpolate G EAP from ACDS on the vector of θ boresight ; Interpolate Rref and T sr from GGADS on the destination coordinates; Compute G rsl from RR and T sr ; Compute G ar from G EAP and G rsl ; Train the interpolator I N using N A R and destination coordinates x; Minimize N and find ∆X; Compute N A Rs ; end end Algorithm 1: Algorithm for shifting annotated thermal noise vectors in the range direction. S denotes sub-swaths, Y denotes destination row coordinates and other notations are the same as in Section III-B.

C. Finding the noise scaling and power balancing coefficients 222
The training of the noise scaling (K ns ) and power balancing 223 (K pb ) coefficients is presented in detail in [1] and is briefly 224 explained here for clarity. The training is performed only 225 once but on a sufficiently large sample of GRD Sentinel-226 1 A and B images collected in two acquisition modes (EW 227 and IW) and in two polarizations (HV and VH) (see Table   228 I). For each combination of the mode and polarization, lists 229 of images from S1A and S1B satellites acquired over areas The individual coefficients k ns i are found under the as-239 sumption that after thermal noise correction and in case of 240 low signal, the signal can be approximated with a linear 241 dependence on incidence angle: σ 0 SN ∼ α i +β i ·x, where α and 242 β are linear regression coefficients, x is a column coordinate. 243 Noise scaling coefficients are found by minimizing the 244 following cost function, which is defined as a difference 245 between a linear polynomial and denoised signal, computed 246 as a difference between the total signal and scaled thermal 247 noise: The minimization of yields not only k ns i but also α and β 249 that are not used further in the computation of K ns . Then an 250 average scaling coefficient is calculated for each sub-swath: where j is sub-swath index, and N is number of sub-blocks 252 from training images.

253
Power balancing aims to remove the discontinuities of 254 signal intensity near sub-swath margins. The discontinuities 255 are a signature of the noise power differences between sub-256 swaths. The training also comprises two major steps. First, 257 values of power balancing coefficients and auxiliary statistics 258 from individual sub-blocks from Sentinel-1 GRD images are 259 computed by minimizing the cost function defined in Eq. 12 260 where k ns is fixed: Individual power balancing coefficients are then computed 262 by finding the difference in signal power approximated at the 263 boundary between two sub-swaths by linear regressions for 264 these two sub-swaths: Finally, the mean power balancing coefficient for all sub-266 swaths is computed: patches in the neighbouring sub-swaths (see Fig.3). RQM is 276 closely related to the Fisher criterion score [11], [12] that 277 measures the discriminability between two samples in one 278 dimension and is given by: where σ 0 M denotes the denoised signal, j is sub-swath index,   The range quality metric (RQM, see Section III-D) was 306 computed for the images listed in Table II. The average 307 RQM for each combination of mode and polarization for the 308 ESA and NERSC methods and the average of the difference 309 between RQM ESA and RQM N ERSC were computed. Fig.4 Fig. 4: Range quality metric (RQM) values for the different imaging mode from the validation dataset of S1-A (upper plot) and S1-B (lower plot) data.
Both these tables and figures illustrate that in most cases 314 the NERSC algorithm outperforms the default denoising pro-315 cedure in almost every combination of mode and polarization. 316 The only mode where the average RQM difference is slightly 317 negative is S1A/IW/HV and mostly in the Desert validation 318  The reason for worse behaviour of the NERSC algorithm  Comparison of the annotated noise with the antenna range 359 gain (G ar , see Eq. 7) before and after the noise shift is applied 360 (grey crosses and black dots on Figures 7, A and D) shows that 361 before the correction, there is a tangible difference between 362 the annotated noise and ARG, and after the correction the 363 difference is negligible. At the same time, if we compare the 364 uncalibrated signal (DN ) with the G ar (see Figures 7, B and 365 E) we can see very little difference for the first image and 366 quite a significant scatter for the second image. The mismatch 367 between G ar and DN on the second image indicates that 368 either the DN were produced at IPF with incorrect antenna 369 pattern gain, or the annotated APG is incorrect. As a result, 370 when we calibrate noise and signal and compare them (see 371 Figures 7, C and F, crosses and dots) in the first case the shift 372 correction improves the NESZ profiles but in the second case 373 the NESZ profiles are deteriorated.

374
In the second case the NESZ profile becomes shifted 375 incorrectly and the noise is overestimated in the near range 376 and underestimated in the far range, leading to the increased 377 discontinuity in the denoised σ 0 on Fig.5(b). The noise scaling 378 procedure is then applied to the incorrectly shifted NESZ vec-379 tor and is not capable of adequately reducing the discontinuity. 380 Another suggestion on the cause of the shift found between 381 range-dependent corrections and noise vectors is the interpola-382 tion of the range noise vector to project it from slant range to 383 ground range geometry. In particular, possibly, the IPF chooses 384 the interpolation starting time incorrectly, taking into account 385 the 'invalid samples' occurring at the beginning of each sub-386 swath. From the end-user perspective, this inconsistency is 387 heavily data-dependent and burst-dependent. It affects both 388 GRDM and GRDH data in all polarizations. Further analysis 389 of mismatch between range spread loss, signal and annotated 390 noise vectors can in future be performed by comparing the 391 radiometric retrievals from GRD products to those from SLC 392 products from the same image acquisition, but goes out of 393 scope of the present publication. The identified problems are 394 accounted for in development of the future IPF versions.  Fig.8 on 405 a Sentinel-1B GRD image acquired in IW mode in VH 406 polarization over the Laptev Sea. The figure shows clear sub-407 swath discontinuities for the ESA algorithm. Adding noise 408 shift and magnitude correction in the new algorithm clearly 409 helps in reducing the discontinuities, Fig.8(b).

410
Another example is for a Sentinel-1A image acquired in EW 411 mode over sea ice in the Weddell Sea as shown in Fig.9, 412 and is a good illustration demonstrating the difference in 413 the thermal noise appearance in the presence of a mixture 414  Therefore, the visual comparison also confirms that the 426 modified NERSC algorithm still has potential for advanced 427 noise removal for IPF versions later than 3.1.  [1] is not possible due to the absence of publicly 431 available data processed with both IPF 2.9 and IPF 3.1. 432 Nevertheless, we evaluated the performance of both algorithms 433 and compared them with the ESA algorithm by processing 434 several Sentinel-1 scenes taken over relatively calm water 435 in the Greenland Sea in June 2019 (IPF 2.9, old algorithm) 436 and June 2021 (IPF 3.2, new algorithm) as shown in Fig.10. 437 The scenes (as shown on Fig.10(a, d)) were acquired from 438 ascending and descending orbits within one week and the 439 time difference between acquisitions did not exceed 2 days. 440 The choice of time and area is dictated by availability of 441 archived images. Profiles of denoised σ 0 were extracted from 442 the processed scenes by a moving average window 100 × 443 100 pix. Understandably, the profiles contain σ 0 values from 444 the same geographical positions but from different incidence 445 angles.

446
Comparison of the profiles clearly shows that both the 447 old and the new versions of the NERSC algorithm enhance 448 the thermal noise removal procedure by reducing the spread 449 of the profiles and improving their shape. For IPF 2.9, see 450 Fig.10(b) and Fig.10(c), the old NERSC algorithm removes 451 clearly wrong high values of σ 0 in the range of 8 -13 • E 452 from profiles 0, 1 and 2. Other NERSC profiles (e.g., 3, 4 453 and 5) show less longitudinal variations (corresponding to the 454 range variations) than the ESA ones. Overall the spread of 455 NERSC profiles in Fig.10(c) is less than in Fig.10(b). 456 The ESA algorithm on data from IPF 3.2, see Fig.10(e) has 457 a clearly better performance than IPF 2.9 Fig.10(b). Neverthe-458 less, a few profiles (e.g. 1 and 5) show strong artificial range 459 in the denoising software [13].

472
The range dependence is corrected by identifying and re- of both signal and noise powers. Ideally there should be no 480 mismatch and noise profiles should be strictly proportional to 481 G ar . In practice, however, a mismatch of 10 -100 pixels is 482 observed in most of the scenes leading to a situation where the 483 noise profile is overestimated on one side of a sub-swath and 484 underestimated on the other side. That, in turn, leads to over-485 or under-correction of thermal noise at inter-swath boundaries 486 and creates discontinuities and obvious stripes on Sentinel-1 487 images after thermal noise correction. The correction of the 488 noise range dependence is very effective in most of the cases 489 in removing these discontinuities. In a few cases, however, the 490 annotated G ar is obviously incorrect -its range dependence 491 does not match with the signal power. In such cases the quality 492 of annotated noise, fitted to the incorrect G ar , deteriorates and 493 the combined noise removal algorithm provides worse results 494 than the standard one.

495
While the magnitude correction is based on the scale/offset 496 coefficients tuned on many images, the range dependence 497 correction is performed individually for each SAR scene. Thus 498 the combined algorithm addresses the problem of universal vs. 499 individual correction raised in [7], but instead of an empirical 500 Python code [13] which provides a convenient Application 519 Programming Interface for performing thermal noise removal 520 from Sentinel-1 data, training the noise scaling and power 521 balancing coefficients, and validation. The repository hosting 522 the code also provides documentation how to install and use 523 the developed software.