A computer-implemented method for calibrating a pair of digital elevation maps

Abstract

The invention refers to a computer-implemented method for calibrating a pair of digital elevation maps (DEMI, DEM2), where the digital elevation maps (DEMI, DEM2) of the pair have been acquired from measurements of a first radar sensor equipment (RE1) and indicate for different time points the height of a predetermined terrain (TR) on a surface of a celestial body, wherein the method processes a time series (ts) of SAR images (TSIM) of the predetermined terrain (TR) acquired from measurements of a second radar sensor equipment (RE2), 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); b) identifying pixel positions as calibration points (cp) within the digital elevation maps (DM1, DM2) of the pair based on the candidate calibration points (ccp); c) determining height value differences (Δ) between the digital elevation maps (DEMI, DEM2) of the pair at the calibration points (cp); d) fitting a two-dimensional curve (CS) to the height value differences (A) at the calibration points (cp), the two-dimensional curve (CS) describing the height value differences (Δ) over the pixel positions representing corresponding points within the predetermined terrain (TR); e) correcting the height values (hvl, hv2) of one digital elevation map (DEMI, DEM2) of the pair by the height value differences (Δ) of the two-dimensional curve (CS).

Description

[0001] A computer-implemented method for calibrating a pair of digital elevation maps

[0002] Description

[0003] The invention refers 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.

[0004] Synthetic aperture radar (commonly abbreviated as SAR) is an established technique for remote sensing of a surface of a celestial body, e.g. the Earth. To do so, a radar equipment on a platform (e. g. on a satellite or an airplane) flying over the ground of the celestial body transmits radar pulses and measures the radio echoes from the pulses scattered on the surface of the celestial body. With well-known SAR processing techniques, SAR images can be derived from the measured radar echoes.

[0005] Those images comprise an amplitude value and a phase value for a plurality of pixels corresponding to points on the surface of the celestial body.

[0006] Digital elevation maps can be derived from pairs of SAR images taken from the same terrain by analyzing the phase differences between the images. A corresponding digital elevation map provides height values for a plurality of pixels at corresponding pixel positions representing points on the surface of the celestial body.

[0007] In order to analyze changes in a terrain on the 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 in order to determine height differences between those digital elevation maps. As a consequence, topographic changes can be monitored and quantified. E.g., the mass balance of glaciers can be monitored, or the amount of lava erupted by a volcano can be quantified. The digital elevation maps used for DEM differencing are derived from different SAR images. Therefore, a precise mutual calibration of the digital elevation maps is a crucial step. Particularly, residual horizontal shifts and vertical offsets and tilts must be compensated, which, if neglected, could lead to errors in the order of several meters.

[0008] Mutual calibration of digital elevation maps is usually achieved by performing a precise spatial co-regi strati on of the geocoded digital elevation maps and by identifying reliable reference calibration points by external measurements, such as GPS or Li- D AR measurements. Document [1] discloses the calibration of digital elevation maps using laser measurements as external reference calibration points.

[0009] The acquisition of the external reference measurements for calibrating digital elevation maps is often expensive and problematic, especially when remote or scarcely accessible regions are considered. Furthermore, common calibration for digital elevation maps solely compensates for an absolute offset between a pair of maps without considering residual tilts which may be caused by a limited accuracy in the estimation of geometric acquisition parameters of the SAR system.

[0010] Document [2] discloses a procedure for identifying permanent scatterers in a temporal series of SAR images. Document [3] describes the use of permanent scatterers in SAR images for improving the generation of a single digital elevation map.

[0011] It is an object of the invention to provide to a computer-implemented method for calibrating a pair of digital elevation maps without the need of external measurements at reference calibration points.

[0012] This object is solved by the method according to claim 1 and the computer program product according to claim 14 and the computer program according to claim 15, respectively. Preferred embodiments of the invention are defined in the dependent claims.

[0013] The invention provides a computer-implemented method for calibrating a pair of digital elevation maps. The digital elevation maps of the pair indicate the height of a predetermined terrain on a surface of a celestial body. Preferably, the celestial body is a planet, and the planet is preferably the Earth. The digital elevation maps each comprise height 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 height of the predetermined terrain earlier in time than the other digital elevation map of the pair.

[0014] Each digital elevation map of the pair is associated with another SAR image pair out of SAR images acquired from measurements of a first radar sensor equipment, where a respective digital elevation map has been derived from phase differences between the associated SAR image pair. The structure of SAR images is well-known for a skilled person. A SAR image is characterized by pixels at a plurality of pixel positions where each pixel in characterized by an amplitude and phase value.

[0015] The method according to the invention processes a time series of SAR images of the predetermined terrain acquired from measurements of a second radar equipment, where each SAR image of the time series comprises pixels at the same pixel positions as the digital elevation maps and as the SAR images acquired from measurements of the first radar equipment. The first radar equipment and the second radar equipment may refer to the same radar equipment. I.e., the first radar equipment may coincide with the second radar equipment. However, the first radar equipment may also be a different radar equipment with respect to the second radar equipment.

[0016] The first radar equipment may be a monostatic radar equipment with one or more radar sensor, where each radar sensor transmits its own radar pulses and exclusively measures the radar echoes of its own radar pulses. Nevertheless, the first radar equipment may also be a bistatic radar equipment with two co-flying radar sensors separated by a predetermined distance, where only one radar sensor transmits radar pulses and both radar sensors measure the radar echoes of the radar pulses. Analogously, the second radar equipment may be a monostatic radar equipment or a bistatic radar equipment. In the detailed description, an embodiment is described where the first radar equipment is a bistatic SAR on satellites of the TanDEM-X mission whereas the second radar equipment is a monostatic SAR on satellites of the Sentinel- 1 mission.

