Dual-Frequency Airborne SAR for Large Scale Mapping of Tidal Flats

: Digital elevation models of tidal ﬂats are a most valuable data source for the water management of coastal areas and need frequent updates to account for changes in sedimentation, erosion and identiﬁcation of damages in building infrastructure. This paper presents the conceptual design, the processing methodology and ﬁrst results of an airborne SAR campaign conducted in July 2019 at the German North Sea coast, showing the potential for accurate monitoring of height changes at decimeter level in mudﬂat areas, as well as indication of vegetation cover and water ﬂooded areas.


Introduction
The water management of coastal areas influenced by tides requires regular high-resolution and accurate digital elevation models (DEMs). Since large areas need to be surveyed, remote sensors are usually employed. Due to their very dynamic behavior, only extremely short time windows corresponding to +/− 1 h around the low tides are available for the remote data acquisition over areas of tidal flats. Hence, airborne sensors are more attractive than spaceborne ones due to their flexibility in terms of acquisition time. Moreover, high resolution airborne SAR systems-like the DLR's F-SAR-have a wider footprint and are less dependent on weather conditions than conventional airborne laser scanners (ALS), which are usually limited to swath widths of <500 m.
The use of SAR and, in particular, of cross-track interferometric SAR (InSAR) data for the monitoring of tidal flats has been the subject of many research studies in the past decades. For example, in [1] the authors successfully used data from the AeS-1 X-band single-pass airborne interferometer to generate a DEM for a intertidal region on the Wadden Sea, in Germany. The generated DEM was obtained using an across-track baseline of 2.4 m, had a resolution of 5 m, and comparisons with ground control points showed differences with a standard deviation of less than 10 cm. In [2], shorelines extracted from ERS-1/2 complex SAR images were used to generate DEMs with around 12.5 m resolution. The authors reported good agreement between the obtained topography map and the aforementioned AeS-1 InSAR DEM. In [3] the authors used backscatter models and coherence analysis to discuss favourable conditions for the use of repeat-pass interferometry for the generation of DEMs over tidal flats. Following this research, results using ERS-1/2 pairs were reported in [4], where the challenges of obtaining high coherence data with spaceborne repeat-pass sensors were highlighted. The authors in [5] addressed the feasibility of monitoring the tidal flats by means of spaceborne repeat-pass InSAR suggesting the use of acquisitions with large across-track baseline and short temporal baseline to cope with the high scene dynamics. The considered temporal baselines interferometers and, additionally, in a repeat-pass (RP) configuration with a nominal vertical baseline of 40 m. For the single-pass interferometry, the half-baseline mode was employed, i.e., one antenna transmits while both antennas receive the SAR signals.
The use of the repeat-pass configuration is necessary due to the short single-pass baseline of the F-SAR (around 1.6 m), which by itself doesn't allow to meet the challenging vertical accuracy requirements. However, the single-pass data can still provide useful information: since master and slave data are acquired simultaneously, the single-pass interferograms tend to be less affected by residual motion errors and atmospheric artifacts, and can be usually used to provide the absolute reference for the repeat-pass data. Moreover, the single-pass interferograms do not suffer from temporal decorrelation. Hence, if desired, DEMs generated with single-pass data can be used to fill voids in the repeat-pass DEMs caused to strong decorrelation (e.g., over flooded areas or vegetation), albeit with reduced vertical accuracy.
The use of different center frequencies is required to solve the processing challenges imposed by the relatively large repeat-pass baseline (which corresponds to heights of ambiguity of less than 2 m for X-band, see Section 4.2). In particular, the phase unwrapping becomes non-trivial, but can be aided with the use of the redundant information in the dual-frequency data-set. Since repeat-pass acquisitions are considered, it is crucial that the data are acquired simultaneously in the different frequency bands, thus ensuring that the corresponding interferograms are affected by residual low-pass artifacts in a consistent manner. If that would not be the case, independent errors in the interferograms (e.g., due to uncompensated residual motion errors or propagation in the atmosphere), which often are of the order of the considered height of ambiguities, could hinder the ambiguity resolving procedure of the multi-channel phases.
The acquisition parameters used for both test data-sets are summarized in Table 1. For the Medem-Rinne test site, six parallel swaths and a perpendicular one were acquired, each one with approximately 3 km × 12 km extension (covering an overall area of around 165 km 2 ). For the Otzumer Balje test site, five parallel swaths and a perpendicular one were acquired, each one with approximately 3 km × 15 km extension (covering an overall area of around 155 km 2 ). The polarimetric image of the scene backscatter (S-band) for the Otzumer Balje test site is given in Figure 1.

