Method for the computer-aided generation of a data-based model for the computer-aided processing of digital SAR data

Abstract

The invention relates to a method for the computer-aided generation of a data-driven model for the computer-aided processing of digital SAR data from a repeat-pass SAR system for determining surface deformations. The SAR data consist of a plurality of interferograms (Γ), each of which is obtained from a first and a second SLC acquisition (s1, s2), as well as a Doppler center frequency (f). dc ), an interferometric phase (ϕ̂) and a coherence (ρ̂), all obtained from a respective pair of first and second SLC acquisitions (s1, s2). The model is configured to produce a zero-Doppler phase (ϕ̂). zd ), an along-track phase (ϕ at ) and an estimate of an associated along-track surface motion (u at ) to determine. The generation of the model comprises the following steps: - Generating a synthetic training dataset (TDS) comprising a multitude of interferograms (Γ), each formed by a pair of a first and a second SLC acquisition (s1, s2), by inducing a model-based and external 3D deformation into a given first SLC acquisition (s1) to obtain the associated second SLC acquisition (s2), which relates to a time point after a deformation event, and by applying a reverse approach defined by a given coherence (ρ) and a noise-free interferometric phase (ϕ). theo ) is specified according to the deformation event; - Training a neural network (NN) with the training data set (TDS), wherein the neural network (NN) has an input for inputting the multitude of interferograms (Γ) and an output for outputting the along-track phase (ϕ). at ), the zero-Doppler phase (ϕ zd ) and the associated along-track surface motion (uat ) for the multitude of interferograms (Γ), where the trained neural network (NN) represents the generated data-driven model.

Description

[0001] The invention relates to a method for the computer-aided generation of a data-based model for the computer-aided processing of digital SAR data, as well as to a method for the computer-aided processing of digital SAR data and a computer program product.

[0002] Terrain Observation by Progressive Scans (TOPS) is a detection mode used in some modern Synthetic Aperture Radar (SAR) imaging systems to increase coverage at the expense of azimuth resolution. The SAR system's antenna uses an electronic beam steering mechanism to direct the transmitted and received signals toward a specific terrain area. The distinguishing feature of the TOPS mode is the detection of bursts by steering the antenna beam in the direction of flight, from back to front.

[0003] Sentinel-1 (S-1) is the European Space Agency's (ESA) satellite constellation for SAR technology, which mainly operates in TOPS detection mode: It continuously observes the Earth's surface with a repeat pass of 6 to 12 days and monitors surface deformations (offsets) using the well-known DInSAR (Differential Interferometric Synthetic Aperture Radar) technique.

[0004] In DInSAR, two SAR images of the same area are acquired at different times. The phase difference between the two images determines the surface motion along the line of sight (LOS), i.e., the direction from the satellite to the imaged surface point. Whenever S-1 acquires data in TOPS mode, the beam swap along the flight path also creates a variable sensitivity of the interferometric system along this direction. Specifically, the sensitivity for motion along the track in the scene increases as the absolute value of the azimuth shift angle increases (at the beginning and end of the burst), so that it has zero sensitivity in the middle of the burst.

[0005] This varying sensitivity is a disadvantage because it generates phase jumps between adjacent bursts, which can disrupt accurate phase detection across the track and potentially also interfere with phase unwrapping. On the other hand, this property can be used to detect offsets along the track near the beginning and end of a burst, provided the phase jump can be estimated correctly. The ability to accurately separate the phase sources generated by the two offset directions allows, firstly, the correct phase unwrapping in the across-track direction and the calculation of the corresponding deformation, and secondly, the provision of a proxy for the inversion of the along-track deformation.

[0006] Two methods have been proposed to solve the problem of source separation in the interferometric TOPS phase.

[0007] In [1], local coregistration is first performed to reduce phase discontinuities and increase coherence. Specifically, this strategy consists of evaluating distance and azimuth shifts using incoherence cross-correlation (ICC) techniques for non-stationary scenarios, followed by rescanning the secondary image in the primary grid and further refinement through spatially configured distance and azimuth shift measurements using spectral diversity (SD) techniques. This solution simplifies phase unwinding, as the interferograms no longer exhibit phase jumps, but it compromises interferometric phase quality because the shift estimation is less accurate and less resolved than the interferometric phase.After phase unwinding, the residual phase, which was implicitly subtracted in the coregistration step, can be added back to the interferogram, thus restoring full interferometric quality for the unwound phase. Therefore, speckle tracking is proposed as a possible solution in reference [1]. This approach only works for extremely large offsets where the interferometric phase information cannot be used. For small offsets, however, it has the disadvantage of potentially increasing the noise of the interferometric phase in low-quality regions (low interferometric coherence).

[0008] Alternatively, [2] proposes a methodology in which the differential interferogram in each overlapping region of the scene is used to initialize an image inpainting algorithm that fills the missing regions, i.e., the gaps between the overlaps, with an offset estimated from a digital elevation model (DEM). The result is the zero-Doppler interferometric phase. Therefore, it is suggested to use overlapping regions of the bursts to estimate the offset only in that region. Although the information present in the overlapping regions is reliable, the fact that the overlapping regions are so small relative to the burst necessitates further processing, such as inpainting / interpolation, to extend the information from the overlap to the entire burst.This technique has several disadvantages: (a) The offset estimation method requires error-prone 2D phase unwrapping, (b) the overlap region is very small compared to the burst, so the offset in the burst is ignored by extending the information from the overlap to the burst via continuity, (c) the computational complexity is higher than that of any other state-of-the-art algorithm. Therefore, while this technique provides better results in the overlap regions than those disclosed in [1], it can still exhibit phase jumps and inaccurate offset estimates.

[0009] In [3], the synergy of orbit directions (descending and ascending) and geometric components (along and perpendicular to the flight path) is proposed for acquiring complete 3D surface offset fields over seismic areas. After global refinements of the azimuth offsets, interferograms are flattened by fitting a bilinear polynomial surface to the difference between external GPS measurements projected onto the corresponding geometry (along or perpendicular to the flight path) and the nearest pixel in the interferogram. This fills gaps between bursts in each of the four interferometric phases and applies phase unfolding to the offset measurement.

