ON THE FEASIBILITY OF APPLYING ORBITAL CORRECTIONS TO SAOCOM-1 DATA WITH FREE OPEN SOURCE SOFTWARE (FOSS) TO GENERATE DIGITAL SURFACE MODELS: A CASE STUDY IN ARGENTINA

In this work we present an orbital correction workflow developed with FOSS tools to compensate for orbital errors present in Synthetic Aperture Radar (SAR) interferograms. The technique is tested in forested areas in Argentina, using full polarimetric images from the argentinean SAR constellation SAOCOM-1 (Satélite Argentino Con Microondas). The results are contrasted with field measurements of canopy height provided by local producers, and the results show that the Root Mean Square Error (RMSE) of the satellite measurements is significantly reduced after the orbital correction. Moreover, forest plantation become more distinguishable in the retrieved Digital Surface Models, especially in those pairs with larger spatial baseline. A section of this article is also dedicated to the discussion on which are the best parameters to run the module, and how different configurations can affect the result.


INTRODUCTION
SAR Interferometry has been widely used for topographic and surface mapping (Zebker and Goldstein, 1986). In the field of forestry applications, the generation of Digital Surface Models (DSMs) is a key input for the estimation of canopy height, and therefore the initial step to the calculation of forest biomass or forest productivity (Soja et al., 2014). In this sense, the launch of the Argentinean SAR constellation SAOCOM-1 has brought new possibilities for the generation of DSMs, especially thanks to the low frequency of its instrument (L-Band) (Giraldez, 2003), which allows the signal to interact directly with tree branches and stems, and is less influenced by temporal decorrelation.
However, single pair interferometry may present challenges in its processing, especially when there are orbital inaccuracies, which can lead to misestimate the height. In fact, a previous work has shown that SAOCOM-1 interferograms can be improved by applying an orbital correction algorithm when they are affected by these errors (Roa et al., 2021).
In this work we present a Python implementation of the orbital correction method proposed by (Pepe et al., 2011), using SAOCOM-1A data acquired over a forest area in Argentina, in order to make it available to the public in the near future. Moreover, this orbital correction module could be integrated within the processing chain of free interferometric software such as GMTSAR and ISCE. * Corresponding author

Study area
In order to asses the effectiveness of the orbital correction module, a case study is presented in the surroundings of the city of Concordia (Entre Ríos, Argentina), located in one of the most important areas in terms of forestry production in Argentina. The approximate footprint and location of the SAR images utilized for this work is presented in figure 1 Entre Rios is among the most important provinces in terms of forest production in Argentina, only surpassed by Misiones and Corrientes (Arturi et al., 2017). In this sense, the retrieval of forest plantations height, which is directly linked to their productivity, is very relevant both for local producers and government agencies.