Baseline Selection and Expected Performance
For the monitoring of tidal flats, digital elevation models with both high spatial resolution and high vertical accuracy are required. In the GeoWAM project, the goal is to demonstrate the use of InSAR for the recovery of DEMs on a 1 m × 1 m grid and with a relative and absolute vertical accuracy (2σ) in the order of decimeters (ideally smaller than 30 cm). With relative vertical accuracy we describe the influence of the statistical errors due to limited signal-to-noise ratio and small scale residual motion errors, while the absolute vertical accuracy requirement is with respect to a large scale bias with respect to the true topographic height. The desired relative vertical accuracy is only theoretically possible if the interferometric baseline is sufficiently large. On the other hand, the spectral shift should be kept to a minimum in order to preserve the range resolution. Moreover, very large baselines correspond to very small heights of ambiguity (e.g., less than 1 m), which are problematic due to the inevitable introduction of phase unwrapping errors. Hence, the baseline cannot be excessively large.
In order to choose the appropriate repeat-pass baseline, a performance analysis has been conducted. The analysis modeled the geometry of the repeat-pass acquisition as defined by the altitude above ground of the master acquisition plus the horizontal (Bh) and vertical (Bv) separation to the slave, as shown in Figure 2. A few options for Bh and Bv were considered (see Table 2), and the corresponding height of ambiguity, spectral shift and, eventually, height sensitivity were assessed. In all cases, we considered an output grid with 1 m × 1 m sampling, an altitude above ground of around 2400 m, and a valid off-nadir range from 25 • to 55 • .

Master
Slave |Bh| |Bv| 2400 m  The sensitivity analysis carried out here assumes that the only sources of decorrelation are the signal-to-noise ratio (SNR) and the temporal decorrelation (for the repeat-pass case only), i.e., we assume that a range-filter has been applied to compensate for the geometrical decorrelation whilst volume decorrelation is considered negligible. This last assumption is reasonable for wet tidal flats, where penetration in X-and S-band is extremely unlikely. However, it is not valid for semi-transparent media. Hence, the sensitivity analysis presented in the following does not apply to vegetated areas or very dry sand. Moreover, since our DFDB data is processed without a priori height information, the applied range spectral filter is not ideal for small features containing non-negligible slopes, e.g., dikes and water branches, in which case, additional decorrelation is expected. For such structures, it is therefore beneficial to use a high resolution input reference DEM in order to preserve fringes and coherence. The employed decorrelation model can be written as while the height accuracy is obtained by numerically solving the known relation [14] where f (φ; γ, L) is the interferometric phase probability density function (pdf) for a certain coherence γ and number of looks L [15], and h 2π is the height of ambiguity. For the considered F-SAR geometry, h 2π can be written as [14,16] with p = 1 for the half-baseline single-pass interferometer and p = 2 for the repeat-pass interferometer. The SNR decorrelation was predicted considering Noise-Equivalent-Sigma-Zero (NESZ) profiles estimated from fully polarimetric data of previous F-SAR campaigns (see Figure 3), and using the relation where σ 0 represents the assumed backscatter coefficient of the scene. According to values reported in the literature [3,10], and from our own observation of data of a previous campaign over the Jade bight reported in [11], the backscatter over the tidal flats can be expected to range from −5 dB to −25 dB, depending on water content, incidence angle and frequency of acquisition. Accordingly, for the repeat-pass baseline selection, backscatter values ranging from −5 to −25 dB were considered. For simplicity, in our analysis the backscatter in (5) is assumed to be constant over the swath. However, as mention above, an incidence-angle-dependent σ 0 is expected in the real data (see also Section 4.3.1) [17]. (b) S-band Figure 3. F-SAR Noise-Equivalent-Sigma-Zero (NESZ) of X-band (a) and S-band (b) for a height above ground of around 2400 m.
A constant value of temporal decorrelation was considered for each frequency band. For the single-pass scenario, γ temp = 1 was assumed, whereas for the repeat-pass one, temporal decorrelation factors of 0.8 and 0.9 were assumed for X-and S-band, respectively. These values were empirically selected considering our previous data acquisition over the Jade bight. For the γ temp selection, a simple analysis was conducted considering an area which mainly contained tidal flats. For this area, the ratio between repeat-pass and single-pass coherences was regarded as a measure of the temporal decorrelation. This assumes that differences in the system parameters/antenna diagram of the two channels of the single-pass interferometer are negligible, that geometric decorrelation is mostly accounted for by the range-spectral filter; and that backscatter changes between master and slave acquisition times can be neglected, at least for the majority of the considered region. The latter assumption, in particular, is not always true due to the high dynamics of the observed area. However, the analysis can still help to grasp the average temporal decorrelation value over tidal flats. The histograms of the repeat-pass to single-pass coherence ratio from the Jade bight data are shown in   Figure 4. Histograms of the ratio between repeat-pass and single-pass coherences for F-SAR DFDB data acquired over the Jade bight. The plot on the left corresponds to X-band, while the one in the right corresponds to S-band Figure 5 shows the expected height of ambiguity and spectral shift for the five repeat-pass cases in Table 2. The results on the top correspond to the X-band repeat-pass data, whereas the results at the bottom correspond to the S-band ones. Given its smaller wavelength, the X-band scenario corresponds to larger spectral shift and smaller height of ambiguity. Hence, we focus on its performance to evaluate which baseline case(s) offers the best compromise for our repeat-pass configuration. The plots at the top of Figure 5 show that the effective baseline of case 1 is too large, resulting in a spectral shift of more than 200 MHz for incidence angles smaller than 35 • (i.e., 2/3 of the available bandwidth), and height of ambiguity smaller than 0.5 m at near range. Case 2 and case 3 present similar metrics, with case 3 having slightly less performance variation across the swath. Both cases have a relative small height of ambiguity at near range (smaller than 1 m up to around 35 • incidence angle). Case 4 and 5 show less variation of the height of ambiguity with range and spectral shifts of 170 MHz and 60 MHz at 25 • . Moreover, they yield heights of ambiguity which are larger than 1 m for the entire swath. Hence, they are considered more appropriate for the reconstruction of the elevation models of the tidal flats, and are taken as candidates for the sensitivity evaluation considered next. Note that for cases 4 and 5, the expected spectral shifts for the S-band repeat-pass interferograms are smaller than 50 MHz (i.e., less than 20% of the available bandwidth) and the height of ambiguity is larger than 3 m for the entire swath. Figure 6 shows the expected repeat-pass height accuracy for cases 4 and 5 (2σ). The results consider that the height standard deviation can be modeled with Equation (2), i.e., only taking into account temporal and SNR decorrelation as error sources. Eventual residual low-pass artifacts or unwrapping errors are not considered at this point. For backscatter coefficients ranging from −5 to −15 dB, case 4 allows for a vertical height accuracy of 2σ < 30 cm for the majority of the swath in X-band. A backscatter of −25 dB will limit the valid swath to around 47 • , if only X-band repeat-pass data is employed for the generation of the final product. Case 5 has worse performance in general, and for the −25 dB backscatter case, the accuracy of X-band repeat-pass goes beyond 50 cm for both near-and far-range. On the positive side, it presents more homogeneous performance within the swath in comparison to case 4. The real flight trajectory deviates from the ideal one, leading to baseline errors up to +− 2 m for the horizontal component and up to −+ 6 m for the vertical component in the case of the F-SAR carrier. As a consequence, the expected DEM performance will also vary along azimuth. Note that baseline errors are more critical for case 5, since its X-band repeat-pass performance is already worse than the assumed threshold. Moreover, since the positive horizontal baseline cause a non-negligible worsening of the performance at near-range, and can, in extreme cases, lead to invalid data due to the alignment of master and slave with the line-of-sight, the optimal baseline choice among the considered ones is the one of case 4, i.e., [0,−40] m.
Finally, note that even for case 4 and considering no additional error sources or baseline deviations, the use of repeat-pass S-band alone is not enough to fulfill the accuracy requirement. Nevertheless, given the negligible penetration expected over the tidal flats [10,11], the S-band height models can be used to generate a final merged XS DEM. In this case, the relative accuracy of the resulting DEM should be at least as good as the X-band one, provided that the merging process accounts for the individual DEM statistics. Moreover, dual-frequency filtering strategies will be employed to further reduce the noise level (see Section 4.2.7).
For completeness, the obtained height accuracy of the single-pass height models are given in Figure 7. As expected, the requirement cannot be fulfilled when using solely single-pass DEMs, since they have an accuracy in the order of meters for typical tidal flats' backscatter values.

