Possibility of using the DInSAR method in the development of vertical crustal movements with Sentinel-1 data

1 Introduction

Several maps of the vertical crustal movements (VCMs) of the Earth’s crust were developed in Poland: Wyrzykowski in 1987 [1], Kowalczyk in 2006 [2], Kontny and Bogusz in 2012 [3], Kowalczyk in 2015 [4], and Kowalczyk and Rapiński in 2017 [5]. For this purpose, data from repeated precision leveling campaigns [6] and data from permanent GNSS stations were used. Both types of sets are characterized by the fact that they are movements determined at specific points [7]. However, they differ in the reference level, method of development, frequency of measurements, etc. In addition to the above measurement techniques, data obtained from aerial and satellite measurement systems are also used to determine changes in the Earth’s surface [8,9]. The development of surface measurement methods such as Airborne laser scanning [10] and satellite data, e.g., synthetic aperture radar (SAR) [11] makes it possible to analyze changes in the surface of the Earth’s crust on large areas at the same time. Mutual measurements of the same area using different techniques allow not only to jointly develop models of displacement or deformation but also to verify the obtained results [12]. The advantage of this solution is the possibility of the simultaneous alignment of several networks created based on data obtained using various techniques called hybrid networks –e.g., leveling and GNSS [13]. For this purpose, the data must be reliable. It cannot differ significantly from the other sets. Common points are identified (for linking networks created based on sets obtained with different measurement techniques) and specific main points which are the reference level with the simultaneous alignment of hybrid networks. A necessary condition for the implementation of this solution is the possession of tie points in these networks [6].

VCM models are developed as point or surface models, relative or absolute. In the case of leveling data, the network is adjusted [14,15,16] and the movements are interpolated [1,2,13] in relation to the adopted point or points and reference level. For GNSS data, time series analyses on permanent stations are used, calculated based on a network solution [17] or in the precise point positioning (PPP) technology [18,19,20,21]. One of the disadvantages of the network solution is the computational difficulty that appears with the increase in the number of processed stations [20] and difficulties in the case of increasing distances between stations. In the case of the PPP strategy, there are no distance constraints, and the estimated position is determined independently for a single station, directly in relation to the established frame of reference [22]. The literature on the development of GNSS time series is very rich and it solves in detail important problems that arise at the stage of determining the trend. The time series of coordinates obtained from PPP technology are not affected by changes taking place at neighboring stations. However, individual stations may be affected by global geophysical processes, mainly due to tidal effects [23,24,25] or local changes such as atmospheric loads [26,27]. This is important when we are modeling VCM using data from GNSS stations developed in PPP technology.

Since SAR data are more commonly used to determine terrain deformation [28,29,30], in this article, we decided to analyze the quality of SAR data prepared by the Differential interferometric synthetic aperture radar (DInSAR) method to validate and be a reference to the obtained GNSS data. We did it on the micro areas (approx. 10 km circles) around four Polish permanent GNSS stations included in the International GNSS service (IGS) as potential main points of the hybrid network. An attempt was also made to assess the quality of the DInSAR data in terms of its applicability to create a vertical movement network that could potentially make another element of such a network.

The main advantage of SAR data is the acquisition rate, spatial resolution, and accuracy. The time of the revisit, i.e., the appearance in the same place above the Earth and the beginning of the next cycle of orbital around it, for one satellite, depending on the mission, ranges from 6 to 12 days. The surveyed areas also vary depending on the mission specifications. The Sentinel-1 satellites selected for this study have different spatial (terrain) resolutions depending on the radar data acquisition mode. In this work, the Interferometric wide (IW) swath mode was used, which is a 240 km belt with an assumed terrain resolution of 5 m × 20 m based on information provided by the European Space Agency (ESA) [31].

The work aimed to acquire basic knowledge in the field of:

• suitability of the DInSAR method for monitoring VCM/displacement of the Earth’s crust,