[0017] The method according to the invention comprises steps a) to e) described in the following. In step a), candidate calibration points are determined from the time series of SAR images based on (amplitude) dispersion measures of respective pixel positions within the SAR images of the time series. The dispersion measure of a respective pixel position describes the statistical dispersion of amplitude values of the pixels at the respective pixel positions over the SAR images of the time series, where each candidate calibration point refers to a pixel position having a dispersion measure below a predetermined threshold. The predetermined threshold may be set in advance. However, in an embodiment described below, the invention also includes a calculation procedure for the predetermined threshold. The dispersion measure may be defined in various ways. It is characterized by the feature that a value of the dispersion measure becomes higher with an increasing statistical dispersion of the amplitude values.

[0018] In step b), pixel positions are identified as calibration points within the digital elevation maps of the pair 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. Nevertheless, a candidate calibration point is not necessarily a calibration point. In other words, a filter may be used for extracting calibration points out of the candidate calibration point. In step c), height value differences between the digital elevation maps of the pair are determined at the calibration points and, in step d), a two-dimensional curve (i.e. a two-dimensional function) is fitted to the height value differences at the calibration points, the two-dimensional curve describing the height value differences over the pixel positions representing corresponding points within the predetermined terrain. I.e., the two-dimensional curve describes the residual height value differences between both digital elevation maps of the pair and thus forms 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 may be used for determining the two-dimensional curve. E.g., a linear interpolation or an interpolation based on high order polynomials may be used to obtain the two-dimensional curve.

[0019] In step e), the height values of one digital elevation map of the pair are corrected by the height value differences of the two-dimensional curve. The corrected digital elevation map and the other digital elevation map (not corrected) form a pair of digital elevation maps which are properly calibrated to each other thanks to the determination of a two-dimensional curve based on suitable calibration points. The mutually calibrated digital elevation maps may then be used for DEM differencing in order to obtain accurate information about the change of the terrain which has occurred during the time span between the two acquisitions of the digital elevation maps.

[0020] The invention is based on the finding that permanent scatterer candidates are selected based on a statistical dispersion measure of a temporal series of SAR images are well suited to become calibration points for a mutual calibration of digital elevation maps.

[0021] 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 being the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position over the SAR images of the time series divided by the mean value of this statistical distribution. The dispersion index as defined above is also used in the method described in document [2], In further preferred embodiments, adequate techniques are used to reduce the number of candidate calibration points determined in step a). In one embodiment, in case that a pixel position within the SAR images of the time series is identified as being affected by geometric distortions, this pixel position is classified as a pixel position not permitted to be a candidate calibration point. Techniques to identify pixel positions as being affected by geometric distortions, e.g. by layover or shadow, are well- known and will not be described in detail herein.

[0022] In another embodiment, in case that a pixel position within the SAR images of the time series is identified as being a point within the predetermined terrain covered by volumetric material, this pixel position is classified as a pixel position not permitted to be a candidate calibration point. Volumetric material has the characteristic that radar waves can penetrate into the volumetric material. The classification of terrain as being covered by volumetric material lies within the expert knowledge of a skilled person. Volumetric material is e.g. snow, ice, vegetation, sand and the like. To identify a pixel position as a point covered by volumetric material, a suitable reference map may be used. This embodiment is based on the finding that regions covered with snow or ice are not stable in time and thus are not suitable to become candidate calibration points.

[0023] In another embodiment, in case that a pixel position within the SAR images of the time series is identified as being a point within the predetermined terrain having a slope higher than a predetermined slope threshold, this pixel position is classified as a pixel position not permitted to be a candidate calibration point. This embodiment avoids candidate calibration points which are not accurately detected.

[0024] Furthermore, appropriate filter techniques may be used in step b) for determining calibration points out of the candidate calibration points. In one embodiment, in case that a pixel position within the digital elevation maps is identified as being affected by geometric distortions, this pixel position is classified as a pixel position not permitted to be a calibration point. As mentioned above, methods for identifying pixel positions as being affected by geometric distortions are well-known for a skilled person.

[0025] In another embodiment, in case that the interferometric coherence for a pixel position of a SAR image pair associated with one of the digital elevation maps of the pair is below a coherence threshold, this pixel position ins classified within the digital elevation maps as a pixel position not permitted to be a calibration point. The interferometric coherence is a well-known quantity which can be estimated pixelwise with respect to image pairs. In the detailed description, the definition of the interferometric coherence is provided.

[0026] In another preferred embodiment, in case that 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, this pixel position is classified within the digital elevation maps as a pixel position not permitted to be a calibration point.

[0027] As mentioned above, the method according to the invention may comprises a calculation procedure for calculating the predetermined threshold for the statistical dispersion measure. This calculation procedure is based on first and second time series of SAR images of a predefined terrain on the surface of the celestial body. The first time series has been acquired from measurements of the first radar sensor equipment and the second time series has been acquired from measurements of the second radar sensor equipment, where all SAR images of the first and second time series comprise pixels at the same pixel positions representing corresponding points within the predefined terrain. The predefined terrain needs not be the same terrain as the predetermined terrain covered by the digital elevation maps.