Data Processing
This section summarizes the SAR focusing, i.e., from raw data to focused and calibrated single-look-complex (SLC) images, as well as the DEM processing, i.e., from SLCs to mosaicked digital elevation models.

SAR Focusing and Calibration
Dedicated acquisitions over the F-SAR calibration site at Kaufbeuren, Germany were used to derive calibration corrections for the processing of SAR data gathered in the field. The calibration site features nine trihedral radar reflectors and one dihedral reflector deployed at off-nadir angles ranging between 20 • and 60 • . The responses of these reference targets were used as input to the calibration procedure described in [18,19], which yields, among other things, precise antenna phase centre baselines for the S-and X-band single-pass interferometers as well as calibration constants for the range delays, the channel gains and the inter-channel phase differences.
In addition, the calibration included the estimation of a residual antenna phase error to correct for small phase inaccuracies in the on-ground characterisation of the antenna patterns. These errors derive from multi-path effects that cannot be faithfully reproduced in the antenna measurement, which only includes the antenna mount but not the entire aircraft, as well as phase noise in the quiet zone of the antenna test range. This antenna phase correction estimate was carried out using interferometric phase differences measured over distributed targets in the range-Doppler domain to yield a 2D phase correction, in terms of off-nadir angle and squint, for each slave antenna of the S-and X-band interferometers. The squint-variant estimate helped to ensure that the correction is applicable for all acquisition geometries encountered in the field. The estimated phase correction for the X-band antenna is illustrated in Figure 8. The phase screen is added to the antenna pattern correction of the single-pass slave channels.

