Linking persistent scatterers to airborne laser scanning points for identifying real objects reflecting SAR signal

3.1 Point scatterer position uncertainty

Persistent Scatterer Interferometry is a widely adopted advanced InSAR technique for displacement monitoring with millimeter-level accuracy [35], 36]. Numerous studies provide a comprehensive discussion of PSI processing details [1], 11], 37], 38]. In this section, we will focus on the methodology involved in estimating the position uncertainty of PSs.

In the radar coordinate system, the position of a scatterer is defined in range (r), azimuth (a), and cross-range (c) coordinates. Range refers to the direction of the radar signal, azimuth to the direction of the satellite’s flight, whereas cross-range is perpendicular to the range and azimuth direction. Range and azimuth can be estimated from the mean of the full stack of SAR images by detecting the effective phase center of the radar scatterer. The units of the estimates used are pixels. The variances (σ ²) of point scatterer (P) position in azimuth and range directions are dependent on the signal-to-clutter ratio (SCR _P ), which describes the uncertainty of the peak position caused by clutter or the presence of more than one dominant scatterer, and oversampling factor Δ _P [23]:

In this study, the SCR _P was calculated in the same manner as in ref. [29] using the normalized amplitude dispersion index. In addition, since no oversampling of the data was applied, an oversampling factor equal to 1 was used to calculate the scatterer variances for the range and azimuth coordinates.

In order to estimate the cross-range variance, InSAR observations are required. According to previous studies [23], 29] the variance of P position along the cross-range direction was calculated by transforming the height variance ( σ h , P 2 ) to the radar geometry (azimuth, range and cross-range) by using the following equation:

where θ is the incidence angle.

The uncertainty of the scatterer position in the radar geometry can then be described using the following variance-covariance matrix (Q _rac):

The Q _rac matrix is diagonal because it is assumed that the range, azimuth and cross-range coordinates are uncorrelated. The PS position in the radar coordinate system is transformed to a geodetic reference frame by a non-linear transformation called geocoding [27]. Then, coordinates can be transformed into a specific coordinate system, e.g., to the system in which LiDAR data are provided. This means that the position uncertainty in the reference coordinate system depends on the range, azimuth, and cross-range uncertainties, as well as the local incidence angle and heading angle of the radar signal that affects the transformation parameters. Considering the variance propagation law and Equation (3), using an S-transformation matrix (R) [23], 39], the variance-covariance matrix in the geodetic reference frame is given by [23]:

where the diagonal and nondiagonal represent the variances and covariances in east (e), north (n), and height (h) directions, respectively. The R matrix is defined based on the radar’s geometry, incorporating the local incidence angle (θ) and heading angle (α):

The geometric representation of the variance-covariance matrix Q _enh is an error ellipsoid with semi-axis lengths equal to the eigenvalues of the Q _enh matrix. These values reflect the random errors in the positions of PSs. Effects of systematic errors causing misalignment of PS positions with reference frame are reduced at the transformation stage using the ICP algorithm (Section 3.3).

3.2 LiDAR point cloud processing

Among the LiDAR techniques mentioned earlier (Section 1), ALS appears to be particularly suitable for the problem discussed because it can cover large areas, and its acquisition geometry is similar to SAR. Compared to PSs, ALS data has a denser point distribution and higher point position accuracy. However, a drawback of ALS is the relatively low number of points on building facades, especially in urban canyons. This limitation poses a challenge when attempting to capture intricate urban structures.

In this study, we assumed that ALS and PS can be the same points, though LiDAR and SAR use different wavelengths and different principles of signal reflection and point creation. C-band and X-band SAR signals can penetrate some objects (e.g., vegetation, snow) to some extent depending on the wavelength (longer wavelengths enable deeper penetration). However, such surfaces generally do not generate PSs because they cannot be classified as stable scatterers. The penetration for more solid surfaces is negligible, considering the accuracy of ALS points. ALS uses typically infrared wavelength that, in general, does not penetrate objects, causing points to be created on the object surfaces with the single laser pulse. Because the laser beam is divergent, parts of the pulse footprint may “fall” on different objects (e.g., different parts of the tree, ground) producing multiple echoes of the single laser pulse. This multiple return capability of the ALS causes an illusory effect of vegetation penetration, but points are always created on the top of the surfaces.