[0010] In [4] and [5], a modified Sentinel-1 interferometer chain is proposed that reduces azimuth coupling to the line-of-sight (LOS) direction through refined coregistration. Specifically, offset matrices are calculated using an external DEM, and the x and y components of the multi-year average ice velocity map are distributed by the Geological Survey of Denmark and Greenland (GEUS). As a result, the use of the updated DEM reduces the phase jumps observed in the interferometric phase and improves coherence in the overlapping regions.

[0011] External references, such as GPS measurements as used in [3], or maps with multi-year average ice velocities as used in [4] and [5], which are usually used to support satellite-based solutions, suffer from the following problems: (a) they are scarce, (b) they are not updated, (c) they are not always freely available, and (d) their accuracy is insufficient in most cases to obtain accurate interferometric measurements.

[0012] EP 3 896 482 A1 describes a computer-implemented method for generating SAR training datasets using a convolutional neural network for SAR interferometry. Interferometric phase images are created from a digital elevation model, where the elevation model contains elevation information from a captured SAR image, and SAR training data containing noise is generated using a neural network.

[0013] CN 1 17 572 420 A describes a method for generating interference-containing training data by fitting a digital elevation model. Based on a U-Net++ network, an extended phase-development network is developed. A measured and an artificially generated dataset are combined into a single dataset used to train the network.

[0014] Reference

[14] describes the detection of deformations, e.g., caused by volcanic activity, in SAR images using neural networks. The neural network is trained using artificially generated training data.

[0015] The object of the present invention is to provide an improved method for interferometric phase source separation using terrain observation in progressive scan mode (TOPS).

[0016] This problem is solved by the method for computer-aided generation of a data-based model according to claim 1, the method for computer-aided processing of digital SAR data according to claim 16, the device for computer-aided processing of digital SAR data according to claim 17, and the computer program product according to claim 18. Further embodiments of the invention are defined in the dependent claims.

[0017] According to a first aspect of the present invention, a method for the computer-aided generation of a data-driven model for the computer-aided processing of digital SAR data from a repeat-pass SAR system is proposed. The SAR system is configured for Terrain Observation by Progressive Scans (TOPS) as its acquisition mode and thus serves to determine surface deformations in a specific area. The data-driven model, in the form of a neural network, is used for interferometric phase source separation. In other words, the data-driven model is intended to process the digital SAR data provided by the SAR system.

[0018] The digital SAR data consist of a multitude of interferograms, each obtained from a first SLC acquisition of the SAR system at a first time point and a second SLC acquisition of the SAR system at a second time point later than the first. As is generally known to those skilled in the art, in single-look complex (SLC) processing, the data are focused in azimuth and range to generate an image. The digital SAR data further consist of a Doppler center frequency, an interferometric phase, and a coherence, all of which are obtained from a respective pair of first and second SLC acquisitions.

[0019] The data-driven model is set up to determine a zero-Doppler phase as the phase contribution associated with the deformation causing a zero-Doppler LOS variation, i.e., the projection of the deformations along the vertical direction, an along-track phase as the phase contribution associated with the deformation in an along-track direction, and an estimate of an associated along-track surface motion.

[0020] The creation of the data-driven model includes the following steps: First, a synthetic training dataset is generated, consisting of a multitude of interferograms. Each interferogram is formed by a pair of first and second SLC acquisitions. This is achieved by inducing a model-based, external 3D deformation into a given first SLC acquisition to obtain the corresponding second SLC acquisition, which relates to a time point after a deformation event. A reverse approach is then applied, specified by a given coherence and a noise-free interferometric phase after the deformation event. The deformation event could be, for example, an earthquake or a glacial flow in the area where the first acquisition occurred.

[0021] After generating the synthetic training dataset, a neural network is trained using this dataset. The neural network includes an input through which the multitude of preprocessed interferograms is fed to obtain the Doppler center frequency, interferometric phase, and coherence. The neural network also includes an output through which the along-track phase, zero-Doppler phase, and associated along-track surface motion for the multitude of interferograms are output. The trained neural network represents the generated data-driven model.

[0022] The invention has the advantage that phase source separation can be performed without the explicit implementation of a rule and without the use of a priori information. For this purpose, a suitably trained neural network is used, which can be flexibly adapted as needed by retraining it when new, pre-existing datasets become available. The invention can preferably be used for SAR systems that employ TOPS as the acquisition mode.

[0023] In general, the generation of the synthetic interferometric acquisition pair is based on the application of a "reverse approach." Given the desired coherence and noise-free interferometric phase after the deformation event, the signal model induces the deformation (also referred to as offset) in the second SLC acquisition, which is assumed after a deformation event. It therefore guarantees data that are specific to the motion. This is achieved by simulating typical behaviors that occur in real InSAR data, which represent the first SLC acquisition.

[0024] According to a preferred embodiment, the first and second SLC acquisitions are defined by the following for the generation of the synthetic training dataset: {s1=A0w1s2=s1ρ*ejϕtheo+A0w21−ρ2 where s1 and s2 are the co-registered first and second SLC acquisitions, A0 is a noise-free amplitude, ρ is the given coherence, w1 and w2 are independent circular complex standard Gaussian variables and ϕ theo the noise-free wrapped-phase sum of the zero-Doppler phase (ϕ) zd ) and the Along-Track phase (ϕ at ), which represent the two deformations.

[0025] The noise-free wrapped-phase sum ϕ theo is determined by: ϕtheo=mod2π(ϕzd+ϕat), where the Along-Track phase ϕ at The phase contribution, which is assigned to the deformation in the along-track direction, is defined by: ϕzd=4πλ[sin(θinc)sin(θN)uE+sin(θinc)cos(θN)uN−cos(θinc)uU]

[0026] In this equation, f dc the Doppler center frequency and v g the speed of the SAR system over ground (ground velocity).