DEM Generation
The DEM generation procedure is based on the dual-frequency/dual-baseline (DFDB) chain introduced in [11,20] and considers a few simplifications in the calibration procedure. The main steps of the chain and main differences to the scheme in [11] are summarized in the following. Moreover, we include details about the mosaicking procedure.

Generation of Interferometric Products
The first step is the generation of the interferometric products. A range-adaptive spectral filter is applied to the repeat-pass data and a mean reference height is considered for the flattening [21].
The multi-squint approach is applied for the correction of residual motion errors in the repeat-pass interferograms (up to constant and linear components) [22]. The estimation is performed using the data from a single frequency of acquisition and applied to both data-sets. The narrower aperture of the X-band antenna in comparison to the S-band one allows for the estimation of higher frequency components of the residual motion errors. Hence, the X-band data-set should ideally be used for the estimation. However, in the results presented in the next section, the estimation was performed based on S-band data due to the large extend of flooded areas, which leads to small residuals on the X-band results. Nevertheless, consistence of the S-and X-band interferograms was found to be sufficient to allow for the phase unwrapping process. Data from next campaigns will be processed considering the X-band as reference.
After the interferograms have been formed, the interferometric phase and coherence are extracted and are interpolated to the geometry of the X-band master image.

Dual-Frequency Phase Unwrapping
The second step is the unwrapping of the interferometric phases and is performed using the dual-channel region-growing approach we suggested in [20]. The approach profits from the redundancy of the dual-frequency data-set, while using contextual information in a region-growing scheme in order to introduce robustness to noise. The algorithm is congruent, i.e., wrapped and unwrapped pixels differ only by 2π multiples, avoiding the introduction of noise induced slopes. As in a conventional single-channel region-growing unwrapping approach [23], the 2π multiple is estimated using the phase gradient computed from available unwrapped pixels in a certain estimation window. Additionally, in our approach, we also enforce the consistence of the gradients in the two channels (after accounting for the different wavelengths and geometries) and the consistence of the underlying height information of the pixel being unwrapped. Moreover, the growing scheme also contains reliability tests based on the statistics of both interferograms, which promotes the choice of the "easiest" unwrapping path, potentially leading to a more robust unwrapping.

Baseline Calibration
After the phase unwrapping, a calibration step is performed. Unlike our previous strategy in [11], the estimation of multi-path artifacts caused by a secondary reflection on the aircraft fuselage is not applied for the GeoWAM data. This is because the F-SAR processing chain now includes the estimation of a residual antenna phase screen which accounts for this effect. This estimation is performed offline using data from calibration campaigns (see Section 4.1), and is considered to be accurate enough for achieving the desired absolute accuracy. The calibration performed at this stage then accounts only for constant and linear baseline errors and a global offset.
The strategy adopted here is a simplified version of the baseline calibration we presented in [11]. This is because within the GeoWAM project, a 10 m resolution ALS DEM acquired in 2014 is available. This information allows for the estimation of the global offset affecting the single-pass interferograms and, consequently, for a better conditioning of the repeat-pass calibration problem [11].
More then 250 ground control points (GCPs) were also measured in situ during the GeoWAM campaign, and are intended to be used for the validation. However, the use of the GCPs for the calibration of the single-pass interferograms is not recommended, since there are not enough points in order to provide an absolute calibration which could lead to a maximum height error in the order of few centimeters. To illustrate this, a Monte-Carlo simulation was performed considering the GeoWAM X-and S-band single-pass parameters. In each realization, a certain number of reference points were simulated with randomly distributed phase values. The statistics of the simulated phases were derived from coherence values typically observed in our single-pass GeoWAM data, and the available number of looks. After the computation of the phase offset using the simulated GCPs, the obtained maximum height error (i.e., error at near-range) corresponding to the offset mismatch is assessed. Figure 9 shows the statistics of this error for varying number of GCPs and coherence values. The plots show, e.g., that over a 100 points with coherence larger than 0.95 are required for obtaining an average error of less than 25 cm in X-band, and over a 1000 points are required to ensure a robust offset estimate leading to errors of less than 10 cm. In the GeoWAM case, the number of GCPs meeting the condition per stripe was typically smaller than 25. In this case, the baseline calibration based on a reference DEM is preferred [24]. Once the single-pass phases are calibrated with the reference DEM, they are used for the calibration of the repeat-pass data as described in [11], i.e., using a conventional baseline error model and employing a joint X-and S-band complex minimization functional to avoid divergence due to either low-pass artifacts or unwrapping errors.