• the impact of the topology of a given area (in this case, the surroundings of active geodetic network-European position determination system (ASG-Eupos) stations) on the results and accuracy of SAR data preparation,

• investigating the relationship between the displacements obtained from GNSS data and those obtained with the radar interferometry technology,

• accuracy and validity of creating hybrid networks from GNSS and SAR data.

With this article, we not only wanted to demonstrate the usefulness of SAR methods in large areas but also to demonstrate the importance of time series and hybrid networks in the VCM study.

2 Data used

In Poland, a geodetic (horizontal and vertical) and gravimetric network are used to monitor land displacements and assist in determining corrections to coordinate systems. In the process of designing the network, the above scientific needs and tectonic structures are taken into account, which results in a different type of stabilization as well as various environmental and anthropogenic conditions around the network points [32].

In the presented study, the analysis used 4 points of the fundamental basic height network as one of the most important main points [32,33]. These points are ASG-EUPOS, EUREF permanent GNSS network (EPN), and IGS stations: Borowiec (BOR1, Poland), Borowa Góra (BOGI, Poland), Lamkówko (LAMA, Poland), and Józefosław (JOZ2, Poland). These points are placed on solid rock, on top of buildings or on a pole mounted on a foundation footing sunk below the ground at freezing level [32]. At each point, there is a GNSS station with a different type of antenna at a different height (on a pole or the roof of the building) (Figure 1).

Figure 1

Locations of antennas at the Borowiec, Borowa Góra, Lamkówko, and Józefosław stations (source: www.asgeupos.pl, acquisition date: 29/12/2020).

Additionally, each GNSS station is located on a different tectonic and geological basis, with different water, anthropogenic, and soil conditions. This allows us to analyze how these factors affect the DInSAR results. Characteristics of the stations are presented in Tables 1 and 2.

Table 2

Characteristics of the vicinity of the stations used: urbanization, geology, hydrology, and DSM and mining areas near four selected ASG-EUPOS stations (source: own study)

Station	Urbanization	Geology	Hydrology	DSM and mining terrain
BOR1
BOGI
LAMA
JOZ2

The vicinity of the BOR1 station is characterized by little urbanization, it is located by the expressway. Most of the surrounding area is covered with sand and gravel of various kinds on deep clay soil. The station and its vicinity are located in the area of two overlapping underground water reservoirs. Most of the terrain is with an open horizon (based on the Digital Surface Model [DSM] analysis). At a distance of about 6 km from the station, there are two large natural gas deposits, marked in grey in Table 2.

The area of the BOGI station is characterized by moderate urbanization, it is located in the suburbs of Warsaw. In its vicinity, there are low-season cultivated plants and a national road. Most of the surrounding area is covered with various kinds of sands. In addition to the proximity of the Narew River, there is also a large underground water reservoir. The analysis of the DSM showed mostly flat terrain, up to the bank of the Narew River (exception: fortifications of the Zegrze Fortress). About 6 km from the station, on the other side of the river, there is a large deposit of quartz sand and sand-lime brick, marked grey in Table 2.

The vicinity of the LAMA station is characterized by almost zero urbanization, but almost complete forestation. Much of the surrounding area is covered with various kinds of clay and peat. In the vicinity of the station, there are two large underground water reservoirs and several larger lakes (Wadąg, Dobrag, and Kiermas). Additionally, there are vast areas of sand and gravel open-cast mines in the vicinity of the station. When analyzing the DSM, large differences in land cover around the station were found (up to 30 m). At a distance of nearly 2 km from the station, there are two large, active natural aggregate mines, marked in grey in Table 2.

The area of the JOZ2 station is highly urbanized, it is located in the suburbs of Warsaw. Most of the surrounding area is covered with sand and clay of various kinds. Despite the proximity of the Vistula River, there are no underground water reservoirs there. When analyzing the DSM, slight differences in the land cover around the station were found. At a distance of about 6 km from the station, there is a large deposit of curative waters (the so-called brines), marked in grey in Table 2.

