Fully Automated Profile-based Calibration Strategy for Airborne and Terrestrial Mobile LiDAR Systems with Spinning Multi-beam Laser Units

LiDAR-based mobile mapping systems (MMS) are rapidly gaining popularity for a multitude of applications due to their ability to provide complete and accurate 3D point clouds for any and every scene of interest. However, an accurate calibration technique for such systems is needed in order to unleash their full potential. In this paper, we propose a fully automated profile-based strategy for the calibration of LiDAR-based MMS. The proposed technique is validated by comparing its accuracy against the expected point positioning accuracy for the point cloud based on the used sensors’ specifications. The proposed strategy was seen to reduce the misalignment between different tracks from approximately 2 to 3 m before calibration down to less than 2 cm after calibration for airborne as well as terrestrial mobile LiDAR mapping systems. In other words, the proposed calibration strategy can converge to correct estimates of mounting parameters, even in cases where the initial estimates are significantly different from the true values. Furthermore, the results from the proposed strategy are also verified by comparing them to those from an existing manually-assisted feature-based calibration strategy. The major contribution of the proposed strategy is its ability to conduct the calibration of airborne and wheel-based mobile systems without any requirement for specially designed targets or features in the surrounding environment. The above claims are validated using experimental results conducted for three different MMS – two airborne and one terrestrial – with one or more LiDAR unit.