Most PSs are related to man-made structures because of their strong scatterer mechanism [40], while ALS data contains points that also belong to other objects. ALS points are usually assigned to real objects thanks to classification into typical point classes defined by ASPRS standard [41]. However, some ALS point classes may contain a mixture of objects suitable and unsuitable for PS forming. Therefore, appropriate ALS data filtering is crucial to minimize the risk of linking PSs to LiDAR points that are not PSs. The selection of appropriate linking candidates from the ALS point cloud is essential at two distinct processing stages: firstly, during the position correction utilizing ICP, and secondly, in the final linking stage (Figure 1). The proposed method involves a few simple filters (Figure 2) that significantly limit the number of wrong PSs candidates in the ALS dataset. In this research, we propose a filtering workflow that, besides the information about existing classification, considers echo number and geometrical features of points that can distinguish some of man-made objects and considers side-looking SAR geometry.

Figure 2:

Methodology LiDAR point cloud filtering for PS linking candidate selection.

The first step of the proposed filtering workflow is ALS point filtering based on the echo information, which is a feature of pulsed scanners exclusively utilized in ALS [42]. In our method, we keep only points being the first echo since PSs cannot be created on the ground surface or other objects under vegetation cover. The second stage contains normal vector and point feature calculation using Principal Component Analysis (PCA). The normal vector is used during ICP registration and for point cloud filtering in specific class similarly as linearity and planarity features. The proposed method can be applied to any laser scanning as well as photogrammetry point cloud that has at least 3 classes: Ground (natural and man-made), Buildings (containing roofs and walls), and Unclassified (all remaining points). Classes Ground and Buildings contain objects that create most of the PSs. The remaining points may be included in the Unclassified class or may belong to other classes. If any other class exists in the ALS dataset, it can be one of three types: (I) completely removed from PS candidate search (e.g., classes: Vegetation, Water, Low points), (II) completely accepted as PS candidates (e.g., Infrastructure class), (III) additionally filtered to search for possible PS candidates similarly as Unclassified class. The class Buildings follows additional filtering to remove points that stay in the shadow due to the SAR satellite line of sight (LOS). Considering the necessity to filter only two classes: Unclassified and Buildings, among the possible PCA-based features, three parameters were used: planarity, linearity, and information about surface orientation (normal vector). For the reliable determination of these features, choosing an appropriate point neighborhood for applying PCA and calculating these features is crucial. As an optimal neighborhood, the authors in ref. [43] defined the largest set of spatially close points that belong to the same object as the point of interest. According to ref. [44], the optimal number of points in the neighborhood that provides the best local description of the point cloud for different classes of objects is 10 points. The neighborhood size should be large enough to enable PCA feature determination for all points of the tested objects while being small enough to prevent the inclusion of points from neighboring objects in the defined local neighborhood. Among the methods for finding neighboring points, a spherical neighborhood was used, with the radius determined experimentally. The tested radius range was 0.5–3.0 m and resulted from the density of the ALS point cloud. The value of the final radius was determined based on tests conducted on ALS points representing several man-made linear objects (such as poles and fences) and planar objects (such as sound barriers and information boards) that can form PS.