The research used the available data from the Sentinel-1 mission, the first of the two satellites (S1-A) launched on April 3, 2014, and the second (S1-B) launched on April 22, 2016. This radar emits waves in the C band with a frequency of 5.3 GHz (a wavelength of approx. 5.6 cm) [31].

For this analysis, 24 radar images from the Sentinel-1A and 1B missions from the ESA [31] were used as Open Access data. The images used in this study are single look complex (SLC) products. The time interval between the acquisition of satellite imagery varies depending on the station, but all pairs were composed based on images from 2020. Four ASG-EUPOS network stations were selected for comparison (BOR1, BOGI, LAMA, and JOZ2) along with the surrounding area, about 10 km in diameter [6,34]. This makes it possible to trace the behavior of the earth’s crust around the station and to check the effect of land cover on the coherence of the photos [35]. At the JOZ2 station, the foundation of the antenna does not coincide with the point of the fundamental network, they are separated from each other by about 84 m, which was taken into account when reading the results.

3 Reference data

GNSS PPP time series have been used as reference data and the vertical velocities of BOR1, BOGI, LAMA, and JOZ2 stations were estimated. Raw GNSS observations from individual stations were obtained as daily data in the compressed RINEX Hatanaka format provided by ASG-EUPOS. RINEX files for these stations were made available for this work by the Central Office of Geodesy and Cartography in Warsaw. The scope of data covered the years 2010–2021 which gives over 4,000 days of observations. The raw GNSS data have been processed in PPP mode with the use of GipsyX Software developed by the Jet Propulsion Laboratory (JPL). To ensure centimeter accuracy, the PPP method requires access to several, high-quality, external corrections to mitigate errors like satellite orbit errors, tropospheric and ionospheric delay, receiver-specific errors, and errors related to the geophysical Earth’s models [20]. In the presented work, the GNSS products (antenna models, final ephemerides clocks, etc.) have been downloaded from the JPL server: https://sideshow.jpl.nasa.gov/pub/JPL_GNSS_Products/Final. Additionally, there are complex models of geometric effects and models of force models for Earth orbits as well as antenna models and Global mapping functions to calculate tropospheric delay implemented in the GipsyX package. Subsequently, the obtained time series have been statistically analyzed to estimate vertical movement. Calculations of daily positions using the PPP method and the development of time series to determine the speed of the station from the GNSS data were performed in the GipsyX program, in the ITRF14 system, which included breaks detection, seasonality determination, removing outliers, and finally trend estimation.

The Nevada Geodetic Laboratory (NGL) reference data were also calculated with the PPP method [36]. Data were obtained from http://geodesy.unr.edu/index.php for BOR1 (1994–2022), BOGI (2001–2022), LAMA (1994–2022), and JOZ2 (2002–2022).

4 Methods

The technology of radar and its use to detect objects date back to the 1870s when James Clerk Maxwell began researching electromagnetic waves and Heinrich Hertz created and emitted the first radio waves. Hertz also discovered the phenomenon of scattering and reflection of radio waves. Thanks to this, in 1903 the first radar was constructed that could detect ships and thus improve safety in maritime transport. Later, radars became an important military technology, specially developed before the outbreak of WW II. The first operational satellite radar device was launched into orbit in 1972 thanks to JPL [37]. The possibility of obtaining the first satellite radar data for civilian studies did not begin until 1978 with the launch of the Seasat satellite. Unfortunately, its work did not last long, because there was a power failure 5 months after the start of the mission. However, this satellite paved the way for further, more elaborate designs for the commercial use of radar interferometry [31].

Radar interferometry is used for various types of environmental monitoring: from determining damage caused by earthquakes [38,39] to forest monitoring [40] and, inter alia, checking ground displacements and deformations [41].