[0028] The calculation procedure is preferably used in case that the first radar equipment is another radar equipment than the second radar equipment so that the threshold being applied to the dispersion measure derived from SAR images of the second radar equipment is adequately chosen for the pair of digital elevation maps being based on measurements of the first radar equipment.

[0029] 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, the dispersion measure of a respective pixel position describing the statistical dispersion of amplitude values of the pixels at the respective pixel position over the SAR images of the first time series.

[0030] 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, the dispersion measure of a respective pixel position describing the statistical dispersion of amplitude values of the pixels at the respective pixel position over the SAR images of the second time series.

[0031] In a step iii) of the calculation procedure, the predetermined threshold is determined based on a preset limit for the first dispersion measures by choosing the predetermined threshold such that both the pixel positions within the respective SAR images of the first time series having a first dispersion measure above the preset limit and the corresponding pixel positions of the respective SAR images of the second time series having a second dispersion measure below the predetermined threshold amount to a preset percentage of 10% or less.

[0032] The first dispersion measures and the second dispersion measures may correspond to the above defined dispersion indices. I.e., the first dispersion measure of a respective pixel position within the SAR images of the first time series may represent or depend on a dispersion index being the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position over the SAR images of the first time series divided by the mean value of this statistical distribution. Analogously, the second dispersion measure of a respective pixel position within the SAR images of the second time series may represent or depend on a dispersion index being the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position over the SAR images of the second time series divided by the mean value of this statistical distribution.

[0033] In a preferred embodiment, the above preset limit has a value between 0.15 and 0.3, preferably 0.25, and / or the above preset percentage has a value of 5% or less, preferably 1%.

[0034] Besides the above method, the invention refers to a computer program product with program code, which is 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.

[0035] Furthermore, the invention refers 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.

[0036] In the following, second, third and fourth aspects of the invention which are independent of the aspects of the invention as described above will be explained.

[0037] The second aspect refers to a method for calibrating digital elevation models (DEMs) of a terranean surface of the Earth, the method comprising:

[0038] (a) obtaining a pair of co-registered DEMs corresponding to a selected region of the terranean surface of the Earth, wherein each of the pair of co-registered DEMs are generated from single-pass synthetic aperture radar (SAR) data;

[0039] (b) obtaining repeat-pass time-series SAR data pertaining to the selected region of the terranean surface;

[0040] (c) selecting a plurality of calibration tie-points (TPs) representative of a portion of the selected region of the terranean surface using both the single-pass SAR data and the repeat-pass time-series SAR data; and (d) generating a pair of mutually-calibrated DEMs corresponding to the selected region of the terranean surface using both the pair of co-registered DEMs and the selected plurality of calibration TPs.

[0041] A preferred embodiment of the second aspect refers to a method, wherein (c) comprises:

[0042] (cl) selecting candidate calibration TPs from the repeat-pass time series data; and (c2) selecting the calibration TPs from the pair of co-registered DEMs using the candidate calibration TPs.

[0043] A preferred embodiment of the second aspect refers to a method, wherein (c) comprises:

[0044] (cl) filtering pixels from the repeat-pass time-series data to provide a plurality of candidate calibration TPs;

[0045] (c2) identifying 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 identified candidate calibration TPs whereby the remaining, unfiltered candidate calibration TPs comprise the selected calibration TPs.

[0046] A preferred embodiment of the second aspect refers to a method, wherein (c) comprises applying a predefined dispersion index amplitude threshold on pixels from the pair of co-registered DEMs.

[0047] A preferred embodiment of the second aspect refers to a method, wherein (c) comprises selecting as the calibration TPs only the calibration TPs meeting both a repeatpass dispersion index amplitude threshold and a single-pass dispersion index amplitude threshold. Preferably, the single-pass dispersion index amplitude threshold is different from the repeat-pass dispersion index amplitude threshold.

[0048] A preferred embodiment of the second aspect refers to a method, wherein (d) comprises: (dl) generating a calibration surface; and

[0049] (d2) applying the calibration surface to one of the pair of co-registered DEMs to generate a calibrated DEM co-registered with the pair of co-registered DEMs, the pair of mutually-calibrated DEMs comprising the other of the pair of co-registered DEMs and the calibrated DEM.

[0050] A preferred embodiment of the second aspect refers to a method, wherein (d) comprises taking the difference between the pair of co-registered DEMs across the calibration TPs selected at (c).

[0051] A preferred embodiment of the second aspect refers to a method, wherein the singlepass SAR data and the repeat-pass time-series SAR data both comprises interferometric synthetic aperture radar (InSAR) data.

[0052] A preferred embodiment of the second aspect refers to a method, wherein the singlepass SAR data is obtained from a first original data source and the repeat-pass timeseries SAR data is obtained from a second original data source that is different from the first original data source.

[0053] A preferred embodiment of the second aspect refers to a method, further comprising: (e) generating a calibrated DEM of difference (DoD) by taking the difference between the mutually-calibrated DEMs.

[0054] The third aspect refers to a method for calibrating digital elevation models (DEMs) of a terranean surface of the Earth, the method comprising:

[0055] (a) obtaining a pair of co-registered DEMs corresponding to a selected region of the terranean surface of the Earth, wherein each of the pair of co-registered DEMs are 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 that is different from the first radar band, the second SAR data pertaining to the selected region of the terranean surface;