The selection of potential PS candidates, being infrastructure objects but located in the Unclassified class, is dependent on planarity or linearity values. Points exceeding predefined S or L thresholds (Figure 2) are accepted as potential PS candidates. For instance, utility poles and cables would exhibit high linearity values, while sound barriers or fences would display elevated planarity values. Natural objects are less likely to meet the specified criteria. Thresholds for S and L were determined empirically by analyzing the abovementioned test objects. The maximum values that ensured all test objects were classified as potential PS candidates were selected as the appropriate thresholds. A set of S and L values from 0.5 to 0.9 with a step of 0.05 were tested. The maximal value of 0.9 was taken based on the literature [45]. The final stage of the proposed PS candidate filtering methodology is intended to remove points from the Buildings class that are in the shadow with respect to the SAR imaging geometry. Due to the significant differences in geometry between SAR and ALS, this step addresses the fact that SAR is a side-looking system, whereas LiDAR captures data from a top-down perspective. This difference can addressed by the SAR signal incidence angle to the object surface (e.g., building façade, roof, etc.). For this purpose, the angle between each point’s normal vector and the vector opposite the LOS of the radar signal is computed. Points with an angle larger than 90° are identified as in shadow and then removed. The LOS-opposite vector is calculated as a unit vector based on angles describing the LOS geometry – the azimuth and incidence angle. The final result of the filtering method is illustrated in Figure 3.

Figure 3:

PS candidate selection methodology. (a) Shows the cross-section of original ALS data and the main classification of the point cloud. (b) Presents the results of the proposed methodology.

3.3 Global transformation parameters estimation using Iterative Closest Point algorithm (ICP)

In radar geometry, the position of PS can be influenced by various factors, which can be categorized into two main groups: systematic factors applicable to the entire dataset and factors specific to a particular PS point. Systematic offsets are caused by, e.g., SAR image timing errors, orbit errors, uncompensated atmospheric signal delays, and a PS reference point height offset. The additional individual PS point position errors are due to uncertainty in the scatterer peak detection (SCR dependent) and the PS height estimate. Using a LiDAR point cloud as well as other datasets such as surface models, it is possible to reduce the effect of errors. In ref. [29] authors proposed an iterative procedure to eliminate systematic height shifts using a gridded surface model obtained by LiDAR data. This method is more precise than single-point correction but does not estimate all shifts independently. In this work, we suggest estimating not only shifts between PSs and LiDAR point clouds but also rotations. We propose to treat the PSs as a sparse point cloud and execute point cloud-to-point cloud registration using the well-known ICP algorithm, which allows the transformation that optimally aligns two or more sets of points. The ICP algorithm has many implementations and modifications aimed at improving its accuracy. During the optimization process in the ICP algorithm, the goal is to minimize the distance between a point in the transformed cloud and a point in the target cloud. This distance minimization can involve point-to-point, point-to-line, or point-to-plane relationships, depending on the specific problem configuration. In this research, we used the point-to-plane distance minimization criteria in the ICP algorithm implementation included in the Open3D Library for Python [46]. Our study focuses on ALS point clouds, where, due to the scanning geometry, low point density on a vertical wall is often a challenge. Therefore, using a point-to-plane approach seems to be an appropriate choice. The point-to-plane distance between two corresponding points is defined as the orthogonal distance of one point to the fitted plane of the other point [47]. The plane for a point from the reference cloud is usually determined based on its neighboring points. In this work, we determined these planes in the previous step (Section 3.2) by calculating the normal vector for each point. Minimizing the distance determined in this strategy proves to be a helpful approach, particularly in situations characterized by areas of low point density. In such cases, the distance to the plane formed by the selected points is expected to be less than the slant distance to the actual point. In this research, as a target point cloud the ALS dataset after the PS candidate selection step was chosen. This is a crucial assumption because, in the ICP algorithm, the closest points are treated as corresponding points. Thus, ALS points that cannot correspond to PSs should be removed.

