A Scalable Method to Improve Large-Scale Lidar Topographic Differencing Results

Abstract

Differencing digital terrain models (DTMs) generated from multitemporal airborne light detection and ranging (lidar) data provide accurate and detailed information about three-dimensional (3D) changes on the Earth. However, noticeable spurious errors along flight paths are often included in the differencing results, hindering the accurate analysis of the topographic changes. This paper proposes a new scalable method to alleviate the problematic systematic errors with a high degree of automation in consideration of the practical limitations raised when processing the rapidly increasing amount of large-scale lidar datasets. The proposed method focused on estimating the displacements caused by vertical positioning errors, which are the most critical error source, and adjusting the DTMs already produced as basic lidar products without access to the point cloud and raw data from the laser scanner. The feasibility and effectiveness of the proposed method were evaluated with experiments with county-level multitemporal airborne lidar datasets in Indiana, USA. The experimental results demonstrated that the proposed method could estimate the vertical displacement reasonably along the flight paths and improve the county-level lidar differencing results by reducing the problematic errors and increasing consistency across the flight paths. The improved differencing results presented in this paper are expected to provide more consistent information about topographic changes in Indiana. In addition, the proposed method can be a feasible solution to upcoming problems induced by rapidly increasing large-scale multitemporal lidar given recent active government-driven lidar data acquisition programs, such as the U.S. Geological Survey (USGS) 3D Elevation Program (3DEP).

1. Introduction

Differencing multitemporal topographic data is a fundamental technique with many valuable applications requiring three-dimensional (3D) Earth surface changes as prerequisite information [1,2,3,4,5,6,7,8,9]. Various platforms with different sensors from space (e.g., synthetic aperture radar (SAR) [2,3]) to the ground (e.g., terrestrial light detection and ranging (lidar) [4,5]) have been used to collect precise topographic information during the last several decades. Airborne lidar is superior to other tools, as it enables efficient and high-quality topography mapping [10,11]. Airborne lidar data have been widely used to monitor 3D changes, providing high accuracy and a high level of detail [12,13,14,15,16].

In order to compute accurate 3D changes between multitemporal airborne lidar data, numerous approaches have been developed [17,18,19]. One of the most frequently used lidar differencing methods is a digital elevation model (DEM) of difference (DoD) process [18] that subtracts a rasterized old elevation product from the corresponding current elevation product or vice versa. Generating rasterized elevation products, such as a digital terrain model (DTM), is a fundamental application of airborne lidar data [20]. A DTM is a basic elevation product derived from airborne lidar data. Therefore, DoD implementation for DTMs is straightforward, and a number of previous studies have demonstrated the DoD technique’s usefulness for landscape monitoring over a coastal area [21], a fluvial area [22], and a volcanic area [13] or for response analysis (e.g., landslide) of extreme precipitation [16,23] and earthquakes [7].

However, despite its prominent advantages, the differencing of airborne lidar data has multiple limitations usually caused by the different data quality between individual data collections [18]. Hence, various sophisticated approaches have been developed to solve the limitations and fully exploit the advantages of airborne lidar. One critical limitation is that inherent errors of airborne lidar data become considerably apparent when differencing multitemporal data [14,16,24]. The inherent errors raise abnormal elevation changes in the differencing results, severely hindering accurate 3D change monitoring. Therefore, previous studies corrected lidar data from the low level to thoroughly solve numerous potential error sources [14,16,25], such as inadequate system calibration, internal flight path misalignments, and the effects of land covers. A standard solution to decrease the error sources is the recalibration of lidar point clouds using flight trajectory files [16]. Some previous studies have used ground control features to calibrate lidar products [23]. However, such supplementary data were often absent, especially from national-scale lidar collection efforts, like the U.S. Geological Survey (USGS) 3D Elevation Program (3DEP) [26], as the program only requires derived 3D point clouds and DTMs. Therefore, an iterative closest point (ICP; [27]) method has frequently been employed to reduce abnormal elevation changes without additional data [28]. These previous approaches enabled us to precisely differentiate the actual elevation changes from the spurious changes. Nevertheless, they were implementable for limited lidar products with supplementary data or where conjugated features for ICP-based algorithms were readily extracted.

Meanwhile, many countries have been actively collecting nationwide 3D elevation data using airborne lidar to address the growing demand for high-quality height information [29,30,31]. As a result, the amount of airborne lidar data has reached an unprecedented level in terms of spatial and temporal coverage. It posed challenges in developing new scalable processing techniques to handle these ‘big data’ efficiently [31,32,33]. Differencing multitemporal topographic data also confronts such challenges such that the previous solutions against the remarkable artificial elevation changes are unscalable and often become impractical at larger spatial scales. The iterative process of large-scale data consumes tremendous computational time, and high-performance computing resources may be required to cover large areas. In addition, the appropriate supplementary data rarely exist for a large area, and the large area may include various land covers where the desirable conjugate features may not exist, which makes it challenging to utilize previously proposed approaches for large-area applications.