[0056] (c) selecting a plurality of calibration tie-points (TPs) representative of a portion of the selected region of the terranean surface using both the first SAR data and the second SAR data; and

[0057] (d) generating a pair of mutually-calibrated DEMs corresponding to the selected region of the terranean surface using both the pair of co-registered DEMs and the selected plurality of calibration TPs.

[0058] A preferred embodiment of the third aspect refers to a method, wherein the first SAR data corresponds to a first radar band and the second SAR data corresponds to a second radar band that is different from the first radar band.

[0059] A preferred embodiment of the third aspect refers to a method, wherein the first SAR data comprises X-band radar data and the second SAR data comprises C-band radar data.

[0060] A preferred embodiment of the third aspect refers to a method, wherein the first SAR data and the second SAR data both comprises interferometric synthetic aperture radar (InSAR) data.

[0061] A preferred embodiment of the third aspect refers to a method, wherein (d) comprises taking the difference between the pair of co-registered DEMs across the calibration TPs selected at (c).

[0062] A preferred embodiment of the third aspect refers to a method, wherein (c) comprises:

[0063] (cl) selecting candidate calibration TPs from the second SAR data; and

[0064] (c2) selecting the calibration TPs from the pair of co-registered DEMs using the candidate calibration TPs. The fourth aspect refers to a system for calibrating digital elevation models (DEMs) of a terranean surface of the Earth, the system comprising: a processor; a non-transitory memory; and one or more applications stored in the non-transitory memory that, when executed by the processor: obtain a pair of co-registered DEMs corresponding to a selected region of the terranean surface of the Earth, wherein each of the pair of co-registered DEMs are generated from single-pass synthetic aperture radar (SAR) data; obtains repeat-pass time-series SAR data pertaining to the selected region of the terranean surface; selects a plurality of calibration tie-points (TPs) representative of a portion of the selected region of the terranean surface using both the single-pass SAR data and the re- peat-pass time-series SAR data; and generates a pair of mutually-calibrated DEMs corresponding to the selected region of the terranean surface using both the pair of co-registered DEMs and the selected plurality of calibration TPs.

[0065] A preferred embodiment of the fourth aspect refers to a system, wherein the one or more applications stored in the non-transitory memory that, when executed by the processor: take the difference between the pair of co-registered DEMs across the selected calibration TPs to generate the pair of mutually calibrated DEMs.

[0066] A preferred embodiment of the fourth aspect refers to a system, wherein the one or more applications stored in the non-transitory memory that, when executed by the processor: select candidate calibration TPs from the repeat-pass time series data; and select the calibration TPs from the pair of co-registered DEMs using the candidate calibration TPs. In the following, embodiments of the invention will be described in detail with respect to the accompanying drawings wherein:

[0067] Fig. 1 shows a flow chart illustrating the steps of an embodiment of a method according to the invention;

[0068] Fig. 2 shows a flow chart illustrating the determination of candidate calibration points as indicated in step STI of Fig. 1;

[0069] Fig. 3 is a diagram illustrating a method for determining the threshold used in step STI of Fig. 1.

[0070] In the following, the invention is described with respect to a mutual calibration of single-pass InSAR digital elevation maps (InSAR = Interferometic SAR). The digital elevation maps are acquired by the bistatic SAR of the TanDEM-X mission (see document [4]). In this bistatic SAR system, two satellites are flying in close formation over the Earth's ground. Each satellite has a radar sensor where the radar sensor of one satellite transmits radar pulses and where the radar sensors of both satellites receive the radar echoes of those radar pulses after reflection on the Earth's ground. As a consequence, SAR image pairs at corresponding time instances are obtained in a single-pass of the satellites over the Earth's ground.

[0071] With a well-known procedure, a digital elevation map can be derived from a corresponding image pair. The digital elevation map comprises height values for a plurality of pixels at pixel positions representing corresponding points on the Earth within the terrain taken by the SAR images. In the embodiment described herein, a pair of digital elevation maps taken at different times from the same terrain on the Earth's ground are mutually calibrated. The mutually calibrated digital elevation maps enable so called DEM differencing where the digital elevation maps are compared in order to identify changes in the terrain, e.g. the melting of glaciers or the deforestation of rain forest.

[0072] As mentioned above, the digital elevation maps considered in the embodiment described herein are derived from a bistatic SAR. Nevertheless, the digital elevation maps may also be acquired by another type of SAR system, e.g. a monostatic SAR overflying in regular intervals the same terrain on the Earth. Based on a pair of SAR images taken from the corresponding terrain at different passes of the monostatic SAR, a digital elevation map can be derived.

[0073] Fig. 1 shows a flow chart of the embodiment described in the following. The method processes a pair of digital elevation maps DEMI and DEM2. Both digital elevation maps cover the same terrain TR on Earth. Digital elevation map DEMI comprises height values hvl for a plurality of pixels representing the terrain TR and digital elevation map DEM2 comprises height values hv2 for the plurality of pixels representing the terrain TR at a different time than map DEMI. Both digital elevation maps DEMI and DEM2 are each derived from a pair IP1 and IP2 of SAR images IM which are acquired by a first radar equipment RE1 which is the bistatic SAR of the TanDEM-X mission.