Identification of the object in the image depends on its actual size, the field resolution of the sensor and whether the object is scattering the radar wave. The traditional SAR interferometry (InSAR) technique is based on the use of interferograms, i.e., raster sets showing the phase difference between two SAR registrations. InSAR uses SAR registrations, performed sequentially with one SAR antenna during successive flights of the satellite over the same area (the so-called Differential InSAR [DInSAR]). Based on the phase differences of the corresponding radar signals from successive SAR images, information is obtained about the relative values of the elevation of the terrain surface or its changes over time [39,42,43]. This method was chosen to test its suitability for monitoring and VCM study in large areas.

SAR technology is based on electromagnetic (radar) waves of different bands (wavelengths) and measurement of the return time of the reflected wave that is needed to travel from the satellite to a specific point on the Earth’s surface. Each wave has a sinusoidal pattern that has explicit minimums and maximums that fall between –π and +π. The signal reflected from the Earth’s surface is registered in a complex form, it maintains its phase and amplitude [37]. The interferometric processing of SAR data calculates the phase differences between the two SAR images. The result of this processing is an image of the phase differences, called an interferogram. Interferograms can be visualized as color sequences corresponding to phase differences, called interferometric fringes [44]. To perform the analysis, at least two photos are needed, based on which the geometric relations between the reference (Master) and the adjusted (Slave) images are examined.

The phase reflected from a single point on the Earth’s surface can be described by the formula (1) [37]

where φ – phase value, R _p  – the range between the radar and a point on the Earth’s surface, λ – radar wavelength, and φ _scat  – the scattering phase contribution which is related to the target’s electrical properties.

The phase differences represented in the interferograms are the sum of many factors. The most important of them are the geometric component (terrain topography), deformation component (vertical and horizontal movements between data acquisition), atmospheric conditions, snow cover, the dielectric constant of the ground, and spatial distance between satellites taking a given image. The value of the change the between photos is made up of relevant factors [37,44,45] (equation (2))

where Δφ – interferometric phase (or phase change between SAR acquisitions), φ _flat  – flat Earth phase, φ _topo  – topographic phase contribution, φ _obrit  – the phase error induced by errors in orbit information, φ _defo  – phase contribution related to ground deformation, φ _tropo  – tropospheric phase contribution, φ _iono  – ionospheric phase contribution, φ _scat  – phase contribution related to the scatterer’s electrical properties, and φ _noise  – the combined noise term.

Depending on what we want to investigate, some of the above values may be considered “noise” or negated as insignificant [44]. This study focuses primarily on the usefulness of the DInSAR method for monitoring vertical ground displacements (φ _defo). Equation (2) can be expanded as follows:

where λ – wavelength, B _║  – value of the parallel distance base for each image pixel, Δz – change in the height of the terrain in relation to the ellipsoid, B _┴  – value of the perpendicular distance base for each image pixel, R – range between target and satellite, Θ – antenna look angle, ΔR _defo  – range change related to terrain deformation, ΔR _tropo  – range change related to the tropospheric delay of the wave, N – number of electrons per square meter of space on a given day.

The quality of the images is most influenced by the topography (φ _topo) and the dielectric constant of the ground (φ _scat). The phenomenon of obscuring objects or their apparent stretching and shortening in the photos depends on the topography of the terrain. The influence of topography must be eliminated with each deformation analysis we perform. The dielectric constant of the substrate determines how high the water content is on the tested soil. This value is defined as a percentage and described as the volumetric water content. Values close to 0% (e.g., dry sand) are as unfavorable as values close to 100% (e.g., water reservoir area). The water content in the ground determines the force with which the wave will be reflected from the Earth’s surface to the satellite, and what part of it will be absorbed or reflected, causing the radar wave to scatter [46,47,48].

DInSAR is a method that was used in this work to learn about its basics and how to use it in practice. Among all measurement methods, DInSAR is widely used to assess deformation in mining areas and tectonically active areas or to monitor the effects of sudden events such as earthquakes or volcanic eruptions [49,50]. Comparisons with measurement results using other techniques, such as precision leveling or GNSS measurements, show that there are differences between these measurement methods and that the obtained results are not unambiguous. This suggests that the SAR data developed using the DInSAR method may provide additional support in the development of VCM models in Poland, but cannot be the only source of such results. The following analyzes were performed to verify this.