Considering the practical difficulties of the previous methods, this paper aims to develop a simple but computationally effective way to alleviate inherent error-induced abnormal topographic changes in large-scale lidar differencing results (i.e., DoD results). Instead of processing the low-level point cloud data, the proposed method estimates and adjusts the abnormal displacements mainly using DTMs, which are usually pre-generated as a primary elevation product derived from airborne lidar data. The proposed method’s effectiveness was assessed using county-level multitemporal airborne lidar datasets of Indiana, USA, where artificial elevation changes were manifested in their differencing results, which were performed using original DTMs from USGS 3DEP. The contributions of this paper can be summarized as (1) proposing a new scalable solution against the erroneous displacements frequently observed in DoD results without raw observation and (2) providing improved DoD results of Indiana 3DEP lidar datasets where the problematic displacements were reported but unsolved [24] with more reasonable terrain change information. This paper is organized as follows: First, Indiana statewide lidar datasets are introduced in Section 2. Then, Section 3 details the developed scalable method. The experimental results using Indiana lidar datasets and discussions are presented in Section 4, and this paper finally concludes with a summary in Section 5.

2. Indiana Statewide Multitemporal Lidar Datasets

Indiana statewide topographic information has been acquired twice using an airborne lidar system from 2011 to 2013 and 2016 to 2020. The statewide lidar data provide broad insights into the challenges in processing and analyzing high-quality and large-scale multitemporal lidar data. They were collected as parts of the 3DEP, a nationwide height information collection project managed by the USGS. The 3DEP aroused expectations to facilitate many valuable downstream applications using elevation data. Several studies have derived meaningful information from 3DEP products [24,30,34] by developing and implementing appropriate approaches to overcome the difficulties of handling these large-scale datasets.

The Indiana statewide multitemporal datasets were collected over the entirety of Indiana (94,324 km²) with an average point density of 1.6 pts/m² for the 2011–2013 dataset [24] and at least 2 pts/m² for the 2016–2020 dataset [35]. Through a complicated and time-intensive data processing workflow, two source products, the lidar point cloud (or LASer file format, LAS) and DTM, were generated from each data collection and became available to the public in 2021 [35]. The DTM provides essential terrain heights (bare-earth, hydro-flattened, gravity-related heights) with spatial resolutions of 1.524 m (5 ft) for the 2011–2013 dataset and 0.762 m (2.5 ft) for the 2016–2020 dataset. The LAS file includes detailed information about all recorded points, such as their exact 3D position, class (e.g., ground or noise), intensity, flight path ID, and global navigation satellite system (GNSS) time. The approximate total data sizes of the LAS and DTM products are 4.9 and 0.5 terabytes for the 2011–2013 data and 13.0 and 0.7 terabytes for the 2016–2020 data. To efficiently handle these voluminous bi-temporal Indiana lidar datasets, entire areas of the individual 92 Indiana counties were divided into square grids (tiles) with a size of 1524 m by 1524 m, and the lidar products (LAS and DTM) were separately archived according to these tiles for both periods [36]. Although the geometric accuracy of each dataset has not been reported publicly in detail, the 2016–2020 dataset should meet USGS quality level 2 [37]. This implies that the absolute vertical accuracy of the root mean square error (RMSE_z) is 0.1 m for non-vegetated areas, which has been announced as twice as accurate as the 2011–2013 dataset by the Indiana Geographic Information Council (IGIC) [35].

Scott et al. [24] first presented Indiana’s statewide elevation changes by designing and implementing an automatic DoD process for these big Indiana datasets. The significant topographic changes over the states, such as river fluvial changes, coastal erosion, and newly constructed roads, were effectively detected based on the DoD results. However, simultaneously, the spurious stripes along north-to-south (NS) and west-to-east (WE) directions were observed in the differenced DTM [38], making analysis of the DoD results challenging. These abnormal elevation changes were concrete evidence of the inherent errors included in Indiana datasets but have not been resolved due to a time- and labor-intensive process using flight trajectory files that have not been distributed to the public.

The striped elevation changes throughout Indiana’s statewide differencing results were likely incurred primarily by inconsistent positioning errors of the flight paths, particularly those of the 2011–2013 dataset [24]. Therefore, correcting the dissimilar 2011–2013 positioning errors can significantly reduce the problematic elevation changes in the DoD result. An ICP-based method has been applied to mitigate relative positioning errors between adjacent flight paths [39]. However, due to its iterative process, the ICP method is inefficient for the massive Indiana dataset. In addition, the ground points within the overlap areas of the 2011–2013 flight paths have not been defined in the provided LAS files. Ground point classification in the overlap areas is thus required prior to implementing the ICP-based method, which is also time-consuming.

The positioning errors can be divided into vertical and horizontal errors, which are assuredly included in the dataset. However, the vertical errors are likely more severe than the horizontal errors throughout Indiana’s DoD results. Unlike vertical errors, horizontal errors do not induce errors in the height values over plains. Contrarily, they bring out the remarkable artificial errors in slopes. Figure 1 shows examples of the county-level DoD results between the 2011–2013 and the 2016–2020 DTMs provided publicly as basic lidar products. Most of the remarkable artificial elevation changes in Figure 1 appear as long, broad stripes irrespective of the terrain slopes. Consequentially, it can be inferred that the vertical errors mainly pose difficulties in accurately analyzing the DoD results. In this respect, this paper developed a novel and practical solution for the vertical positioning errors of the 2011–2013 datasets to improve the quality of their differencing results.

3. Methodology