The ICP algorithm implementation used in this work has a few parameters that need to be provided: initial coarse transformation parameters, distance threshold, and convergence criteria. Initial coarse transformation parameters include rotation and translation parameters that describe the initial alignment of point clouds. Because PSs are transformed earlier from radar geometry to the geodetic frame of the LiDAR point cloud, both point clouds are in the same coordinate system. Thus, coarse transformation parameters are equal to 0. The distance threshold is the maximum point-to-plane distance accepted. The computations do not include points for which the distance threshold is exceeded. In this work, the distance threshold is set to equal the maximum pixel size in the radar geometry of the SAR image, i.e., 14 m and 2 m for Sentinel-1 and TerraSAR-X data, respectively. The convergence criteria decide when the iteration process ends. In this work, two parameters of convergence criteria were used: a maximum number of iterations of 100 iterations and a change in RMSE values between the two last iterations of smaller than 0.001 m. The RMSE is calculated from point-to-plane distances for inlier points that met the ICP threshold criteria. The iteration process ends when any of these criteria are met. The results of the executed ICP algorithm are the transformation parameters from PS to LiDAR point cloud and ICP quality assessment values. Estimated transformation parameters were applied to PS the point cloud, and thus, the corrected PS point cloud was prepared for final linking.

3.4 Final linking

Using the ICP method, it is possible to improve the accuracy of the position of PSs, but they are still not linked to real objects because the same translation and rotation value is applied to all points. Therefore, an individual approach to each point is required. For the final PS linking, we adopted the methodology proposed in ref. [22] involving a nearest neighbor search process concerning the radar geometry that links the scatterers to their most likely point in the LiDAR point cloud. Searching is performed inside the error ellipsoid defined by the variance-covariance matrix which was described in detail in Section 3.1. An error ellipsoid at the 2-sigma confidence level was used, meaning that the true reflection point is expected to be within the ellipsoid in approximately 95 % of the cases.

Since the original error ellipsoids based on Q _rac are typically non-spherical due to differences in the variances in the azimuth, range and cross-range directions, this also applies to the error ellipsoids after transformation. Consequently, finding the closest LiDAR point and, thereby, the most likely PS scatterer location is not straightforward. Therefore, a whitening transformation [48] is adapted to decorrelate the original data, providing a set of uncorrelated data with unit variances. The transformation matrix (W) is formed using the eigenvalues (E) and eigenvectors (D) of the covariance matrix of PS in the geodetic reference frame (Q _enh ):

A transformation was performed on the ALS and PS point clouds. After this transformation, searching the nearest neighbor in the LiDAR point cloud to each PS can be performed in the Euclidian space. The PS point was assigned with the selected LiDAR point coordinates (corrected position) and within the confidence ellipsoid. Otherwise, its PS position was not corrected. If the confidence ellipsoid contained more than one point, the selection of the single linking point was based on the point class and the distance to the ellipsoid center. The priority of links was put on the class Buildings, then Ground and Infrastructure (if it exists), and in the last order, the newly created Other class. The prioritization of classes was applied according to object type suitable to form PS and to prevent links to noise points that were not completely removed from the Unclassified class during filtering but are closer to the ellipsoid center than points from other classes. If more than one point of the same class was included in the confidence ellipsoid, the nearest neighbor to the ellipsoid center was selected as the PS link.

The final output of the entire linking procedure is the original PS dataset enriched with corrected PSs coordinates obtained during final linking, the LiDAR-based class assigned to the PSs, and information regarding the determined variances in the radar geometry. PS points not linked to LiDAR points did not get corrected coordinates and class.

3.5 Validation

3.5.2 Qualitative evaluation – 3D visualization

Typically, PSI processing results are visualized in 2D maps as a point with information about the displacement. This type of visualization does not give information about the PS location in 3D space. For instance, PS visualized on a bridge may belong to the bridge deck, bridge support, or another element. By performing 3D visualization, we can more precisely determine from which object or its specific part the PS originates. Regardless of the dimensionality of the visualization, linking PSs to the LiDAR point cloud, as described in the workflow above, leads to new coordinates for numerous PSs. This can pose a challenge in visually verifying the accuracy of the proposed method. Visualizing PS’s location before and after linking data, even in 3D, may not give adequate insights if the number of points is very large. In our contribution, we propose another way to visualize the linking results that enable showing PS’s original location (P) as well as its corrected location after ICP registration (P′) and the location of the linked LiDAR point (P″) together with the linking vector (Figure 4). In addition, error ellipsoids are visualized to show the uncertainty of PS locations (Figure 4). This approach enables visual verification of the linking process and facilitates the interpretation of the PSI results.