[0074] Each of the SAR images IM is characterized by a plurality of pixels covering the terrain TR, where each pixel has an amplitude value a in the phase value ph. The digital elevation map DEMI is derived from the image pair IP1 of the SAR images IM simultaneously acquired by the bistatic SAR at a corresponding point in time and the digital elevation map DEM2 is derived from the other image pair IP2 of the SAR images IM taken simultaneously by the bistatic SAR at a later point in time. The two digital elevation maps DEMI and DEM2 form the input of the calibration method described herein.

[0075] In a first step STI of this method, calibration points cp are determined by using a so called amplitude dispersion index D and a preset threshold th for this dispersion index. A calibration point cp refers to a permanent scatterer candidate within the predetermined terrain TR, i.e. to a point which has not changed from one digital elevation map to the other and thus is a good point for calibrating the digital elevation maps DEMI and DEM2. Each calibration point is represented by the same pixel position within the two digital elevation maps DEMI and DEM2.

[0076] For determining the dispersion index D, the method described herein processes a time series ts of SAR images TSIM of the terrain TR acquired by a second radar equipment RE2 being different from the first radar equipment RE1. Typically, the time series ts has at least 20 images TSIM. As for the SAR images IM, each of the SAR images TSIM includes a plurality of pixels covering the terrain TR where each pixel is characterized by an amplitude value a and a phase value ph. The second radar equipment RE2 is the monostatic SAR of the Sentinel-1 mission with two satellites flying on the same nominal orbit allowing for a repeat-pass acquisition of the same terrain on Earth. In this mission, the nominal acquisition plan foresees a revisit time of six days over Europe and selected test areas and of twelve days for the rest of the world. This allows for continuous acquisitions of InSAR time series at global scale.

[0077] Differently, though having the uniqueness of bistatic InSAR allowing for the generation of high-quality digital elevation maps, the TanDEM-X mission does not allow for the constant acquisition of long InSAR time series due to limited onboard resources. Therefore, in the embodiment described herein, the time series ts for determining the dispersion index D is taken from the Sentinel-1 mission. Nevertheless, in case that the first radar equipment RE1 provides besides the digital elevation maps a long time series of SAR images, this time series may also used for determining the dispersion index D. In this case, the first radar equipment RE1 coincides with the second radar equipment RE2.

[0078] In the following, step STI of the method of Fig. 1 is described in detail with respect to Fig. 2. Step STI comprises several sub-steps ST101 to ST107. In sub-steps ST101 to ST 104, the time series ts of SAR images TSIM acquired by the second radar equipment RE2 is processed. As a result of this processing, candidate calibration points ccp being candidates for becoming calibration points cp within the terrain TR are determined. Before starting the steps of Fig. 2, both TanDEM-X and Sentinel- 1 data are mapped to the same coordinate system and properly interpolated on the same grid. This is a well-known pre-processing step which is omitted in Fig. 2 for the sake of simplicity.

[0079] In steps ST101 to ST103, a pre-filtering of the pixel positions within the images TSIM is performed. According to this pre-filtering, pixel positions are identified which are not allowed to become candidate calibration points so that those pixel positions are disregarded in subsequent step ST104.

[0080] In step ST101, pixel positions within the SAR images TSIM are identified which correspond to points on the Earth's ground affected by geometric distortions GD. Geometric distortions are well-known effects in SAR data acquisition and result in an improper positioning of scatterers in the SAR image. Geometric distortions are caused by the side-looking nature of the SAR and comprise inter alia the well-known effects of layover and shadow. Corresponding techniques of classifying pixel positions in SAR images as points affected by geometric distortions are well-known for a skilled person and, thus, will not be described in detail. Step ST101 results in a reduced set of pixel positions no longer including points affected by geometric distortions.

[0081] In step ST 102, the reduced number of pixels is further reduced by the use of an external reference land cover map LCM indicating those pixel positions belonging to areas covered by volumetric material such as snow or ice. Points in those areas may vary in height and may be affected by radar penetration effects. Therefore, these pixel positions are discarded, resulting in a further reduced number of pixels. In step ST103, the reduced set of pixels resulting from step ST102 is reduced again by considering the slope si of the terrain at the corresponding pixel positions. Those pixel positions are discarded which have a slope si higher than a predetermined slope threshold sit. Pixel positions with a high slope result in inaccuracies and shall not be used as calibration points. The slope of the pixel positions is derived by using an external reference digital elevation map (e.g. the global TanDEM-X DEM product).

[0082] The set of pixel positions remaining within the images TSIM after having performed step ST103 are then analyzed based on the dispersion index D. The dispersion index is e.g. known from document [2] and is given by the following expression:

[0083] D = - (1) n

[0084] The dispersion index D is calculated for each pixel position remaining after step ST 103. For a respective pixel position, p represents the temporal mean of the amplitude values of the pixels at the respective pixel position over the images TSIM of the time series ts. Furthermore, o corresponds to the standard deviation of the statistical distribution of the amplitude values of the pixels at the respective pixel position over the images TSIM of the time series ts. Points on earth are considered to be stable over time if the corresponding pixel position has a dispersion index D lower than the preset threshold th. An adequate value for this threshold th may be chosen based on the specialist knowledge of a skilled person. E.g., a value of 0.25 may be used for the threshold. However, in an optional calculation procedure, the threshold th may also be automatically calculated, as will be described below with respect to Fig. 3.