Freeware SNAP v8.0.9 from ESA along with SNAPHU Unwrapping v8.0.0 and the Linux Ubuntu 18.04.6 LTS operating system were used for the analysis and processing of the images.

The methodology of the study was based on the “Sentinel-1 Toolbox Interferometry Tutorial” [51] and our own experiences. Due to the use of SLC photos from the Sentinel-1 mission for this study, the activities shown in Figure 2 were performed.

Figure 2

Diagram of the methodology and research process based on the “Sentinel-1 Toolbox Interferometry Tutorial” and own experiences (source: own elaboration).

After analyzing the data availability and performing image compatibility tests, identifiable time intervals were determined and pairs of photos were collected for 4 stations along with the marking of ascending and descending orbits (Table A1 in Appendix nr 1).

In the next step, orbital corrections were taken and automatic photo correction was performed (full images). Interferograms were made with the removal of the influence of the terrain topography on the results of the studies and with the calculation of the coherence of the images. An automatic subtract topographic phase was performed with the aid of the auto-sampled SRTM 1Sec HGT. Following equation (4), the correlation coefficient of the interferometric phase of a pair of SAR images was calculated – coherence.

Coherence is an important indicator of the reliability of the photos. It has values from 0 to 1, where 1 is the value of the perfect fit of both photos. A coherence of 0 indicates a decorrelation of photos. The coherence value consists of

where γ _total – total correlation (interferometric coherence), γ _spatial – spatial baseline decorrelation, γ _Doppler – Doppler centroid decorrelation, γ _temporal – temporal decorrelation and γ _thermal – thermal decorrelation.

The spatial base correlation is related to the horizontal spacing between the two orbits of the satellites. Doppler centroid decorrelation occurs when the orientation of a satellite changes (e.g., off course, slope, or orbital altitude) between the Master and Slave images. Temporal correlation occurs when the physical properties of scatterers on the Earth’s surface change. This form of decorrelation prevails in heavily forested areas. Thermal decorrelation, also called noise, is in most cases negligible when using interferometry [44]. Then, the gaps in the bands were removed (using the TOPSAR-Deburst tool) and saved in the source data and the images (Subset) were trimmed in the area of about 10 km around the station (Table 3).

Table 4

Results of unpacking the phase from the ascending and descending orbit, at 4 ASG-EUPOS stations, with a time base of 24 days

Time base 24 days	BOR1	BOGI	LAMA	JOZ2
Ascending
Descending

Cropped images were subjected to noise reduction (Multilooking) and filtering (Goldstein phase filtering). The same parameters were used for all pairs of images.

“Unpacking,” also known as “unfolding” a phase, was done with the SNAP add-on, “SNAPHU Unwrapping,” version 8.0.0. Phase unpacking is a process whose final effect is influenced by several parameters (number of columns and rows, coverage between them, and unpacking method). The selection of these parameters is crucial; therefore, some additional tests and trials were performed before the final phase unpacking was done. Independent tests (different sets of parameters) were performed for each image. Ultimately, it was decided to use common settings (Statistical cost mode as Deformation mode [DEFO]) in all images. Both Minimum Cost Flow (MCF) and Minimum Spanning Tree algorithms [52] were tested as an unpacking method. The number of rows and columns and their mutual coverage were tested in numerous variants: 1 × 1_0_0; 2 × 2_0_0; 2 × 2_200_200; 10 × 10_0_0; 10 × 10_100_100; 20 × 20_0_0; 20 × 20_50_50; A × B_C_D, where A is the number of columns, B is the number of rows, C is the mutual coverage of rows, and D is the mutual coverage columns.