Figure 4:

Scheme of 3D visualization of PS-LiDAR linking results for one point. PS point location: original – P (orange), after ICP registration P′ (red), linked LiDAR point P″ (violet) and the error ellipsoid (black). The shape of the ellipsoid is described by the semi-axes (green, blue, orange) and oriented according to the radar geometry. P′ is the center of the ellipsoid. The linking vector (yellow) to the LiDAR point is within the ellipsoid.

The visualization module was also implemented in Python, using packages such as PyVista (https://docs.pyvista.org/) and LasPy (https://laspy.readthedocs.io/). Depending on the analysis, certain visualization elements, such as ellipsoid shapes, vectors, or PSs, may be excluded from the selected stage to improve the interpretation of the results.

5.1 PS position uncertainties

The standard PSI approach was used with two software programs to obtain PS displacements in the presented areas of interest: DePSI and SarScape were employed for the Amsterdam and Ruda Śląska test sites, respectively. Based on accuracy metrics provided during PSI processing, we estimated each PS’s position uncertainty in range, azimuth, and cross-range directions. Figure 6 presents boxplots of the obtained uncertainty in each dimension for both case studies and the different SAR datasets.

Figure 6:

Uncertainties of the PS point positions in r – range, a – azimuth and c – cross–range for the Amsterdam (a–c) and Ruda Śląska (d–f) test sites.

The major uncertainty was encountered in both sites for S1 PSI products. The larger is directly related to the pixel spacing of the SAR data, which is 2.3 m × 14.1 m and 1.5 m × 1.8 m for the S1 and TSX missions, respectively, in the range and azimuth directions. The mean uncertainties for S1 in Amsterdam, where DePSI software for PSI calculations was used are 0.9 m, 5.4 m, and 4.4 m for range, azimuth, and cross-range, respectively, and correspond to the results presented in ref. [50]. If we compare solutions from two different software, it is clearly visible that the uncertainty values, as well as their variations, are lower for the DePSI solution (Figure 6a–c) than for SarScape (Figure 6d–f). It should be highlighted that DePSI calculates the coordinates of points in the radar coordinate system as sub-pixel positions, while SarScape lacks this capability, which will significantly impact the obtained results. The results confirmed the need to improve the accuracy of PS localization with other methods apart from the standard PSI technique processing, as described in the following sections of this paper.

5.2 LiDAR point cloud processing results

During the initial PS candidate selection stage, any LiDAR points with return numbers greater than 1 were removed. Then geometric parameters such as planarity, linearity, and normal vectors were estimated based on PCA (see Figure 2). Classes with suitable candidates for linking were found for the Dutch ALS dataset from the Ground (class 2), Building (class 6), Infrastructure (class 26), and Unclassified (class 1) classes. For the Polish ALS dataset, the infrastructure class was not specifically identified, resulting in the selection of solely the Ground (class 2), Building (class 6), and Unclassified (class 0) classes. It is important to note that the Dutch point cloud lacks the vegetation class, unlike the Polish ALS point cloud, causing more challenges in searching for man-made structures in the Unclassified class. Geometrical filtering in class 1 (Dutch dataset) or 0 (Polish dataset) is based on the planarity or linearity values. The test that was executed on optimal hyperparameter selection resulted in the radius of the local neighborhood being equal to 2 m, planarity threshold S = 0.7, and linearity threshold L = 0.6 for both ALS datasets. Datasets with similar characteristics, especially point cloud densities, are expected to result in similar values of the above hyperparameters. However, for point clouds with notably higher or lower densities and accuracies (e.g., those obtained through different techniques), experimental verification of radius, planarity, and linearity values should be executed. Figure 7 presents the quantitative results of ALS point cloud filtering for ICP registration and further linking with the Sentinel-1 ascending dataset (removed points are marked with hatches). According to the proposed workflow, the filtering applied for selecting potential PS candidates reduced the ALS point clouds to 51 % and 46 % (Figure 7, sum of green values) of their original sizes for the Amsterdam and Ruda Śląska datasets, respectively. Depending on the geometry of the SAR data with which the ALS point cloud is linked, this value changes because the number of points rejected from the Buildings class due to SAR geometry varies. These differences range between 1 and 2 %.

Figure 7:

Filtering results of ALS point clouds for linking with ascending Sentinel-1 mission data: (a) Amsterdam, (b) Ruda Śląska. The pie charts illustrate point classes and the filtering process applied to ALS point clouds: (1) colors represent individual classes existing in ALS data. (2) Hatches indicate points removed at each filtering stage based on specific criteria, including class-based filtering, return number filtering, LOS geometry filtering, and geometric feature filtering. (3) Labels show the point percentages of each part of pie chart relative to the total number of points in the dataset; green labels indicate percentages of points used for the final PS-ALS linking, while black labels indicate percentage of points removed during filtering.

The results of filtering the Unclassified class were added to a new class – Other (marked in solid gray in Figure 7). Depending on the area of interest, it contained 0.03 % and 1 % of all ALS points for Amsterdam and Ruda Śląska datasets, respectively. Thanks to filtering based on geometrical features, it was possible to find suitable PS linking candidates in the Unclassified class and belonging to linear objects (Figure 8a), road barriers, or structural elements of bridges (Figure 8b). Additionally, it incorporated points from flat surfaces, specifically some misclassified building points, enabling them for further analysis and linking.

Figure 8:

Examples of the Dutch LiDAR point cloud before (top) and after (bottom) PS candidate filtering. Colors mean point classes: grey – class 1 (unclassified), brown – class 2 (ground), orange – class 6 (building), blue – class 8 (water), yellow – class 26 infrastructure), pink – class 27 (other). (a) Power line mast and surroundings. (b) Bridge and surroundings.