[0085] Step ST 104 results in an amount of pixel positions being candidate calibration points ccp. From those 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 DEMI and DEM2 are derived. In step ST105, pixel positions being affected by geometric distortions GD are identified within the image pairs IP1 and IP2 in the same way as for the images TSIM in step ST101. This step assures that, neither in the Sentinel-1 acquisition geometry nor in the TanDEM-X acquisition geometry, finally selected elevation points are affected by geometric distortions. The pixel positions affected by geometric distortions are discarded from the candidate calibration points ccp.

[0086] In the step ST106, a threshold ct is set on the bistatic interferometric SAR coherence of both TanDEM-X bistatic acquisitions (i.e. both images IP1 and IP2) in order to assure a good quality of the TanDEM-X digital elevation map in correspondence of the selected targets. The interferometric coherence is a well-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 herein 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: and u2represent the images of the corresponding image pair IP1 or IP2. E is the expectation value and * is the complex conjugate operator. The interferometric coherence ytotquantifies the amount of noise in the interferogram and represents therefore a key quantity for assessing the quality of each digital elevation map. The interferometric coherence can be estimated for each pixel position over a number of adjacent pixels by well-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 larger, particularly a value between 0.6 and 0.7. Those pixel positions which have estimated interferometric coherences below the predefined threshold ct are discarded from the candidate calibration points ccp. The reduced number of candidate calibration points ccp resulting from step ST 106 is processed by another filtering step ST 107. In this step, the signal-to-noise ratio SNR of each pixel of the images of the image pairs IP1, IP2 is computed as follows:

[0087] The quantities in the above formula are well-known. (3° identifies the radar brightness and / 3nis the noise equivalent beta nought representing the SAR system noise floor. A candidate calibration point is discarded in case that the corresponding pixel position of the candidate calibration point belongs to a pixel within the image pairs having a signal-to-noise ratio SNR below a predefined threshold st. This assures that targets are used as calibration points which not only have sufficient amplitude stability but also sufficient phase stability.

[0088] The candidate calibration points ccp remaining after the filtering of steps ST 105 to ST 107 form the calibration points cp which are processed in the subsequent steps ST2 to ST4 of Fig. 1. In step ST2 of Fig. 1, height value differences A between the digital elevation maps DEMI and DEM2 are calculated for the calibration points cp. Those height value differences can be expressed by the following formula:

[0089] A = DEMl(cp) - DEM2(cp) (4)

[0090] Thereafter, in step ST3, the height value differences at the calibration points are used as data points for fitting a two-dimensional curve CS to the height value differences. The two-dimensional curve CS is referred as to the calibration surface in the following. Any known interpolation technique may be used for fitting the two-dimensional curve. The interpolation may be a linear interpolation or any other interpolation of higher order. As a result, the calibration surface CS describes a function providing a height value difference for each of the pixels within the digital elevation maps DEMI, DEM2. The calibration surface CS is used in the final step ST4 of Fig. 1 for correcting the digital elevation map DEM2 where calibration DEMI is left unchanged. The corrected digital elevation map DEM2' is calculated per pixel position within DEM2 based on the following formula:

[0091] DEM2' = DEM2 + Acs(5)

[0092] In the above formula, Acs- refers 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.

[0093] For the above correction, the digital elevation map DEMI is regarded as the reference digital elevation map. Nevertheless, digital elevation map DEM2 may also be regarded as the reference elevation map. In this case, the indices for the digital elevation maps in the above formulae (4) and (5) are interchanged.

[0094] The method of Fig. 1 results in two properly mutually calibrated digital elevation maps DEMI and DEM2' thanks to the use of permanent scatterer candidates as the calibration points. The calibrated digital elevation maps may be used for DEM differencing by identifying height differences between the digital elevation maps. This enables the determination of regions where the height of the terrain has changed over time, e.g. in order to detect the melting of glaciers or the deforestation of rain forest.

[0095] In an optional modification, the threshold th for the dispersion index D used in step STI may be calculated in a separate calculation procedure performed before the actual calibration method. This calculation procedure is illustrated in Fig. 3. The calculation procedure is used in case that the first radar equipment RE1 is another radar equipment than the second radar equipment RE2. This is the case in the embodiment described herein as the radar equipment RE1 refers to the TanDEM-X mission whereas the radar equipment RE2 refers to the Sentinel- 1 mission.

[0096] The method of Fig. 3 processes a first time series tsl of SAR images TSIM1 acquired by the radar equipment RE1 of the TanDEM-X mission. Furthermore, the method processes a second time series ts2 of SAR images TSIM2 acquired by the radar equipment RE2 of the Sentinel- 1 mission. All images cover the same predefined terrain TR' on the Earth's surface with the same resolution where each point in the terrain TR' corresponds to the same pixel position within the images. The predefined terrain TR' needs not be the same terrain as the terrain TR considered in the method of Fig. 1.

[0097] The method of Fig. 3 is based on the finding that the threshold th applied to the Sentinel-1 data shall be made compatible with the TanDEM-X data from which the digital elevation maps being calibrated are derived. Indeed, the data acquisitions of the Sentinel- 1 mission and the TanDEM-X mission have substantial differences. The Sentinel- 1 repeat-pass data are acquired at C-band whereas the bistatic TanDEM-X data are acquired at X-band. When considering radar sensors at different frequencies, a different behavior in terms of backscatter levels as well as interferometric performance, i.e. coherence, is expected.

[0098] The calculation procedure described in the following is an empirical approach developed for setting a proper threshold th on the dispersion index D computed for the Sentinel- 1 repeat-pass time series, which will assure the selection of stable calibration points when considering the same area acquired at X-band using TanDEM-X bistatic acquisitions.

[0099] In the method of Fig. 3, as first dispersion indices DI, the dispersion indices of the respective pixel positions within the SAR images TSIM1 of the first time series tsl referring to the TanDEM-X mission are determined based on the above formula (1). Analogously, as second dispersion indices D2, the dispersion indices of the respective pixel positions within the SAR images TSIM2 of the second time series ts2 referring to the Sentinel- 1 data are calculated based on the above formula (1). In a next step, the histogram HI depending on both dispersion indices DI and D2 is calculated. 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') having the values of dispersion indices DI and D2 equivalent 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 refers to a specific value of DI and the specific value of D2 where the histogram indicates the number of pixel positions having the specific value of DI and the specific value of D2.

[0100] For the first dispersion indices DI, a specific limit li is preset. In the embodiment described herein, this limit is chosen to be 0.25, analogously to document [2], Other values may be used for the limit li as well. Preferably, the limit lies in the range between 0.15 and 0.3. The limit li is indicated by a dash-dotted line in Fig. 3. In a next step, the threshold th for the dispersion indices D2 indicated by a dashed line in Fig.

[0101] 3 is chosen such that the number of pixel positions (i.e. points in the terrain TR') in the right lower quadrant limited by the lines of the limit li and the threshold th corresponds to a preset percentage pc of pixel positions (10% or less) within the total number of pixel positions in the respective images TSIM1, TSIM2 where this total number is the same for all images TSIM1 and TSIM2 and corresponds to the total number of points within the terrain TR'. In other words, the threshold th is chosen such that both the pixel positions within the respective SAR images TSIM1 of the first time series tsl having a first dispersion measure DI above the preset limit 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 th amount to the preset percentage.

[0102] The threshold th is chosen as a trade-off between the maximum acceptable rate of false alarms and the minimum number of selected targets to be maintained for the calibration. False alarms correspond to targets which, if selected, would not satisfy the stability requirement for the TanDEM-X data (first dispersion index DI above limit li) but would satisfy the stability requirement for the Sentinel-1 data (dispersion index D2 smaller than threshold th). In the embodiment described herein, the percentage pc is chosen to be 1%. This results in a value of the threshold th at approximately 0.19. Preferably, the percentage pc is chosen to be 5% or less.

[0103] The method as described in the foregoing has been evaluated by the inventors based on corresponding data of the TanDEM-X and the Sentinel- 1 missions. The terrain being considered for the digital elevation maps refers to a sandy / rocky area in the Sahara Desert, Egypt. The height differences between the uncalibrated digital elevation maps and the calibrated digital elevation maps were determined, and it was found that residual tilts are correctly compensated after applying the calibration method whereas the residual tilts are still visible without mutual calibration.

[0104] The invention as described in the foregoing has several advantages. Particularly, suitable calibration points for a pair of digital elevation maps acquired at different time points are extracted from SAR data without the use of external reference measurements. The calibration points are chosen to be permanent scatterer candidates by considering the amplitude dispersion index of corresponding pixel positions. For determining the dispersion index, the method may use data from SAR acquisitions other than the SAR acquisition from which the digital elevation maps to be calibrated are derived. In an optional variant, the method also provides a calculation procedure for the threshold to be applied to the dispersion index for selecting candidate calibration points.

[0105] List of references

[0106] [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), p. 39-49.

[0107] [2] A. Ferretti, C. Prati, F. Rocca: "Permanent Scatterers in SAR Interferometry", IEEE Transactions on Geoscience and Remote Sensing, vol. 39, no. 1, Jan. 2001.

[0108] [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, vol. 15, no. 1, p.954-974.

[0109] [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, vol. 45, no. 11, pp. 3317-3341, Oct 2007.

Claims

Claims1. A computer-implemented method for calibrating a pair of digital elevation maps (DEMI, DEM2), where the digital elevation maps (DEMI, DEM2) of the pair indicate the height of a predetermined terrain (TR) on a surface of a celestial body and comprise height values (hvl, hv2) for a plurality of pixels at pixel positions representing corresponding points within the predetermined terrain (TR), where one digital elevation map (DEMI, DEM2) of the pair indicates the height of the predetermined terrain (TR) earlier in time than the other digital elevation map (DEMI, DEM2) of the pair, where each digital elevation map (DEMI, DEM2) of the pair is associated with another SAR image pair (IP1, IP2) out of SAR images (IM) acquired from measurements of a first radar sensor equipment (RE1), where a respective digital elevation map (DEMI, DEM2) has been derived from phase differences between the associated SAR image pair (IP1, IP2), wherein the method processes a time series (ts) of SAR images (TSIM) of the predetermined terrain (TR) acquired from measurements of a second radar sensor equipment (RE2), where each SAR image (TSIM) of the time series (ts) comprises pixels at the same pixel positions as the digital elevation maps (DEMI, DEM2) and as the SAR images (IM) acquired from measurements of the first radar sensor equipment (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), the dispersion measure (D) of a respective pixel position describing the statistical dispersion of amplitude values (a) of the pixels at the respective pixel position over the SAR images (TSIM) of the time series (ts), where each candidate calibration point (ccp) refers to a pixel position having a dispersion measure (D) below a predetermined threshold (th); b) identifying pixel positions as calibration points (cp) within the digital elevation maps (DEMI, DEM2) of the pair based on the candidate calibrationpoints (ccp), where each calibration point (cp) corresponds to a pixel position of a candidate calibration point (ccp); c) determining height value differences (A) between the digital elevation maps (DEMI, DEM2) of the pair at the calibration points (cp); d) fitting a two-dimensional curve (CS) to the height value differences (A) at the calibration points (cp), the two-dimensional curve (CS) describing the height value differences (A) over the pixel positions representing corresponding points within the predetermined terrain (TR); e) correcting the height values (hvl, hv2) of one digital elevation map (DEMI, DEM2) of the pair by the height value differences (A) of the two- dimensional curve (CS).