Correction of Unwrapping Errors
After the trends and global offsets have been calibrated, unwrapping errors are detected and corrected using the DFDB approach we presented in [11]. This step is necessary mainly due to the presence of medium and large scale unwrapping errors, which relate to the fact that different areas of the image might be completely isolated from each other, preventing the use of contextual information. In the case of tidal flats, this occurs mainly for sand banks surrounded by flooded areas and water branches. The identification of the unwrapping errors is divided in two steps. The first step is a rough detection of suspicious pixels using the difference between the dual-baseline products (i.e, single-and repeat-pass of a certain frequency band). The evaluation of this difference leads to a mask indicating possible pixels affected by unwrapping errors and the corresponding 2π multiples necessary for their correction. However, given the large difference between single-and repeat-pass heights of ambiguity, this first estimation is severely impacted by noise. Moreover, it is impacted by low-frequency artifacts which affect mostly the repeat-pass products. The second step is the refinement of the estimated error map. This is accomplished by including the difference between the dual-frequency products (i.e., X-and S-band) in the scheme and considering a regularized active-contour-based algorithm. When ensuring the consistence between the repeat-pass products we avoid biases caused by low-frequency artifacts, and when searching for regularized borders of relevant size for the affected areas, we introduce robustness to noise.

Phase-to-Height Conversion
In the next step, the calibrated phases are individually transformed to height maps, still in slant-range geometry, and the expected height error maps based on the interferometric coherences are calculated. At this point in the chain, we have the height maps (and height error maps) corresponding to single-pass X-band, single-pass S-band, repeat-pass X-band and repeat-pass S-band.

Residual Calibration
At this stage, the repeat-pass data-sets are still impacted by uncompensated low-frequency artifacts. These might originate from atmospheric disturbances, but can also be related to limitations of the residual motion compensation algorithm (e.g., due to the presence of large incoherent areas spanning through the whole imaged swath, which is typical for tidal flat regions). Hence, a residual calibration step is performed after the height conversion. As proposed in [11], we consider here that Xand S-band height maps are affected by the same errors. Moreover, we consider that the errors have a limited low-pass spectrum and that they can be estimated using a simple low-pass filter. For the work presented here, a frequency domain Gaussian filter was selected. Hence, the filter frequency response is governed by a single-parameter, σ f , which is estimated by exploiting the redundancy provided by the DFDB data-set. Namely, we consider the two available dual-baseline errors maps, i.e., the X-band single-to repeat-pass height difference and the S-band one. The filter parameter is then estimated by minimizing the difference between the unfiltered error map with smallest expected noise level (h diff,L ) and the filtered error map with largest expected noise level (h diff,H ), i.e., where G is the filter frequency response and F () −1 represents the inverse Fourier Transform. As discussed in [11], the use of wider filters causes the out-of-band noise in h diff,H to increase the objective function in Equation (6). On the other hand, the use of narrower responses cause information about the low-frequency error spectrum to be missed, also increasing Equation (6). Hence, the used minimization scheme offers a compromise. The estimation is performed block-wise using blocks of around 1 km × 1 km with an overlap of approximately 200 m, thus allowing for a certain spatial adaptability which is relevant specially in the case of atmospheric artifacts. The main assumption of the residual calibration strategy is that the single-pass height maps do not contain low-frequency artifacts. This assumes that (1) the antenna correction discussed in Section 4.1 is enough to lead to negligible multipath artifacts; and (2) that the single-pass DEMs are not affected by residual motion errors. If that is not the case, and no correction is applied, the repeat-pass data will inherit the single-pass artifacts. If the first assumption is not true, the proposed data-dependent multipath solution suggested in [11] can be used. In case the second assumption cannot be verified, the available low-resolution reference laser DEM can be used to perform a parametric estimation/removal of residual baseline errors prior to the estimation of the repeat-pass errors, e.g., in a similar way as proposed in [25]. In this case, we recommend to pre-process the residual height map by • considering the repeat-pass coherence to filter out flooded regions and vegetated areas, which are not consistent in the airborne SAR and reference DEMs. • consider an outlier removal in order to account for differences in the topography, which might have happened in between acquisitions.
If the pre-processing leads to invalid azimuth bins, the estimated azimuth profile is interpolated assuming the smoothness of the residual motion error.

Generation of the Combined Repeat-Pass Height Map
At this stage, the X-band and S-band repeat-pass interferograms are merged into the final elevation model considering a wavelet-based noise mitigation strategy, as described in [11]. For that, the multi-resolution wavelet decomposition of X-and S-band height maps is performed. The approach then considers that most of the energy from the topography information is compressed into large magnitude coefficients at higher (or coarser) scales, while the wavelet coefficients of fine features, e.g., small terrain fissures, can be spread throughout the lower scales along noise. Since for tidal flats the useful topography content is assumed to be the same in S-and X-band height maps, their wavelet coefficients correlate at both higher and lower scales. On the other hand, noise contributions are independent, and no correlation of the coefficients is expected. Hence, by given more weight to wavelet coefficients with higher XS correlation prior to the inverse discrete wavelet transform, an efficient regularization is achieved. Moreover, a combined XS product is obtained by considering in addition a weighted average of the wavelet coefficients of the X-and S-band decompositions using data-derived noise statistics.
In the GeoWAM project we do not merge the information of SP and RP interferograms, i.e., the combined repeat-pass DEM will be marked as invalid for flooded areas, whereas the single-pass one might contain the information of the water level (depending on the roughness of the water surface).