The proposed scalable method estimates the erroneous displacements mainly induced by dissimilar vertical positioning errors along the flight paths by comparing bi-temporal DTMs provided as source products of Indiana lidar datasets rather than time-intensive point cloud adjustment. The estimated displacement values were then utilized to adjust the 2011–2013 flight paths. Using the existing DTMs only, we developed a simple but efficient and practical method to alleviate the problematic errors in large-scale Indiana DoD results with a high degree of automation. The proposed method consists of three main steps: (1) extracting additional information for effective comparison of bi-temporal DTMs, (2) estimating the erroneous displacements between the 2011–2013 and the 2016–2020 DTMs, and (3) generating improved differencing results by adjusting the 2011–2013 DTMs, as shown in Figure 2.

This paper used the 2016–2020 DTMs as references, and the 2011–2013 DTMs were aligned with them. The advances in lidar systems have increased lidar position accuracy [40], and it was reported that the 3DEP products’ vertical accuracy have consistently increased [30]. As a result, the 2011–2013 dataset has lower vertical accuracy than the 2016–2020 dataset, and its vertical positioning errors can be reduced once aligned with the 2016–2020 dataset. Furthermore, the height values of the bi-temporal DTMs are mostly identical if there is no vertical positioning error. The Earth’s surface rarely changes evenly over a specific and broad area, especially like the striped areas of the individual 2011–2013 flight paths. Accordingly, the proposed method in this paper compared the height values of the bi-temporal DTMs using pairs of the 2011–2013 flight paths and their corresponding 2016–2020 paths, which were spatially overlapped, to examine whether they have certain artificial displacements. If the displacements were detected, they were regarded as induced by the vertical errors of the 2011–2013 flight path, and the proposed algorithm described in the following section alleviated them from the 2011–2013 DTMs. For brevity, the ‘target’ and the ‘reference’ datasets refer to the 2011–2013 and the 2016–2020 datasets for the rest of this paper. The entire process described in the following sections, except data download, was implemented using various packages with Python programming languages, such as the Geospatial Data Abstraction Library (GDAL), laspy, OpenCV, pyproj, and SciPy.

3.1. Extracting Additional Information

Because we observed a systematic striping pattern from the original DTM differencing results, we hypothesized that the striping pattern originated from issues within the individual flying strips. For this reason, we extracted additional information to effectively compare the height values from the target (2011–2013) and the reference (2016–2020) DTMs in this paper. Two different types of information were extracted from other data except for the DTMs prior to the following comparison step. One was water masks to mask out the water bodies, as lidar measurements are less reliable than ones over the ground. Hence, DTM values over the water body might show significant elevation differences between the two DTMs. The other was the boundary information of each flight path composing Indiana statewide datasets. The flight path boundary information would help independently adjust the vertical displacements of flight strips.

3.1.1. Water Mask

By the nature of the Earth’s surface, the values of the bi-temporal DTMs should be similar on a macroscale throughout the individual counties. However, some regions, like water bodies, unavoidably show significant and distinct differences between the target and the reference DTMs. As the area of the water bodies accounted for a considerable proportion of the entirety of Indiana, the water bodies could severely disrupt the comparison of the bi-temporal DTMs. Therefore, the water bodies, such as lakes, ponds, streams, and rivers, needed to be masked out during the comparison.

The water mask in each county was generated using National Hydrography Dataset Plus High Resolution (NHDPlus HR) available to the public [41]. NHDPlus HR provides valuable information on the U.S. inland waters with diverse feature classes [42]. Two classes related to the water surface, NHDArea and NHDWaterbody, were used to determine the water surface. Most NHDPlus HR products have been published on different dates according to hydrologic second-level units (or Hydrologic Unit (HU)-4 Subregion) with a scale of 1:24,000 or better. Entirety of Indiana comprised eight different HU-4 Subregions, of which the products were published in 2019 and 2020. The eight products were downloaded via USGS national map download v2.0 [41]. Indiana statewide water masks were generated by merging the two feature classes of the eight downloaded NHDPlus HR products and dividing them into the same tiles of the DTMs.

3.1.2. Flight Path Boundary Information

Because the bi-temporal datasets were compared using the flight path pairs, the spatial boundaries of the individual flight paths needed to be defined. Unfortunately, the flight path boundary has not been published as a separate product. However, instead, information about the flight paths was obtainable from the LAS files because the flight paths of both datasets were indexed by number, and the numbers were stored as point_source_id in the LAS files. The LAS files provide various information on individual points with their precise 3D coordinates, unlike the DTMs based on two-dimensional (2D) image coordinates. Therefore, the following procedure was required to define the flight path boundary by aggregating the individual multiple points’ point_source_ids of the LAS file with the DTM’s pixels.

A tiled LAS file was first imported as input data. While the imported LAS file contained numerous points from various objects, only points reflected from the ground were used to define flight path boundaries in this paper. The ground points were thus extracted based on the class information in the LAS file produced by data vendors before 3DEP data publishing. Then, an empty flight boundary data file was created as an output raster file with the same coordinate system of the tiled target DTM (1.524 m spatial resolution and 1000 pixels by 1000 pixels size), which corresponded to the imported LAS file. The output file’s pixel values were determined using the point_source_id (index flight path number) of the LAS file’s ground points located within individual pixels. The corresponding ground points for each pixel were identified by converting their horizontal 2D coordinates, x and y, into the 2D image coordinates (i.e., pixel locations), m and n, of the output file.

During this procedure, the pixels on the water bodies were assigned an invalid value using the water masks generated in the previous step. When two different flight paths occupied one pixel because of the side overlap, the most frequent index number of the existing points was assigned as the index number. Contrarily, when there was a pixel with no assigned point_source_id, the nearest neighboring pixel’s ID was assigned as the index number.