Considering the selection of possible candidates based on their position in relation to the acquisition geometry, we determined the angles between the LiDAR normal vectors and the PS LOS vector Figure 9 presents the results for selection based on the TSX ascending (Figure 9b) and descending (Figure 9c) geometry for the same building. This approach makes it possible to reject LiDAR points obscured by other points (Figure 9b) from further analysis.

Figure 9:

Example of building in Dutch LiDAR point cloud before (a) and after PS candidate filtering depending on the LOS direction for ascending (b) and descending (c) orbit. Colors mean point classes: grey – class 1 (unclassified), brown – class 2 (ground), orange – class 6 (building), blue – class 8 (water), yellow – class 26 (infrastructure), pink – class 27 (other).

5.3 PS global transformation results

The next step is the estimation of global transformation parameters for the PSs. The ICP algorithm was adopted to calculate global transformation parameters between PS and LiDAR point clouds. The LiDAR point cloud, which contains only points after PS linking candidate selection (Section 3.2), was chosen as the reference cloud. Utilizing the transformation obtained from the ICP method with the Open3D library, the effectiveness of this process is assessed using two metrics, including point fitness, describing the percentage of PSs that are inliers, and RMSE, determining the accuracy of alignment through the root mean square error for correctly matched points. RMSE is calculated based on point-to-plane distances for inlier points. Global transformation parameters and ICP registration performance measures are listed in Table 3.

Despite the significant difference in density between the ALS and PS point clouds, the ICP algorithm found correspondences to 78.1 %–98.3 % PSs in the LiDAR point cloud that met the point-to-plane maximal distance criteria. Obviously, these correspondences are not necessarily linking points (P″) (Figure 4) since they may be outside the ellipsoid error. The S1 data showed over 90 % matching for both case studies, potentially due to the higher acceptable search radius. Considering the obtained RMSE values presented in Table 3, it can be observed that for the TSX data sets the value is much lower and is about 0.9 m, while for S1 it is between 3.4 and 4.1 m, which is also coincident with the spatial resolution of S1 and TSX data. The process of determining global parameters using ICP was successful. The transformation parameters were determined based on over 75 % of all PSs, and the average distances between them (based on RMSE) are smaller than the resolution of SAR images.

