Computer-implemented method for calibrating a pair of digital elevation maps
The method calibrates digital height maps using SAR image time series to identify stable points and fit a two-dimensional curve, addressing errors in existing methods by compensating for residual slopes and eliminating the need for external measurements, thus enhancing terrain change monitoring accuracy.
Patent Information
- Application Number
- DE102023123310
- Authority / Receiving Office
- DE · DE
- Patent Type
- Patents
- Current Assignee / Owner
- Filing Date
- 2023-08-30
- Publication Date
- 2024-12-19
- Estimated Expiration
- 2043-08-30
AI Technical Summary
Existing methods for calibrating digital height maps in synthetic aperture radar (SAR) require expensive external reference measurements and fail to account for residual slopes, leading to errors of several meters.
A computer-implemented method that calibrates digital height maps using a time series of SAR images to identify candidate calibration points based on amplitude dispersion measures, fitting a two-dimensional curve to height value differences, and correcting the maps without external measurements.
Achieves accurate calibration of digital height maps, reducing errors and enabling precise monitoring of terrain changes by compensating for residual slopes, without the need for costly external reference measurements.
Smart Images
Figure 00000000_0000_ABST
Abstract
Description
[0001] The invention relates to a computer-implemented method for calibrating a pair of digital elevation maps, as well as to a corresponding computer program product and a corresponding computer program.
[0002] Synthetic aperture radar (commonly abbreviated to SAR) is a common technique for remotely sensing the surface of a celestial body, such as Earth. A radar system on a platform (such as a satellite or an aircraft) flying over the surface of the celestial body transmits radar pulses and measures the radio echoes from the pulses scattered on the surface of the celestial body. Using well-known SAR processing techniques, SAR images can be derived from the measured radar echoes. These images include an amplitude value and a phase value for a large number of pixels corresponding to points on the surface of the celestial body.
[0003] Digital elevation maps can be derived from pairs of SAR images taken of the same terrain by analyzing the phase differences between the images. A corresponding digital elevation map provides elevation values for a large number of pixels at corresponding pixel positions, representing points on the surface of the celestial body.
[0004] To analyze changes in terrain on Earth or another celestial body, the technique of DEM differencing (DEM = Digital Elevation Map) is used. Based on this technique, a pair of digital elevation maps taken at different times from the same terrain are compared to determine elevation differences between these digital elevation maps. This allows topographic changes to be monitored and quantified. For example, the mass balance of glaciers can be monitored, or the amount of lava ejected by a volcano can be quantified.
[0005] The digital elevation maps used for DEM differentiation are derived from various SAR images. Therefore, precise mutual calibration of the digital elevation maps is a crucial step. In particular, residual horizontal shifts and vertical offsets and tilts must be compensated for, which, if neglected, can lead to errors on the order of several meters.
[0006] Mutual calibration of digital elevation maps is typically achieved by performing precise spatial co-registration of the geocoded digital elevation maps and determining reliable reference calibration points through external measurements, such as GPS or LiDAR measurements. Document [1] discloses the calibration of digital elevation maps using laser measurements as external reference calibration points.
[0007] Obtaining external reference measurements for calibrating digital elevation maps is often expensive and problematic, especially in remote or difficult-to-access areas. Furthermore, conventional calibration of digital elevation maps only compensates for an absolute offset between a map pair, without considering residual tilts that may be caused by limited accuracy in estimating the SAR system's geometric acquisition parameters.
[0008] Document [2] discloses a method for detecting permanent scatterers in a time series of SAR images. Document [3] describes the use of permanent scatterers in SAR images to improve the creation of a single digital elevation map.
[0009] Document [5] discloses a method for the mutual calibration of a pair of digital elevation maps. In the method, calibration points within the digital elevation maps are determined using a time series of SAR images. Among other things, a dispersion index is used, which is the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position across the SAR images of the time series divided by the mean of this statistical distribution. Pixel positions with a dispersion index below a threshold value are identified as candidates for calibration points. A calibration curve is fitted to the elevation value differences between the two digital elevation maps at the calibration points, using which the elevation values of one of the digital elevation maps are corrected.
[0010] It is an object of the invention to provide a computer-implemented method for calibrating a pair of digital elevation maps without the need for external measurements at reference calibration points.
[0011] This object is achieved by the method according to claim 1 and the computer program product according to claim 13 or the computer program according to claim 14. Preferred embodiments of the invention are defined in the dependent claims.
[0012] The invention provides a computer-implemented method for calibrating a pair of digital elevation maps. The digital elevation maps of the pair indicate the elevation of a predetermined terrain on a surface of a celestial body. Preferably, the celestial body is a planet, and the planet is preferably Earth. The digital elevation maps each include elevation values for a plurality of pixels at pixel positions representing corresponding points within the predetermined terrain. One digital elevation map of the pair indicates the elevation of the predetermined terrain at an earlier time than the other digital elevation map of the pair.
[0013] Each digital elevation map of the pair is linked (associated) with another SAR image pair obtained from measurements of a first radar sensor system, with each digital elevation map derived from phase differences between the associated SAR image pair. The structure of SAR images is known to those skilled in the art. A SAR image is characterized by pixels at a plurality of pixel positions, each pixel characterized by an amplitude and phase value.
[0014] The method according to the invention processes a time series of SAR images of the predetermined terrain obtained from measurements by a second radar system, wherein each SAR image of the time series comprises pixels at the same pixel positions as the digital elevation maps and as the SAR images obtained from measurements by the first radar system. The first radar system and the second radar system can refer to the same radar system. This means that the first radar system can be identical to the second radar system. However, the first radar system can also be a different radar system than the second radar system.
[0015] The first radar system can be a monostatic radar system with one or more radar sensors, where each radar sensor transmits its own radar pulses and measures exclusively the radar returns of its own radar pulses. However, the first radar system can also be a bistatic radar system with two radar sensors flying at a predetermined distance from each other, where only one radar sensor transmits radar pulses and both radar sensors measure the radar returns of the radar pulses. Analogously, the second radar system can be a monostatic radar system or a bistatic radar system. The detailed description describes an embodiment in which the first radar system is a bistatic SAR system on satellites of the TanDEM-X mission, while the second radar system is a monostatic SAR system on satellites of the Sentinel-1 mission.
[0016] The method according to the invention comprises the steps a) to e) described below. In step a), candidate calibration points are determined from the time series of SAR images based on (amplitude) dispersion measures of the respective pixel positions within the SAR images of the time series. The dispersion measure of a respective pixel position describes the statistical scatter (dispersion) of the amplitude values of the pixels at the respective pixel positions across the SAR images of the time series, with each candidate calibration point referring to a pixel position that has a dispersion measure below a predetermined threshold value. The predetermined threshold value is calculated according to the invention based on a calculation procedure that is explained in more detail below. The dispersion measure can be defined in various ways. It is characterized in that a value of the dispersion measure increases with increasing statistical scatter of the amplitude values.
[0017] In step b), pixel positions are determined as calibration points within the pair's digital height maps based on the candidate calibration points, where each calibration point corresponds to a pixel position of a candidate calibration point. In other words, each calibration point is located at a position of a candidate calibration point. However, a candidate calibration point is not necessarily a calibration point. In other words, a filter can be used to extract calibration points from the candidate calibration point.
[0018] In step c), elevation differences between the pair's digital elevation maps at the calibration points are determined, and in step d), a two-dimensional curve (i.e., a two-dimensional function) is fitted to the elevation differences at the calibration points. The two-dimensional curve describes the elevation differences over the pixel positions representing corresponding points within the predetermined terrain. This means that the two-dimensional curve describes the remaining elevation differences between the two pair's digital elevation maps, thus forming a mutual calibration surface. According to step c), a two-dimensional curve is interpolated based on the calibration points. Any known fitting or interpolation technique can be used to determine the two-dimensional curve.For example, linear interpolation or interpolation based on high-order polynomials can be used to obtain the two-dimensional curve.
[0019] In step e), the elevation values of one digital elevation map in the pair are corrected using the elevation value differences of the two-dimensional curve. The corrected digital elevation map and the other (uncorrected) digital elevation map form a pair of digital elevation maps that are correctly calibrated to each other thanks to the determination of a two-dimensional curve based on appropriate calibration points. The mutually calibrated digital elevation maps can then be used for DEM differencing to obtain accurate information about the terrain changes that have occurred in the time period between the two digital elevation map acquisitions.
[0020] The computational procedure for calculating the predetermined threshold for the statistical dispersion measure is based on a first and a second time series of SAR images of a predefined terrain on the surface of the celestial body. The first time series was obtained from measurements of the first radar sensor system, and the second time series was obtained from measurements of the second radar sensor system. All SAR images of the first and second time series include pixels at the same pixel positions, representing corresponding points within the predefined terrain. The predefined terrain does not have to be the same as the predetermined terrain covered by the digital elevation maps.
[0021] The calculation procedure is particularly advantageous when the first radar is a different radar than the second radar, so that the threshold applied to the dispersion measure derived from the SAR images of the second radar is appropriately selected for the pair of digital elevation maps based on measurements of the first radar.
[0022] In a step i) of the calculation procedure, first dispersion measures are determined, which are dispersion measures of respective pixel positions within the SAR images of the first time series, wherein the dispersion measure of a respective pixel position describes the statistical scatter of amplitude values of the pixels at the respective pixel position across the SAR images of the first time series.
[0023] In a step ii) of the calculation procedure, second dispersion measures are determined, which are dispersion measures of respective pixel positions within the SAR images of the second time series, wherein the dispersion measure of a respective pixel position describes the statistical scatter of amplitude values of the pixels at the respective pixel position across the SAR images of the second time series.
[0024] In a step iii) of the calculation procedure, the predetermined threshold is determined based on a predetermined limit value for the first dispersion measures by selecting the predetermined threshold such that both the pixel positions within the respective SAR images of the first time series that have a first dispersion measure above the predetermined limit value and the corresponding pixel positions of the respective SAR images of the second time series that have a second dispersion measure below the predetermined threshold value are a predetermined percentage of 10% or less.
[0025] The invention is based on the finding that permanent scatterer candidates based on a statistical dispersion measure of a time series of SAR images are well suited to become calibration points for a mutual calibration of digital elevation maps.
[0026] In a particularly preferred embodiment, the dispersion measure of a respective pixel position within the SAR images of the time series represents or depends on a dispersion index, which is the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position across the SAR images of the time series divided by the mean of this statistical distribution. The dispersion index, as defined above, is also used in the method described in document [2].
[0027] In further preferred embodiments, suitable techniques are used to reduce the number of candidate calibration points determined in step a). In one embodiment, if a pixel position in the SAR images of the time series is determined to be affected by geometric distortions, this pixel position is classified as a pixel position that is not accepted as a candidate calibration point. Techniques for determining pixel positions as affected by geometric distortions, for example, by layover or shading, are known and will not be described in detail here.
[0028] In a further embodiment, in the event that a pixel position within the SAR images of the time series is considered to be a point within the predetermined
[0029] If terrain covered by volumetric material is determined, this pixel position is classified as a pixel position that is not permitted as a candidate calibration point. Volumetric material has the property that radar waves can penetrate the volumetric material. The classification of terrain as covered by volumetric material is within the expertise of a person skilled in the art. Volumetric material is, for example, snow, ice, vegetation, sand, and the like. In order to determine a pixel position as a point covered by volumetric material, a suitable reference map can be used. This embodiment is based on the recognition that areas covered with snow or ice are not stable over time and are therefore not suitable as calibration points.
[0030] In another embodiment, if a pixel position within the SAR images of the time series is determined to be a point within the predetermined terrain that has a slope greater than a predetermined slope threshold, this pixel position is classified as a pixel position that is not accepted as a candidate calibration point. This embodiment prevents candidate calibration points from being inaccurately acquired.
[0031] Furthermore, in step b), suitable filtering techniques can be used to determine calibration points from the candidate calibration points. In one embodiment, if a pixel position in the digital height maps is determined to be affected by geometric distortions, this pixel position is classified as a pixel position that is not accepted as a calibration point. As mentioned above, methods for determining pixel positions as affected by geometric distortions are known to those skilled in the art.
[0032] In another embodiment, if the interferometric coherence for a pixel position of a SAR image pair associated with one of the pair's digital elevation maps is below a coherence threshold, that pixel position within the digital elevation maps is classified as a pixel position that is not allowed as a calibration point. Interferometric coherence is a known quantity that can be estimated pixel by pixel with respect to image pairs. The definition of interferometric coherence can be found in the detailed description.
[0033] In a further preferred embodiment, if the signal-to-noise ratio of a pixel at a pixel position within at least one SAR image of the SAR image pairs associated with the digital elevation maps of the pair is below a signal-to-noise ratio threshold, that pixel position within the digital elevation maps is classified as a pixel position that is not accepted as a calibration point.
[0034] The first dispersion measures and the second dispersion measures can correspond to the dispersion indices defined above. This means that the first dispersion measure of a respective pixel position within the SAR images of the first time series can represent or depend on a dispersion index that is the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position across the SAR images of the first time series divided by the mean of this statistical distribution. Analogously, the second dispersion measure of a respective pixel position within the SAR images of the second time series can represent or depend on a dispersion index that is the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position across the SAR images of the second time series divided by the mean of this statistical distribution.
[0035] In a preferred embodiment, the above predetermined limit value has a value between 0.15 and 0.3, preferably 0.25, and / or the above predetermined percentage has a value of 5% or less, preferably 1%.
[0036] In addition to the above method, the invention relates to a computer program product with program code stored on a machine-readable carrier for carrying out a method according to the invention or one or more preferred embodiments thereof when the program code is executed on a computer.
[0037] Furthermore, the invention relates to a computer program with program code for carrying out a method according to the invention or one or more preferred embodiments thereof when the program code is executed on a computer.
[0038] In the following, a second, third and fourth aspect are explained, which are independent of the aspects of the invention described above.
[0039] The second aspect relates to a method for calibrating digital elevation models (DEMs) of an Earth's surface (ground surface of the Earth), the method comprising: (a) obtaining a pair of co-registered DEMs corresponding to a selected area of the Earth's surface, each DEM of the pair of co-registered DEMs being generated from single-pass synthetic aperture radar (SAR) data; (b) obtaining multi-pass time series SAR data relating to the selected area of the Earth's surface; (c) selecting a plurality of calibration tie points (TPs) representing a portion of the selected area of the Earth's surface using both the single-pass SAR data and the multiple-pass time series SAR data; and (d) generating a pair of mutually calibrated DEMs corresponding to the selected area of the Earth's surface using both the pair of co-registered DEMs and the selected plurality of calibration TPs.
[0040] A preferred embodiment of the second aspect relates to a method in which (c) comprises: (c1) selecting candidate calibration TPs from the multi-overflight time series data; and (c2) Selecting the calibration TPs from the pair of co-registered DEMs using the candidate calibration TPs.
[0041] A preferred embodiment of the second aspect relates to a method in which (c) comprises: (c1) filtering pixels from the multi-pass time series data to provide a plurality of candidate calibration TPs; (c2) determining the candidate calibration TPs in the pair of co-registered DEMs; and (c3) filtering at least some of a plurality of pixels forming the pair of co-registered DEMs corresponding to the determined candidate calibration TPs, wherein the remaining unfiltered candidate calibration TPs comprise the selected calibration TPs.
[0042] A preferred embodiment of the second aspect relates to a method wherein (c) comprises applying a predefined dispersion index amplitude threshold to pixels from the pair of co-registered DEMs.
[0043] A preferred embodiment of the second aspect relates to a method wherein (c) comprises selecting as calibration TPs only those calibration TPs that satisfy both a multi-pass dispersion index amplitude threshold and a single-pass dispersion index amplitude threshold. Preferably, the single-pass dispersion index amplitude threshold differs from the multiple-pass dispersion index amplitude threshold.
[0044] A preferred embodiment of the second aspect relates to a method in which (d) comprises: (d1) generating a calibration surface; and (d2) applying the calibration surface to one DEM of the pair of co-registered DEMs to produce a calibrated DEM co-registered with the pair of co-registered DEMs, wherein the pair of mutually calibrated DEMs comprises the other DEM of the pair of co-registered DEMs and the calibrated DEM.
[0045] A preferred embodiment of the second aspect relates to a method in which (d) comprises determining the difference between the pair of co-registered DEMs over the calibration TPs selected in (c).
[0046] A preferred embodiment of the second aspect relates to a method in which both the single-pass SAR data and the multiple-pass time series SAR data comprise interferometric synthetic aperture radar (InSAR) data.
[0047] A preferred embodiment of the second aspect relates to a method in which the single-pass SAR data is obtained from a first original data source and the multiple-pass time series SAR data is obtained from a second original data source that is different from the first original data source.
[0048] A preferred embodiment of the second aspect relates to a method further comprising: (e) Generating a calibrated DEM of difference (DoD) by determining the difference between the mutually calibrated DEMs.
[0049] The third aspect relates to a method for calibrating digital elevation models (DEMs) of an Earth's surface, the method comprising: (a) obtaining a pair of co-registered DEMs corresponding to a selected area of the Earth's surface, each DEM of the pair of co-registered DEMs being generated from first synthetic aperture radar (SAR) data corresponding to a first radar band; (b) obtaining second SAR data corresponding to a second radar band different from the first radar band, the second SAR data relating to the selected area of the Earth's surface; (c) selecting a plurality of calibration tie points (TPs) representing a portion of the selected area of the Earth's surface using both the first SAR data and the second SAR data; and (d) generating a pair of mutually calibrated DEMs corresponding to the selected area of the Earth's surface using both the pair of co-registered DEMs and the selected plurality of calibration TPs.
[0050] A preferred embodiment of the third aspect relates to a method in which the first SAR data corresponds to a first radar band and the second SAR data corresponds to a second radar band different from the first radar band.
[0051] A preferred embodiment of the third aspect relates to a method in which the first SAR data comprises X-band radar data and the second SAR data comprises C-band radar data.
[0052] A preferred embodiment of the third aspect relates to a method in which both the first SAR data and the second SAR data comprise interferometric synthetic aperture radar (InSAR) data.
[0053] A preferred embodiment of the third aspect relates to a method wherein (d) comprises determining the difference between the pair of co-registered DEMs over the calibration TPs selected in (c).
[0054] A preferred embodiment of the third aspect relates to a method in which (c) comprises: (c1) selecting candidate calibration TPs from the second SAR data; and (c2) Selecting the calibration TPs from the pair of co-registered DEMs using the candidate calibration TPs.
[0055] The fourth aspect relates to a system for calibrating digital elevation models (DEMs) of an Earth's surface, the system comprising: a processor; a non-transitory memory; and one or more applications stored in non-transitory memory which, when executed by the processor: obtain a pair of co-registered DEMs corresponding to a selected area of the Earth's surface, each of the pair of co-registered DEMs being generated from single-pass synthetic aperture radar (SAR) data; Obtain multi-flyover time series SAR data covering the selected area of the Earth's surface; select a plurality of calibration tie points (TPs) representing a portion of the selected area of the Earth's surface using both the single-pass SAR data and the multi-pass time series SAR data; and generate a pair of mutually calibrated DEMs corresponding to the selected area of the Earth's surface using both the pair of co-registered DEMs and the selected plurality of calibration TPs.
[0056] A preferred embodiment of the fourth aspect relates to a system wherein when the one or more applications stored in the non-transitory memory are executed by the processor: determine the difference between the pair of co-registered DEMs over the selected calibration TPs to generate the pair of mutually calibrated DEMs.
[0057] A preferred embodiment of the fourth aspect relates to a system wherein the one or more applications stored in the non-transitory memory, when executed by the processor: select the candidate calibration TPs from the multiple-overflight time series data; and select the calibration TPs from the pair of co-registered DEMs using the candidate calibration TPs.
[0058] In the following, embodiments of the invention are described in detail with reference to the accompanying drawings, in which: Fig. 1 shows a flowchart illustrating the steps of an embodiment of a method according to the invention; Fig. 2 shows a flowchart illustrating the determination of the candidate calibration points as in step ST1 of Fig. 1 indicated, illustrated; Fig. 3 is a diagram illustrating a method for determining the value of the Fig. 1 threshold used.
[0059] In the following, the invention is described with reference to a mutual calibration of digital single-pass InSAR elevation maps (InSAR = interferometric SAR). The digital elevation maps are obtained using the bistatic SAR of the TanDEM-X mission (see document [4]). In this bistatic SAR system, two satellites fly in close formation over the Earth's surface. Each satellite has a radar sensor, with the radar sensor of one satellite transmitting radar pulses, and the radar sensors of both satellites receiving the radar returns of these radar pulses after they are reflected off the Earth's surface. Thus, in a single pass of the satellites over the Earth's surface, SAR image pairs are obtained at corresponding times.
[0060] Using a known procedure, a digital elevation map can be derived from a corresponding image pair. The digital elevation map includes elevation values for a plurality of pixels at pixel positions representing corresponding points on the Earth within the terrain captured by the SAR images. In the embodiment described here, a pair of digital elevation maps captured at different times of the same terrain on the Earth's surface are mutually calibrated. The mutually calibrated digital elevation maps enable so-called DEM differencing, in which the digital elevation maps are compared to determine changes in the terrain, such as glacial melting or deforestation.
[0061] As mentioned above, the digital elevation maps in the embodiment described here are derived from a bistatic SAR. However, the digital elevation maps can also be obtained from another SAR system, for example, a monostatic SAR system that flies over the same terrain on Earth at regular intervals. A digital elevation map can be derived based on a pair of SAR images acquired of the corresponding terrain during different overflights of the monostatic SAR system.
[0062] Fig. 1 shows a flowchart of the embodiment described below. The method processes a pair of digital elevation maps DEM1 and DEM2. Both digital elevation maps cover the same terrain TR on Earth. The digital elevation map DEM1 includes elevation values hv1 for a plurality of pixels representing the terrain TR, and the digital elevation map DEM2 includes elevation values hv2 for the plurality of pixels representing the terrain TR at a different time than the map DEM1. The two digital elevation maps DEM1 and DEM2 are each derived from a pair IP1 and IP2 of SAR images IM obtained by a first radar installation RE1, the bistatic SAR of the TanDEM-X mission.
[0063] Each of the SAR images IM is characterized by a plurality of pixels covering the terrain TR, with each pixel having an amplitude value a at a phase value ph. The digital elevation map DEM1 is derived from the image pair IP1 of the SAR images IM, acquired simultaneously at a corresponding time by the bistatic SAR, and the digital elevation map DEM2 is derived from the other image pair IP2 of the SAR images IM, acquired simultaneously at a later time by the bistatic SAR. The two digital elevation maps DEM1 and DEM2 form the input for the calibration procedure described here.
[0064] In a first step ST1 of this method, the calibration points cp are determined using a so-called amplitude dispersion index D and a predetermined threshold value th for this dispersion index. A calibration point cp refers to a permanent scatterer candidate within the predetermined terrain TR, i.e., a point that has not changed from one digital elevation map to another and is thus a good point for calibrating the digital elevation maps DEM1 and DEM2. Each calibration point is represented by the same pixel position within both digital elevation maps DEM1 and DEM2.
[0065] To determine the dispersion index D, the method described here processes a time series ts of SAR images TSIM of the terrain TR acquired by a second radar station RE2, different from the first radar station RE1. Typically, the time series ts comprises at least 20 TSIM images. As with the SAR images IM, each of the SAR images TSIM contains a plurality of pixels covering the terrain TR, each pixel characterized by an amplitude value a and a phase value ph. The second radar station RE2 is the monostatic SAR of the Sentinel-1 mission, in which two satellites fly in the same nominal orbit, allowing the same terrain on Earth to be surveyed in multiple overpasses. For this mission, the nominal acquisition plan provides for a revisit time of six days over Europe and selected test areas and twelve days for the rest of the world.This enables the continuous recording of InSAR time series on a global scale.
[0066] Although the TanDEM-X mission has the unique feature of a bistatic InSAR, which enables the generation of high-quality digital elevation maps, the limited on-board resources make it impossible to continuously acquire long InSAR time series. Therefore, in the embodiment described here, the time series ts for determining the dispersion index D is acquired by the Sentinel-1 mission. However, if the first radar station RE1 also provides a long time series of SAR images in addition to the digital elevation maps, this time series can also be used to determine the dispersion index D. In this case, the first radar station RE1 and the second radar station RE2 coincide.
[0067] The following describes step ST1 of the procedure for Fig. 1 with reference to Fig. 2. Step ST1 comprises several substeps ST101 to ST107. In substeps ST101 to ST104, the time series ts of the SAR images TSIM acquired by the second radar system RE2 is processed. As a result of this processing, the candidate calibration points ccp are determined, which are candidates for calibration points cp within the terrain TR. Before proceeding with the steps in Fig. 2, the data from TanDEM-X and Sentinel-1 are mapped to the same coordinate system and appropriately interpolated on the same grid. This is a well-known preprocessing step that is used in Fig. 2 is omitted for simplicity.
[0068] In steps ST101 to ST103, a pre-filtering of the pixel positions within the TSIM images is performed. This pre-filtering identifies pixel positions that are not accepted as candidate calibration points, so that these pixel positions are not considered in the subsequent step ST104.
[0069] In step ST101, pixel positions within the SAR images TSIM are determined that correspond to points on the ground affected by geometric distortions GD. Geometric distortions are known effects in SAR data acquisition and lead to incorrect positioning of scatterers in the SAR image. Geometric distortions are caused by the side-looking nature of SAR and include, among other things, the well-known effects of layover and shadowing. Corresponding techniques for classifying pixel positions in SAR images as points affected by geometric distortions are known to those skilled in the art and are therefore not described in detail. Step ST101 results in a reduced set of pixel positions that no longer contain any points affected by geometric distortions.
[0070] In step ST102, the reduced number of pixels is further reduced by using an external reference land cover map (LCM), which indicates those pixel positions corresponding to areas covered by volumetric material such as snow or ice. The points in these areas may vary in elevation and be affected by radar penetration effects. Therefore, these pixel positions are discarded, resulting in a further reduced number of pixels.
[0071] In step ST103, the reduced set of pixels resulting from step ST102 is reduced again by taking into account the slope sl of the terrain at the corresponding pixel positions. Those pixel positions with a slope sl greater than a predetermined slope threshold slt are discarded. Pixel positions with a large slope lead to inaccuracies and should not be used as calibration points. The slope of the pixel positions is derived using an external digital reference elevation map (e.g., the global TanDEM-X DEM product).
[0072] The set of pixel positions remaining in the images TSIM after step ST103 has been performed are then analyzed based on the dispersion index D. The dispersion index is known, for example, from document [2] and is given by the following expression: D=σμ
[0073] The dispersion index D is calculated for each pixel position remaining after step ST103. For a respective pixel position, µ represents the temporal mean of the amplitude values of the pixels at the respective pixel position across the images TSIM of the time series ts. Furthermore, σ corresponds to the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position across the images TSIM of the time series ts. Points on the earth are considered to be temporally stable if the corresponding pixel position has a dispersion index D that is smaller than the predetermined threshold value th. The threshold value th is calculated automatically, as described below with reference to Fig. 3 is described.
[0074] Step ST104 results in a number of pixel positions that are candidate calibration points (ccp). From these candidate calibration points (ccp), suitable calibration points (cp) are selected based on the SAR image pairs IP1 and IP2, from which the corresponding digital elevation maps DEM1 and DEM2 are derived. In step ST105, the pixel positions affected by geometric distortions (GD) within the image pairs IP1 and IP2 are determined in the same way as for the TSIM images in step ST101. This step ensures that the ultimately selected elevation points are not affected by geometric distortions in either the Sentinel-1 acquisition geometry or the TanDEM-X acquisition geometry. The pixel positions affected by geometric distortions are excluded from the candidate calibration points (ccp).
[0075] In step ST106, a threshold value ct is set for the bistatic interferometric SAR coherence of the two bistatic TanDEM-X images (i.e., the two images IP1 and IP2) to ensure good quality of the digital TanDEM-X elevation map in accordance with the selected targets. The interferometric coherence is a known quantity and is defined as the absolute value of the normalized cross-correlation coefficient between the images of an interferometric image pair. The embodiment described here processes the interferometric coherence between the images of the pair IP1 and the interferometric coherence between the images of the pair IP2. The interferometric coherence is given by the following expression: γtot=|E[u1u2]|E[|u1|2|u2*|2] u1 and u2 represent the images of the corresponding image pair IP1 or IP2. E is the expectation value and * is the complex conjugate operator. The interferometric coherence γ tot quantifies the amount of noise in the interferogram and therefore represents a key parameter for assessing the quality of any digital height map. The interferometric coherence can be estimated for each pixel position across a number of neighboring pixels using known algorithms. The estimated interferometric coherence for each pixel position in both image pairs IP1 and IP2 is compared with a predefined threshold ct. Preferably, the predefined threshold has a value of 0.6 or greater, in particular a value between 0.6 and 0.7. Those pixel positions that have estimated interferometric coherences below the predefined threshold ct are excluded from the candidate calibration points ccp.
[0076] The reduced number of candidate calibration points ccp resulting from step ST106 is processed in a further filtering step ST107. In this step, the signal-to-noise ratio (SNR) of each pixel of the images of the image pairs IP1, IP2 is calculated as follows: SNR=(β0−βn)βn
[0077] The quantities in the above formula are known. β 0 indicates the radar brightness and β nBeta-zero is the noise equivalent, which represents the noise floor of the SAR system. A candidate calibration point is discarded if the corresponding pixel position of the candidate calibration point belongs to a pixel within the image pairs that has a signal-to-noise ratio (SNR) below a predefined threshold (st). This ensures that the calibration points used are targets that exhibit not only sufficient amplitude stability but also sufficient phase stability.
[0078] The candidate calibration points ccp remaining after filtering in steps ST105 to ST107 form the calibration points cp, which are calculated in the subsequent steps ST2 to ST4 of the Fig. 1. In step ST2 of the Fig. 1, elevation value differences Δ between the digital elevation maps DEM1 and DEM2 are calculated for the calibration points cp. These elevation value differences can be expressed by the following formula: Δ=DEM1(cp)−DEM2(cp)
[0079] Then, in step ST3, the height value differences at the calibration points are used as data points to fit a two-dimensional curve CS to the height value differences. The two-dimensional curve CS is referred to below as the calibration surface. Any known interpolation technique can be used to fit the two-dimensional curve. The interpolation can be a linear interpolation or another higher-order interpolation. As a result, the calibration surface CS describes a function that provides a height value difference for each of the pixels within the digital elevation maps DEM1, DEM2. The calibration surface CS is defined in the final step ST4 of the Fig. 1 is used to correct the digital elevation map DEM2, leaving the calibration DEM1 unchanged. The corrected digital elevation map DEM2' is calculated for each pixel position within DEM2 based on the following formula: DEM2'=DEM2+Δcs
[0080] In the above formula, Δ CS to the height value difference given by the calibration surface CS at the respective pixel position, and DEM2 is the height value of DEM2 at the respective pixel position.
[0081] For the above correction, the digital elevation map DEM1 is considered the digital reference elevation map. However, the digital elevation map DEM2 can also be considered the reference elevation map. In this case, the indices for the digital elevation maps in the above formulas (4) and (5) are swapped.
[0082] The procedure of Fig. 1 results in two appropriately mutually calibrated digital elevation maps, DEM1 and DEM2, thanks to the use of permanent scatterer candidates as calibration points. The calibrated digital elevation maps can be used for DEM differencing by determining the elevation differences between the digital elevation maps. This allows the identification of areas where the elevation of the terrain has changed over time, for example, to capture the melting of glaciers or deforestation.
[0083] The threshold value th for the dispersion index D used in step ST1 is calculated in a separate calculation procedure that is performed before the actual calibration procedure. This calculation procedure is described in Fig. 3. The calculation procedure is particularly advantageous when the first radar system RE1 is a different radar system than the second radar system RE2. This is the case in the embodiment described here, since the radar system RE1 relates to the TanDEM-X mission, while the radar system RE2 relates to the Sentinel-1 mission.
[0084] The procedure of Fig. 3 processes a first time series ts1 of SAR images TSIM1 acquired by the radar RE1 of the TanDEM-X mission. Furthermore, the method processes a second time series ts2 of SAR images TSIM2 acquired by the radar RE2 of the Sentinel-1 mission. All images cover the same predefined terrain TR' on the Earth's surface with the same resolution, with each point in the terrain TR' corresponding to the same pixel position within the images. The predefined terrain TR' does not have to be the same terrain as the terrain TR acquired in the method of Fig. 1 is considered.
[0085] The procedure of Fig. 3 is based on the realization that the threshold th applied to the Sentinel-1 data must be made compatible with the TanDEM-X data from which the digital elevation maps to be calibrated are derived. The data acquisitions of the Sentinel-1 mission and the TanDEM-X mission do indeed exhibit significant differences. The Sentinel-1 multi-pass data are acquired in C-band, while the bistatic TanDEM-X data are acquired in X-band. When considering radar sensors at different frequencies, one expects different behavior in terms of backscatter levels and interferometric performance, i.e., coherence.
[0086] The computational procedure described below is an empirical approach developed to establish a suitable threshold th for the dispersion index D calculated for the Sentinel-1 multi-pass time series, ensuring the selection of stable calibration points when considering the same area acquired in X-band using bistatic TanDEM-X images.
[0087] In the proceedings of the Fig. 3, the dispersion indices of the respective pixel positions within the SAR images TSIM1 of the first time series ts1, which relate to the TanDEM-X mission, are determined as first dispersion indices D1 based on the above formula (1). Analogously, the dispersion indices of the respective pixel positions within the SAR images TSIM2 of the second time series ts2, which relate to the Sentinel-1 data, are calculated as second dispersion indices D2 based on the above formula (1). In a next step, the histogram HI is calculated as a function of the two dispersion indices D1 and D2. The brightness of a point within the histogram HI encodes the number of pixel positions within the images TSIM1, TSIM2 (i.e., the number of corresponding points within the terrain TR') that have the values of the dispersion indices D1 and D2 that correspond to the position of the point in the histogram.The higher the brightness, the higher the number of pixel positions. In other words, each point in the histogram corresponds to a specific value of D1 and a specific value of D2, where the histogram indicates the number of pixel positions that have the specific value of D1 and the specific value of D2.
[0088] For the first dispersion indices D1, a specific limit value li is specified. In the embodiment described here, this limit is set to 0.25, analogous to document [2]. However, other values for the limit value li can also be used. Preferably, the limit value lies in the range between 0.15 and 0.3. The limit value li is in Fig. 3 is indicated by a dotted line. In a next step, the threshold value th for the dispersion indices D2, which is Fig.3 is indicated by a dashed line, is chosen so that the number of pixel positions (ie points on the terrain TR') in the lower right quadrant bounded by the lines of the limit value li and the threshold value th corresponds to a preset percentage pc of the pixel positions (10% or less) within the total number of pixel positions in the corresponding images TSIM1, TSIM2, this total being the same for all images TSIM1 and TSIM2 and corresponding to the total number of points within the terrain TR'.In other words, the threshold value th is selected such that both the pixel positions within the respective SAR images TSIM1 of the first time series ts1 which have a first dispersion measure D1 above the predetermined limit value li and the corresponding pixel positions of the respective SAR images TSIM2 of the second time series ts2 which have a second dispersion measure D2 below the predetermined threshold value th amount to the predetermined percentage.
[0089] The threshold th is chosen as a compromise between the maximum acceptable false alarm rate and the minimum number of selected targets that must be retained for calibration. False alarms correspond to targets that, if selected, would not meet the stability requirement for the TanDEM-X data (first dispersion index D1 above the threshold li), but would meet the stability requirement for the Sentinel-1 data (dispersion index D2 less than the threshold th). In the embodiment described here, the percentage pc is set to 1%. This results in a value for the threshold th of approximately 0.19. Preferably, the percentage pc is set to 5% or less.
[0090] The method described above was evaluated by the inventors based on corresponding data from the TanDEM-X and Sentinel-1 missions. The terrain used for the digital elevation maps refers to a sandy / rocky area in the Sahara Desert, Egypt. The elevation differences between the uncalibrated digital elevation maps and the calibrated digital elevation maps were determined, and it was found that the remaining inclinations are correctly compensated after applying the calibration procedure, while the remaining inclinations remain visible without mutual calibration.
[0091] The invention, as described above, has several advantages. In particular, suitable calibration points for a pair of digital elevation maps acquired at different times are extracted from SAR data without using external reference measurements. The calibration points are selected by considering the amplitude dispersion index of the corresponding pixel positions as candidates for permanent scatterers. To determine the dispersion index, the method can use data from SAR images other than the SAR image from which the digital elevation maps to be calibrated are derived. Furthermore, the method also provides a calculation procedure for the threshold to be applied to the dispersion index to select candidate calibration points. List of references [1] A. Gruber, B. Wessel, M. Huber, A. Roth: „Operational TanDEM-X DEM calibration and first validation result“, ISPRS Journal of Photogrammetry and Remote Sensing 73 (2012), S. 39-49. [2] A. Ferretti, C. Prati, F. Rocca: „Permanent Scatterers in SAR Interferometry“, IEEE Transactions on Geoscience and Remote Sensing, Bd. 39, Nr. 1, Januar 2001. [3] Yongjiu Feng, Yilun Zhou, Yanling Chen, Pengshuo Lia, Mengrong Xi and Xiaohua Tang, "Automatic selection of permanent scatterers-based GCPs for refinement and reflattening in InSAR DEM generation", International Journal of Digital Earth, 2022, Bd. 15, Nr. 1, S. 954-974. [4] G. Krieger, A. Moreira, H. Fiedler, 1. Hajnsek, M. Werner, M. Younis, and M. Zink, „TanDEM-X: A satellite formation for high-resolution SAR interferometry“, IEEE Transactions on Geoscience and Remote Sensing, Bd. 45, Nr. 11, S. 3317-3341, Oktober 2007. [5] RIZZOLI, Paola [et al.], Sentinel-1-Aided Calibration of TanDEM-X DEMs for Estimating Winter Snow Accumulation over Alpine Glaciers. 24.05.2022. URL: https: / / earth.esa.int / livingplanet-symposium-2022-presentations / 24.05.Tuesday / Nairobi_3-4 / 1540-1720 / 05_Rizzoli.pdf
Claims
[1] A computer-implemented method for calibrating a pair of digital elevation maps (DEM1, DEM2), wherein the digital elevation maps (DEM1, DEM2) of the pair indicate the elevation of a predetermined terrain (TR) on a surface of a celestial body and comprise elevation values (hv1, hv2) for a plurality of pixels at pixel positions representing corresponding points within the predetermined terrain (TR), wherein one digital elevation map (DEM1, DEM2) of the pair indicates the elevation of the predetermined terrain (TR) at an earlier time than the other digital elevation map (DEM1, DEM2) of the pair, wherein each digital elevation map (DEM1, DEM2) of the pair is linked to a different SAR image pair (IP1, IP2) from SAR images (IM) acquired from measurements of a first radar sensor system (RE1), wherein a respective digital elevation map (DEM1, DEM2) is derived from phase differences between the linked SAR image pair (IP1, IP2) has been derived,wherein the method processes a time series (ts) of SAR images (TSIM) of the predetermined terrain (TR) obtained from measurements of a second radar sensor system (RE2), each SAR image (TSIM) of the time series (ts) comprising pixels at the same pixel positions as the digital elevation maps (DEM1, DEM2) and as the SAR images (IM) obtained from measurements of the first radar sensor system (RE1), the method comprising the following steps: a) determining candidate calibration points (ccp) from the time series (ts) of SAR images (TSIM) based on dispersion measures (D) of respective pixel positions within the SAR images (TSIM) of the time series (ts), wherein the dispersion measure (D) of a respective pixel position describes the statistical dispersion of amplitude values (a) of the pixels at the respective pixel position across the SAR images (TSIM) of the time series (ts), wherein each candidate calibration point (ccp) refers to a pixel position having a dispersion measure (D) below a predetermined threshold value (th); b) determining pixel positions as calibration points (cp) within the digital height maps (DEM1, DEM2) of the pair based on the candidate calibration points (ccp), wherein each calibration point (cp) corresponds to a pixel position of a candidate calibration point (ccp); c) determining height value differences (Δ) between the digital height maps (DEM1, DEM2) of the pair at the calibration points (cp); d) fitting a two-dimensional curve (CS) to the height value differences (Δ) at the calibration points (cp), wherein the two-dimensional curve (CS) describes the height value differences (Δ) over the pixel positions representing corresponding points within the predetermined terrain (TR); e) correcting the elevation values (hv1, hv2) of a digital elevation map (DEM1, DEM2) of the pair using the elevation value differences (Δ) of the two-dimensional curve (CS); wherein the method comprises a computational procedure for calculating the predetermined threshold value (th) based on a first and a second time series (ts1, ts2) of SAR images (TSIM1, TSIM2) of a predefined terrain (TR') on the surface of the celestial body, wherein the first time series (ts1) has been obtained from measurements of the first radar sensor system (RE1) and the second time series (ts2) has been obtained from measurements of the second radar sensor system (RE2), wherein all SAR images of the first and second time series (ts1, ts2) comprise pixels at the same pixel positions, representing corresponding points within the predefined terrain (TR'), wherein the computational procedure comprises the following steps: i) determining dispersion measures of respective pixel positions within the SAR images (TSIM1) of the first time series (ts1) as first dispersion measures (D1), wherein the dispersion measure of a respective pixel position describes the statistical scatter of amplitude values (a) of the pixels at the respective pixel position over the SAR images (TSIM1) of the first time series (ts); ii) determining dispersion measures of respective pixel positions within the SAR images (TSIM2) of the second time series (ts2) as second dispersion measures (D2), wherein the dispersion measure of a respective pixel position describes the statistical scatter of amplitude values (a) of the pixels at the respective pixel position over the SAR images (TSIM2) of the second time series (ts2); iii) determining the predetermined threshold value (th) based on a predetermined limit value (li) for the first dispersion measures (D1) by selecting the predetermined threshold value (th) such that both the pixel positions within the respective SAR images (TSIM1) of the first time series (ts1) having a first dispersion measure (D1) above the predetermined limit value (li) and the corresponding pixel positions of the respective SAR images (TSIM2) of the second time series (ts2) having a second dispersion measure (D2) below the predetermined threshold value (th) amount to a predetermined percentage of 10% or less. [2] Method according to claim 1, wherein the dispersion measure (D) of a respective pixel position within the SAR images (TSIM) of the time series (ts) represents or depends on a dispersion index which is the standard deviation of the statistical distribution of the amplitude values (a) of the pixels at the respective pixel position over the SAR images (TSIM) of the time series (ts) divided by the mean of this statistical distribution. [3] Method according to claim 1 or 2, wherein in step a), in the event that a pixel position within the SAR images (TSIM) of the time series (ts) is determined to be affected by geometric distortions (GD), this pixel position is classified as a pixel position which is not admitted as a candidate calibration point (ccp). [4] Method according to one of the preceding claims, wherein in step a), in the event that a pixel position within the SAR images (TSIM) of the time series (ts) is determined to be a point within the predetermined terrain (TR) covered by snow or ice, this pixel position is classified as a pixel position which is not admitted as a candidate calibration point (ccp). [5] Method according to one of the preceding claims, wherein in step a), in the event that a pixel position within the SAR images (TSIM) of the time series (ts) is identified as a point within the predetermined terrain (TR) having a slope (sl) greater than a predetermined slope threshold (slt), this pixel position is classified as a pixel position that is not admitted as a candidate calibration point (ccp). [6] Method according to one of the preceding claims, wherein in step b), in the event that a pixel position within the digital height maps (DEM1, DEM2) is determined to be affected by geometric distortions (GD), this pixel position is classified as a pixel position which is not admitted as a calibration point (cp). [7] Method according to one of the preceding claims, wherein in step b) in the event that the interferometric coherence (γ tot ) for a pixel position of a SAR image pair (IP1, IP2) associated with one of the digital elevation maps (DEM1, DEM2) of the pair is below a coherence threshold (ct), that pixel position is classified as a pixel position not eligible as a calibration point (cp). [8] Method according to one of the preceding claims, wherein in step b), if the signal-to-noise ratio (SNR) of a pixel at a pixel position within at least one SAR image (IM) of the SAR image pairs (IP1, IP2) associated with the digital elevation maps (DEM1, DEM2) of the pair is below a signal-to-noise ratio threshold (st), this pixel position is classified as a pixel position that is not admitted as a calibration point (cp). [9] Method according to one of the preceding claims, wherein the first radar system (RE1) is a different radar system than the second radar system (RE2). [10] Method according to one of the preceding claims, wherein the first radar system (RE1) is a monostatic radar system with one or more radar sensors, in which each radar sensor transmits its own radar pulses and measures exclusively the radar echoes of its own radar pulses, or a bistatic radar system with two radar sensors separated by a predetermined distance, in which only one radar sensor transmits radar pulses and both radar sensors measure the radar echoes of the radar pulses, and / or wherein the second radar system (RE2) is a monostatic radar system with one or more radar sensors, in which each radar sensor transmits its own radar pulses and measures exclusively the radar echoes of its own radar pulses, or a bistatic radar system with two radar sensors separated by a predetermined distance, in which only one radar sensor transmits radar pulses and both radar sensors measure the radar echoes of the radar pulses. [11] Method according to one of the preceding claims, wherein the first dispersion measure (D1) of a respective pixel position within the SAR images (TSIM1) of the first time series (ts1) represents or depends on a dispersion index which is the standard deviation of the statistical distribution of the amplitude values (a) of the pixels at the respective pixel position across the SAR images (TSIM1) of the first time series (ts1) divided by the mean value of this statistical distribution, and the second dispersion measure (D2) of a respective pixel position within the SAR images (TSIM2) of the second time series (ts2) represents or depends on a dispersion index which is the standard deviation of the statistical distribution of the amplitude values (a) of the pixels at the respective pixel position across the SAR images (TSIM2) of the second time series (ts2) divided by the mean value of this statistical distribution. [12] Method according to one of the preceding claims, wherein the predetermined limit value (li) has a value between 0.15 and 0.3, preferably 0.25, and / or the predetermined percentage (pc) has a value of 5% or less, preferably 1%. [13] Computer program product comprising program code stored on a machine-readable carrier for carrying out a method according to any one of the preceding claims when the program code is executed on a computer. [14] Computer program with program code for carrying out a method according to one of claims 1 to 12, when the program code is executed on a computer.