[0027] Starting from the estimated interferometric phase, ϕ̂, which is calculated from the pair of first and second SLC acquisitions, s1, s2, by ϕ=∠E{s1s2*} The zero Doppler phase (ϕ) will be zd ) and the Along-Track phase (ϕ αt ) separated from each other, where E{·} is a mathematical expectation, * specifies the complex conjugation operator and ∠(·) denotes an angle operator.

[0028] The model-based and external 3D deformation is provided in an East-North-Up (ENU) coordinate system, as used by the SAR system, from external and model-based sources representing general surface deformation areas in both land ice (e.g., glacial flows) and land-based scenarios (e.g., earthquakes), with the synthetic surface deformations being converted into the zero-Doppler phase ϕ zd and the Along-Track phase ϕ atare induced. The 3D maps of the deformation in the East-North-Up (ENU) coordinate system are referred to as u E , u N and u U The scenarios for land ice can be represented by the ice velocity maps produced, for example, under the Greenland Ice Sheet Monitoring Programme (PROMICE) [6]. In particular, a time series of annual winter ice velocity mosaics based on Sentinel-1 SAR offset tracking was published as part of the project.

[0029] As a result, a training dataset from multiple sources is provided, in which the synthetic surface deformations are incorporated into the zero-Doppler phase ϕ. zd and the Along-Track phase ϕ at can be induced by applying the equations given above.

[0030] No offset refinement, such as enhanced spectral diversity (ESD) and incoherent cross-correlation (ICC), is required to generate the interferometric phase from the synthetic pair. Interferogram generation can be performed burst by burst by applying standard multilooking analysis through a maximum likelihood estimator and using a moving average with a square window Ω. Γ=∑Ω{s1s2*}=|Γ|ejϕ

[0031] The interferometric phase ϕ̂ and coherence ρ̂ can then be retrieved as follows: ϕ^=∠Γ ρ^|Γ|√(∑Ω{|s1|2}∑Ω{|s2|2})

[0032] The next step is a burst merging operation, which involves debursting or mosaicking along the azimuth and creates a continuous image with respect to azimuth time by removing redundant lines between the bursts (overlaps). As a result, the amplitudes and interferometric coherence ρ appear continuous along the azimuth, while the interferometric phase shows the phase jumps at the overlaps if there is along-track motion in the given pixel.

[0033] It should be noted that the desired coherence ρ addresses the noise in interferometric products: the lower its value, the higher the noise contribution in the estimated quantities retrieved in the equations above. Using a ρ-variant dataset is advantageous for noise reduction filtering in the neural network.

[0034] After the dataset is generated, patches are preferably created to define the inputs for the proposed neural network. The extraction of these patches is generally based on the geometric properties of the TOPS system, such as the position of the connection between two merged TOPS bursts and the size of each burst. In particular, the size of the patches is dimensioned to ensure centering at the phase shift while simultaneously guaranteeing the final aggregation of the predicted patches. According to a preferred embodiment, the synthetic training dataset is generated following the step of creating patches with a dimension defined by a selection of the (square) window Ω when multilooking is applied to the interferogram.

[0035] The neural network is based on a supervised training strategy to learn the aforementioned phase separation problem. The unique aspect of this strategy lies in its multitask learning approach, which simultaneously learns two phases to be separated (in the along-track and across-track directions) and the offset in the along-track direction. To achieve this, the neural network is a multitask learning (MTL) convolutional neural network (CNN). Preferably, the MTL CNN is based on the well-known U-Net model, which uses three layers in both the encoder and decoder.

[0036] According to a preferred embodiment, the encoder is followed by a bridge layer and two branches, each comprising a three-layer decoding path with skip connections.

[0037] A first branch initially predicts the real and imaginary parts of the along-track phase and then uses the predicted (reconstructed) along-track phase and the interferometric phase, which are assigned to the real and imaginary parts of a first spur link, to predict the zero-Doppler phase and its real and imaginary parts.

[0038] The normalized Doppler center frequency is used as input for the second branch to apply an RMSE (Root Mean Squared Error) regularization to the along-track phase estimated in the first branch and the same quantity estimated from the along-track deformation in the second branch.

[0039] The predicted (reconstructed) along-track phase is compared with the along-track phase resulting from the output of the second branch to determine the along-track deformation.

[0040] According to a second aspect of the present invention, a method for the computer-aided processing of digital SAR data, consisting of a plurality of interferograms, is proposed. Each interferogram is obtained from a first SLC acquisition of a SAR system at a first time point and a second SLC acquisition of the SAR system at a second time point later than the first time point, as well as a Doppler center frequency, an interferometric phase, and a coherence, all of which are obtained from a respective pair of first and second SLC acquisitions.A zero-Doppler phase is the phase contribution associated with the deformation that causes a zero-Doppler LOS variation, an along-track phase is the phase contribution associated with the deformation in an along-track direction, and an estimate of an associated along-track surface motion is generated using a data-driven model generated by a method according to one or more embodiments of the invention.

[0041] According to a third aspect of the present invention, a device for the computer-aided processing of digital SAR data is proposed. The digital SAR data consist of a plurality of interferograms, each obtained from a first SLC acquisition of the SAR system at a first time point and a second SLC acquisition of the SAR system at a second time point later than the first, as well as a Doppler center frequency, an interferometric phase, and a coherence, all of which are obtained from a respective pair of first and second SLC acquisitions. The device comprises a data-driven model generated by a method according to one or more embodiments of the invention.

[0042] According to a fourth aspect of the present invention, a computer program product is proposed which comprises program code stored on a machine-readable medium for carrying out a method according to one or more embodiments of the invention when the program code is executed on a computer.

[0043] The invention is explained in more detail with reference to the attached figures in the drawing. Fig. Figure 1 shows a general overview of the supervised convolutional neural network for the interferometric phase source separation problem according to the invention. Fig. Figure 2 shows a schematic diagram illustrating patch extraction based on the geometric properties of the TOPS system to be simulated. Fig. Figure 3 shows a schematic diagram illustrating the evolution of the Doppler center frequency f. dcillustrated along the azimuth time t within a single burst. Fig. Figure 4 shows the network architecture of the convolutional neural network according to the invention. Fig. Figure 5 shows a structure of residual blocks used in the proposed architecture. Fig. Figure 6 shows an exemplary test on a real interferometric pair over Greenland, where different rows refer to a different TOPS substrip, while each column indicates a feature. Fig. Figure 7 shows a time-frequency diagram of a TOPS system describing two successive bursts, where a temporal overlap corresponds to a spectral separation.