In the next stage of data processing, the phase values were converted into displacement values. A layer with coherence values has been added to the images. Range-Doppler terrain correction was used to compensate for distortions caused by the incidence angles of radio waves other than the nadir angle, terrain topography, and satellite tilt [53]

Low coherence may disturb the interpretation of the results and show erroneous values [35]. Therefore, a filter mask for points with a coherence ≥0.2 [37] was added to exclude values with the lowest likelihood. The next process was to read the displacement values of the terrain around and in the station in the direction of the satellite’s Line of Sight (LOS). Table 1 shows the values of the station coordinates from which the displacement values were taken. Due to the ambiguity and diversity of the displacement values in space, it was decided to use an additional proof test. The test consisted in determining the value of the displacement from several neighboring pixels around the coordinates of the station (a square of 5 × 5 pixels – area of about 25 m × 100 m). The displacement values were determined using the weighted average

where defo_station – deformation value in a given period of time at a given station, pix₁ – deformation value of the first-left top pixel, coh₁ – first-left top pixel coherence value, pix _n  – deformation value of the n pixel, coh _n  – coherence value of the n pixel, Σcoh – the sum of the coherence values in the 5 × 5 pixel area around the station.

To obtain the value of the vertical displacement, instead of the displacement in the LOS direction, the obtained values were converted [54]

where d _los  – deformation value in the LOS direction, θ – the angle of the “view” of the satellite, Δα – the difference in the value of the direction of the satellite movement between the ascending and descending orbit, d _up – vertical displacement value, d _hald – projection of the horizontal deformation decreasing towards the azimuth.

The obtained results were saved in geo-referenced files in the .kmz and .geotiff formats.

5 Results

Figures 3 and 4 show the examples of identified differences in the unpacked phase for photos around the BOR1 station.

Figure 3

Differences between the parameters for unpacking photos of the area around the BOR1 station. On the left photo, 1 × 1_0_0 and on the right, 10 × 10_0_0 (source: own elaboration).

Figure 4

Differences between the parameters for unpacking photos of the area around the BOR1 station. On the left a photo, 20 × 20_0_0 and on the right, 20 × 20_50_50 (source: own elaboration).

Finally, the following variant was used for all photos: DEFO MCF 20 × 20_0_0. In Figure 5a and b, sample results of unpacking images are shown:

Figure 5

(a) Interferogram and the unpacked aASCBOR1 phase and (b) interferogram and the unpacked aASCLAMA phase.

The phase values converted into displacement for the variant of 24 days interval, without the coherence mask (coherence ≥ 0.2), are shown in Table 4.

VV (Vertical-Vertical) polarization, IW swath sensor mode and product type, and SLC were used for all pairs of satellite images. The results of comparison and parameters for unpacking individual pairs of photos and measurements at individual stations are presented in Tables 6–8, Table A1 and Figures 6–9. Additionally, Table 9 presents comparison of the results from the DInSAR method and the DInSAR 5 × 5 pixels (an area of about 25 m × 100 m in the field), with the GNSS results prepared using the PPP method from own calculations (UWM in Olsztyn) and the results from the NGL service.

Figure 6

Results of analyzes at the BOR1 station in the direction of LOS – graphic representation.

Figure 7

Results of analyzes at the BOGI station in the direction of LOS – graphic representation.

Figure 8

Analysis results at the LAMA station in the direction of LOS – graphic representation.

Figure 9

Results of analyzes at the JOZ2 station in the direction of LOS – graphic representation.

From UWM PPP, the calculated average annual vertical velocities and the accuracy of their determination, respectively, are as follows- BOR1: 0.25 mm/year with ±0.01 mm/year error, BOGI: 0.31 mm/year with ±0.02 mm/year error, LAMA 0.53 mm/year with ±0.02 mm/year error, and JOZ2: −0.47 mm/year with ±0.02 mm/year error.