Geocoding
The retrieved height maps can now be interpolated back to their master geometry (when necessary). The last step of the DFDB chain for a single stripe is the geolocation of the obtained SP, RP and combined elevation models, i.e., their transformation from radar geometry to UTM coordinates.

Mosaicking
Since several stripes are available, the final stage of the "DEM Generation" is the mosaicking of the available products. In this study we choose to handle each product independently, i.e., we have at the end: mosaicked X-band SP, S-band SP, X-band RP, S-band RP and XS-band RP DEMs. The procedure takes into consideration the expected relative height errors derived in the "Phase-to-height conversion" step and, additionally, performs feathering in order to avoid border artifacts. The mosaicked height can be simply written as where h k is the k-th acquisition (of N available acquisitions) of the height value for a certain pixel and w k are the mosaicking weights. In the GeoWAM case, for two adjacent DEM stripes we have either near-range/near-range overlap (around 200 m) or far-range/far-range overlap (around 800 m). Considering additionally the overlap with the orthogonal stripe, the number of available height measurements for a certain pixel varies from 0 to 3. The feathering is approached using the same philosophy used for the mosaicking of the TanDEM-X global product [26,27], i.e., the weights are artificially decreased towards the borders, We start with the explanation of the denominator of the previous equation: σ h,k is the expected height standard deviation and σ border,k is the artificial standard deviation introduced to promote smooth border transitions. Similarly to what is suggested in [27], we define this standard deviation as a function of the distance to the border of the valid polygon, where ω feat is the square feathering window with length N feat , δ x and δ y are the samplings in both directions, and M polygon,k is a validity mask which assumes 1 if the pixel is inside the valid geocoded polygon and 0 otherwise. Note that this is not the data validity mask (M val,k ), which might also assume 0 due to invalid unwrapping or due to a expected σ h,k higher than a certain threshold. This distinction is necessary to ensure that the border handling is only applied to the actual borders of the valid polygon. The definition of σ border,k in (9) ensures that pixels close to this border are assigned an artificial standard deviation which is much larger than the expected σ h,k . On the other hand, for pixels which are closer to the center of the data stripe, no artificial standard deviation is introduced and the mosaicking weight is entirely defined by the expected height standard deviation. Note that starting from a pixel at the border and moving towards the outside of the valid polygon, σ border,k in (9) decreases, reaching 0 if M polygon,k = 0 for the entire feathering window. However, this does not impact the mosaicking since a pixel outside the valid polygon is invalid (i.e., M val,k = 0) and, according to (8), it is not considered for the merge operation. The term m k in (8) is included as a consistence check and intends to mask out pixels which considerably differ among the overlapping products. This can happen, e.g., due to extreme noise which is underestimated by the data derived statistics, unwrapping errors or residual low-pass artifacts. This term also promotes a smooth transition between adjacent stripes and avoids that bad pixels have relatively larger weights due to feathering. It is calculated as where h k,min corresponds to the height measurement of the overlapping pixel with the minimum expected standard deviation, and σ h,min to its associated height standard deviation.

Auxiliary Masks for Data Interpretation
In addition to the derived DEMs, the monitoring of the tidal flats using InSAR can also profit from the availability of the SAR amplitude and coherence maps, and from the fully polarimetric information available in S-band.
In this section, we describe the generation of an interferometry and a polarimetry mask which aid on the interpretation of the land cover. For example, the auxiliary masks can be used alongside freely available land use and land cover products such as OpenStreetMap [28] and CORINE land cover [29] to produce reliable fused training data sets that are suitable as input data for machine learning classification approaches, an on-going research topic of the GeoWAM university partners [30,31]. Moreover, information retrieved from the interferometry mask can be used to mask the GeoWAM DEM products at a post-processing step.