[0044] Fig. Figure 1 shows a general overview of the training of a data-driven model for processing digital SAR data from a repeat-pass SAR system using the known TOPS acquisition mode for the interferometric phase source separation problem according to the invention. The data-driven model is based on a convolutional neural network (CNN) architecture and a supervised multitask learning strategy specifically designed for interferometric TOPS data.

[0045] With reference to Fig. 1. A general overview of the method according to the invention is given. In the left part of Fig. Figure 1 shows the training process (“Training” 110) of the data-driven model, while the right part of Fig. 1 The processing of real SAR data from the SAR system (“Prediction” 120) using a trained data-driven model is illustrated to determine the desired phase information.

[0046] To train a data-driven model for processing digital SAR data from a repeat-pass SAR system, a synthetic training dataset (TDS) is generated (Block 111). The synthetic training dataset TDS is formed from a variety of interferograms, each consisting of a pair of first and second single-look complex (SLC) acquisitions. This is achieved by inducing a model-based and external 3D deformation to a given first SLC acquisition to obtain the associated second SLC acquisition, which is referenced to a time point after a deformation event. A reverse approach is then applied, specified by a given coherence and a noise-free interferometric phase after the deformation event.The process of generating the synthetic training dataset TDS comprises a first step of generating the multitude of interferograms that produce patches (Block 111G), and a second step of generating (i.e., extracting) patches (Block 111P) from the interferograms of Step 111G. As will be explained later, patch extraction is generally based on the geometric properties of the TOPS system, such as the position of the connection between two merged TOPS bursts and the size of each burst. In particular, the patch size is tailored to ensure centering at the phase shift while also guaranteeing the final aggregation of the predicted patches.

[0047] After generating a large number of patches, the training dataset TDS is generated. The training dataset TDS is divided into a first part 113 (“training dataset”), i.e., data used for training a data-driven model 115, and a second part 114 (“validation dataset”), which is used for validating the quality of the trained data-driven model 115. As is known to those skilled in the art, the data-driven model 115 comprises a model 116 (“CNN”) to be trained and an evaluation part 117 (“evaluation”) for evaluating the quality of the model 116. The model 116 uses a CNN architecture and a supervised multitask learning strategy specifically designed for interferometric TOPS data, as will be explained later.

[0048] Model 116 (hereinafter also referred to as CNN), which was trained with training dataset 113 and validated with validation dataset 114, represents the data-driven model 121 (“Model”) for the computational processing of digital SAR data (“Real Data”) in the “Prediction” part 120 of Fig. 1. The digital SAR data (“Real Data”) are provided in a database 122 to predict specific information (“TOPS AT / ZD Products”) 123 that is suitable for phase source separation. This specific information can be stored in a database 123 for further evaluation.

[0049] The input data transferred to the data-driven model 115 from the training dataset 113 and the validation dataset 114 comprise the input data X train , as they could be retrieved from database 122, and the output data Y train, which correspond to the specific information provided by CNN 121 as the output to be calculated.

[0050] The following section describes in detail the generation of synthetic datasets 111 (111G and 111P).

[0051] For the generation of synthetic 111G datasets, synthetic interferometric acquisition pairs are created by applying a reverse approach. Given the desired coherence and noise-free interferometric phase after an event, the signal model induces the offset in the second acquisition (i.e., after the event). The term "event" is to be understood as a deformation. The deformation could be, for example, an earthquake or a glacial flow. The proposed method guarantees data that are specific to the movement.

[0052] Offset maps, which provide the necessary information for inducing offsets, can be selected from various sources. For example, maps of external surface deformations, such as ice velocity maps produced under the Greenland Ice Sheet Monitoring Programme (PROMICE) [6], and more sophisticated mathematical models, such as the one proposed by Okada to describe seismic events [7], can be used. As a result, a set of synthetic TOPS interferograms is provided for patch-building step 111P, which is responsible for splitting the generated data into the training dataset 113 and the validation dataset 114.