For all analyzed datasets, the transformation primarily involves translation, with the rotation angles being close to 0. The determined translation parameter values are higher for Sentinel-1 data, that can be attributed to the lower PS position accuracy estimated during processing (Figure 6). Additionally, the highest values are observed in the vertical component (Tz). For the TSX datasets, the difference between the translations for the horizontal (Tx, Ty) and vertical components is much lower. When the transformation parameters were calculated, the point cloud of the PSs was transformed, resulting in their corrected localization (P′) (Figure 4).

5.4 Final linking – results validation

The final linking procedure for PS within the 2-sigma confidence level ellipsoid was conducted in two study areas for points with corrected localization (P′, Figure 4): Amsterdam and Ruda Śląska, using three datasets for each location. The obtained linking results, categorized into individual point classes according to the ALS data classification, as a percentage distribution is depicted in Figure 10a. In addition, the results obtained for processing, excluding one of two key steps of the proposed method, are also included. Figure 10b shows linking percentages for processing without systematic error removal by applying global transformation using parameters determined with the ICP algorithm (Section 3.3), and Figure 10c shows linking percentages for processing without ALS point cloud filtering (Section 3.2).

Figure 10:

Percentage of linked PSs categorized into point classes according to the ALS data classification and PSs that could not be linked: (a) according to the proposed method, (b) excluding global transformation utilizing ICP, (c) excluding ALS data filtering. The infrastructure class is limited to the ALS dataset from Amsterdam, where it is labeled.

Thanks to the proposed approach, 80 % of the PS points from the Sentinel-1 ascending orbit and 72 % and 71 % of the PS from TerraSAR-X ascending and descending orbits, respectively, were linked to real objects in the Amsterdam case study (Figure 10a). The results were comparable for the Ruda Śląska case study, with 88 %, 90 %, and 65 % linked PSs determined from the Sentinel-1 descending, Sentinel-1 ascending, and TerraSAR-X ascending orbits, respectively (Figure 10a).

The percentage distribution between different classes varied depending on the analyzed study area. For Amsterdam, the area analyzed is located in the city’s center, where high-density housing occurs. For the Ruda Śląska area, single-family housing dominates. In addition, there is a highway crossing the Ruda Śląska test site, and points are classified as Ground, making this class the majority of linked PS (50 %). Additionally, the Other class (Figure 10a) was extracted from the Unclassified class utilizing PCA planarity and linearity features. It contains linking points for infrastructure objects, such as high-voltage poles (Figure 11a), thus increasing the percentage of linked PSs. This is particularly true for Sentinel-1 data, where the percentages of the Other class for Amsterdam are 19 % (S1 ASC), 8 % (S1 DSC), and 9 % (S1 ASC) for Ruda Ślaska.

Figure 11:

Linking results for selected objects: (a) high-voltage pole and S1 ASC PSs with error ellipsoids (green) for each point – Amsterdam test site, (b) TSX PS linking results with error ellipsoids for each point from both ASC (green ellipsoids) and DSC (orange ellipsoids) orbits – Amsterdam test site. PSs colors: original – red, after ICP registration – green, linked with LiDAR – yellow. The arrows connect the points after the ICP with the linked position, and their color describes the distance between the points.