Interferometry Mask
The interferometry mask uses the information of the S-band VV amplitude of master and slave, as well as the S-band repeat-pass coherence, and it can help, e.g., with the separation of flooded areas from non-flooded ones.
The repeat-pass coherence is useful to distinguish between water and ground, since the former completely decorrelates in repeat-pass interferograms. However, semi-transparent media (e.g., forests or crops) and areas with steep or discontinued topography (e.g., buildings) also tend to present very low coherence in the repeat-pass interferograms. On the other hand, contrary to water, such targets usually present high amplitudes, which are stable in master and slave acquisition passes.
Tidal mudflats usually correspond to medium-low to high repeat-pass coherences depending on the water dynamics and moisture. They tend to present low amplitude, which can even be lower than the one of water if the acquisition is performed during strong wind conditions. Moreover, due to the tidal dynamics, the water level might vary from master to slave pass, causing the amplitude to change. Hence, for the data interpretation it is important to consider both master and slave amplitudes.
The generation of the interferometry mask is conducted in three steps. First, a post-processing amplitude correction is performed in order to remove the incident angle dependent behavior observed over water and mudflats (see, e.g., brighter areas in between stripes in the mosaicked amplitude in Figure 1). The second step is the computation of coherence and amplitude thresholds which will be used to differentiate between coherent and incoherent regions, and between incoherent regions with low and high amplitudes. The final step is the actual generation of the mask according to the classes defined in Table 3. After the assignment of the classes, a morphological opening is applied to mitigate the effects of estimation biases in low coherence areas. The focusing of the F-SAR data includes a default radiometric calibration which yields amplitude images corresponding to the radar brightness (β 0 ). Normalization to the ground plan (σ 0 ) or to the wave-front plan (γ 0 ) can be done a-posteriori, considering the incidence angles [32]. However, observation of the GeoWAM data shows that none of these corrections is enough to remove the incidence-angle dependence of the resulting amplitudes over water and tidal flats. As a consequence, the resulting masks would be affected by incidence-angle dependent artifacts. The incidence dependent scattering of water and rough surfaces has been widely researched by the scientific community [17,33] and the findings show that water and mudflats can present different behaviors and should, in principle, be characterized independently. However, in this study, we implemented an empiric evaluation of the available data to ensure amplitude normalization across the swath using a simple scalar threshold. For that, a stripe which is known to contain mostly water and tidal mudflats is first selected by the user. A profile along range is generated using only areas of low coherence, i.e., presumably corresponding to water. A polynomial is then fitted to the profile (as a function of the incidence angle) and the coefficients are used to create a correction for all stripes of the mosaic, which flattens the amplitude over water in an optimum and those of tidal flats in a sufficiently approximate way. This is considered to provide sufficiently reliable input data for subsequent thresholding. Figure 10 shows, in the left, the obtained amplitude profile of the chosen reference stripe in the Otzumer Balje test site. The plot shows, additionally, the fitted polynomial, and the σ 0 and γ 0 correction terms. As it can be observed, unlike the fitted polynomial, σ 0 and γ 0 corrections are not enough to flatten the amplitude. The plot on the right show the same type of curves. However, in this case, the polynomial fit is obtained from the coefficients estimated with the reference track. It is possible to see that removing the polynomial correction will also lead here to a fairly good amplitude flattening.
For the computation of the coherence threshold, the histogram of the repeat-pass coherence of the entire acquired mosaic is computed. Typically, two main local maxima are identified: the first one corresponds mainly to flooded areas, vegetation, areas of layover and shadow, urban features with sharp height changes, and tidal flats presenting low coherence. The second main lobe, at higher coherence values, generally includes tidal flats which remain dry in both master and slave acquisitions, bare soil and grassland. The threshold is selected by searching the local minimum in between these two main histogram lobes. In the left, Figure 11 shows the obtained coherence histogram for the Otzumer Balje test site, and indicates the estimated threshold.  For the computation of the amplitude threshold, only areas of low coherence are considered (identified by the previously calculated coherence threshold). The main goal is to separate areas of low coherence with bright amplitude from those with low amplitude. A typical amplitude histogram of low coherent repeat-pass InSAR areas for the GeoWAM campaign will contain at least two lobes at lower and medium amplitudes: the main one containing flooded areas (i.e., water) and a smaller one including tidal flats of very low backscatter leading to SNR decorrelation. As wind conditions are affecting the backscatter from the water areas, the position of the two lobes may be exchanged. Moreover, given the nature of the imaged test-sites, an even smaller peak should be present at higher amplitudes, corresponding to vegetation and urban areas, both causing low coherence due to volume decorrelation. The amplitude threshold is set between the lobe corresponding to higher amplitudes and the remaining ones. In the right, Figure 11 shows the obtained amplitude histograms of the low coherent areas for the Otzumer Balje test site, and indicates the estimated threshold. The square root of the occurrences is shown to promote the visualization of the lobe at higher amplitudes. As we generally do not expect amplitude differences between master and slave (except for few tidal channels), we use only the master image amplitude for the computation of the threshold.  Figure 12 shows the interferometry mask computed for the Otzumer Balje test site. Water can be seen in dark blue, while the majority of the tidal flats appears in pink.

Polarimetry Mask
In the case of the polarimetry mask, the standard entropy-alpha unsupervised classification scheme was performed [34]. The entropy-alpha classifier utilizes the polarimetric information that is not considered in the interferometry mask. This allows for the removal of inconsistencies in the interferometric mask due to difficulties in separating vegetation from low-coherence mudflats.
A small modification on the ordering of the original entropy-alpha classes was performed, going from low to high entropy within the same type of reflection (see Figure 15 and Table 4). This reordering step is necessary since feathering on overlapping classes is not defined. For the mosaicking of labels, a ceiling approach (taking the highest class value on overlapping pixels) is applied instead of border feathering, prioritizing the preservation of high entropy classes. Nonetheless, some areas of the polarimetric mask are still affected by different water levels and incidence angles leading to residual border effects during the mosaicking. Figure 13 shows the polarimetry mask computed for the Otzumer Balje test site. Figure 14 shows a small crop where the polarimetry information separates an area detected as class 4 in the interferometric mask (shown as pink) into surface and short vegetation scattering (blue and green in the polarimetry mask, corresponding to classes [0,1 and 4] in Figure 15).