Introduction
Laser scanning has gained widespread popularity over recent years owing to the continuous improvement in the performance, size, and cost of available sensors. Recent advancements in direct georeferencing technology using the Global Navigation Satellite System/Inertial Navigation System (GNSS/INS) have further boosted the ability to collect accurately georeferenced data when combined with laser units onboard mobile systems. Mobile mapping systems (MMS) equipped with one or more laser scanners are heavily relied on for a multitude of applications, including urban modeling [1-3], transportation corridor mapping [4][5][6], precision agriculture [7][8][9][10], infrastructure monitoring [11,12], shoreline monitoring [13][14][15], archaeological mapping [15][16][17], and digital documentation of cultural heritage sites [18][19][20]. Each application determines the necessary architectural requirements for the MMS depending on the accessibility of the survey site, desired level of accuracy and details, and so on. The design of an MMS entails decisions in terms of the platform (airborne or terrestrial), onboard GNSS/INS unit, number and type of mounted LiDAR units, and any other types of sensor (RGB,

System Description and Calibration Test Field
As mentioned earlier, this research proposes a fully automated approach for calibrating airborne and terrestrial MLMS consisting of spinning multi-beam laser scanners. In this study, we worked with three such mobile mapping systems-two airborne and one terrestrial-equipped with one or more spinning multi-beam LiDAR units along with a GNSS/INS unit for direct georeferencing. The following subsections describe the system architecture and calibration test fields for each of these systems in detail.

Airborne MLMS
The two UAV-based LiDAR systems used in this study are shown in Figure 1 and are denoted henceforth as Airborne MLMS 1 and 2, respectively. Each system has one LiDAR unit integrated with a GNSS/INS unit for direct georeferencing. The airborne MLMS 1 is equipped with a Velodyne VLP32C LiDAR unit, which consists of 32 radially oriented laser rangefinders, whereas the airborne MLMS 2 has a Velodyne VLP16 Puck Lite LiDAR unit, which has 16 radially oriented laser beams. The manufacturer specifications for each of these LiDAR units are listed in Table 1 [34,35]. Each of the systems consists of an APX-15 UAV V3 GNSS/INS unit, which has a post-processing positional accuracy of 2 to 5 cm, 0.025 • accuracy for the roll/pitch angles, and 0.08 • accuracy for the heading [36]. In order to generate a point cloud from these airborne systems, the maximum reconstruction range was set to 70 m and the FOV for reconstruction across the flying direction was ±70 • . Based on the aforementioned individual sensor accuracies and reconstruction parameters, we could estimate the resultant accuracy of point positioning for such a system using the LiDAR Error Propagation Calculator Remote Sens. 2020, 12, 401 4 of 39 developed by Habib et al. [37]. For a flying height of 50 m, the point positioning accuracy was found to be around 5 to 6 cm in the horizontal and vertical directions at nadir. At the edge of the swath, the horizontal accuracy would be about 8 to 9 cm and the vertical accuracy would be 5 to 6 cm.  600,000 points per second 300,000 points per second Apart from the system components for the two airborne MLMS, another difference lies in the relative orientation in which the LiDAR unit is mounted with respect to the GNSS/IMU body frame for the two systems. The nominal boresight angles relating the LiDAR unit and GNSS/IMU body frame in airborne MLMS 1 are 90°, 90°, and 0° for the roll, pitch, and heading angles, respectively, and those for airborne MLMS 2 are 90°, −90°, and 0°, respectively. One should note that these nominal values indicate that both systems will have a gimbal lock during boresight calibration due to a secondary rotation of ±90°. As suggested by Ravi et al. [33], such a gimbal lock problem can be mitigated by introducing a virtual laser unit coordinate frame ( ) approximately aligned with the IMU body frame and thereafter estimating the boresight angles relating the virtual frame and IMU body frame during calibration. The discussed coordinate systems (original and virtual laser unit coordinate frames and IMU body frame) are shown in Figure 1.
A 3D point cloud of the calibration test field used for airborne MLMS in this study is shown in Figure 2a (colored by height), along with the flight lines configuration (also colored by height with increasing height from black to white). The test field consisted of 16 highly reflective sign board targets laid on the ground along with 5 hut-shaped targets. The east side of the test field has a building, which is about 12 m high with a gable roof, while the west side is covered by vegetation. Figure 2b also shows zoomed-in images of a few of the reflective and hut-shaped targets deployed in the field. These targets are not required for the proposed profile-based calibration study. However, they are deployed in the field in order to be used to conduct the manually-assisted feature-based calibration strategy proposed by Ravi et al. [33], the results from which were to validate the accuracy  600,000 points per second 300,000 points per second Apart from the system components for the two airborne MLMS, another difference lies in the relative orientation in which the LiDAR unit is mounted with respect to the GNSS/IMU body frame for the two systems. The nominal boresight angles relating the LiDAR unit and GNSS/IMU body frame in airborne MLMS 1 are 90 • , 90 • , and 0 • for the roll, pitch, and heading angles, respectively, and those for airborne MLMS 2 are 90 • , −90 • , and 0 • , respectively. One should note that these nominal values indicate that both systems will have a gimbal lock during boresight calibration due to a secondary rotation of ±90 • . As suggested by Ravi et al. [33], such a gimbal lock problem can be mitigated by introducing a virtual laser unit coordinate frame (Lu ) approximately aligned with the IMU body frame and thereafter estimating the boresight angles relating the virtual frame and IMU body frame during calibration. The discussed coordinate systems (original and virtual laser unit coordinate frames and IMU body frame) are shown in Figure 1.
A 3D point cloud of the calibration test field used for airborne MLMS in this study is shown in Figure 2a (colored by height), along with the flight lines configuration (also colored by height with increasing height from black to white). The test field consisted of 16 highly reflective sign board targets laid on the ground along with 5 hut-shaped targets. The east side of the test field has a building, which is about 12 m high with a gable roof, while the west side is covered by vegetation. Figure 2b also shows zoomed-in images of a few of the reflective and hut-shaped targets deployed in the field. These targets are not required for the proposed profile-based calibration study. However, they are deployed in the field in order to be used to conduct the manually-assisted feature-based calibration strategy proposed by Ravi et al. [33], the results from which were to validate the accuracy of fully automated profile-based calibration strategy. One should note that the test field was similar for both UAV systems, with minor variations in the distance between the deployed targets. As shown in Figure 2c, the flight configuration for both systems consisted of 18 tracks at three different flying heights and different lateral separation. There were six tracks flown at 15 m height, another six tracks at 25 m, and six more at 45 m height. Each set of six tracks at a given flying height was comprised of three laterally separated pairs of flight lines in opposite directions.
Remote Sens. 2020, 12, 401 5 of 39 of fully automated profile-based calibration strategy. One should note that the test field was similar for both UAV systems, with minor variations in the distance between the deployed targets. As shown in Figure

Terrestrial MLMS
The terrestrial mobile mapping system used in this research (shown in Figure 3) consisted of four LiDAR units (three Velodyne HDL32E and one Velodyne VLP-16 Puck Hi-Res) along with an

Terrestrial MLMS
The terrestrial mobile mapping system used in this research (shown in Figure 3) consisted of four LiDAR units (three Velodyne HDL32E and one Velodyne VLP-16 Puck Hi-Res) along with an Applanix POSLV 220 GNSS/INS unit for direct georeferencing. The individual sensor specifications for the onboard LiDAR units are listed in Table 2 [38,39]. Following GNSS/INS post-processing, the Remote Sens. 2020, 12, 401 6 of 39 POSLV 220 can attain an accuracy of less than 2 cm in position, and an accuracy of 0.02 • and 0.025 • in the roll/pitch and heading, respectively [40]. Based on the manufacturer's specifications for sensor accuracies, we were able to derive the expected accuracy for the computed mapping frame coordinates using the LiDAR Error Propagation calculator developed by Habib et al. [37]. The calculator suggested an expected accuracy of about 2-4 cm at a range of 30 m.  Table 2 [38,39]. Following GNSS/INS post-processing, the POSLV 220 can attain an accuracy of less than 2 cm in position, and an accuracy of 0.02° and 0.025° in the roll/pitch and heading, respectively [40]. Based on the manufacturer's specifications for sensor accuracies, we were able to derive the expected accuracy for the computed mapping frame coordinates using the LiDAR Error Propagation calculator developed by Habib et al. [37]. The   Figure 4 (a) and (b), respectively. The test field was a residential area with multiple houses lined up on both sides of the street. The drive-run configuration is overlaid on Figure 4 (a) in white showing the four tracks-two opposite drive-runs oriented in the N-S direction and two more with the E-W orientation.   Figure 4a,b, respectively. The test field was a residential area with multiple houses lined up on both sides of the street. The drive-run configuration is overlaid on Figure 4a in white showing the four tracks-two opposite drive-runs oriented in the N-S direction and two more with the E-W orientation.

Methodology
As mentioned earlier, this study proposes a calibration technique that can accurately estimate the mounting parameters, i.e., the lever-arm and boresight angles, relating the onboard LiDAR units to the GNSS/IMU position and orientation system. Inaccurate estimates of mounting parameters will result in a discrepancy between point clouds for the same area captured from different sensors along different drive-runs/flight lines. The conceptual basis of calibrating an MLMS is to minimize the above-mentioned discrepancy between point clouds. One should note that, hereafter in this paper, the point cloud captured by a single LiDAR unit in a given drive-run/flight line is regarded as an individual track. So, given the number of onboard sensors ( ) and number of drive-runs/flight lines , there will be a total of * tracks, which will be used in the system calibration strategy. There are two major components associated with the development of any calibration strategy: (1) definition of an optimal set of calibration primitives and drive-run/flight line configuration, and (2) development of a point-pairing strategy along with an optimization function for calibration. Ravi et al. [32,33] conducted a theoretical bias impact analysis for terrestrial and airborne MLMS to suggest an optimal track and calibration primitive configuration. Their study recommended to include calibration primitives that provide variability in three-dimensional coordinates of constituent points

Methodology
As mentioned earlier, this study proposes a calibration technique that can accurately estimate the mounting parameters, i.e., the lever-arm and boresight angles, relating the onboard LiDAR units to the GNSS/IMU position and orientation system. Inaccurate estimates of mounting parameters will result in a discrepancy between point clouds for the same area captured from different sensors along different drive-runs/flight lines. The conceptual basis of calibrating an MLMS is to minimize the above-mentioned discrepancy between point clouds. One should note that, hereafter in this paper, the point cloud captured by a single LiDAR unit in a given drive-run/flight line is regarded as an individual track. So, given the number of onboard sensors (n sensors ) and number of drive-runs/flight , there will be a total of tracks, which will be used in the system calibration strategy.
There are two major components associated with the development of any calibration strategy: (1) definition of an optimal set of calibration primitives and drive-run/flight line configuration, and (2) development of a point-pairing strategy along with an optimization function for calibration. Ravi et al. [32,33] conducted a theoretical bias impact analysis for terrestrial and airborne MLMS to suggest an optimal track and calibration primitive configuration. Their study recommended to include calibration primitives that provide variability in three-dimensional coordinates of constituent points with respect Remote Sens. 2020, 12, 401 8 of 39 to the tracks capturing these primitives to ensure sufficient control for an accurate calibration. In this study, we propose an automated targetless calibration procedure using thinly sliced profiles oriented along and across the flight/driving direction as the calibration primitives. A flowchart of the proposed methodology is shown in Figure 5.
Remote Sens. 2020, 12, 401 8 of 39 with respect to the tracks capturing these primitives to ensure sufficient control for an accurate calibration. In this study, we propose an automated targetless calibration procedure using thinly sliced profiles oriented along and across the flight/driving direction as the calibration primitives. A flowchart of the proposed methodology is shown in Figure 5. Here, the first challenge waas to develop a strategy to extract an optimal set of profiles that would ensure sufficient control along the X, Y, and Z directions for an accurate estimation of the mounting parameters. In keeping with the findings of the aforementioned work, we developed an automated two-step profile extraction strategy, shown in Block 1 of the flowchart in Figure 5. 1. Template Profile Selection: We inspected all the tracks to identify and extract template profiles, i.e., the ones that were comprised of multiple linear segments with sufficient angular variability of the fitted lines. 2. Matching Profile Identification: Next, the selected template profiles were used to automatically identify their corresponding matches in the remaining tracks (drive-runs/flight lines and sensors) using a height map correlation-based strategy. Finally, a point-pairing scheme was proposed for the extracted template and matched profiles from different tracks and an optimization function was developed to estimate the system calibration parameters-i.e., lever-arm and boresight angles for all the LiDAR units onboard the mobile system-that minimized the discrepancy between the generated point pairs, as shown in Block 2 of Here, the first challenge waas to develop a strategy to extract an optimal set of profiles that would ensure sufficient control along the X, Y, and Z directions for an accurate estimation of the mounting parameters. In keeping with the findings of the aforementioned work, we developed an automated two-step profile extraction strategy, shown in Block 1 of the flowchart in Figure 5.

1.
Template Profile Selection: We inspected all the tracks to identify and extract template profiles, i.e., the ones that were comprised of multiple linear segments with sufficient angular variability of the fitted lines.

2.
Matching Profile Identification: Next, the selected template profiles were used to automatically identify their corresponding matches in the remaining tracks (drive-runs/flight lines and sensors) using a height map correlation-based strategy.
Finally, a point-pairing scheme was proposed for the extracted template and matched profiles from different tracks and an optimization function was developed to estimate the system calibration parameters-i.e., lever-arm and boresight angles for all the LiDAR units onboard the mobile Remote Sens. 2020, 12, 401 9 of 39 system-that minimized the discrepancy between the generated point pairs, as shown in Block 2 of the flowchart. In this section, we discuss the profile-based calibration strategy in detail, starting with the introduction of the mathematical model for point positioning in systems equipped with spinning multi-beam LiDAR units. Then, we describe the automated profile extraction strategy consisting of template profile selection and profile matching, followed by the point pairing scheme and optimization function that constitute the profile-based calibration strategy for a generic MLMS equipped with several spinning multi-beam LiDAR units.

Mathematical Model for LiDAR Point Positioning
A given point, I, acquired from an MLMS can be reconstructed in the mapping coordinate system using Equation (1), which is graphically illustrated in Figure 6. The reconstruction is done by applying a coordinate transformation based on the relationship between the laser unit frame, IMU body frame, and the mapping frame. For the laser unit frame, the origin is defined at the laser beams firing point and the z-axis is along the axis of rotation of the laser unit. For a spinning multi-beam LiDAR unit, each laser beam is fired at a fixed vertical angle, β; the horizontal angle, α, is determined based on the rotation of the unit; and the range, ρ, is defined by the distance between the firing point and its footprint. So, the coordinates of the 3D point, I, captured by the j th laser beam relative to the laser unit coordinate system, r Lu(t) I ( j), is defined by Equation (2). The laser unit frame is related to the IMU body frame by a rigidly defined lever arm, r b Lu , and boresight matrix, R b Lu . The GNSS/INS integration provides the time dependent position, r m b(t) , and rotation, R m b(t) , relating the IMU body frame and mapping coordinate systems.
Remote Sens. 2020, 12, 401 9 of 39 the flowchart. In this section, we discuss the profile-based calibration strategy in detail, starting with the introduction of the mathematical model for point positioning in systems equipped with spinning multi-beam LiDAR units. Then, we describe the automated profile extraction strategy consisting of template profile selection and profile matching, followed by the point pairing scheme and optimization function that constitute the profile-based calibration strategy for a generic MLMS equipped with several spinning multi-beam LiDAR units.

Mathematical Model for LiDAR Point Positioning
A given point, I, acquired from an MLMS can be reconstructed in the mapping coordinate system using Equation (1), which is graphically illustrated in Figure 6. The reconstruction is done by applying a coordinate transformation based on the relationship between the laser unit frame, IMU body frame, and the mapping frame. For the laser unit frame, the origin is defined at the laser beams firing point and the z-axis is along the axis of rotation of the laser unit. For a spinning multi-beam LiDAR unit, each laser beam is fired at a fixed vertical angle, ; the horizontal angle, , is determined based on the rotation of the unit; and the range, , is defined by the distance between the firing point and its footprint. So, the coordinates of the 3D point, I, captured by the laser beam relative to the laser unit coordinate system, ( ) ( ), is defined by Equation (2). The laser unit frame is related to the IMU body frame by a rigidly defined lever arm, , and boresight matrix, . The GNSS/INS integration provides the time dependent position, ( ) , and rotation, ( ) , relating the IMU body frame and mapping coordinate systems.

Automated Profile Extraction
In this section, we present an automated approach to extract an optimal set of thinly sliced profiles from the different tracks to be used for calibration. The conceptual basis for profile-based calibration is to use the thinly sliced profiles extracted from different tracks and treat them as two

Automated Profile Extraction
In this section, we present an automated approach to extract an optimal set of thinly sliced profiles from the different tracks to be used for calibration. The conceptual basis for profile-based calibration is to use the thinly sliced profiles extracted from different tracks and treat them as two dimensional entities to minimize the discrepancy between them in the vertical and along profile directions. In other words, each profile will contribute to discrepancy minimization in two directions. Having stated the overview of the contribution of a profile towards calibration, we proceed to state the desired characteristics of a profile that would ensure its contribution towards establishing sufficient control in three-dimensional space for an accurate calibration.

1.
Since a profile contributes towards discrepancy minimization in the vertical and profile-length directions, it should have a unique definition along these directions. This implies that profiles which are mainly constituted of points with planar neighborhoods should not be used for calibration. Rather, profiles with a majority of linear neighborhoods would facilitate an accurate point pairing for discrepancy minimization in the two aforementioned directions.

2.
In order to ensure discrepancy removal in the vertical and profile-length directions, a profile should not be monotonous, i.e., there should be variability within the points constituting the profile. This criterion is ensured by selecting profiles that are comprised of multiple linear segments with sufficient angular variability of the fitted line segments.
The above two conditions constitute the basis for the development of an automated algorithm to identify and extract template profiles from one of the tracks covering the calibration test field that would be used for calibration. In the following subsections, we first develop an algorithm for automated template profile selection and extraction, followed by proposing a strategy for automated profile matching within the remaining tracks for calibration.

Template Profile Selection
As discussed earlier, template profiles are selected by analyzing the angular variability of individual line segments constituting the profile. The selection of template profiles from the given point clouds captured from different tracks is conducted heuristically. In other words, the entire point cloud is parsed to analyze the characteristics of candidate profiles and decide whether each profile qualifies as a template for calibration. One should note that each of the steps involved in the proposed template profile selection algorithm is illustrated using a sample 3D point cloud illustrated in Figure 7.
In order to select template profiles for calibration, we first determine the bounding box that contains the union of the point clouds captured from all the tracks, denoted henceforth as the parent bounding box, shown as the black rectangle encompassing the 3D point cloud in Figure 7. The parent bounding box is tiled using pre-defined dimensions (L Tile ), as shown in Figure 7 as white grid lines. The desired length L p and depth D p of profiles is also specified by the user and these parameters are used as the basis for extracting the profiles that will be assessed for their validity to be used as a template profile. Once the point clouds are tiled collectively, each tile corresponding to each track is processed separately to extract the best candidate template profiles along as well as across flight/driving direction. The best profile in each direction is designated to be the one with the highest angular variability of the fitted line segments within the profile. Within each tile, a seed point for profile extraction is chosen at one of the corners of the tile to extract two profiles with the designated profile width-one along the flight/driving direction and another across the flight/driving direction-extending on either sides of the seed point to result in the desired profile length. One should note that the extracted profile points are allowed to extend beyond the tile boundaries when the seed points are located at the edges. This is done to ensure that the algorithm does not miss the extraction of a qualifying template profile that is located at the edge between two adjacent tiles. This profile extraction step is repeated by shifting the seed points throughout the tile bounds in the X and Y directions with a step size of half the designated profile length, as shown in Figure 7 by the red circles inside the top left tile. Once all the profiles within the tile are extracted, each profile is treated as a two-dimensional entity, so all the points within the profile are transformed to obtain the corresponding two-dimensional coordinates which are aligned along the vertical and profile-length directions. So, the variability of the points along the profile depth (which should be small enough to treat the profile as a 2D entity) is completely eliminated. Next, the transformed 2D profile points undergo the first screening criterion (as stated in Section 3.2), which targets the removal of any profile that is predominantly planar, such as profiles extracted along a building wall face, since point pairings along such profiles cannot be used to minimize discrepancies in the vertical and profile-length directions. Each 2D point within the profile is labeled as linear or planar by conducting a Principal Component Analysis (PCA) [41] of its neighborhood (with a user-defined search radius determined based on the local point spacing [42]) and any profile with a majority of planar points is immediately rejected. A sample of profiles-one accepted and one rejected-using this linearity criterion is shown in Figure 8, where the blue points denote planar neighborhoods and red points denote linear neighborhoods. The retained profiles undergo further processing wherein the identified linear neighborhood points within the profile are used to conduct a 2D segmentation to cluster the points into different line segments. One should note that the line segments are tested on their lengths and the number of points to ensure that too short or sparsely populated line segments are not considered for further assessment of angular variability. One such profile is shown in Figure 9, where the black points denote the points classified as a part of planar neighborhoods (thus, not used for 2D segmentation), white points denote the points that were classified as belonging to linear neighborhoods but not included as part of any of the segmented lines due to either the sparse nature or short length of resultant line segments from such points. The rest of the points are colored according to the individual line segments they belong to. The angles subtended by the fitted line segments and the axis denoting the direction along the profile length are also computed. For the profile displayed in Figure 9, the line segments have angles of 3.08 • , 89.03 • , and 0.32 • for blue, green, and red colored clusters, respectively. Finally, these recorded angles of each qualifying line segment within a profile are used to derive the angular variability within the profile. All the profiles oriented along the flight/driving direction within a tile are sorted based on the angular variability and the one with the highest angular variability that also exceeds a user-defined variability threshold is selected and stored as the template profile within that tile. The same is repeated to select the best profile across the flight/driving direction within the tile. One should note that within each tile, the point cloud from each of the tracks/sensors is analyzed individually and the best profile within the tile is extracted by comparing the angular variability of the profiles coming from all the tracks covering this tile. The same approach is repeated for each tile within the parent bounding box. So, we will have a maximum of two qualifying template profiles from each tile, as shown in Figure 7 as red lines within each tile with a central solid circle denoting the corresponding seed point. One should note that there might be tiles with no template profiles due to none of the candidate profiles within the tile satisfying the criteria for template profile selection. Given the heuristic approach that exhaustively searches through the entire point cloud, the resultant set of template profiles is bound to have the best possible distribution of profile orientation along and across the flight/driving direction for a given test field, as required for an accurate system calibration.
As discussed earlier, template profiles are selected by analyzing the angular variability of individual line segments constituting the profile. The selection of template profiles from the given point clouds captured from different tracks is conducted heuristically. In other words, the entire point cloud is parsed to analyze the characteristics of candidate profiles and decide whether each profile qualifies as a template for calibration. One should note that each of the steps involved in the proposed template profile selection algorithm is illustrated using a sample 3D point cloud illustrated in Figure  7.  of two qualifying template profiles from each tile, as shown in Figure 7 as red lines within each tile with a central solid circle denoting the corresponding seed point. One should note that there might be tiles with no template profiles due to none of the candidate profiles within the tile satisfying the criteria for template profile selection. Given the heuristic approach that exhaustively searches through the entire point cloud, the resultant set of template profiles is bound to have the best possible distribution of profile orientation along and across the flight/driving direction for a given test field, as required for an accurate system calibration.

Matching Profile Identification
Once the template profiles are extracted with each template coming from a specific track (denoted henceforth as the reference track), the next step is to identify the corresponding matching profiles in the remaining (non-reference) tracks. In this research, we propose an automated profile matching technique based on height map correlation. Using the seed point for each extracted template profile as the center, a template height map is generated for points belonging to a square window around the seed point with dimensions equal to the extracted profile length . The

Matching Profile Identification
Once the template profiles are extracted with each template coming from a specific track (denoted henceforth as the reference track), the next step is to identify the corresponding matching profiles in the remaining (non-reference) tracks. In this research, we propose an automated profile matching technique based on height map correlation. Using the seed point for each extracted template profile as the center, a template height map is generated for points belonging to a square window around the seed point with dimensions equal to the extracted profile length L p . The generated height map is gridded into smaller cells whose dimensions are determined based on the local point spacing of the point cloud. Each cell within the template height map is assigned the 95th percentile height of all points inside the cell. Next, a larger search window with the tile dimensions used for template profile selection (L Tile ) is created around the seed point and the points from the non-reference tracks covering this window are isolated to generate a gridded search height map. Now, a moving window of size L p is used to find the location within the search height map that has the highest correlation (more than 90%) with the previously established template height map. A highly correlated 95th percentile height map indicates the similarity of the two neighborhoods being compared. So, the central point for this highest correlation location is designated as the seed point for the matching profile within the non-reference track and is used to extract the matching profile with the same profile length and depth as the template profile. This procedure is applied to each non-reference track and for each extracted template profile to find and extract the matching profiles for calibration. A sample template profile (colored by its individual line segments) is shown in Figure 10a and the same profile along with its corresponding matches in 17 other tracks having discrepancies ranging up to 2 m are illustrated in Figure 10b, where a unique color (blue to red) is assigned to each non-reference track and the points belonging to the reference track are colored in black.
Remote Sens. 2020, 12, 401 13 of 39 generated height map is gridded into smaller cells whose dimensions are determined based on the local point spacing of the point cloud. Each cell within the template height map is assigned the 95 th percentile height of all points inside the cell. Next, a larger search window with the tile dimensions used for template profile selection ( ) is created around the seed point and the points from the non-reference tracks covering this window are isolated to generate a gridded search height map. Now, a moving window of size is used to find the location within the search height map that has the highest correlation (more than 90%) with the previously established template height map. A highly correlated 95 th percentile height map indicates the similarity of the two neighborhoods being compared. So, the central point for this highest correlation location is designated as the seed point for the matching profile within the non-reference track and is used to extract the matching profile with the same profile length and depth as the template profile. This procedure is applied to each nonreference track and for each extracted template profile to find and extract the matching profiles for calibration. A sample template profile (colored by its individual line segments) is shown in Figure  10a and the same profile along with its corresponding matches in 17 other tracks having discrepancies ranging up to 2 m are illustrated in Figure 10b, where a unique color (blue to red) is assigned to

Point Pairing Scheme and Optimization Function for Calibration
Once the template and matching profiles are extracted, the next step is to proceed with the calibration and establish a suitable point pairing scheme along with an appropriate objective function for the estimation of mounting parameters for all the LiDAR units onboard the mobile system. In order to conduct a profile-based calibration, the objective is to minimize the discrepancy between the mapping frame coordinates for all the point pairs formed between different tracks for each extracted ( )

Point Pairing Scheme and Optimization Function for Calibration
Once the template and matching profiles are extracted, the next step is to proceed with the calibration and establish a suitable point pairing scheme along with an appropriate objective function for the estimation of mounting parameters for all the LiDAR units onboard the mobile system. In order to conduct a profile-based calibration, the objective is to minimize the discrepancy between the mapping frame coordinates for all the point pairs formed between different tracks for each extracted profile. Each pairing between conjugate points will result in a random misclosure vector → e , as given in Equation (3). The point pairs for each profile are formed by firstly sorting all the tracks for this profile in a decreasing order in terms of the number of points. Starting from the track with the second most number of points, each point within this track is used to find its closest counterpart in the first track, which is also closer than a pre-defined distance threshold. A similar procedure is adopted for all the tracks to pair them with the closest point found in the preceding tracks, starting with the reference one. Once the point pairs are formed for all the extracted profiles, the last step is to estimate the mounting parameters by minimizing the discrepancy between the point pair coordinates. However, one should note that the profile matching strategy proposed in Section 3.2.2 would ensure the profile uniqueness only along the profile length but not along its depth. For instance, for the hut-shaped target (approximately 0.60 m wide) shown in Figure 11a, the selected template profile is shown in red and the corresponding matched profile is shown in blue in Figure 11b,c. As evident from such illustration, the profile uniqueness is guaranteed in the vertical and profile-length directions, whereas the uniqueness cannot be ensured along the profile depth since the hut will have very similar variation for all the profiles extracted anywhere within the 0.60 m wide area of the hut-shaped target. So, in this case, non-conjugate point pairs will be formed, which will result in an additional non-random component → D in the misclosure vector, given in Equation (4). One should note that this non-random component is aligned in the profile-depth direction. So, as mentioned earlier, the objective function for calibration would aim to minimize the discrepancy between profile point pairs only in the vertical and corresponding profile-length directions. Mathematically, this could be achieved by defining a modified weight matrix (P ), which would nullify the non-random component of the misclosure vector, → D, as given in Equation (5) [43]. From the automated profile extraction strategy, the orientation of each profile is known, denoted henceforth by θ p , which is the angle between the profile length direction and the X-axis of the mapping frame coordinate system. Then, the rotation matrix given in Equation (6) will relate the mapping frame coordinate system (XYZ) to the profile coordinate system (UVW). Here, the U-axis is oriented along the profile length, V-axis is oriented along the profile depth, and W-axis is in the vertical direction, as shown in Figure 11b,c. The weight matrix, P XYZ (which depends on the point cloud accuracy), in the mapping coordinate system is transformed to a weight matrix, P UVW , in the profile coordinate system according to the law of error propagation (Equation (7)). Next, the weight matrix, P UVW , is modified by assigning zero weights to the elements corresponding to the direction along the profile depth (Equation (8)) and the modified weight matrix, P UVW , is transformed back into the mapping frame to obtain P XYZ using Equation (9). Finally, the obtained modified weight matrix is applied to the condition in Equation (4) to only retain the discrepancy between non-conjugate point pairs in the vertical and profile-length directions. The mounting parameters are estimated with the objective to minimize the resultant discrepancies between the non-conjugate point pairs. After each iteration of mounting parameters estimation, the point pairing is conducted again by identifying the closest points using the newly reconstructed point coordinates from the revised estimates of mounting parameters. The next section demonstrates the results of the proposed profile-based calibration strategy for several mobile mapping systems-both airborne and terrestrial-equipped with one or more LiDAR units.
Remote Sens. 2020, 12, 401 14 of 39 note that this non-random component is aligned in the profile-depth direction. So, as mentioned earlier, the objective function for calibration would aim to minimize the discrepancy between profile point pairs only in the vertical and corresponding profile-length directions. Mathematically, this could be achieved by defining a modified weight matrix ( ), which would nullify the non-random component of the misclosure vector, ⃗ , as given in Equation (5) [43]. From the automated profile extraction strategy, the orientation of each profile is known, denoted henceforth by , which is the angle between the profile length direction and the X-axis of the mapping frame coordinate system. Then, the rotation matrix given in Equation (6) will relate the mapping frame coordinate system ( ) to the profile coordinate system ( ). Here, the -axis is oriented along the profile length, -axis is oriented along the profile depth, and -axis is in the vertical direction, as shown in Figure  11 (b) and (c). The weight matrix, (which depends on the point cloud accuracy), in the mapping coordinate system is transformed to a weight matrix, , in the profile coordinate system according to the law of error propagation (Equation (7)). Next, the weight matrix, , is modified by assigning zero weights to the elements corresponding to the direction along the profile depth (Equation (8)) and the modified weight matrix, , is transformed back into the mapping frame to obtain using Equation (9). Finally, the obtained modified weight matrix is applied to the condition in Equation (4) to only retain the discrepancy between non-conjugate point pairs in the vertical and profile-length directions. The mounting parameters are estimated with the objective to minimize the resultant discrepancies between the non-conjugate point pairs. After each iteration of mounting parameters estimation, the point pairing is conducted again by identifying the closest points using the newly reconstructed point coordinates from the revised estimates of mounting parameters. The next section demonstrates the results of the proposed profile-based calibration strategy for several mobile mapping systems-both airborne and terrestrial-equipped with one or more LiDAR units.

Experimental Results and Discussion
In this section, we present the calibration results for two different UAV-based LiDAR systems (equipped with one LiDAR unit each) and one terrestrial mobile LiDAR system (equipped with four LiDAR units). For each of the systems, we evaluated the performance of the two components of the proposed profile-based calibration -(1) automated profile selection and matching results, and (2) resultant mounting parameter estimates after calibration. The automated profile selection algorithm was validated by checking the fitted line segments and their angular variation within each of the selected profiles. Furthermore, the profile matching strategy was verified by a qualitative inspection of the matched profiles as well as the report for all the tracks where each template profile is matched. The parameters used for automated template profile selection and profile matching for the three MLMS are listed in Table 3. Owing to the similarity of the calibration test fields and flight configuration for both airborne MLMS, the same values are assigned for the different parameters used. The tile dimensions for airborne MLMS were set to 7 m since the average spacing between the different huts ranged between 6 m to 9 m; thus, a tile dimension of 7 m ensured that no two huts lie within the same tile. However, in the case of terrestrial MLMS where the test site covered a large area with a vast collection of similar-looking entities (in this case, houses lined up on either side of the road), the tile size was increased to 20 m to avoid redundant profiles for calibration. The profile length was set to 2 m for both airborne MLMS, which was chosen based on the calibration test field such that there would be sufficient angular variability within the total length of the extracted profile. In the case of terrestrial MLMS, the profile length was increased to 10 m, which would result in profiles spanning the road surface to building roof. The profile depth should be a small value in order to ensure that the resultant extracted profile is thin enough to be treated as a two-dimensional entity for calibration. Another criterion for choosing profile depth is to ensure a unique definition of the profile in the vertical and profile-length directions and minimum variability in the profile-depth direction. In this research, a profile depth of 0.10 m was chosen for airborne and terrestrial MLMS. The chosen value for profile depth applied to all the MLMS, even though the VLP16 has a higher local point spacing than VLP32C/HDL32E. This is due to the fact that the local point spacing was less than 0.10 m for all the sensors, thus ensuring that there were sufficient points within the extracted profile to define it uniquely in the vertical and profile-length directions. Moreover, the terrain of the calibration test fields ensured that there was negligible variability within 0.10 m along the profile-depth direction. One should note that the calibration results are not sensitive to the profile depth or profile length thresholds and these can be approximated based on the characteristics of the calibration test field and local point spacing. Finally, any profile with an angular variance greater than 45 • was extracted as a template. For profile matching, the template and search height map cell size was set to 0.20 m based on the local point spacing and the height map correlation threshold was set to 90%, i.e., any profile was deemed as a match to the corresponding template profile if the correlation between template and search height maps was more than 90%. The extracted template and matched profiles were used for calibrating the different MLMS using the proposed profile-based calibration strategy, whose results are validated using the following criteria:

1.
Square root of aposteriori variance factor (σ 0 ) after calibration, which was used as a representative of the discrepancy between the established point pairs.

2.
Standard deviation of the estimated parameters, which was expected to indicate the respective confidence intervals.

3.
Qualitative assessment based on the alignment between profiles from different tracks (drive-runs/flight lines and sensors) before and after calibration, which would also be reflected by the aposteriori variance factor.
A previous study by Ravi et al. [32,33] proposed a manually assisted feature-based calibration strategy for MLMS and proved its ability to attain the best possible accuracy in keeping with the manufacturer's specifications for the onboard sensors. Therefore, in addition to conducting a profile-based calibration, all the systems were also calibrated using the existing manually assisted feature-based calibration technique. Several planar features, such as highly reflective sign boards, hut facades, building facades, roof patches, and ground patches, were manually extracted within each of the calibration test fields in order to conduct the feature-based calibration. A comparative analysis of the two calibration approaches wasconducted based on the following criteria:

1.
Similarity ofσ 0 and standard deviation of the estimated parameters from the two approaches.

2.
Similarity of estimated mounting parameters from the two calibration strategies by assessing their impact on the reconstructed profiles while considering the expected point positioning accuracy. This similarity was qualitatively and quantitatively evaluated as follows: • Qualitative evaluation: visual assessment of the alignment between the reconstructed profile points using parameters from the two calibration strategies. • Quantitative evaluation: difference in 3D mapping frame coordinates (using Equation (1)) of points computed based on the different calibration parameters.
In the following subsections, we provide detailed results for the extracted template profiles, profile matching, and qualitative and quantitative evaluation of calibration results for the three MLMS.

Airborne MLMS 1
The automated template profile selection resulted in a total of 23 template profiles being extracted in the calibration test field for the airborne MLMS 1 and these are shown in Figure 12, where each template profile is colored by its individual fitted line segments. The resultant angular variation of the fitted line segments within the template profiles are listed in Table 4. Proceeding further, the template profiles were matched to find corresponding matching profiles in the remaining tracks, which were used for point pairing and subsequent system calibration. The overall layout of all the automatically extracted template and matched profiles in the calibration test field is shown in Figure 13. The figure displays the profile IDs along with the track in which the corresponding template profile was selected. The matched profiles (including the template) for the 23 extracted profiles are shown individually in Figure 14, where each profile is colored according to the track in which it was captured. One should note that the observed discrepancy between different tracks using the manually measured initial approximations of the mounting parameters ranged up to 2 m in the horizontal directions and up to 0.50 m in the vertical direction. Table 5 shows a list of all tracks in which each profile was extracted, where the track corresponding to template profile is highlighted in yellow. The automatically extracted template and matched profiles are used to perform a system calibration using the proposed point pairing scheme to minimize the discrepancy between the profiles captured from different flight lines. The manually measured initial approximations, profile-based calibration results, and feature-based calibration results for the mounting parameters are listed in Table 6 along with the standard deviations for each of the parameters as obtained from profile-based and feature-based calibrations. horizontal directions and up to 0.50 m in the vertical direction. Table 5 shows a list of all tracks in which each profile was extracted, where the track corresponding to template profile is highlighted in yellow. The automatically extracted template and matched profiles are used to perform a system calibration using the proposed point pairing scheme to minimize the discrepancy between the profiles captured from different flight lines. The manually measured initial approximations, profilebased calibration results, and feature-based calibration results for the mounting parameters are listed in Table 6 along with the standard deviations for each of the parameters as obtained from profilebased and feature-based calibrations.          The calibration results were evaluated qualitatively and quantitatively based on the previously mentioned criteria. Figure 15 shows the reconstructed profiles before and after profile-based calibration, where each of the profiles are colored by the track in which they were captured. It can be seen that the alignment between different profiles exhibits significant improvement after calibration for all the profiles, which is also reflected by theσ 0 value of 1.77 cm listed in Table 6. The ability of the proposed profile-based calibration technique to reduce the misalignment from about 2 m to less than 2 cm proves its feasibility for cases where the initial approximations of the mounting parameters are significantly different from the true values. Theσ 0 value for manually-assisted feature-based calibration is 2.20 cm and the standard deviations of the estimated parameters from the two calibration approaches are similar, which indicates that the two approaches perform at par in terms of the resultant alignment after calibration as well as the confidence interval of the estimated parameters. The profiles are reconstructed using the mounting parameters from the two different calibration approaches and are shown in Figure 16, where the points in blue and red are reconstructed using profile-based and feature-based calibration results, respectively. Figure 16 indicates that the mapping frame coordinates derived from the two sets of mounting parameters are similar for all the profiles. Table 7 shows the resultant mean, standard deviation, and root mean squared error values for the difference between mapping frame coordinates (X, Y, and Z) derived using the estimated mounting parameters from both calibration strategies. The differences in mapping frame coordinates can be observed to lie in the range of 1 to 5 cm, which is less than the expected point positioning accuracy of ±8-9 cm based on the individual accuracies of the involved sensors and reconstruction parameters (such as range and horizontal FOV of reconstruction). Since the accuracy of feature-based calibration strategy is already proven in prior work, the difference of 1 to 5 cm in the mapping frame coordinates derived using profile-based and feature-based calibration strategies validates the accuracy of the proposed calibration procedure.

P0
Remote Sens. 2020, 12, 401 22 of 39 calibration approaches and are shown in Figure 16, where the points in blue and red are reconstructed using profile-based and feature-based calibration results, respectively. Figure 16 indicates that the mapping frame coordinates derived from the two sets of mounting parameters are similar for all the profiles. Table 7 shows the resultant mean, standard deviation, and root mean squared error values for the difference between mapping frame coordinates (X, Y, and Z) derived using the estimated mounting parameters from both calibration strategies. The differences in mapping frame coordinates can be observed to lie in the range of 1 to 5 cm, which is less than the expected point positioning accuracy of ±8-9 cm based on the individual accuracies of the involved sensors and reconstruction parameters (such as range and horizontal FOV of reconstruction). Since the accuracy of feature-based calibration strategy is already proven in prior work, the difference of 1 to 5 cm in the mapping frame coordinates derived using profile-based and feature-based calibration strategies validates the accuracy of the proposed calibration procedure.

Airborne MLMS 2
There were a total of 19 template profiles extracted in the calibration test field for the airborne MLMS 2. Since the test field was similar to the one used for airborne MLMS 1, only four out of the nineteen profiles were used for qualitative assessment of the results from each step within profilebased calibration. The chosen subset of profiles includes one along the wall of a building, two profiles over hut-shaped targets (one oriented along and another across the flying direction), and the last profile along a pole in the test area. These four profiles are shown individually in Figure 17, where each template profile is colored by its fitted line segments. One should note that the point cloud from airborne MLMS 2 (VLP16 Puck Lite) has a lower point density compared to that from airborne MLMS

Airborne MLMS 2
There were a total of 19 template profiles extracted in the calibration test field for the airborne MLMS 2. Since the test field was similar to the one used for airborne MLMS 1, only four out of the nineteen profiles were used for qualitative assessment of the results from each step within profile-based calibration. The chosen subset of profiles includes one along the wall of a building, two profiles over hut-shaped targets (one oriented along and another across the flying direction), and the last profile along a pole in the test area. These four profiles are shown individually in Figure 17, where each template profile is colored by its fitted line segments. One should note that the point cloud from airborne MLMS 2 (VLP16 Puck Lite) has a lower point density compared to that from airborne MLMS 1 (VLP32C), but as stated earlier, the chosen profile depth of 0.10 m is sufficient to extract points that define the profile uniquely in the vertical and profile-length directions. The resultant angular variation of the fitted line segments within each profile is listed in Table 8. Next, the template profiles were matched to find corresponding profiles in the remaining tracks. The overall layout of all the automatically extracted template and matched profiles in the calibration test field is shown in Figure 18 along with the profile ID and reference track ID for each profile. The matched profiles (including the template) for the subset of extracted profiles are shown individually in Figure 19, where each profile is colored by the track in which it was captured. The misalignment between the different tracks due to initial approximations of the mounting parameters ranges up to 3 m in the horizontal directions and about 0.50 m in the vertical direction. Table 9 shows all the tracks where each profile has been extracted, with the track corresponding to the template profile highlighted in yellow. The automatically extracted template and matched profiles were used to perform a system calibration. The manually measured initial approximations, profile-based calibration results, and feature-based calibration results for the mounting parameters are listed in Table 10, along with their respective standard deviations.
Remote Sens. 2020, 12, 401 26 of 39 1 (VLP32C), but as stated earlier, the chosen profile depth of 0.10 m is sufficient to extract points that define the profile uniquely in the vertical and profile-length directions. The resultant angular variation of the fitted line segments within each profile is listed in Table 8. Next, the template profiles were matched to find corresponding profiles in the remaining tracks. The overall layout of all the automatically extracted template and matched profiles in the calibration test field is shown in Figure  18 along with the profile ID and reference track ID for each profile. The matched profiles (including the template) for the subset of extracted profiles are shown individually in Figure 19, where each profile is colored by the track in which it was captured. The misalignment between the different tracks due to initial approximations of the mounting parameters ranges up to 3 m in the horizontal directions and about 0.50 m in the vertical direction. Table 9 shows all the tracks where each profile has been extracted, with the track corresponding to the template profile highlighted in yellow. The automatically extracted template and matched profiles were used to perform a system calibration. The manually measured initial approximations, profile-based calibration results, and feature-based calibration results for the mounting parameters are listed in Table 10, along with their respective standard deviations.       Profile ID   Profile ID Figure 19. Airborne MLMS 2: automatically extracted template and matched sample profiles, where each profile is colored by the track in which it was captured (discrepancy between tracks ≈ 3 m). Table 9. Airborne MLMS 2: list of tracks (flight lines or sensors) in which each profile has been extracted.

VLP16 Puck Lite T-1 T-2 T-3 T-4 T-5 T-6 T-7 T-8 T-9 T-10 T-11 T-12 T-13 T-14 T-15 T-16 T-17 T-18
The track corresponding to the template profile highlighted in yellow.  The resultant sample profiles before and after calibration are depicted in Figure 20, indicating a significant improvement in the alignment between different tracks for each of the profiles. The quality of alignment is also reflected by theσ 0 value of 1.98 cm listed in Table 10. In other words, the proposed calibration strategy reduces the misalignment from approximately 3 m to less than 2 cm. Theσ 0 value of 2.33 cm for feature-based calibration and the similarity in standard deviations of the estimated parameters from the two calibration approaches indicate the compliance of the results from the two approaches in terms of the alignment after calibration as well as the confidence interval of estimated parameters. Figure 21 shows the qualitative comparison for each of the sample profiles reconstructed using the mounting parameters from the two different calibration strategies (blue denotes the reconstruction after profile-based calibration and red denotes the same for feature-based calibration). Table 11 lists the mean, standard deviation and RMSE of differences between the mapping frame coordinates of all the profile points as reconstructed using the estimated mounting parameters from the two calibration strategies. Again, the differences are observed to lie in the range of 0.5 to 3 cm, which is less than the expected point positioning accuracy of ±8-9 cm based on the individual accuracies of the involved sensors and reconstruction parameters. This again validates the accuracy of the profile-based calibration based on its comparison with the previously established feature-based calibration strategy.

Terrestrial MLMS
There were a total of 11 template profiles extracted in the calibration test field for this system, as shown individually in Figure 22, where each template profile is colored by its fitted line segments. The resultant angular variation of the fitted line segments within each profile is listed in Table 12. Next, the template profiles were matched to find corresponding profiles in the remaining tracks. One should note that there were four LiDAR units mounted on the terrestrial MLMS, each capturing point clouds from four drive-runs, so there was a total of sixteen tracks used for calibration. The overall layout of all the extracted template and matched profiles is shown in Figure 23. The matched profiles (including the template) for a sample of extracted profiles are shown in Figure 24, where each profile is colored by the track in which it was captured. The initial misalignment between different tracks before calibration is observed to be about 1 m in the horizontal and vertical directions. Table 13 depicts all the tracks in which each profile was identified and extracted with the track corresponding to template profile highlighted in yellow. The extracted template and matched profiles were used to perform a system calibration and the manually measured initial approximations, profile-based calibration results, and feature-based calibration results for the mounting parameters are listed in Table 14

Terrestrial MLMS
There were a total of 11 template profiles extracted in the calibration test field for this system, as shown individually in Figure 22, where each template profile is colored by its fitted line segments. The resultant angular variation of the fitted line segments within each profile is listed in Table 12. Next, the template profiles were matched to find corresponding profiles in the remaining tracks. One should note that there were four LiDAR units mounted on the terrestrial MLMS, each capturing point clouds from four drive-runs, so there was a total of sixteen tracks used for calibration. The overall layout of all the extracted template and matched profiles is shown in Figure 23. The matched profiles (including the template) for a sample of extracted profiles are shown in Figure 24, where each profile is colored by the track in which it was captured. The initial misalignment between different tracks before calibration is observed to be about 1 m in the horizontal and vertical directions. Table 13 depicts all the tracks in which each profile was identified and extracted with the track corresponding to template profile highlighted in yellow. The extracted template and matched profiles were used to perform a system calibration and the manually measured initial approximations, profile-based calibration results, and feature-based calibration results for the mounting parameters are listed in Table 14 along with the respective standard deviations.           The resultant profiles before and after calibration are depicted in Figure 25, which indicates a significant improvement in the alignment between different tracks for each of the profiles, also reflected by theσ 0 value of 1.72 cm. In other words, the profile-based calibration succeeds in reducing the misalignment from 1 m to less than 2 cm. Again, the similarity in the standard deviations of the estimated parameters from the two calibration strategies along with theσ 0 value of 1.57 cm for feature-based calibration indicate that the two approaches attain similar resultant alignment and confidence intervals for the estimated parameters after calibration. The qualitative comparison between the mounting parameters estimated from profile-based and feature-based calibration strategies are shown for some sample profiles in Figure 26, wherein the points in blue are reconstructed using the parameters from profile-based calibration and those in red are the same points reconstructed using the mounting parameters derived from feature-based calibration. Table 15 lists the mean, standard deviation and RMSE of differences between the mapping frame coordinates of all the profile points as reconstructed from the mounting parameters derived using the two different calibration strategies. Here, the differences between the mapping coordinates derived from the two sets of mounting parameters is observed to be less than 2 cm, thus in coherence with the expected point positioning accuracy of 2-4 cm based on the manufacturer's specifications for the sensor accuracies. As observed earlier for the airborne mobile mapping systems, it can be seen that the results from profile-based calibration for terrestrial systems with multiple LiDAR units are also consistent with the corresponding feature-based calibration results.

T-1 T-2 T-3 T-4 T-1 T-2 T-3 T-4 T-1 T-2 T-3 T-4 T-1 T-2 T-3 T-4
Remote Sens. 2020, 12, 401 33 of 39 of the estimated parameters from the two calibration strategies along with the value of 1.57 cm for feature-based calibration indicate that the two approaches attain similar resultant alignment and confidence intervals for the estimated parameters after calibration. The qualitative comparison between the mounting parameters estimated from profile-based and feature-based calibration strategies are shown for some sample profiles in Figure 26, wherein the points in blue are reconstructed using the parameters from profile-based calibration and those in red are the same points reconstructed using the mounting parameters derived from feature-based calibration. Table  15 lists the mean, standard deviation and RMSE of differences between the mapping frame coordinates of all the profile points as reconstructed from the mounting parameters derived using the two different calibration strategies. Here, the differences between the mapping coordinates derived from the two sets of mounting parameters is observed to be less than 2 cm, thus in coherence with the expected point positioning accuracy of 2-4 cm based on the manufacturer's specifications for the sensor accuracies. As observed earlier for the airborne mobile mapping systems, it can be seen that the results from profile-based calibration for terrestrial systems with multiple LiDAR units are also consistent with the corresponding feature-based calibration results.

Conclusions and Recommendations for Future Research
In this paper, we presented a novel approach for a fully automated system calibration procedure that does not rely on specially designed targets or a manual feature extraction procedure. Instead, we propose a profile-based calibration technique wherein thin profiles are automatically selected and extracted from the point clouds based on a two-dimensional linear segmentation within the profile and the resultant angular variability of the fitted lines. Moreover, the selected template profiles are used to find their corresponding profiles in the remaining tracks using a height map correlationbased matching strategy. The extracted template and matched profiles are then used for performing the LiDAR system calibration by pairing points within the profile coming from different driveruns/flight lines and/or sensors and minimizing the discrepancy between the point pair coordinates in the vertical and profile-length directions. Finally, the proposed fully automated calibration technique was experimentally validated using two different airborne systems with a single LiDAR unit and a terrestrial system with four LiDAR units.
The proposed strategy was seen to reduce the misalignment between different tracks from approximately 2 to 3 m before calibration down to less than 2 cm after calibration for airborne as well as terrestrial MLMS. These results emphasize the achievable relative accuracy of the proposed strategy capability along with its ability to derive accurate estimates of mounting parameters even on starting from significantly different initial approximations. Furthermore, the accuracy of the obtained results was also verified by comparing them with results from a previously proven

Conclusions and Recommendations for Future Research
In this paper, we presented a novel approach for a fully automated system calibration procedure that does not rely on specially designed targets or a manual feature extraction procedure. Instead, we propose a profile-based calibration technique wherein thin profiles are automatically selected and extracted from the point clouds based on a two-dimensional linear segmentation within the profile and the resultant angular variability of the fitted lines. Moreover, the selected template profiles are used to find their corresponding profiles in the remaining tracks using a height map correlation-based matching strategy. The extracted template and matched profiles are then used for performing the LiDAR system calibration by pairing points within the profile coming from different drive-runs/flight lines and/or sensors and minimizing the discrepancy between the point pair coordinates in the vertical and profile-length directions. Finally, the proposed fully automated calibration technique was experimentally validated using two different airborne systems with a single LiDAR unit and a terrestrial system with four LiDAR units.
The proposed strategy was seen to reduce the misalignment between different tracks from approximately 2 to 3 m before calibration down to less than 2 cm after calibration for airborne as well as terrestrial MLMS. These results emphasize the achievable relative accuracy of the proposed strategy capability along with its ability to derive accurate estimates of mounting parameters even on starting from significantly different initial approximations. Furthermore, the accuracy of the obtained results was also verified by comparing them with results from a previously proven manually assisted feature-based calibration strategy. The quantitative difference between the two approaches derived on the basis of the comparison of mapping frame coordinates was found to be in the range of 1 to 5 cm for airborne MLMS and 0.5 to 3 cm for terrestrial MLMS. Both of these values indicate that the two calibration approaches have a similar performance within the range of expected point positioning accuracy derived based on the sensor specifications. Therefore, the fully automated targetless profile-based calibration strategy proposed in this paper is proven in terms of its feasibility and achieved accuracy for airborne as well as terrestrial mobile mapping systems equipped with one or more 3D LiDAR units.
Future research in this area will focus on evaluating the applicability of the proposed profile-based calibration technique for systems equipped with 2D LiDAR units. Moreover, the proposed approach will be extended to suggest a similar fully automated calibration technique for the simultaneous calibration of cameras and LiDAR units onboard a mobile mapping system.