3.2. Estimating Erroneous Displacements of Target (2011–2013) DTMs

3.2.1. Gathering Height Values

Indiana’s lidar data were collected with numerous flight paths during both periods, and a single flight path covered large areas across multiple separate tiles. Therefore, the height values according to each flight path pair had to be gathered first from the multiple different tiles and then compared. Figure 3 depicts the gathering process used in this paper. Once a tiled input dataset (DTMs and flight path boundaries of two periods) was imported, the flight path pairs within the tile were defined using the flight path boundary data. For example, as illustrated in Figure 3, if there were two paths in the target data (#8–9) and three paths in the reference data (#20–22) in a tile, the four flight pairs, of which flight paths were overlapping, could be defined within the tile. These pairs were indexed as Pair (Target id, Reference id), like Pairs (#8, #20), (#8, #21), (#9, #21), and (#9, #22). According to each indexed pair, the height values of the valid pixels (i.e., non-water pixels) in target DTM and a resampled reference DTM were individually extracted and reshaped as a two-row matrix like the right of Figure 3. During the gathering process, the reference DTM was resampled using bilinear interpolation because its spatial resolution was half that of the target DTM. All the tiled DTMs in each Indiana county were sequentially processed in the same manner, and the rearranged height values from the different tiles were accumulated according to their flight path pair.

3.2.2. Histogram-Based Comparison

The long and wide spurious stripes in the original DoD results, as shown in Figure 1, revealed constant offsets between the target and the reference DTMs caused by vertical positioning errors. The offsets between the bi-temporal DTM values were estimated using the gathered height values according to the flight path pairs indexed in the previous process. In this paper, the target and the reference values were compared through their histograms to estimate offsets rather than compared pixel-by-pixel considering inevitable noises in the DTMs [43,44,45]. A histogram is a simple and effective way to describe data and has been used for many applications, including comparing multiple data, like in change detection [46] or image registration [47,48].

Figure 4 illustrates the histogram-based comparison implemented in this paper. As the area of the single flight path pair is extensive, the area where the terrain has changed generally is not broad enough to affect a histogram at a macroscale. Thus, the two histograms from the bi-temporal height values within one flight path pair are expected to be similar but have a displacement in the elevation domain induced by a constant offset due to vertical errors, as depicted with a red solid line and a black solid line in Figure 4a. To estimate this offset, the target histogram was gradually shifted by adding increasing values from −1 m to 1 m with a fixed interval of about 0.015 m (0.05 ft) to their height values, like the red histogram’s movement from left to right in Figure 4a. The similarity between the shifted target histogram and the reference histogram was calculated at every movement. When the target histograms shifted as the same offset amount, the two histograms would almost correspond, and their similarity would become the maximum value. Hence, the offset between the bi-temporal histograms could finally be determined as the shift with the maximum similarity (Figure 4b).

Numerous similarity measures exist between two histograms [49], each with advantages or disadvantages. Choosing the most appropriate measures for the compared histograms’ properties is thus essential; however, the properties of the histograms (such as their shapes) from the numerous flight path pairs over a large area are diverse depending on land covers or terrain slopes. Accordingly, we utilized five representative similarity measures between two histograms and aggregated the offsets individually determined with the five measures. The selected five measures were correlation d₁, intersection d₂ [50], Bhattacharyya distance d₃ [51], chi-squared distance d₄, and Kolmogorov–Smirnov (KS)-test d₅ as follows:

$$d_1=\frac{n\left[\sum _{i=1}^n{h}_1\left(i\right)h_2\left(i\right)\right]-\sum _{i=1}^n{h}_1\left(i\right)\sum _{i=1}^n{h}_2\left(i\right)}{\sqrt{\left[n\sum _{i=1}^n{h}_1^2\left(i\right)-{\left(\sum _{i=1}^n{h}_1\left(i\right)\right)}^2\right]\left[n\sum _{i=1}^n{h}_2^2\left(i\right)-{\left(\sum _{i=1}^n{h}_2\left(i\right)\right)}^2\right]}},$$

(1)

$$d_2=\frac{\sum _{i=0}^{n-1}\mathrm{m}\mathrm{i}\mathrm{n}\left[h_1\left(i\right), h_2\left(i\right)\right]}{\sum _{i=0}^{n-1}{h}_2\left(i\right)},$$

(2)

$${\left(d_3\right)}^2=1-{\sum }_{i=1}^n\sqrt{h_1\left(i\right) h_2\left(i\right)},$$

(3)

$${\left(d_4\right)}^2=\frac{1}{2}{\sum }_{i=1}^n\frac{{\left[h_1\left(i\right)-h_2\left(i\right)\right]}^2}{h_1\left(i\right)+h_2\left(i\right)},$$

(4)

$$d_5=\mathrm{m}\mathrm{a}\mathrm{x}\left|g_1\left(i\right)-g_2\left(i\right)\right|$$

(5)

where h₁ and h₂ indicate the two relative frequencies in bin i of the target and the reference elevation values, respectively. g₁ and g₂ are the cumulative relative frequencies of the two values. The number of bins for the histograms is n. In this paper, the two histograms were generated using the common bins, of which the entire range was between the minimum and the maximum of all the compared values, and the width was 0.15 m (0.5 ft), which was ten times the interval used for shifting the target values.

Correlation d₁ presents a linear relationship between two data and ranges from −1 to 1. Intersection d₂ calculates an overlapped rate between two histograms and ranges between 0 and 1. Since both methods quantify the similarity between two histograms, the two histograms correspond most when the calculated intersection value is closest to 1. On the contrary, the other three methods (Bhattacharyya distance d₃, chi-squared distance d₄, and KS test d₅) measure the dissimilarity between two histograms, and the low calculated values indicate that the two histograms are similar. Bhattacharyya distance d₃ ranging from 0 to 1 is related to the amount of overlap between two data, like intersection. Chi-square distance d₄ and KS test d₅ come from the statistics comparing the observed and the expected distributions. Chi-square distance ranges from 0 to infinity and considers the differences between small bins more important than large bins [52], as the differences are divided by the summed frequencies according to Equation (4). KS test distance is a standard nonparametric statistical test [53], ranging from 0 to 1. Accordingly, we first extracted five individual shifts with the maximum values for correlation and intersection then with the minimum values for the other measures as the offset candidates of each flight path pair. Afterward, the final offset for each pair was determined as the majority (i.e., mode) of the extracted five shifts.

3.3. Generating Improved Differencing Results

The estimated offsets for the flight path pairs were regarded as induced by the vertical errors of the target (2011–2013) flight paths and adjusted in this paper. All target DTMs with their corresponding flight path boundary data were subsequently loaded, and the adjusted target DTMs were generated by subtracting the estimated offsets from the loaded DTMs according to each flight path pair boundary.

DoD is a straightforward method compared to other lidar differencing methods once rasterized elevation products are generated from the lidar [54]. Since Indiana’s target and reference DTMs were divided into common tiles, DoD was simply implemented by subtracting the heights of the adjusted target DTMs, $z_{tar}^{adj}$, from the corresponding reference DTMs, $z_{ref}^{res}$, which were resampled as conducted in the previous comparison step:

$$z_{dif}=z_{ref}^{res}-z_{tar}^{adj}$$

(6)

The elevation difference $z_{dif}$ presents the terrain changes between the two periods with reduced erroneous displacements according to vertical positioning errors. The final differencing results were generated by merging all the tiled DoD results within each Indiana county.

4. Experimental Results

4.1. Study Sites

As Indiana’s statewide lidar products were collected and provided by counties, the proposed method was implemented on a county scale, and we assessed its feasibility and effectiveness based on the generated DoD results. A total of 8 out of 92 Indiana counties were selected as the experimental sites, representing different levels and appearances of the errors in the original DoD results (which were generated using the as-provided DTMs), as shown in Figure 1: Boone, Brown, Carroll, Decatur, Kosciusko, Monroe, Starke, and Tippecanoe. The lidar products of the eight counties were downloaded through the Integrated Digital Forestry Initiative at Purdue website [36]. The constant abnormal changes in the original DoD results clearly appear along the 2011–2013 flight paths over the selected counties. Most errors seem within the expected range of ±0.3 m based on the RMSE_z values announced by the IGIC [7]. Regardless of the severity of the errors, the existence of inherent errors is manifested over the counties, which disrupts a macroscale topographic change analysis. The selected counties’ detailed specifications and locations are presented in

Table 1. Specification of selected counties.

County	The Number of Tiles	Area (km2)	The Number of Flight Path Pairs	Acquisition Months
				2011–2013	2016–2020
Boone	520	808.1	1057	2011/03, 2011/04, 2011/09	2018/03, 2018/04
Brown	414	1147.1	183	2011/03	2017/03, 2017/04, 2018/03
Carroll	450	964.0	799	2011/03, 2011/04	2018/03, 2018/04
Decatur	469	965.0	956	2012/03, 2012/04, 2012/12	2017/03, 2017/04, 2018/03
Kosciusko	651	1376.3	342	2011/03, 2011/04, 2012/03	2017/03, 2017/04
Monroe	525	1021.8	2241	2011/03, 2011/09	2017/04, 2018/03, 2018/04, 2018/12, 2019/03
Starke	399	800.6	316	2011/03, 2011/04	2018/03, 2018/04
Tippecanoe	648	1294.5	284	2013/03, 2013/04	2018/03, 2018/04
Total	4076	8377.4	6178	-	-

and Figure 5, respectively. The proposed method processed 4076 tiled datasets (bi-temporal DTMs and flight boundary data with water masks), and the height values from 6178 flight path pairs were gathered and compared throughout the 8377 km² study sites, which are sufficient for evaluating the proposed method’s scalability (i.e., feasibility and effectiveness) for a large-scale 3DEP dataset.

4.2. Effectiveness of Proposed Method

This paper developed a new method to reduce the artificial displacements along the flight paths observed throughout the differencing results. However, the appropriate reference data, which enabled quantifying the individual flight path’s accuracy over Indiana, was absent. As an alternative, the proposed method’s results, (1) the displacements estimated with the histogram-based comparison and (2) the adjusted DoD results, were compared with the original DoD results in Figure 1.

4.2.1. Estimated Displacements by Vertical Positioning Errors

The original DoD results in Figure 1 revealed the existence of the inherent vertical errors in the Indiana datasets, but the severity of the errors could not be quantitatively measured. In this respect, we investigated whether the estimated displacements from the proposed histogram-based comparison could quantitively represent the errors displayed in Figure 1.

Table 2. The estimated displacement ranges based on the 5 and 95 percentiles. The minimum and maximum values are bold (unit: meter).

County	Boone	Brown	Carroll	Decatur
5th	−0.131	−0.345	−0.101	−0.131
95th	0.006	0.052	0.143	0.021
County	Kosciusko	Monroe	Starke	Tippecanoe
5th	−0.070	−0.162	−0.040	−0.147
95th	0.143	0.006	0.147	0.021

summarizes the range of each county’s estimated displacements based on its 5 and 95 percentile values. The county-level estimated displacements range between −0.345 m and 0.147 m. Specifically, 39.0% of the eight counties’ total flight path pairs were determined to include minimal displacements with absolute values of less than 0.025 m according to Figure 6a, which is the relative frequency histogram of the entire flight path pairs’ estimated displacements. The additional 57.8% of pairs contained absolute vertical displacements less than 0.15 m. These estimated values agreed with the expected range based on the RMSE_z announced by the IGIC [35] as well as the visually observed ranges based on the original DoD results in Figure 1.

The estimated errors of each county are also demonstrated to be consistent with the problematic elevation changes found in the original DoD results (Figure 1). The relative frequency histograms of each county’s estimated displacements are presented in Figure 6b,c. For instance, according to Boone County’s histogram in Figure 6b, over 70% of the displacements were estimated between −0.10 and −0.05 m. These estimations were considered reasonable compared to Boone County’s original DoD results, doubtfully indicating that the terrain commonly increased throughout the county. Similarly, the displacements of the other counties, where the elevations appeared generally and evenly increasing or decreasing in their DoD results, seemed as adequately estimated by the proposed method according to their histograms.

One noteworthy county is Brown County, which had the most severe errors in the original differencing results among the eight counties, as demonstrated in Figure 1. The more significant number of Brown County’s flight path pairs had displacements around −0.3 m (Figure 6b), which appeared as the distinct blue regions in Figure 1, compared to those of other counties. Furthermore, its histogram is demonstrated to consist of two distinguished distributions concentrated around −0.3 m and 0 m, which could also be visually identified in the original DoD results. In contrast, Tippecanoe County was presented to have the least errors according to the DoD results (Figure 1) and the estimated displacements by the proposed method. Figure 6c exhibits that a large number (72.9%) of Tippecanoe County’s flight path pairs have minimal displacements with absolute values of less than 0.025 m. Such consistency in the estimated displacements’ ranges (Table 1) and distributions (Figure 6) with the visual inspection results of the original DoD results (Figure 1) demonstrates that the proposed histogram-based comparison could quantitatively measure the vertical-error-induced displacements in the DoD results. Consequentially, it can be expected that the problematic striped errors will be reduced significantly in the differencing results presented in the following section by adjusting these properly estimated erroneous displacements.

4.2.2. Improved County-Level DoD Results

• Comparison on a county-level scale

The target (2011–2013) DTMs were adjusted using the estimated displacements, and the adjusted DoD results were finally generated. Figure 7 presents the thumbnail images of the adjusted county-level DoD results. The results demonstrate a higher consistency throughout the counties with fewer abnormal striped changes than the original DoD results in Figure 1. Even though the vertical error within each flight path pair was modeled as one constant value (i.e., offset) in this paper, the DoD results in Figure 7 reveal that the proposed method is effective overall in the eight counties regardless of the errors’ severities or patterns.

To quantitatively compare the original and the adjusted county-level DoD results, we calculated and compared the median values of the individual 2011–2013 flight paths for both DoD results. A common 2011–2013 flight path area is a stripe of several tens of kilometers by more than 1 km, including up to twenty million pixels of the DoD results. Even though the terrain heights changed within several regions during the data collection periods, the changes likely randomly occurred over small areas, and thus the terrain height could be regarded as stable overall in this massive flight path area. As a result, the calculated medians of the DoD values within each 2011–2013 flight path should be near zero unless there are systematic vertical errors within the paths, and we can measure the level of abnormal elevation changes in the DoD results by analyzing the median values.

The calculated medians within each county were analyzed using their absolute mean and standard deviations, as summarized in

Table 3. Statistical comparison of the original and improved DoD medians of individual 2011–2013 flight paths: the absolute means (Mean), the standard deviations (STD), and the improvement ratios (R). Original and adjusted mean of the original and the adjusted DoD results, respectively.

County		Boone		Brown		Carroll		Decatur
DoD		Original	Adjusted	Original	Adjusted	Original	Adjusted	Original	Adjusted
Mean	m	0.0629	0.0025	0.1701	0.0200	0.0546	0.0023	0.0562	0.0026
R	%	-	96.0	-	88.3	-	95.9	-	95.5
STD	m	0.0336	0.0030	0.1462	0.0240	0.0680	0.0031	0.0424	0.0037
County		Kosciusko		Monroe		Starke		Tippecanoe
DoD		Original	Adjusted	Original	Adjusted	Original	Adjusted	Original	Adjusted
Mean	m	0.0441	0.0195	0.0740	0.0210	0.0624	0.0084	0.0175	0.0065
R	%	-	55.7	-	71.6	-	86.5	-	62.8
STD	m	0.0618	0.0223	0.0545	0.0343	0.0499	0.0117	0.0170	0.0108

. It is noticed that the absolute means and the level of the abnormal elevation changes are related by comparing the absolute means in Table 3 and the county-level DoD results in Figure 1 and Figure 7; as the absolute mean value increased, the more abnormal changes were observed in the DoD results. For example, Carroll County, with the lowest absolute mean (0.0023 m) among the values in Table 3, appeared as the adjusted DoD result with the least abnormal elevation changes in Figure 7. Boone and Decatur Counties, which have similar absolute mean values to that of Carroll County, also present comparably fewer abnormal elevation changes in their adjusted DoD results. In contrast, the counties with high levels of systematic errors in their adjusted DoD results (Figure 7), such as Brown and Monroe Counties, show about ten times higher absolute mean values of 0.0240 m and 0.0210 m (Table 3), respectively, than that of Carroll County. As a result, based on this relationship, we can quantitatively assess the proposed method’s effectiveness by calculating a relative decrease in the absolute means between the original and adjusted DoD results. Equation (7) shows an improvement ratio R used to measure the relative decrease in this paper.

$$\mathrm{R}=\raisebox{1ex}{$(m_1-m_2)$}\!\left/ \!\raisebox{-1ex}{$m_1$}\right.\times 100$$

(7)

where m₁ and m₂ are the absolute means of the original and the adjusted DoD medians, respectively.

The R values in Table 3 illustrate the effectiveness of the proposed method in this paper to reduce the abnormal elevation changes in the original DoD results. The high R values of more than 50% in Table 3 illustrate that the medians of the DoD values within the individual flight paths became close to zero as the proposed method mitigated the possible abnormal elevation changes. In addition, the decrease in the standard deviation of the median values in the adjusted DoD results for each county (Table 3) can be interpreted as the consistency between the adjacent flight paths increasing, as visually identified in Figure 7. Meanwhile, the different R values of the eight counties imply that the proposed method was inevitably affected by inherent problems in the large datasets, such as irregularly overlapped areas of the 2011–2013 and the 2016–2020 flight paths, errors in the water masks induced by the temporal differences between the NHDPlus HR datasets and bi-temporal Indiana datasets, and errors in the point_source_id of the provided LAS files. Nevertheless, the decreases in their absolute means and standard deviations in Table 3 quantitatively present the proposed method’s effects in reducing the errors along the flight paths observed in the original DoD results, as shown in Figure 7.

• Detailed comparison with examples

For further investigation of the improvements due to the proposed method, the original and the adjusted DoD results were compared based on their profiles (or cross-sections). The profiles of the representative three counties’ DoD results were extracted, as shown in Figure 8. Brown, Carroll, and Tippecanoe Counties were selected as the representatives for the large errors, the moderate errors, and the small errors, respectively, considering the diversity of the error levels and patterns. The locations of the profiles were randomly selected, and their directions were determined perpendicular to the 2011–2013 flight paths. Due to the inherent random noises of the DTMs, the profile values were extracted by averaging 10 pixels by 10 pixels (i.e., 15.24 m or 50 ft).

Figure 8. Profiles (cross-section) of the three representative counties’ DoD results: (a) Brown County in the west-to-east (WE) direction; (b) Carroll County in the north-to-south (NS) direction; and (c) Tippecanoe County in the WE direction. The trace for each profile is presented as red lines in Figure 1 and Figure 7. The profiles were extracted based on the average values of every 10 pixels by 10 pixels (i.e., 15.24 m or 50 ft) area. The black and red lines indicate the values of the original and the adjusted DoD results, respectively. The blue dotted lines indicate the boundaries between the adjacent 2011–2013 flight paths, and the grey lines present theoretically ideal elevation changes at the macroscale (i.e., where the DoD value is 0). Some conspicuous values were actual elevation changes or the water bodies, which are not removed with NHDPlus HR.

Figure 8. Profiles (cross-section) of the three representative counties’ DoD results: (a) Brown County in the west-to-east (WE) direction; (b) Carroll County in the north-to-south (NS) direction; and (c) Tippecanoe County in the WE direction. The trace for each profile is presented as red lines in Figure 1 and Figure 7. The profiles were extracted based on the average values of every 10 pixels by 10 pixels (i.e., 15.24 m or 50 ft) area. The black and red lines indicate the values of the original and the adjusted DoD results, respectively. The blue dotted lines indicate the boundaries between the adjacent 2011–2013 flight paths, and the grey lines present theoretically ideal elevation changes at the macroscale (i.e., where the DoD value is 0). Some conspicuous values were actual elevation changes or the water bodies, which are not removed with NHDPlus HR.

The profiles again revealed that the proposed method reduced the abnormal elevation changes due to the inherent vertical errors along the flight paths. The original DoD results’ profiles (the black lines in Figure 8) illustrate the two typical problematic elevation changes visually identified in Figure 1. One is the consistently positive or negative average DoD values throughout the multiple flight paths presented over several kilometers across the flight paths in Brown and Carroll Counties’ profiles (Figure 8a,b). The most obvious identified erroneous changes are observed in the western and eastern regions of Brown County, and their widths are about 8 km. The other problematic change is a discontinuity between adjacent flight paths, such as sudden dramatic rises or falls between the adjacent flight paths, as seen in Figure 8a. Such problems possibly result in misinterpretation of the actual county-level topography changes.

These problems were mitigated in the adjusted DoD results, as demonstrated with the red profiles in Figure 8. As the vertical erroneous displacements decreased with the proposed method, the profiles along several tens of kilometers were aligned to 0 m, which were considered reasonable elevation changes on a macroscale. In addition, the profiles became more consistent along the WE or NS directions across the flight paths. Simultaneously, it was found that the variations of the elevation changes within the individual flight paths were unaffected by the proposed method. This indicates that meaningful elevation changes can still be easily observed in the adjusted DoD results.

The findings from the profiles are also observed in the tile-size DoD results (1524 m by 1524 m), as shown in the examples in Figure 9. Figure 9a shows the unrealistic elevation changes, as the entire area had evenly risen over 0.3 m between 2011–2013 and 2016–2020. Figure 9c illustrates a distinct artificial line in the middle of the tile induced by the discontinuity between the adjacent flight paths. Through the examples, it is predictable that such problematic changes in the original DoD results will bring an inaccurate terrain change analysis. The proposed method alleviated these abnormal changes, and thus more reasonable DoD results were created, as shown in Figure 9b,d. Furthermore, the detailed elevation changes in the overall areas are maintained in the adjusted DoD results. Therefore, the meaningful elevation changes in Figure 9d are expected to be analyzed more accurately as having less discrepancy between the adjacent flight paths.

Figure 9. Examples of the tile-size DoD results. (a,c): the original tiled DoD results. (b,d): the adjusted tiled DoD results. (a,b) are mountainous areas in Brown County, and (c,d) are near a quarry in Monroe County. The exact locations are marked as red boxes in Figure 1 and Figure 7.

Figure 9. Examples of the tile-size DoD results. (a,c): the original tiled DoD results. (b,d): the adjusted tiled DoD results. (a,b) are mountainous areas in Brown County, and (c,d) are near a quarry in Monroe County. The exact locations are marked as red boxes in Figure 1 and Figure 7.

To complement the missing appropriate reference data for assessment, we comprehensively compared the original and our adjusted DoD results using the flight path median values, the county-level profiles, and the tile-level examples. Through three different assessments, it was consistently demonstrated that the proposed method could considerably reduce the apparent problems in the original lidar differencing results, and the adjusted differencing results generated in this paper could provide more reasonable information about the elevation changes in eight county-level study sites. Given that its effectiveness was validated over large 8377 km² study sites composed of 6178 flight path pairs (Table 1), the proposed method is expected to be scalable as well as universally applicable for various levels and patterns of errors over different land covers. Furthermore, the proposed method can be implemented automatically once additional data regarding the water bodies is prepared, which is often available to the public [55,56], including the NHDPlus HR. Although few parameters were required for histogram-based comparison, the advanced quality differencing results were generated without changing them throughout the experiments in this study. Consequentially, we expect the proposed method to be a feasible solution to mitigate the apparent errors in multitemporal lidar datasets expectedly collected given the recent active government-driven lidar data acquisition without a time-intensive process.

4.3. Remaining Horizontal Positioning Errors

Although the proposed method was effective against critical vertical positioning errors, minor abnormal changes caused by the other error sources inevitably remained. One representative remaining error source is the horizontal positioning error. The horizontal error does not cause an apparent problem in plains. Still, it may become a problem in slopes because terrain aspects and slopes determine the sign and magnitude of the abnormal changes caused by horizontal errors. Figure 9b demonstrates the inherent horizontal errors over a mountainous area of Brown County. The opposite elevation differences in the opposite aspect, as shown in Figure 9b, are the visual evidence of the horizontal positioning errors in the dataset. In this respect, the counties with steep hills or mountains, such as Brown and Monroe Counties, have prominent remaining abnormal elevation changes in the adjusted DoD results (Figure 7) compared to the other relatively plain counties.

It would be ideal that the horizontal errors could also be reduced with simple implementation. However, it seems nearly impossible to adjust the horizontal errors using the provided DTM, like the proposed method in this paper, because the horizontal errors in Indiana datasets are predicted to be less than the provided DTM’s spatial resolution of 1.524 m (5 ft). Precise processing of low-level point cloud data, like a conventional method appropriate for a microscale, is inevitably required to achieve an accurate terrain change analysis over the mountain areas. Therefore, solving the horizontal errors is out of the scope of this study, which was developing a feasible method for a massive lidar dataset, and we would like to handle this problem in a future study.

5. Conclusions

Considering the previous methods’ practical limitations, this paper has presented a scalable method to mitigate the spurious elevation changes in large-scale multitemporal lidar differencing results. The proposed method focused on resolving critical vertical positioning errors by comparing bi-temporal DTMs readily provided as basic lidar products. By using the existing DTMs, implementing the proposed method became simple and feasible for a massive lidar dataset. The proposed method was applied to reduce the reported but unsolved abnormal changes in the DoD results of massive Indiana multitemporal 3DEP datasets [24]. The improved county-level DoD results from the proposed method were comprehensively compared with the original DoD results on macro- and microscales to complement the absent reference data. The proposed method’s effectiveness was demonstrated as (1) significantly reducing the vertical-error-induced abnormal elevation changes, (2) improving the consistency across the flight paths, and (3) simultaneously preserving details of the actual elevation changes within the individual flight paths. The remaining horizontal-induced displacements in the DoD results, which are out of the scope of this study, will be handled in a future study. The proposed method can be automatically implemented and universally applied with easily accessible water body information. Moreover, the improved Indiana lidar topographic differencing results produced in this paper are expected to facilitate more accurate and efficient downstream applications. Given the continuously increasing amount of high-quality lidar data, the proposed method can be a practical solution against the apparent errors in lidar differencing results that disturb the accurate and systematic analysis of macroscale elevation changes.