Results
In this section, preliminary results from the Otzumer Balje test site are given. Figure 16 shows the retrieved X-band single-pass DEM, while Figure 17 shows the fused XS repeat-pass DEM. No post-processing masking was applied. As expected, flooded areas are naturally marked as invalid in the repeat-pass DEM, due to their low coherence. On the other hand, the single-pass DEM contain valid values for the water surface. For the results shown here, the low-resolution ALS DEM acquired in 2014 (10 m resolution) was used as reference to correct for residual low-pass artifacts in the single-pass DEMs, as explained in Section 4.2.6. However, residual artifacts are observed for areas which contain water for the majority of the swath. This is either due to interpolation errors of the estimated residual baseline errors (when no data is available) or poor conditioning of the minimization problem (when only a small portion of the swath has valid data).
The differences between the obtained SAR DEMs and the available full-resolution ALS DEM acquired in 2014 (1 m resolution) are shown in Figures 18 and 19, scaled from −1 m to 1 m. Overall, a good consistency is observed between the products. Strong changes in the coastal areas in the north are also visible (larger than 1 m), as well as smaller changes within the tidal flats. Changes in vegetation height are also clearly visible in the fields in the coastal mainland areas in the south. Figure 20 shows a zoom of around 1.5 km × 1.5 km located at the north-west portion of the image for the different SAR products: (top left) X-band SP, (top right) S-band RP, (bottom left) X-band RP and (bottom right) XS RP. As expected, the noise level improves from the single-pass DEM to the fused repeat-pass one. This can be also seen in Figure 21, which shows the normalized histogram of the differences between laser reference and SAR products for the region of interest. Note that the differences encompass not only noise, but also residual low-pass artifacts and topographic changes. Nevertheless, we obtain standard deviations of around 26 cm, 16 cm, 15 cm and 13 cm for the differences, attesting for the good quality of the SAR products.

Summary and Discussion
This paper presented an overview of the GeoWAM project which aims for the monitoring of tidal flats by means of SAR interferometry and polarimetry. Whilst the project is ongoing, we presented the adopted SAR processing approach and first results demonstrating the potential for deriving accurate high resolution DEMs as well as supplementary information like indication of water flooded areas and vegetation cover. The presented results focused on the test-site of Otzumer Balje, which is a tidal basin of 118 km 2 south of the island of Spiekeroog. The airborne F-SAR data have been acquired simultaneously in X-and S-band within the tidal window of +/− 1 h around low tide. Five partially overlapping parallel stripes were flown in repeat-pass configuration with minimum time lag of approximately 10 min between the overpasses and an optimized vertical baseline of 40 m. The results we presented indicate the feasibility to generate large area DEMs with high spatial sampling of 1 m and relative and absolute height accuracies in the order of 1-3 dm. To ensure the bias-free height reconstruction, we relied on a heavily sub-sampled reference DEM generated 5 years earlier. Its high resolution version served as input for a first comparison and for quick identification of changes, some of which we reported. We also computed mosaics of interferometric and polarimetric masks from the radar data to support the further processing and evaluation of the DEMs. The same methodology is currently being applied for data of the Medem-Rinne test-site, located at the estuary of the Elbe river and covering 150 km 2 . These have been collected along 6 parallel stripes within a second low tide window.
Future work includes the tuning of the processing chain to improve the quality of the retrieved DEMs, especially regarding the presence of residual low-pass artifacts in the single-pass interferograms. The present computation of the interferometry and polarimetry masks is partially deteriorated by incidence angle dependent properties of the backscattering. Our first attempt to normalize the amplitudes across the swath for the interferometric mask may be considered sufficient to discriminate water from non-water surfaces, but partially failed for other specific land cover types. Therefore, a class dependent normalization will be targeted in cooperation with the project partners dealing with the advanced polarimetric classification of the acquired data. However, the most important subsequent step is the validation of the obtained DEMs with laser models acquired by an ALS sensor a few weeks after the F-SAR acquisition, a task which is to be performed primarily by the federal authorities with support from the project partners. The obtained SAR DEMs as well as the classification maps generated by the partners will subsequently be used to create a final digital bathymetry model of the complete basin, which includes as additional input profile measurements from ship-borne acoustic echo sounding. For this purpose, the SAR interferometry mask turns out to be an essential input for efficient planning of the ship routes to the few remaining water covered areas at low tide. Finally, a further campaign is planned for the summer of 2020, and the newly generated results will be used to assess the annual changes.