From NGL PPP the calculated average annual vertical velocities and the accuracy of their determination, respectively, are as follows: BOR1: −0.48 mm/year with ±0.42 mm/year error, BOGI: 0.64 mm/year with ±0.61 mm/year error, LAMA −0.28 mm/year with ±0.48 mm/year error, and JOZ2: −0.54 mm/year with ±0.52 mm/year error.

The differences in the vertical displacement determined by the three methods are very diverse (Table 9). With a time base of 336 days, the differences between SAR and PPP (both UWM and NGL) range from about 20 mm (BOGI) to about 170 mm (BOR1). With the 24 days time base, these differences are smaller, from approx. −25 mm (LAMA) to approx. 60 mm (BOGI). With the 12 days time base, the values are even smaller: from about 7 mm (BOR1) to about 55 mm (BOGI) with one exception – JOZ2 with a value of 161 mm. Statistically smaller differences occur with shorter time bases of images, which results from the better matching of photos and smaller differences that may have occurred in a given area.

The averaged values on the 5 × 5 pixel area around the station (about 25 m × 100 m in the field) do not differ significantly from the pixel value at the station location. With a time base of 336 days, the differences are the highest. They range from approx. 2 mm (BOGI and JOZ2) to approx. 8 mm (BOR1 and LAMA). With the 24 days time base, the values of the differences are significantly smaller: from approx. 0.2 mm (BOR1) to approx. 0.4 mm (LAMA and JOZ2). With the 12 days time base, the values of the differences range from approx. 0.2 mm (BOGI) to approx. 5 mm (LAMA).

The overall coherence between the satellite imagery at co-registration ranged from 0.62 (aJOZ2 and aBOGI) to 0.98 (cBOR1, cBOGI, and cJOZ2). The coherence of photos is closely correlated with the time base. The longer the time base between views, the lower the coherence. Pixel coherence ranged from 0.24 (cLAMA) to 0.82 (aBOGI). The coherence in the 5 × 5 pixel square is lower in most cases, with slight improvement for cBOGI (from 0.80 to 0.83) and bLAMA (from 0.38 to 0.42). The coherence for a 5 × 5 square ranges from 0.21 (aLAMA) to 0.83 (cBOGI).

The biggest differences in the DInSAR method can be observed between the ascending and descending orbits, about 270 mm (aBOR1ASC and aBOR1DSC). In addition, the values of the calculated displacements usually have reverse signs, which is also confirmed by the results at other stations (Tables 5–8 and Figures 6–9). This may be because the ascending and descending orbits in some cases came from a Sentinel 1A or 1B. The lowest value of the spatial base is on aBOR1ASC and was 5.87 m, and the highest bBOGIDSC was −140.32 m.

There were also differences between the PPP results calculated by UWM and NGL. UWM calculations were made for the 2010–2021 time base with high accuracy, and NGL data were from different, longer time bases but with higher errors. Both PPP methods turned out to be very useful when verifying the results of the DInSAR method.

6 Discussion