Data inputs
For this case study, three SAOCOM-1A images were selected according to the needs for interferometric processing (Table 1), and three interferometric pairs were formed according to Table  2. SAOCOM-1A is part of a constellation of two satellites launched by the argentinean space agency -Comisión Nacional de Actividades Espaciales (CONAE) -in 2018 and 2020 respectively. The SAR instrument mounted on these spaceborne platforms works in L-band (Λ 23cm) and its use in forestry applications is of great interest due to the way in which it interacts with trees, as showed in figure 2.
The selected scenes for this study correspond to the Stripmap imaging mode (around 10 m. spatial resolution) and all of them The red frame is the approximate footprint of the SAR images, the blue-dotted rectangle indicates the areas for which field measurements were provided. Figure 2. Common interaction of SAR L-and C-bands with the forest structure (Omar et al., 2017) are full polarimetric, which means the information is available for the VV, VH, HV and HH polarization channels.
On the other hand, field measurements of tree heights were provided by local producers 1 in order to assess the accuracy of the obtained DSMs, and the open-to-public Shuttle Radar Topography Mission elevation model (https://www2.jpl.nasa. gov/srtm/) was used to refer these measurements to the ellipsoidal height.

Orbital corrections in SAR interferometry
The principle of SAR interferometry (InSAR) technique relies on the measurement of the phase difference between two or more SAR images, acquired from different positions and/or at different times. (Pepe and Calò, 2017). The interferometric phase, related to any point seen within the illumination footprint of the SAR sensor, can be expressed as: where RM and RS are the sensor-target distance for the master and slave acquisitions respectively, λ is the radar wavelength, and B, ϑ and α are the baseline distance, the look angle, and the tilt angle respectively.
Considering the Taylor expansion of equation 1 around the scene center (range position R0 and the incidence angle ϑ0), which is a valid approximation for spaceborne systems, the interferometric phase can be expressed as: 2019-08-20 2020-02-28 320 192 Table 2. Interferograms formed with the available images. Bp stands for perpendicular baseline, and Btemp stands for temporal baseline.
Where B and B ⊥ are the parallel and perpendicular components of the baseline, respectively.
Then, the range phase gradient is given by: Where ΩR is the local terrain slope in range direction. Equation  3 shows that there exists a phase gradient in the range direction which is modulated by the topography. It also shows that the phase gradient can be affected by an error in the perpendicular baseline value.
On the other side, non-parallel orbits may produce baseline variations in azimuth direction also introducing a phase gradient. Expressing equation 2 as a function of azimuth position, both the parallel and perpendicular components of the baseline depend on the azimuth spatial coordinate as: Then the phase gradient in azimuth direction is given by: From equation 5 and following Pepe et al. (2011) we can infer from the first term that an azimuth phase ramp artifact may be generated by errors on the azimuth rate of the parallel baseline. The second term introduces a variability of the azimuth slope with range. Moreover this term is significantly smaller than the first one for satellite geometries, and can be neglected as well the third term.
The parameters B ⊥ and ∂B ∂x significantly affect the interferometric phase as we can see from equations 3 and 5. Moreover these parameters can be estimated using the algorithm proposed by Pepe et al. (2011). In the algorithm the range and azimuth phase gradient obtained by processing the data with the nominal orbital values are compared with the ones computed by using an external DEM.
Then the correction value for the perpendicular (∆B ⊥ ) baseline is given by the following equation: are the estimated and nominal perpendicular baseline values and ∂φm ∂R and ∂φ sint ∂R are the measured and synthetic interferometric phase range gradients. It is clear that should the nominal orbit information be precise and in the absence of deformation, atmospheric and residual topographic contributions, the differential gradient would be zero and the baseline correction would vanish (Pepe et al., 2011). The latter means that the measurement and the synthetic phase range gradients are equal as we can see from equation 6.
In the same way, to compensate for incorrectly estimated nonparallel orbits we can compute the parallel baseline rate correction value (∆ ∂B ∂x ) using equation 5 and considering only the first term as follow: Equation 6 and 7, which relate range and azimuth phase gradients and orbital parameters, are approximations. In a real case, when the exact equations should be used, the the differential range gradient is given by: Where f (·) is a function that maps the perpendicular baseline value into the range phase gradient. Then by using an iterative numerical approach, the problem is reduced to find the value To compute the parallel baseline rate in azimuth a similar approach is used. The differential azimuth gradient can be expressed by the following equation: In this case g(·) represents a function that maps the parallel baseline azimuth rate into the azimuth phase gradient. In the same way, the problem is reduced to find the value of ∆ ∂B ∂x that makes g For sake of simplicity, from now on the terms ∆ ∂B ∂x and ∆ ∂B ∂x are simply referred to as ∆B ⊥ and ∆B respectively, and ∆B will refer to the generic case of the latter two.
As stated before, the algorithm makes use of iterative numeric methods to firstly get a coarse estimation of ∆B ⊥ and ∆B . This is achieved by calculating the P G for two equally distanced ∆B values from 0 (i.e. +10 and -10), and therefore getting two points in the ∆B-P G plane. These two pairs are used to find the initial coarse estimation of ∆B by applying the secant method ( figure 3).
Based on this initial approximation, the P G is re-computed for N values of ∆B around the coarse estimation, creating N pairs of (P G,∆B). Then, the fine estimation of ∆B is obtained by fitting a linear function to these points and finding its root, which makes the P G zero ( figure 4).  The computation of P G is a key step in the correction process. This will be readdressed in section 4, and the processing workflow is detailed in the Appendix.

InSAR processing and orbital correction
The interferometric pairs were formed and processed both with GMTSAR and ISCE up to the interferogram formation. At this point, the outputs from this step was ingested by the orbital correction module in order to compensate for the orbital inaccuracies. In order to mask out the non-coherent points, a coherence threshold (γ) of 0.25 was applied. The result of this rectification was re-inserted in the processing chain to continue with the unwrapping stage and the final conversion of phase to elevation.
To evaluate the impact of this correction, the same process was ran without making use of the orbital module, and the results were compared by contrasting with the field data measurements (see figure 5). As a mean of visually illustrating the effect of the orbital corrections on the interferograms, some examples are presented in figure 6. Both the differential and topographic interferograms are shown, since in the latter the orbital fringes are blurred by the topographic ones and the effect of the correction is not so evident.

Results
Although the focus of this work are the topographic interferograms, the orbital ramps are not so clearly visible when displaying them, and that is why the differential ones are also shown. When looking at the original differential interferogram (A in figure 6), there is a clear orbital fringe pattern, mainly in the range direction. The corrected interferofram in B does not present these phase patterns, so the correction seems to have worked effectively.
The effect of the orbital correction on the final DSM is evident especially in forested areas. If applied, forest plantations are clearly differentiated from their flat surroundings owing to their height. On the other hand, if no orbital correction is applied the forest canopy is not so clearly distinguishable from the rest (see figure 7). These results can be also quantified by computing the Root Mean Square Error (RMSE) of the differences between the estimaed elevation and the actual measured tree height. Table 3 depicts this calculation for the four processed polarimetric channels. Figure 6. Results of the orbital correction for Pair 1 -HH Polarization. A-B are the differential interferograms before and after correction, and C-D are the topographic interferograms before and after correction, respectively (Note: These products are mapped in slant-range coordinates).  Table 3. Summary of Root Mean Square Errors calculated for the four polarimetric channels, in meters. RMSE corresponds to processing workflow B (no orbital correction applied), and RMSE* corresponds to workflow A (orbital correction applied)

DISCUSSION
It is interesting to note that the pairs presented in Table 2 are different not only in terms of their spatial and temporal baselines, but also in their interferometric coherence. This parameter affects directly the phase gradient on which the estimations of ∆B ⊥ and ∆B are based. Given a coherence threshold (γ) used to mask out non-coherent pixels, in a noisier interferogram a smaller number of pixels will we available to compute the phase gradient, compared to a more coherent one (see figures 8, 9 and 10). This may influence the selection of the pairs of (∆B, P G) that define the coarse and refined estimations of the orbital correction, as discussed in 2.3. This is exemplified in figure 11, where HH and HV data yield different estimations of fine ∆B ⊥ . Figures 11 and 12 show how the estimations for ∆B ⊥ converge to a more similar result when changing the γ threshold. The effect of this convergence is not neutral: In the first case (γ = 0.25) the module does not achieve a proper correction for the HV channel (figure 13), and when changing γ to 0.10 the correction on this channel seems to improve (figure 14). Figure 8. A subset of the coherence maps for pair #2. HV is much noiser than HH and this has an impact on the orbital correction Figure 9. Histogram of coherence for the subsets of HH and HV of figure 8, after applying a γ threshold of 0.25. Figure 10. Histogram of coherence for the subsets of HH and HV of figure 8, after applying a γ threshold of 0.1. The number of points to be processed in both channels is now more equal. Figure 11. Fine estimation of ∆B ⊥ from HH and HV data in pair #2 The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/W2-2021 FOSS4G 2021 -Academic Track, 27 September-2 October 2021, Buenos Aires, Argentina Figure 12. Fine estimation of ∆B ⊥ from HH and HV data in pair #2 Figure 13. A: Original HV differential interferogram, B: Corrected HV interferogram after using the ∆B ⊥ shown in figure 11. The orbital fringes are not properly corrected; indeed, they get worse. Figure 14. A: Original HV differential interferogram, B: Corrected HV interferogram after using the ∆B ⊥ shown in figure 12.
The correction in this case seems to have improved significantly.
In this sense, it is important to provide the users of this module with flexibility in terms of the parameters to run the correction. Apart from the γ threshold used, the output of the correction also seems to be sensitive to the range and azimuth sampling factors applied to the input interferograms. Different combinations of these parameters lead to diverse results depending on the interferogram used as input.
However, what is the best solution? Which of the different combinations of ∆B ⊥ and ∆B that give (visually) good results is the one that will lead to the best correction? It is very difficult to give an answer that will satisfy all the cases, since the input interferometric data will vary in every case. As an example, the RMSE for the HH and HV channel of pair #2 was re-computed after processing the data with the coherence threshold and multilooking factors used in this section.  Table 4. RMSE obtained after running the orbital correction with different combinatios of Multilook and γ thresholds. In case of N/A, this means the correction could not be achieved (similar to the case of figure 13).
These figures and partial results suggest that the fact of working with a higher coherence threshold, and having a large number of points (when working at higher resolution) guarantee better results, but this is subject to further discussion and cannot be established as a rule of thumb. Figure 15 plots the resulting RMSE after running the correction with different coherence thresholds and the same multilook factor, and it is clear that using very high γ values may not necessarily give better outputs.

CONCLUSIONS
The results presented in this article show that the orbital correction algorithm proposed by Pepe et al. (2011) can be implemented in Python and integrated with the processing chain of FOSS interferometric software. We have been able to show that the compensation for the orbital ramps in the interferograms can significantly reduce the height error in forested areas, where the interaction of SAOCOM L-Band with volumetric targets such as trees is fully exploited.
However, it is still necessary to fine tune the algorithm in order to define the proper parameters to use with each data set.
The tests were carried out on a study area with field measurements, but where the availability of SAOCOM imagery was not the ideal in terms of temporal separability and variability in spatial baselines. It would be desirable to perform this correction on larger datasets of SAOCOM scenes, being able to experiment with different interferometric configurations and shorter temporal baselines to prevent temporal decorrelation. While these issues remain unsolved, the user of this module must be very conscious of the effect of coherence and multilooking in the estimation of the phase gradient, a key step in the correction process.
A good approach would be to compute the ∆B ⊥ and ∆B for the most coherent polarimetric channel within the same pair, and apply them to the rest of the channels, avoiding the divergence problems exposed in the previous section. For forested areas, the more coherent channels are typically the VV and HH ones (Neumann et al., 2010).
Regarding the potentiality of SAOCOM images to map forest heights, a second stage of this work should include the integration of this Orbital Correction Module with processing techniques that make use of full polarimetric data to optimize the interferometric phase and get better results of the height, such as Polarimetric SAR Interferomtry (PolInSAR). These kind of approches require short temporal baselines to mitigate decorrelation, so the main challenge for future publications will be to get access to SAOCOM scenes that satisfy these constrains. APPENDIX Figure 16 shows the general workflow of the orbital corrector module described and implemented in this work.  Figure 17 shows more in detail the steps for the coarse and fine estimations presented in the general workflow. As explained in section 2.3 the algorithm generates a number of ∆B values and uses them to introduce a perturbation in the orbits of the images, re-compute the interferogram and generate a P G map for each iteration (2 iterations for the coarse estimation, and 6 for the fine one). Figure 17. Detailed workflow to describe the coarse and fine step present in the general workflow ( Figure 16).
In each step, the orbit perturbation takes as argument both ∆B ⊥ and ∆B . It is important to note that when estimating ∆B ⊥ the initial values for both parameters are zero. Instead, when estimating ∆B , the fine estimation of ∆B ⊥ is passed as a fixed parameter, as illustrated by figure 16.