The results of processing that excludes one of the processing steps (Figure 10b and c) demonstrate that systematic bias removal by utilizing ICP for the estimation of global transformation parameters is a crucial step. Omitting the global transformation executed with ICP leads to a decrease in the percentage of linked points in all case studies, except for the Sentinel-1 dataset in Amsterdam, where a slight increase of 1 % is observed. However, in this case, the proportion of points linked to the building and ground classes changes (Figure 10b). Moreover, the execution of global transformation also ensures that the percentages of PSs assigned to each class are similar for both DSC and ASCdatasets of the same SAR system. For example, without the global transformation step, in the S1 ASC and S1 DSC datasets for Ruda Śląska, a notable discrepancy existed in the number of PSs classified for each class, e.g., 34 % and 61 % for Ground class for ASC and DSC datasets, respectively (Figure 10b). Applying global transformation using parameters estimated with ICP allowed for balancing these proportions, e.g., 51 % and 56 % for the Ground class for ASC and DSC datasets, respectively (Figure 10a). This reflects better the fact that both datasets represent the same area and product, while small differences occur due to differences in LOS geometry. Moreover, for the S1 DSC dataset for Ruda Śląska, the positive value of vertical translation (Table 3) corrected misclassifications of PSs (initially below the terrain), linking them correctly to buildings instead of incorrectly to the ground (compare percentages for Buildings and Ground classes in Figure 10a and b, Ruda Śląska S1 DSC). In contrast, for the ASC dataset, where PSs were initially above the terrain, causing that most of them to be linked to buildings (Figure 10b, Ruda Śląska S1 ASC), the global transformation allowed to improve their linkage to ground points (Figure 10a, Ruda Śląska S1 ASC). The global transformation utilizing ICP effectively mitigated systematic errors, enhancing the reliability and accuracy of the linking process.

ALS data filtering step was similarly vital for improving linking accuracy. Filtering reduces the likelihood of linking PSs to classes unlikely to generate persistent scatterers, such as vegetation, water, or unclassified objects, e.g., vehicles or temporary structures (Figure 10c, green and grey bars). While skipping this step increases the overall percentage of linked points, it results in decreased linking quality, with PSs more frequently misclassified.

The quality of the lining process may be asses also by the distribution of the linking vector lengths between the corrected PS (P′, Figure 4) and the linked PS (P″, Figure 4). Histograms of linking vector lengths (Figure 12) demonstrate that the magnitude of the vectors corresponds well to the data resolution. The 3D vector lengths mostly do not exceed 2 m for TSX and 10 m for S1 in both case studies.

Figure 12:

Histograms of the linking vectors (P′P″) for both TSX and S1 datasets in Amsterdam (top row) and Ruda Ślaska (bottom row), respectively.

Despite the higher initial accuracy in positioning PS for TerraSAR-X data (Figure 6), finally, fewer percentage of points were linked to real objects compared to Sentinel-1 data (Figure 10a). This is particularly evident in the case of the Upper Silesia region. Visual analysis of the obtained results reveals that the positioning of PS from TerraSAR-X often occurs near the vertical walls of buildings (Figure 13, green dots). One of the key factors contributing to the observed difference in linking percentages between S1 and TSX is the limitation of ALS data, particularly its scan geometry. Gaps in the ALS data frequently occur on building walls or in shadows, especially in densely built-up areas, as clearly shown in Figure 13. When PSs form near buildings, a double-bounce reflection often occurs, i.e., wall-ground reflections, where the actual PS position lies at the intersection of these two surfaces. Since these vertical surfaces are typically absent from ALS data, the correct links cannot be created. For individual buildings, TSX, due to its higher data resolution compared to Sentinel-1, often generates more points corresponding to building facades. As a result, fewer percentage of points are linked in these areas, making it more challenging to establish reliable correspondences. The smaller positional uncertainty of TSX PSs, represented by their error ellipsoids, restricts the spatial area within which a match can be identified. Conversely, the larger positional uncertainty of Sentinel-1 PSs facilitates linkage, as the broader error ellipsoid increases the likelihood of finding a match within the available ALS data. The absence of links to PSs located on walls supports the validity of the proposed method. It was designed to identify only reliable links rather than to maximize their number, which would lead to numerous erroneous links.

Figure 13:

Example of buildings with missing LiDAR points on the walls in Ruda Śląska test site. The green color indicates PSs obtained from TSX ASC orbit before linking (P′), and the yellow color indicates linked LiDAR point (P″). Arrows connect the new positions of the points, and their color describes the distance between them.