When determining the VCM from GNSS data, the elimination of noise is an important element [55]. As a rule, noise is the result of a large discrepancy between neighboring values. As shown by the results, the discrepancy in displacement values differs between the SAR and PPP results. The displacement values from the DInSAR method are sometimes several times greater than the PPP methods. Taking into account the time difference between the pairs of images allows us to determine the vertical movement at the station using the approximate method. As shown by the results [2,3,20], movements at the four selected stations should be at the level of a few millimeters per year. The obtained values of displacements only at the stations JOZ2 and BOGI allow us to assume that millimeter values can be obtained. For the rest, it is up from a few to several centimeters. Such differentiation does not give credibility to the obtained results from the point of view of VCM modeling. This will only work for areas with proven large vertical changes in the ground surface at short intervals (earthquakes, landslides, mining subsidence, volcanic eruptions, etc.) [38,49,56,57]. The changes that follow such events on the Earth’s surface are significant and very easily marked. On the other hand, monitoring the fine vertical movements of the earth’s crust with this method is more difficult and would require the use of a time series of images to monitor these changes with the desired accuracy. Difficulties with the selection of a pair of images (poor coherence) show that for the areas of Poland, obtaining more reliable SAR data with a fixed time interval for selected areas (e.g., over 100 locations of the Polish ASG-EUPOS network) is very labor-intensive. It may be assumed that in some locations it is impossible. The best coherence and the smallest difference in displacement values were obtained at the BOGI station. This confirms the stability of this station and the selection of the right location. The applied PPP method [20] and the use of two independent datasets proved to be an element of verification of displacements determined from SAR data. The determined vertical displacements in the micro-areas of 10 km around the stations are used to show that these changes are not the same, both in the directions of displacement and distance from the GNSS stations. Therefore, for this type of analysis, methods based on the generation of time series that are simultaneously easily controllable by supporting the GNSS results are more suitable. Such methods include, for example, persistent scatterer InSAR (PSInSAR) [58,59] or short baseline subsets [60,61,62]. Creating hybrid networks for this kind of analysis seems to be recommended due to the possibility of simultaneous verification of each method. It is also recommended because it can give us a wider spectrum of factors that contribute to final VCM values.

7 Conclusion

Taking into account the above results and the entire data processing, the following conclusions were obtained. The DInSAR method is used primarily to detect and monitor changes occurring in short time intervals, therefore the research aimed to investigate the possibility of using it to develop vertical movements of the Earth’s crust. VCM are developed using geological methods (long time periods) or geodetic methods (short time periods of several dozen years). Geodetic methods require much greater precision and accuracy.

Using the DInSAR method, many times greater displacements were found than those determined from the GNSS data. From the stations used, the best intermediate results and the final result were obtained for the BOGI station. The displacements determined from the GNSS data developed by the two centers (UWM and NGL) are convergent. Using more pixels around the station did not improve the results much, but the pixel coherence decreased by averaging it. LOS displacements determined from the ascending and descending satellites differ, sometimes significantly, which reduces the reliability of the results. The areas around the station were characterized by a different geological structure, water relations, urbanization process, and land cover. These factors significantly influenced the coherence of the images. The presence of afforestation, high underground water levels, and anthropogenic changes (e.g., excavation of aggregate and sand, roads, and building constructions) negatively affect the coherence of the photos. The dielectric constant also changes during aggregate extraction, which affects the reflection of the signal. High coherence occurs in images with a large number of anthropogenic structures. Forest and agricultural areas show high coherence only in a short period of time; – therefore, the SAR data should be taken with the shortest possible time base (6 days and 12 days). Phase development must be preceded by separate tests for each micro-area around GNSS stations. The time interval between the SAR satellite images affects the reliability of the values obtained from the calculations. There is no unequivocal way to verify the results obtained from DInSAR. Due to the much greater displacement values from the SAR data than from the GNSS data, it can be assumed that building the time series of the photos, with a shorter time base, would have a significantly greater amplitude than the GNSS time series. This would translate into an increase in the share of noise in the data. It was found that micro-areas around permanent Polish GNSS stations show divergent vertical displacements. Taking into account the above assumptions, the obtained results, and the entire data processing process, it was found that the DInSAR method with such a large time bases is not appropriate for the determination of VCM in Poland. However, it could be a significant support in selecting areas requiring detailed monitoring of vertical movements of the Earth’s crust or support the process of selecting the location for new ASG-EUPOS stations. Therefore, to bypass the basic limitations of the DInSAR method, it is developed toward time series analysis (MTI – Multi-Temporal InSAR). This is possible thanks to the use of point methods, based on the selection of pixels that maintain coherence over time.

This article is part of a dissertation on VCM monitoring using different datasets. Further research is planned to test the PSInSAR method for determining VCM as a method based on constant scatters, which allows for better control of the obtained results. It is planned to designate VCM in micro-areas around several dozen Polish ASG-EPOS stations. The first trial carried out gave satisfactory results.