[0053] The problem of generating the training set is solved by simulating typical behaviors that occur in real-world InSAR data. Specifically, the following signal model is used: {s1=A0w1s2=s1ρ*ejϕtheo+A0w21−ρ2 where s1 and s2 are the two SLC acquisitions of the interferometric pair (already coregistered in the geometric sense), A0 is the noise-free amplitude, ρ is the desired coherence, w1 and w2 are independent circular complex standard Gaussian variables, and ϕ is the phase. theo is the noise-free wrapped sum of the two deformations, the zero-Doppler phase ϕ zd and the Along-Track phase ϕ at : ϕtheo=mod2π(ϕzd+ϕat)

[0054] The zero-Doppler phase ϕ zd The phase contribution associated with the deformation that causes a zero-Doppler LOS variation is defined as: the projection of the deformations along the vertical direction and their orthogonal plane onto the perpendicular LOS between target and sensor. ϕzd=4πλ[sin(θinc)sin(θN)uE+sin(θinc)cos(θN)uN−cos(θinc)uU]

[0055] Here, λ is the carrier wavelength, θ incthe angle of incidence, θ N the north angle, which corresponds to a 90° rotation of the TOPS heading angle, while u E , u N and u U The east, north, and upward components of the deformation are respectively (ENU coordinate system used by the SAR system). The along-track phase ϕ at is the phase contribution assigned to the deformation in the along-track direction: ϕat=2πfdcvg[−cos(θN)uE+sin(θN)uN] where f dc the Doppler center frequency and v g The ground speed is.

[0056] The aim is to separate the two quantities in equations (3) and (4) starting from the estimated interferometric phase ϕ̂, which is calculated from the SLC pair defined in equations (1) and (2): ϕ^=∠E{s1s2*} where E{·} is the mathematical expectation value, represents the complex conjugation operator, and ∠(·) denotes the angle operator.

[0057] The 3D maps of the deformation in the East-North-Up (ENU) coordinate system, specified as u E , u N and u U , are generated from the above-mentioned external and model-based sources, which represent general areas of surface deformation in both scenarios with land ice (e.g., glacial rivers) and scenarios with solid earth (e.g., earthquakes).

[0058] The scenarios for the land ice are represented by ice velocity maps produced under the Greenland Ice Sheet Monitoring Programme (PROMICE) [6]. In particular, a time series of annual winter ice velocity mosaics based on Sentinel-1 (S-1) SAR offset tracking was published as part of the project. Each geospatial product covers the period from 1 November 2016 to 28 February of the following year, with a nominal grid spacing of 500 m. The mosaic includes a time-averaged cross-correlation between the individual reflectivity maps and incorporates all 6- and 12-day pairs recorded during this period.

[0059] The differences with respect to the proposed methodology lie in (1) the low resolution and (2) the estimation errors due to the cross-correlation between the amplitude of the SLCs (no interferometry is applied) and the (3) time average, which covers the entire winter season and does not allow for a direct comparison with the actual data. To correct the errors in the products, a Gaussian filter can be applied to the geocoded ice velocity maps using a window size of 11 × 11 pixels. Subsequently, the data can be re-geocoded into the slant-range raster of the SLC. The resulting ice velocity components in meters per day are scaled to an offset that the repeat-pass SAR system, here S-1, can observe by applying the following formula: uENU=vENUΔT where ΔT is the repetition time in days.

[0060] Scenarios for solid Earth layers are treated with a more sophisticated mathematical model proposed by Okada to describe seismic events (see reference [7]). Once the geometric parameters of the fault are determined, the surface deformation is calculated in Cartesian coordinates (XYZ) derived from a Green's function solution for the elastic half-space problem. Since the local SAR geometry is known from the state vector information, the deformation provided by the model can be locally rotated into the radar geometry (ENU) using a suitable rotation matrix (see reference [8]).

[0061] As a result, a training dataset from multiple sources is provided, in which the synthetic surface deformations are incorporated into the zero-Doppler phase ϕ. zd and the Along-Track phase ϕ atcan be induced by applying equation (3) or equation (4).

[0062] For the generation of the interferometric phase ϕ̂ from the synthetic pair, no offset refinement such as Enhanced Spectral Diversity (ESD) and Incoherent Cross-Correlation (ICC) is performed, and the interferogram generation can be carried out directly burst by burst by applying standard multilooking through a maximum likelihood estimator, as proposed in reference [9], and using a moving average with square window Ω: Γ=∑Ω{s1s2*}=|Γ|ejϕ

[0063] The interferometric phase ϕ̂ and coherence ρ̂ are obtained from Γ as follows: ϕ^=∠Γ and: ρ^=|Γ|√(∑Ω{|s1|2}∑Ω{|s2|2})

[0064] Finally, a burst merging operation (debursting or mosaicking along the azimuth) is applied, which produces a continuous image with respect to the azimuth time, as in Fig. Figure 2 shows that redundant lines between bursts B1, B2, B3, B4 (overlaps) are removed. As a result, the amplitudes and interferometric coherence ρ appear continuous along the azimuth (az-time), while the interferometric phase exhibits phase jumps PJ1, PJ2, PJ3 at the overlaps when there is motion along the track in the respective pixel.

[0065] It is important to note that the desired coherence ρ represented in equation (1) addresses the noise in interferometric products: the lower its value, the higher the noise contribution in the estimated quantities obtained in equations (8) and (9). Using a ρ-variant dataset is advantageous for noise reduction filtering in the neural network.

[0066] Next, the generation of patches according to block 111P will be carried out. Fig. As described in section 1, patch formation for the TOPS phase source separation problem is generally system-based, meaning its dimension and selection depend on the geometric properties of the TOPS data. In particular, the selection of patches leads to the selection of a specific window Ω when multilooking is applied to the interferogram.

[0067] For example, in the case of S-1, since the window size is chosen according to the azimuth and ground area resolution of S-1 Interferometric Wide (IW) mode data, 14 m and 3.7 m respectively

[10] , the sampling at the output is set to 100 m, which corresponds to a 7 × 27 averaging window size. Fig. Figure 2 shows the main advantage of using this assumption.

[0068] Fig. Figure 2 shows an interferogram Γ for a substrip 200 of the SAR system, in which redundant lines between bursts B1, B2, B3, B4 (overlaps) have already been removed. The dashed lines labeled MB refer to the midbursts of bursts B1, B2, B3, B4. From the perspective of the substrip (subswath), the division into patches 210 (usually on the left side of Fig. (2 shown) a center overlap. In the present example, two patches 211, 212 selected from the interferogram are shown. The person skilled in the art will understand that the patches 211, 212 (generally: 210) are selected to cover the entire interferogram. The patches 210 are used to determine where phase jumps are located in the interferogram. The patches 210 are arranged such that a corresponding phase jump PJ1, PJ2, PJ3 is located in the center of a corresponding patch 210. In the present embodiment, the phase jump PJ1 is located in the center of patch 211, and the phase jump PJ2 is located in the center of patch 210.

[0069] Since the burst duration in an S-1 IW product is approximately 3 s and the azimuth time interval is approximately Δt = 2 ms (reference

[10] ), the nominal burst size along the azimuth at full resolution is always approximately 1500 lines. After burst removal and multilooking, the effective burst size along the azimuth is nazML, approximately 214 pixels, which is suitable for a patch size of M = 256 and ensures an overlap between the patches.

[0070] The proposed CNN will take into account the real and imaginary parts of the interferometric phase ϕ, which is denoted as ξ Re = Re{e jϕ} or ξ Im = Im{e jϕ} is specified. In addition, a map describing the sensitivity of the first TOPS image along the azimuth direction is considered. This image is specifically enhanced by adding a linear azimuth dependence to each Doppler center of gravity distance variant vector using the vector calculated in

[11] and in Fig. 3 spectral separation Δf was obtained.

[0071] Since the convergence of a CNN is generally faster when the average of each input variable over the training set is close to zero

[12] , the resulting Doppler center-of-mass matrix is ​​normalized to 2.8 kHz, a value not exceeded in S-1 TOPS data, thus guaranteeing an input feature between -1 and 1. In this document, this quantity is referred to as fdcnorm designated.

[0072] The proposed CNN 115 and its phase separation training are described in detail below using the following examples: Fig. 4 described. The CNN 400 in Fig. 4, which CNN 115 in Fig. The approach described in point 1 is based on a supervised training strategy to learn the phase separation problem explained above. The distinctive feature of the proposed strategy lies in how the network is trained using a multi-task learning (MTL) strategy; that is, it simultaneously estimates two phases to be separated (in the along-track and across-track directions) and the offset in the along-track direction.

[0073] The main goal of the proposed MTL-CNN 400 is to separate the along-track phase ϕ. at from the zero Doppler phase ϕ zd to optimize and the associated along-track surface movement u at to appreciate. The MTL-CNN 400 is based on the well-known U-Net model, which uses three layers in both the encoder and the decoder. In the schematic representation of Fig. The CNN comprises an encoder 410 and two branches 420 and 430, each containing a three-layer decoding path with jump links 440 and 450, respectively. The first branch 420 includes a phase encoder 421. The second branch 430 includes an AT deformation decoder 431.

[0074] The in Fig. The four graphic elements used have the following meaning. Entrance OUT output FM Feature Maps CFM Linked Feature Maps REG Regularization Arrow P1 Residual block (RB) Arrow P2 Maxpooling + RB Arrow P3 RB + Upsampling Arrow P4 1×1 fold CSC Copy-Skip Connection SC standard connection |(·)| Modulus ∠(·) angle (·) 2π Wrapping R / I separation Real / Imaginary part T at {·} Along-track (AT) deformation by phase AT Along-Track ZD Zero Doppler

[0075] Both the encoder 410 and the decoders 421 and 431 have three layers, each connected to a residual block (RB) architecture, as shown in Fig. 5 shown.

[0076] At each level, instead of a classic cascade of folds, the following occurs: Fig. Figure 5 shows a residual block architecture 500. Here, the input IN is connected to the output OUT via a short-circuit connection (1 × 1 convolution) to match the dimensions of the output to those of the input. More precisely, the residual block architecture (RB) comprises a first 9x9 convolution block 501, followed by a first batch normalization block 502, a first ReLU activation function 503, a second 9x9 convolution block 504, and a second batch normalization block 505. A third 1x1 convolution block 506 branches from the input IN, and its output is added to the output of the second batch normalization block 505 in an addition block 507. The output of the addition block 507 is forwarded to a second ReLU activation function 508, which provides the output for the next layer.

[0077] Replacing standard convolutional layers with residual layers is a proposed solution to reduce computation time and limit overfitting during training

[13] . The input form is [M × M × C], where M is the patch size and C is the number of channels, i.e., the features selected as input for the proposed neural network. The patch size M is set to 256 and the number of input features C to 4, and C consists of the real and imaginary parts of the complex coefficient ξ associated with the interferometric phase. Re or ξ Im (Phasor), in addition to the coherence ρ and the sensitivity along the azimuth, which are shown in the normalized Doppler center frequency map fdcnorm are hidden. The number of filter sizes at the output is set to a multiple of P = 32 filters for all layers.

[0078] Therefore, in the MTL-CNN architecture of the CNN, one layer was removed from the original U-Net, and a three-layer implementation was used. The reason for this choice lies in the fact that small input patches (256 x 256 pixels) are processed, and that a flatter architecture is easier to train. In fact, it has fewer hyperparameters to optimize and therefore requires fewer input samples (training patches).

[0079] As already mentioned, the encoder 410 is followed by a bridge layer and two branches 420, 430, each comprising a three-layer decoding path with jump links 440, 450: The upper branch 440 first predicts the real and imaginary parts of the along-track phase. {ξazRe,ξazIm} and then uses the reconstructed along-track phase together with the interferometric phase, which corresponds to the real and imaginary parts {ξ Re , ξ Im} is assigned (network inputs) for predicting the zero-Doppler phase and, in particular, the real and imaginary parts {ξzdRe,ξzdIm} used. The same reconstructed along-track phase is compared with the same phase measured at the output of the lower branch 450, i.e., the along-track deformation u. at .

[0080] In particular, it is important to consider the normalized Doppler center frequency. fdcnorm The input used for the network allows RMSE regularization to be applied to the along-track phase estimated in the upper branch 440 and to the same quantity estimated from the along-track deformation in the lower branch 450. In this way, the network can determine the weights on both branches 440 and 450 to also satisfy the relationship between the outputs. Furthermore, the proposed neural network guarantees regularization of the unit modulus in the along-track phase and regularization of the continuity of the zero-Doppler phase.

[0081] The method described above was successfully demonstrated by the inventors through the execution of a series of simulations.

[0082] In Fig. Figure 6 shows the results for three sub-strips iw1, iw2, iw3 of a real interferogram. Each row in Fig. Column 6 is assigned to a different substrip. Columns (a) and (b) represent the coherence ρ and interferometric phase ϕ of the S-1 pair, respectively; these quantities are used indirectly as input for CNN 121 and 500, respectively. Columns (c) and (d) represent the reconstructed zero-Doppler phase ϕ. zd and the Along-Track phase ϕ at , while column (e) is the along-track deformation u at in meters, which is assigned to the latter size.

[0083] First, let's consider the first line (iw1) of Fig. 6. Thus, it can be observed that the trained network is able to denoise the interferometric phase when the coherence is low. This is motivated by the random coherence variation applied to the trained patches. An additional check to evaluate the validity of the deformation along the track can be performed by comparing the enveloping phase reconstructed from the estimated deformation with the phase assigned to the differential interferogram within the overlapping region. The former can be derived from equation (4) as: ϕatovli=mod2π(ϕazF,ovli−ϕazB,ovli)=mod2π(2πΔfdcovlivguat) where ϕazF,ovli and ϕazB,ovli The along-track phases are those that correspond to the overlap area of ​​two successive bursts (forward and backward). Equation (9) can be solved by considering the scheme in Fig. 7 can be justified, which represents the time-frequency diagram of two successive bursts of the same sub-strip in a specific area.

[0084] In particular, it shows Fig. 7, that Δfdcovli the difference between the forward and backward Doppler center frequencies at the overlap ovl i is and relates to the previously set maximum Doppler center frequency of 2.8 kHz fdcmax max can be approximated to a slightly varying term α between -1 and -2.

[0085] The proposed approach enables the use of existing systems in TOPS mode to estimate deformation along the flight path, which is typically the direction in which SAR systems in LEO orbit are less sensitive. Furthermore, this technique paves the way for future systems, such as ESA's upcoming Harmony mission, to better estimate deformation maps along the flight path. The technique can be used to monitor glaciers in the Arctic, e.g., over Greenland, using Sentinel-1 data; to measure large deformations over ice land surfaces, e.g., glacial rivers; to measure large deformations over solid Earth, e.g., earthquakes; and, as a generalization of the ESD technique, to mitigate the phase jumps present in TOPS interferograms when motion is present along the track in the scene. References [1] R. Scheiber, M. Jäger, P. Prats-Iraola, F. De Zan, and D. Geudtner, „Speckle Tracking and Interferometric Processing of TerraSAR-X TOPS Data for Mapping Nonstationary Scenarios,“ IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 8, no. 4, pp. 1709-1720, 2015. [2] P. Prats-Iraola, M. Nannini, N. Yague-Martinez, M. Pinheiro, J.-S. Kim, F. Vecchioli, F. Minati, M. Costantini, S. Borgstrom, P. De Martino, V. Siniscalchi, M. Foumelis, and Y.-L. Desnos, „Interferometric investigations with the Sentinel-1 constellation,“ in 2017 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), 2017, pp. 5537-5540. [3] R. Grandin, E. Klein, M. Metois, and C. Vigny, „Three-dimensional displacement field of the 2015 Mw8.3 Illapel earthquake (Chile) from across- and along-track Sentinel-1 TOPS interferometry,“ Geophysical Research Letters, vol. 43, no. 6, pp. 2552-2561, 2016. [4] J. K. Andersen, A. Kusk, J. P. Merryman Boncori, C. S. Hvidberg, and A. Grinsted, „Improved Ice Velocity Measurements with Sentinel-1 TOPS Interferometry,“ Remote Sensing, vol. 12, no. 12, 2020. [5] A. Kusk, J. K. Andersen, and J. P. Merryman Boncori, „Sentinel-1 TOPS Interferometry for Ice Flow Mapping,“ in EUSAR 2021; 13th European Conference on Synthetic Aperture Radar, 2021, pp. 1-5. [6] A. Solgaard, A. Kusk, J. P. Merryman Boncori, J. Dall, K. D. Mankoff, A. P. Ahlstrøm, S. B. Andersen, M. Citterio, N. B. Karlsson, K. K. Kjeldsen, N. J. Korsgaard, S. H. Larsen, and R. S. Fausto, „Greenland ice velocity maps from the PROMICE project,“ Earth System Science Data, vol. 13, no. 7, pp. 3491-3512, 2021. [7] Y. Okada, „Surface deformation due to shear and tensile faults in a halfspace,“ Bulletin of the Seismological Society of America, vol. 75, no. 4, pp. 1135-1154, Aug. 1985. [8] M. Grewal, L. Weill, and A. Andrews, Global Positioning Systems, Inertial Navigation, and Integration: Second Edition. John Wiley Sons, Ltd, Feb. 2007, ch. Appendix C: Coordinate Transformations, pp. 456-501. [9] J.-S. Lee, K. Papathanassiou, T. Ainsworth, M. Grunes, and A. Reigber, „A new technique for noise filtering of SAR interferometric phase images,“ IEEE Transactions on Geoscience and Remote Sensing, vol. 36, no. 5, pp. 1456-1465, 1998.

[10] M. Bourbigot, H. Johnsen, and R. Piantanida, „Sentinel Product Definition,“ European Space Agency (ESA), Tech. Rep. S1-RSMDA-52-7440, March 2016. [Online]. Available: https: / / sentinel.esa.int / documents / 247904 / 1877131 / Sentinel-1-Product-Definition

[11] N. Yagüe-Martinez, P. Prats-Iraola, F. Rodriguez Gonzalez, R. Brcic, R. Shau, D. Geudtner, M. Eineder, and R. Bamler, „Interferometric Processing of Sentinel-1 TOPS Data,“ IEEE Transactions on Geoscience and Remote Sensing, vol. 54, no. 4, pp. 2220-2234, 2016.

[12] Y. A. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller, Neural Networks: Tricks of the Trade: Second Edition. Springer Berlin Heidelberg, 2012, ch. Efficient BackProp, pp. 9-48.

[13] K. He, X. Zhang, S. Ren, and J. Sun, „Deep residual learning for image recognition,“ in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016, pp. 770-778.

[14] ANANTRASIRICHAI N. [et al.]: A deep learning approach to detecting volcano deformation from satellite imagery using synthetic datasets. In: Remote Sensing of Environment, Vol. 230, 2019. S. 1-11. - ISSN 0034-4257

Claims

[1] A method for the computer-aided generation of a data-based model for the computer-aided processing of digital SAR data of a repeat-pass SAR system configured to determine surface deformations, wherein the digital SAR data consists of a plurality of interferograms (Γ) each obtained from a first SLC acquisition (s1) of the SAR system at a first time and a second SLC acquisition (s2) of the SAR system at a second time later than the first time, and having a Doppler center of gravity frequency (f dc ), an interferometric phase (ϕ̂) and a coherence (ρ̂), all obtained from a respective pair of first and second SLC acquisitions (s1, s2), and wherein the data-driven model is arranged to provide a zero-Doppler phase (ϕ zd), which represents the phase contribution associated with the deformation causing a zero-Doppler LOS variation, an along-track phase (ϕ αt ) representing the phase contribution associated with deformation in an along-track direction and an estimate of an associated along-track surface motion (u at ), whereby the generation of the data-driven model includes the following steps: - generating a synthetic training data set (TDS) comprising a plurality of interferograms (Γ), each formed by a pair of a first and a second SLC acquisition (s1, s2), by inducing a model-based and external 3D deformation into a given first SLC acquisition (s1) to obtain the associated second SLC acquisition (s2) related to a time point after a deformation event, and by applying an inverse approach defined by a given coherence (ρ) and a noise-free interferometric phase (ϕ theo ) is specified after the deformation event; - Training a neural network (NN) with the training data set (TDS), wherein the neural network (NN) has an input for inputting the plurality of interferograms (Γ) and an output for outputting the along-track phase (ϕ at ), the zero-Doppler phase (ϕ zd ) and the associated along-track surface movement (uat ) for the plurality of interferograms (Γ), wherein the trained neural network (NN) represents the generated data-based model. [2] Method according to claim 1, wherein for generating the synthetic training data set (TDS) the first and second SLC acquisitions (s1, s2) are defined by: {s1=A0w1s2=s1ρ*ejϕtheo+A0w21−ρ2 where s1 and s2 are the co-registered first and second SLC acquisitions, A0 is a noise-free amplitude, ρ is the given coherence, w1 and w2 are independent circular complex standard Gaussian variables and ϕ theo the noise-free wrapped-phase sum of the zero-Doppler phase (ϕ zd ) and the along-track phase (ϕ at ), which represent two deformations. [3] The method of claim 2, wherein the noise-free wrapped phase sum (ϕ theo ) is determined by: ϕtheo=mod2π(ϕzd+ϕat). [4] Method according to claim 3, wherein the along-track phase (ϕ at ) is the phase contribution associated with the deformation in the along-track direction, defined by: ϕzd=4πλ[sin(θinc)sin(θN)uE+sin(θN)cos(θN)uN−cos(θinc)uU] where f dc is the Doppler center frequency and v g a ground speed of the SAR system. [5] Method according to claim 4, wherein starting from the estimated interferometric phase (ϕ̂) calculated from the pair of first and second SLC acquisitions (s1, s2), by ϕ^=∠E{s1s2*} the zero-Doppler phase (ϕ zd ) and the along-track phase (ϕ at ) are separated from each other, whereby E{·} is a mathematical expectation, specifies the complex conjugation operator and ∠(·) denotes an angle operator. [6] A method according to any one of the preceding claims, wherein the model-based and external 3D deformation in an East-North-Up (ENU) coordinate system as used by the SAR system are provided by external and model-based sources representing common surface deformation areas in both land ice and continent scenarios, wherein the synthetic surface deformations are converted to the zero-Doppler phase (ϕ zd ) and the along-track phase (ϕ αt ) can be induced. [7] A method according to any one of the preceding claims, wherein the generation of the interferograms (Γ) is carried out burst by burst by applying standard multilooking by a maximum likelihood estimator and using a moving average with a square window Ω. [8] A method according to claim 7, wherein a burst merging operation is performed which creates a continuous image with respect to azimuth time by removing redundant lines between bursts (overlaps). [9] A method according to claim 7 or 8, wherein the step of generating the synthetic training data set (TDS) is followed by the generation of patches having a dimension defined by a selection of the window (Ω) when multilooking is applied to the interferogram (Γ). [10] Method according to one of the preceding claims, wherein the neural network (NN) is a convolutional neural network (500) with multitask learning, MTL-CNN. [11] The method of claim 10, wherein the MTL-CNN (500) is based on the U-net model, which uses three layers in both the encoder (510) and the decoder (521, 531). [12] Method according to claim 10 or 11, wherein the encoder (510) is followed by a bridge layer (514) and two branches, each comprising a three-layer decoding path (520, 530) with jump connections (540, 550). [13] The method of claim 12, wherein a first branch (520) first calculates the real and imaginary parts of the along-track phase ({ξzdRe,ξzdIm}) predicted and then the predicted along-track phase and the interferometric phase corresponding to the real and imaginary parts ({ξzdRe,ξzdIm}) a first jump connection (540) are used to determine the zero-Doppler phase and its real and imaginary parts ({ξzdRe,ξzdIm}) to predict. [14] A method according to claim 12 or 13, wherein the normalized Doppler centroid frequency (fdcnorm) is used as input to the second branch (530) to apply an RMSE regularization to the along-track phase estimated in the first branch (520) and to the same quantity estimated from the along-track deformation in the second branch (530). [15] Method according to one of claims 12 to 14, wherein the predicted along-track phase is compared with the same phase resulting from the output of the second branch (530) to determine the along-track deformation (u at ) to determine. [16] Method for the computer-aided processing of digital SAR data consisting of a plurality of interferograms (Γ), each of which is obtained from a first SLC acquisition (s1) of a SAR system at a first time and a second SLC acquisition (s2) of the SAR system at a second time later than the first time, and a Doppler centroid frequency (f dc), an interferometric phase (ϕ̂) and a coherence (ρ̂), all obtained from a respective pair of first and second SLC acquisitions (s1, s2), where a zero-Doppler phase (ϕ zd ) is the phase contribution associated with the deformation causing a zero-Doppler LOS variation, an along-track phase (ϕ at ) is the phase contribution associated with deformation in an along-track direction, and an estimate of an associated along-track surface motion (u at ) is generated with a data-driven model generated with a method according to one of claims 1 to 15. [17] Apparatus for the computer-aided processing of digital SAR data consisting of a plurality of interferograms (Γ), each of which is obtained from a first SLC acquisition (s1) of a SAR system at a first time and a second SLC acquisition (s2) of the SAR system at a second time later than the first time, and a Doppler centroid frequency (f dc ), an interferometric phase (ϕ̂) and a coherence (ρ̂), all obtained from a respective pair of first and second SLC acquisitions (s1, s2), the apparatus comprising a data-driven model generated by a method according to any one of claims 1 to 15. [18] A computer program product comprising program code stored on a machine-readable medium for carrying out a method according to any one of claims 1 to 15 when the program code is executed on a computer.

Images