2. The 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 being 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 value of this statistical distribution.

3. The method according to claim 1 or 2, wherein in step a), in case that a pixel position within the SAR images (TSIM) of the time series (ts) is identified as being affected by geometric distortions (GD), this pixel position is classified as a pixel position not permitted to be a candidate calibration point (ccp).

4. The method according to one of the preceding claims, wherein in step a), in case that a pixel position within the SAR images (TSIM) of the time series (ts) is identified as being a point within the predetermined terrain (TR) covered by snow or ice, this pixel position is classified as a pixel position not permitted to be a candidate calibration point (ccp).

5. The method according to one of the preceding claims, wherein in step a), in case that a pixel position within the SAR images (TSIM) of the time series (ts) is identified as being a point within the predetermined terrain (TR) having a slope (si) higher than a predetermined slope threshold (sit), this pixel position is classified as a pixel position not permitted to be a candidate calibration point (cep).

6. The method according to one of the preceding claims, wherein in step b), in case that a pixel position within the digital elevation maps (DEMI, DEM2) is identified as being affected by geometric distortions (GD), this pixel position is classified as a pixel position not permitted to be a calibration point (cp).

7. The method according to one of the preceding claims, wherein in step b), in case that the interferometric coherence (ytot) forapixel position of a SAR image pair (IP1, IP2) associated with one of the digital elevation maps (DEMI, DEM2) of the pair is below a coherence threshold (ct), this pixel position is classified as a pixel position not permitted to be a calibration point (cp).

8. The method according to one of the preceding claims, wherein in step b), in case that that 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 (DEMI, DEM2) of the pair is below a sig- nal-to-noise ratio threshold (st), this pixel position is classified as a pixel position not permitted to be a calibration point (cp).

9. The method according to one of the preceding claims, wherein the first radar equipment (RE1) is another radar equipment than the second radar equipment (RE2).

10. The method according to one of the preceding claims, wherein the first radar equipment (RE1) is a monostatic radar equipment with one or more radarsensors, where each radar sensor transmits its own radar pulses and exclusively measures the radar echoes of its own radar pulses, or a bistatic radar equipment with two radar sensors separated by a predetermined distance, where 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 equipment (RE2) is a monostatic radar equipment with one or more radar sensors, where each radar sensor transmits its own radar pulses and exclusively measures the radar echoes of its own radar pulses, or a bistatic radar equipment with two radar sensors separated by a predetermined distance, where only one radar sensor transmits radar pulses and both radar sensors measure the radar echoes of the radar pulses.

11. The method according to one of the preceding claims, wherein the method comprises a calculation procedure for calculating the predetermined threshold (th) based on first and second time series (tsl, ts2) of SAR images (TSIM1, TSIM2) of a predefined terrain (TR1) on the surface of the celestial body, where the first time series (tsl) has been acquired from measurements of the first radar sensor equipment (RE1) and the second time series (ts2) has been acquired from measurements of the second radar sensor equipment (RE2), where all SAR images of the first and second time series (tsl, ts2) comprise pixels at the same pixel positions representing corresponding points within the predefined terrain (TR1), where the first radar equipment (RE1) preferably is another radar equipment than the second radar equipment (RE2), wherein the calculation procedure comprises the following steps: i) determining, as first dispersion measures (DI), dispersion measures of respective pixel positions within the SAR images (TSIM1) of the first time series (tsl), the dispersion measure of a respective pixel position describing the statistical dispersion 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, as second dispersion measures (D2), dispersion measures of respective pixel positions within the SAR images (TSIM2) of the second time series (ts2), the dispersion measure of a respective pixel position describing the statistical dispersion 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 (th) based on a preset limit (li) for the first dispersion measures (DI) by choosing the predetermined threshold (th) such that both the pixel positions within the respective SAR images (TSIM1) of the first time series (tsl) having a first dispersion measure (DI) above the preset limit (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 (th) amount to a preset percentage of 10% or less.

12. The method according to claim 11, wherein the first dispersion measure (DI) of a respective pixel position within the SAR images (TSIM1) of the first time series (tsl) represents or depends on a dispersion index being the standard deviation of the statistical distribution of the amplitude values (a) of the pixels at the respective pixel position over the SAR images (TSIM1) of the first time series (tsl) 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 being the standard deviation of the statistical distribution of the amplitude values (a) of the pixels at the respective pixel position over the SAR images (TSIM2) of the second time series (ts2) divided by the mean value of this statistical distribution.

13. The method according to claim 11 or 12, wherein the preset limit (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, preferable 1%.

14. A computer program product with program code, which is stored on a machine-readable carrier, for carrying out a method according to one of the preceding claims when the program code is executed on a computer.

15. A computer program with program code for carrying out a method according to one of claims 1 to 13 when the program code is executed on a computer.

Images