Extraction of Detailed 3D Coseismic Displacements in the 2024 Noto Peninsula Earthquake from Airborne LiDAR Data

1. Introduction

At 07:10 UTC on January 1, 2024, a moment magnitude (Mw) 7.5 earthquake occurred near the northeastern tip of the Noto Peninsula in central Japan [1,2,3]. The earthquake caused significant crustal deformation across the entire peninsula, with the seafloor emerging above sea level over a wide area along the northern coast [4,5]. Seismic waveforms and permanent displacements [6,7], together with tsunami waveforms [8,9] and tsunami run-up heights [10], provide fundamental information for constructing earthquake source-fault models. Regional crustal deformation and localized ground failures both contribute to three-dimensional (3D) displacement of the ground surface. Earthquake-induced crustal deformation also changes the positions of surveying reference points, such as triangulation and leveling points [11], and can shift road and residential boundary lines [12]. Because remeasuring these survey reference points takes time, about 1,300 real-time electronic reference stations (GEONET) based on satellite positioning (GNSS) had already been deployed throughout Japan [13]. But only six or seven stations were placed in the northern Noto region [14,15]. Temporary GNSS stations had also been installed by universities after the active swarm that began in November 2020 near the northeastern tip of the peninsula (Suzu City) [16,17], and additional stations were operated by SoftBank Corp. [17,18], but these datasets were not openly available.

Evaluating crustal deformation is a major challenge in assessing the impact of the 2024 Noto Peninsula earthquake. Wide-area estimates have been obtained using differential interferometric SAR (DInSAR) and pixel-offset (tracking) analyses of satellite SAR data [7,17,19]. However, coherence was low because of severe ground failures in the northern peninsula, and phase-unwrapped satellite line-of-sight (LOS) displacements could not be obtained from the DInSAR analyses. A three-dimensional topographic differencing (3DTD) analysis using airborne LiDAR data has also been reported [20], but further verification is still needed. Coseismic displacements at strong-motion observation points can also be calculated by integrating seismometer records [21,22], but the observation network on the Noto Peninsula was not dense enough to estimate the detailed distribution of crustal deformation.

In Wajima and Suzu cities in the northern Noto Peninsula, where crustal deformation was large and slope failures were frequent, surveying control points such as triangulation and leveling points remained unusable for an extended period after the earthquake [23]. However, after on-site conditions had been confirmed and GNSS observations had been completed, the triangulation points in the northern peninsula resumed their role as surveying control points at the end of May 2025 [23]. In this study, GNSS survey results obtained before and after the earthquake at more than 150 surveying control points in the northern Noto region, provided by the Geospatial Information Authority of Japan (GSI), are used as reference data for the crustal deformation caused by the earthquake.

This study aims to derive detailed, continuous 3D ground-surface displacements from pre- and post-earthquake airborne LiDAR data, thereby linking GNSS observations at survey control points with wide-area satellite-image analysis results. Using point-cloud data from pre-earthquake airborne LiDAR measurements conducted by Ishikawa Prefecture [24] and digital elevation models (DEMs) from post-earthquake airborne LiDAR data compiled by the Forestry Agency of Japan [25], we perform a 3D topographic differencing (3DTD) analysis based on the Iterative Closest Point (ICP) algorithm [26,27,28]. Because the pre- and post-earthquake DEMs are generated and compared on 0.5-m grids, this method is expected to provide horizontal accuracy of a fraction of the grid spacing [29,30,31]. However, because the results of the ICP method depend on the size and arrangement of the calculation window (tile), this study proposes a moving-average method using a window whose side length is an odd-number multiple of the 50-m square tile, which is considered sufficiently small for the ICP analysis.

In this study, we also attempt to separate ground-surface displacement into two components: crustal movement and local soil movement. The moving-average results are compared with GNSS observations at survey control points (triangulation points and GEONET stations) to verify the broad-area coverage but accurate 3D crustal deformation in the northern Noto Peninsula caused by the 2024 earthquake.

2. Materials and Methods

2.1. The 2024 Noto Peninsula Earthquake and Crustal Deformation

Figure 1 shows a satellite map of the study area, the northern part of the Noto Peninsula, Ishikawa Prefecture, central Japan. In the figure, the epicenter is indicated by a red symbol [32], and the submarine fault segments that moved during the 2024 Noto Peninsula earthquake are shown by yellow lines [33]. The large crustal deformation was caused by sequential slips on northeast-southwest-striking, southeast-dipping submarine faults along the northern coast of the peninsula, following the continuing swarm since late 2020 [34,35]. Strong acceleration records exceeding 1 g were observed at several K-NET and KiK-net stations operated by the National Research Institute for Earth Science and Disaster Resilience (NIED) [36] and at one station operated by the Japan Meteorological Agency (JMA) [37]. Due to the intense shaking, numerous landslides and soil movements were observed in mountainous areas [38], as also plotted in the figure. Coseismic displacements were recorded at GEONET CORS stations shown by pink triangles and at triangulation points shown in blue (Table A1), with their ranks indicated by Roman numerals. Of the six GEONET stations, real-time coordinates and elevations were released for four stations (Wajima, Suzu, Noto, and Anamizu) soon after the earthquake, showing that huge crustal movements had occurred on the peninsula. In contrast, two stations (Wajima-2a in Machino and Wajima-3 in Monzen) were deployed after the event.

GNSS surveys were conducted by the GSI at 150 triangulation points before the earthquake (April 2014) and after the earthquake (from April to December 2024). The post-event data were available online [39], and the pre-event data were provided by the GSI in June 2025. Note that the Japanese geodetic datum was updated from the “Japanese Geodetic Datum 2011 (JGD2011)” to the latest “Japanese Geodetic Datum 2024 (JGD2024),” effective April 1, 2025 [40]. This update aims to eliminate differences among survey results caused by the long-term accumulation of crustal deformation and to enable more rapid determination of orthometric height by GNSS observations using the newly implemented “Geoid 2024: Japan and its surroundings.” Therefore, to obtain the change in orthometric height at the 150 triangulation points, the post-event data were adjusted using the API [41]. In the study area, Geoid 2024 is 17–22 cm lower than Geoid 2011; therefore, these values should be added to convert JGD2024 data to JGD2011 data.

Figure 2 and Table A1 show the location and the horizontal and vertical coseismic displacements at the 4 GEONET CORS and 150 triangulation points in the study area. Vertical displacements were upward along the northern coast and almost zero along the southern coast. The maximum uplift (4.11 m) was recorded at the Isuzu point (TR35636051901), followed by 3.64 m at the Miyamaruyama point (TR35536755901), both in the Monzen district of Wajima City. Westward displacements dominated the horizontal field, with a maximum value of 2.22 m at the Kuroshima point (TR45536753801) in Monzen. Although the density of triangulation points is relatively high in the southern part of the peninsula, sparse or blank areas exist along the northern coastline. Thus, spatial interpolation of the observed displacements is not straightforward.

Coseismic displacements can also be estimated from satellite SAR data. Figure 3 shows the 2.5-dimensional results from pixel-offset (tracking) analyses of multi-path ALOS-2 intensity data [2,42]. Although the 2.5D analysis provides only quasi-UD and quasi-EW displacements, it captures the overall trend of the horizontal displacement field because the NS displacement is relatively small, as shown in Figure 2. The maximum quasi-UD displacement was about 4 m in the Monzen area, consistent with the GNSS observations. Uplift of about 2 m was estimated along the northeastern coastline of Suzu City, where triangulation points are sparse. The UD displacement decreases from the northern coastline toward the southern coastline. In contrast, the quasi-EW displacement is distributed more uniformly across the peninsula, from about 2 m along the northern coast to about 1 m along the southern coast. To perform ICP analysis of airborne LiDAR data, square regions excluding the sea must be selected. Therefore, ten square regions measuring 6.0 km (EW) × 4.5 km (NS) were selected, as shown in Figure 2, taking into account the distribution of GNSS survey points and coseismic displacements.

2.2. Airborne LiDAR Data Used in the Study

The areas covered by the pre-event and post-event airborne LiDAR data used in this study are shown in Figure 4 and Table 1. The pre-event data were acquired by the Department of Agriculture, Forestry and Fisheries of the Ishikawa Prefectural Government as part of a digital-transformation project for forestry information in fiscal years 2020 and 2022. The dataset covers the northern part of the peninsula and consists of high-density original point-cloud data (4.0 points/m²). Since February 2024, the dataset has been available as open data for non-commercial use to support recovery and reconstruction activities [24].

The post-event airborne LiDAR data were acquired through a joint survey project of the Geospatial Information Authority of Japan (GSI) and the Forestry Agency of Japan (FA). The dataset covers the entire peninsula at a density of 4.0 points/m². The original point-cloud data were processed by the FA into digital elevation (terrain) models (DEMs/DTMs) with a 0.5-m grid by removing trees and buildings, and the DEMs were released as open data in April 2025 [25]. By comparing these pre- and post-event airborne LiDAR datasets, we carried out a 3D topographic differencing (3DTD) analysis based on the Iterative Closest Point (ICP) algorithm [26,27,28].

To examine the validity of the ICP algorithm for extracting coseismic displacements associated with the 2024 Noto Peninsula earthquake, the Monzen district in the western part of Wajima City was selected because the largest crustal movement was observed there. Figure 5 shows the projected map of the pre-event LiDAR data for Monzen. Considering the computational cost of the ICP analysis, a sample area of 6.0 km (EW) × 4.5 km (NS) was selected; this area is composed of 6 × 6 (=36) standard airborne LiDAR data units (1.0 km × 0.75 km) in Japan. Using ENVI LiDAR software [43], a DEM with 0.5-m spacing in GeoTIFF format was created from the original point-cloud (LAS) data.

The post-event LiDAR data were released in much larger units, as shown in Figure 6. The processed DEM data with 0.5-m spacing cover a rectangle measuring 20 km (EW) × 15 km (NS). An area of the same size and location as the pre-event data in Monzen was extracted in GeoTIFF format using ENVI software [43], as also shown in the figure. Both the pre-event and post-event DEM data were then converted to point-cloud data (LAS format) using CloudCompare freeware [44] to conduct the ICP analysis.

2.3. The Iterative Closest Point (ICP) Method

In our previous study on extracting crustal movements from airborne LiDAR data for the 2016 Kumamoto earthquake [45], the pixel-offset method was used to obtain horizontal movements, and vertical differencing was then conducted to derive vertical displacements in Mashiki Town, Kumamoto Prefecture, Japan. In recent years, a 3D topographic differencing analysis tool based on the Iterative Closest Point (ICP) algorithm has become available online [46], and Python and MATLAB programs have also been released [47]. We tested the Python program using digital surface model (DSM) data for Mashiki Town and obtained results almost identical to those from our pixel-offset analysis. Therefore, in this study we use the ICP method available through OpenTopography [47,48] to extract crustal movements associated with the 2024 Noto Peninsula earthquake from airborne LiDAR data.

Figure 7 shows a schematic flow of the ICP method available on the OpenTopography (OT) website. First, core points representing the centers of square tiles (windows) are assigned for the pre- and post-event point clouds. Then, a 3D rigid-body transformation with translation and rotation is iteratively searched to bring the two sets of core points as close as possible in 3D space. For a pre-event (comparison) window width x, the post-event (reference) window width is set to x + buffer, allowing for a possible offset corresponding to the maximum correlation outside the common tiles. This method has been successfully applied to the 2008 Iwate-Miyagi (Mw 6.9) earthquake [27], the 2011 Fukushima-Hamadori (Mw 7.1) earthquake [27], and the 2016 Kumamoto earthquake [28].

2.4. Moving Average of the ICP Results

In the above-mentioned ICP method, the most important parameter is the common tile (window) size for the pre- and post-event LiDAR point-cloud data. Scott et al. [48] suggested that the proper window size is a function of LiDAR point-cloud density: 45 m for all original points and 32 m for ground points. We tested several window sizes for the Noto LiDAR datasets. In the northern part of the peninsula, where coseismic displacements are large, and numerous landslides and soil movements were observed (Figure 1). If a window includes part of a landslide-affected zone, the ICP method may still find a high-correlation 3D displacement, but that displacement may not represent the regional crustal movement. Therefore, a moving-average scheme is introduced to reduce the effects of ground-surface irregularity and obtain smooth coseismic displacements. Another reason for introducing the moving-average scheme is to reduce the effect of the grid layout of the window tiles. A large window is suitable for representing crustal movement over a wide area, but it increases the influence of the grid layout. Thus, the moving-average scheme is introduced to represent wide-area crustal movement while minimizing grid-layout effects.

The 3D coseismic ground-surface displacement, d^T = [d_x, d_y, d_z], is expressed as the sum of crustal movement (d_crustal) and local ground movement (d_local) as

$$\mathbf{d}={\mathbf{d}}_{crustal}+{\mathbf{d}}_{local}$$

(1)

The wavelength of the crustal movement is considered to be sufficiently long, on the order of a few hundred meters in the horizontal (x, y) plane. In contrast, local ground movement has a much smaller scale, on the order of several tens of meters in 2D space. The two-dimensional (2D) simple moving average (SMA) for a window centered at (n, m) with size (2k+1)(2l+1) is written as

$${\mathbf{d}}_{crustal}\left(n,m\right)={\displaystyle \frac{1}{\left(2k+1\right)\left(2l+1\right)}}\sum _{i=n-k}^{n+k}\sum _{j=m-l}^{m+l}\mathbf{d}(i,j)$$

(2)

where d_crustal(n, m) is the averaged (filtered) value at row n and column m, d(I, j) is the original input data value at row I and column j, and k and l are the radii of the moving-average window in the y and x directions, respectively. The total number of tiles in the neighborhood (window size) is (2k+1)(2l+1). A schematic of the 2D moving average for the case k = l = 2 is shown in Figure 8. If the moving window extends beyond the pre-event data range and the central unit tile of the window is still within the post-event data range, the averaging operation is applied using only the valid tiles within the data range.

The original input tile size was selected as 50 m × 50 m after several trials in this study. The local ground movement at row n and column m is then obtained as

$${\mathbf{d}}_{local}\left(n,m\right)=\mathbf{d}\left(n,m\right)-{\mathbf{d}}_{crustal}\left(n,m\right)$$

(3)

where d_local(n, m) is the residual coseismic displacement after removal of the crustal-movement (trend) component.

3. Results

3.1. 3D Crustal Movement in Monzen District, Wajima City

For the prepared pre- and post-event DEM data for the Monzen district in the western part of Wajima City, the ICP analysis was carried out using the Python code [47]. For the selected study area of 6.0 km (EW) × 4.5 km (NS), the analysis result for 50-m windows (tiles) is shown in Figure 9(a), where horizontal coseismic displacements are plotted as arrows at 100-m intervals and vertical displacements are shown by colors for 50-m square tiles. To examine the causes of variations in the coseismic displacements, the landslide polygons and surface scarp lines [38], both visually extracted from the same LiDAR datasets, are shown in Figure 9(b) together with five triangulation points in this study area. Considerable parts of the ground surface, especially forested areas, were clearly affected by slope failures (orange color). Disturbances in the coseismic displacements are co-located with the slope failures; for example, at the largest landslide location, the vertical displacement shows blue (negative) or gray (no-data) colors. The orientation of the horizontal displacement also becomes unstable at these soil-failure locations.

Using the results for the 50-m tiles (windows), the 2D simple moving-average calculation was carried out for various window sizes (k = l = 1, 2, …, 19) in Equation (2). Figure 10 and Figure 11 show the results for a 250-m square (5 × 5) moving window and a 550-m square (11 × 11) moving window, respectively, for the crustal movement (a) and local ground movement (b) components. As the window size increases, the crustal movement becomes smoother for both the horizontal (H) and vertical (V) components. Pixels (50-m tiles) without a solution diminish, and only one sinkhole-like settlement corresponding to the largest landslide remains in the V component for the 250-m moving-window result.

The crustal movement becomes even smoother for the 550-m moving window, and the sinkhole in the V component finally disappears. In contrast, the local soil-movement component looks almost the same as that for the 250-m moving-window case because it consists only of the local residual movements after removal of the mean trend (crustal-movement) component, and the trend is dominant in Monzen. Because the objective of this study is to extract crustal deformation from the pre- and post-event LiDAR data, further discussion of local ground movement is deferred to future work.

The moving-average results for various window sizes (k = l = 1, 2, …, 19) were evaluated using the root mean square error (RMSE) between the GNSS-observed displacement at point I (d_i^obs) and the displacement estimated by the ICP analysis (d_i^est), as

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(d_i^{est}-d_i^{obs}\right)}^2}$$

(4)

where n is the number of data points.

Figure 12 plots the relationship between moving-window size and RMSE at the five GNSS survey points in Monzen for each direction and for the sum of the three components. The RMSE is smallest for the 250-m window for the three-component sum because uplift is dominant, and the RMSE for the vertical component also reaches its minimum at the 250-m window. In contrast, the RMSE values for the two horizontal components do not show minimum values within this window range (150–950 m) because they are less affected by local conditions in the study area. The RMSE for the three-component sum shows the second-smallest value at the 550-m window; therefore, the 550-m window is used as another representative size together with the 250-m window hereafter.

The results for the original 50-m tiles and the moving averages for the 250-m and 550-m windows are compared with the observed displacements at the five survey points in Monzen, as shown in Figure 13 and Table A2 and Table A3. For both window sizes, the two horizontal displacement components agree well with the observed values (RMSE < 0.2 m, as also shown in Figure 12), whereas the 50-m tile results exhibit much greater scatter. For the vertical component, however, good agreement is seen at points MNZ2 and MNZ3 for both the 250-m and 550-m moving windows, whereas overestimation of more than 0.5 m is observed at the other three points. This discrepancy is discussed in the Discussion section.

3.2. 3D Crustal Movement in Anamizu Town

The ICP analysis of pre- and post-event LiDAR DEM data for 50-m tiles and the moving-average analysis based on those 50-m tile results were also tested in Anamizu Town, where the vertical coseismic displacements (uplift) are much smaller than those in Monzen but the horizontal displacements exceed 1 m westward. There are 11 survey control points (one GEONET station and 10 triangulation points) in the study area shown in Figure 2, the largest number of GNSS reference points among the 10 selected study areas.

Figure 14(a) shows the ICP analysis result for 50-m DEM tiles in the Anamizu study area (6.0 km × 4.5 km): vertical displacements are shown by the colors of the 50-m tiles, and horizontal displacements are shown by arrows at 100-m spacing. The vertical displacements are clearly much smaller than those in Monzen, with a maximum value of about 0.2 m. The horizontal displacements are more uniform in both amplitude and orientation than those in Monzen. The landslide polygons [38] are shown in Figure 14(b), together with the 11 GNSS survey control points in the study area. Some landslide-affected zones also exist in Anamizu along roadways and in forests, but they are much fewer than those in Monzen.

Figure 15 shows the moving-average results for the 250-m square (5 × 5) moving window and the 550-m square (11 × 11) moving window. Both the horizontal and vertical displacements become more stable as the window size increases, as expected from Figure 14.

The moving-window size was also examined in Anamizu, as shown in Figure 16, which plots the variation in RMSE for the 11 GNSS observation points with respect to window size. The RMSE values are much smaller for each direction and for the sum of the three components. The RMSE for the three-component sum reaches its minimum at the 750-m moving window, but it changes very little and remains below 0.1 m for window sizes larger than 350 m. This observation confirms the validity of the moving-average scheme for ICP analysis with a smaller tile size.

The results for the original 50-m tiles and the moving averages for the 250-m and 550-m windows are compared with the observed displacements at the 11 GNSS reference points in Anamizu, as shown in Figure 17. The 50-m tile results exhibit errors greater than 0.5 m in the EW and UD directions. The EW components still show some errors (AMZ4 and AMZ11) for the 250-m window, but these errors become much smaller for the 550-m window. For the 550-m moving window, the maximum error is 0.21 m for the NS component at AMZ4 and 0.27 m for the vertical component at AMZ7. These results further support the effectiveness of the moving-average scheme for estimating crustal movements from airborne LiDAR data.

3.3. 3D Crustal Movement in Machino District, Wajima City

As the third study area for extracting coseismic displacements from airborne LiDAR data, the Machino district in the eastern part of Wajima City was selected because crustal deformation was the second largest after Monzen and slope failures were frequent, as shown in Figure 1 and Figure 2.

The result of the ICP analysis for 50-m DEM tiles in the Machino study area (6.0 km × 4.5 km), together with an optical satellite image (Google Earth) showing ground-failure locations and three survey control points, is presented in Figure 18. The vertical displacement (uplift), shown by color, is smaller than that in Monzen but much larger than that in Anamizu. There are some gray pixels (50-m tiles) without solutions and some white/blue pixels indicating settlement. These irregular pixels correspond to landslide-related polygons in forested areas [38]. Two additional types of ground-surface irregularity data are also shown: scarp lines [38], indicated by pink lines, and surface cracks [49], indicated by yellow lines. The scarp lines are mainly observed in cropland, whereas the cracks occur along river and road embankments.

Figure 19 shows the moving-average results for the 250-m and 550-m windows. Owing to the effect of moving averaging, the irregular areas in the original 50-m tile results become smoother as the window size increases. However, in this study area, nonuniform crustal movements still remain even within this window range.

The results for the original 50-m tiles and the 250-m and 550-m moving windows are compared with the observed displacements at the three reference points in Machino, as shown in Figure 20 and Table A2 and Table A3. The ICP analysis for the original 50-m tiles could not provide meaningful results (less than 5.0 m for each horizontal direction) at MCN2 and MCN3, probably because of soil-surface irregularity around these points. However, the moving-average calculation yielded meaningful results for larger window sizes. The estimated vertical displacements almost perfectly match the GNSS observations after moving averaging. The deviations in the estimated horizontal displacements for the 250-m window become smaller for the 550-m window. In spite of the ground irregularity near the control points, the proposed calculation scheme is considered effective for reducing errors.

3.4. 3D Crustal Movement in Noroshi District, Suzu City

As the fourth study area of crustal deformation, the Noroshi district at the northeastern tip of the Noto Peninsula in Suzu City was selected because seafloor uplift was large there and the epicenter was located nearby. The earthquake swarm had continued in this area from late 2020 until the main shock occurred at 16:10:22.5 JST on January 1, 2024 [32].

The result of the ICP analysis for 50-m DEM tiles in the Noroshi study area, together with an optical satellite image showing ground-failure locations [38,49] and four GNSS survey points, is shown in Figure 21. To select a square area of 6.0 km × 4.5 km that includes as many GNSS points as possible, the current layout was adopted even though it includes some sea areas. These sea areas correspond to no-data (gray) or zero-vertical-displacement (white) pixels. A group of blue pixels in the southeastern port of the study area corresponds to a dam reservoir; therefore, the negative values there indicate a reduction in water level rather than ground settlement due to the earthquake.

Figure 22 shows the moving-average results for the 250-m and 550-m windows. The irregular areas in the original 50-m tile results become smoother as the window size increases. The effects of the sea and reservoir are reduced by the moving-average calculation, but they still remain even for the 550-m window case.

The results for the original 50-m tiles and the moving averages for the 250-m and 550-m windows are compared with the observed displacements at the four reference points in Noroshi, as shown in Figure 23 and Table A2 and Table A3. Scattered values for the NS displacement in the original 50-m tiles become closer to the reference data after moving averaging.

However, errors still remain in the NS component at NRS2 and NRS4, both of which are located close to sea cliffs. Slope failures were observed within 50 m from NRS2, which might generate errors in the ICP analysis and moving average calculation. Seafloor uplift occurred along the coast line in the study area. NRS4 is located 150-m south from the coast, and this position might affect the moving average calculation.

3.5. 3D Crustal Movement for the 10 Study Areas with 59 GNSS Observation Points

In the same manner as for the four study areas described above, 3D coseismic displacements were also calculated for six additional study areas shown in Figure 2. The results for the 10 study areas, including a total of 59 GNSS survey points provided by the GSI, are summarized in Figure 24 and Table A2 and Table A3. The corresponding RMSE values for the 59 survey points are listed in Table 2 for the original 50-m tiles and the 250-m and 550-m moving averages. The RMSE for the sum of the three components decreases from 0.55 m for the 50-m tiles to 0.19 m for the 250-m moving average and 0.16 m for the 550-m moving average. The coseismic displacements were much smaller and slope failures were less frequent in the study areas facing Toyama Bay and Nanao Bay.

The ICP and moving-average results for the remaining six study areas (Wajima City center, Noto Airport, Notocho Town center, Ogi district, Horyu district, Suzu City center) shown in Figure 2 are provided in Supplementary Materials.

4. Discussion

4.1. Data Type of Airborne LiDAR for the Extraction of Coseismic Displacements

Airborne LiDAR data are available in several forms, as shown in Figure 25. Raw observations are recorded as original point clouds in xyz coordinates, including all objects on the ground surface, such as trees, buildings, bare ground, and poles. Ground (terrain) point clouds are then obtained by cleaning and filtering nonground objects from the raw data. Digital surface models (DSMs) and digital terrain (elevation) models (DTMs/DEMs) are generated from original point clouds and ground point clouds, respectively. Both DSMs and DTMs (DEMs) can be produced either as triangular irregular networks (TINs) or in raster (grid) form [50].

In our case study of the 2024 Noto Peninsula earthquake, the pre-event data were available as original point clouds, whereas the post-event data were available as DTMs with a 0.5-m grid. Therefore, we first produced ground point clouds from the original point clouds, and then generated the DTMs using ENVI LiDAR software. To perform the ICP analysis with the Python program [47], the pre- and post-event DTMs were converted back to LAS-format point clouds after selecting the common study area. Because of this data limitation, we used only DTMs (ground point clouds). However, if the original point-cloud data for both epochs were available, ICP analysis could also be performed on DSMs (or the original point clouds).

An important question is which data type – DSMs or DTMs – is more suitable for obtaining coseismic displacements from a pair of pre- and post-event LiDAR datasets. In our experience with LiDAR differencing for the 2018 Hokkaido-Iburi earthquake [51], DTMs represented coseismic ground displacements better than DSMs because DSMs included tree growth and forestry activities that occurred during the time interval between the two datasets.

Fallen trees and collapsed buildings caused by earthquakes also introduced errors in change detection because they might be recognized as ground rather than as above-ground objects, as illustrated schematically in Figure 26. Landslides and associated fallen trees may cause overestimation of ground-surface elevation [52]. However, when there is little or no time lag between the two LiDAR acquisitions, DSMs can provide better coseismic displacement estimates, especially in urban areas [45], because buildings and urban infrastructure serve as good markers for cross-correlation between the DSMs.

4.2. Ground Failures and Errors in the ICP Analysis

As shown in Figure 24, the proposed method in this study – ICP analysis for 50-m tiles combined with moving-average windows – provided accurate results for the 59 reference points in the 10 study areas on the Noto Peninsula. Some results for the 50-m tiles, however, contain considerable errors, especially in the two horizontal components. Almost all of these errors diminish after the moving-average operation. For the vertical component, however, errors of 0.7–0.8 m still remain in the Monzen area, where the upward crustal movements were largest during the event. To investigate the cause of the vertical errors, optical images of the surrounding 1.5-km-square areas around the five GNSS reference points in Monzen are presented in Figure 27, together with visible ground failures (landslide polygons and scarp lines [38], and cracks [49]) identified from the LiDAR data and aerial photographs. At MNZ1, MNZ4, and MNZ5, where the vertical-component errors are large, landslide zones are recognized within 150 m of the triangulation points. As discussed before, landslides and fallen trees may cause overestimation of DTMs. The errors in the vertical component at these three sites can be explained by this mechanism.

4.3. Spatial Interpolation of Observed Coseismic Displacements

In this study, the pre- and post-earthquake high-density airborne LiDAR data were used to estimate 3D coseismic displacements associated with the 2024 Noto Peninsula earthquake. GNSS observations at 59 locations were used to validate the ICP analysis results. However, in the upper part of the peninsula, 154 survey reference points with pre- and post-event GNSS observations are available (Figure 2), although data from 150 triangulation points were acquired and released more than one year after the earthquake. Because the density of observation points is relatively high, spatial interpolation of the 154 GNSS displacement measurements for each (EW, NS, UD) component can be performed directly without using physical models.

Spatial interpolation was carried out using the Kriging tool of ArcGIS Pro [53]. Kriging is a well-known geostatistical method that estimates values at unsampled locations based on spatial autocorrelation among observed data points. It assumes that the distance and relative position between sample points reflect spatial dependence, which is quantified using a semivariogram model [54,55]. Ordinary Kriging was used in this study, in which the empirical semivariogram—constructed from the average squared differences between paired observations at given lag distances—is fitted with a spherical model. The spherical model is commonly used and represents a gradual increase in semivariance with distance until reaching a range beyond which spatial autocorrelation becomes negligible [56].

For interpolation, a variable search radius was adopted, allowing the number and distribution of neighboring points to vary depending on the local data density. The number of neighboring points used in each estimation was set to 12, following the default setting in ArcGIS Pro. The predicted value at each location is computed as a weighted sum of surrounding observations, where the weights are derived from the fitted semivariogram and the spatial configuration of the data points [57].

Ordinary Kriging was performed separately for the three-components over the upper peninsula using a 50 m × 50 m grid. After interpolation, non-land areas were masked out using boundary data for Ishikawa Prefecture obtained from the National Land Numerical Information download service [58]. The obtained results (interpolated displacement surfaces (mean estimators) and their variances) are presented in Figure 28.

For the UD component, the interpolated displacement surface appears similar to the 2.5D pixel offset result derived from ALOS-2 data (Figure 3). While Kriging reproduces the observed values at measuring points, the estimation uncertainty (variance) increases as the distance between an estimation point and nearby observations becomes larger. Along the northern coastline of the peninsula, there are no observation points between Wajima and Machino, or between Machino and Noroshi. Consequently, the variances (σ²) are large (maximum 0.222 m²; σ = 0.47 m) in these areas, and Kriging does not provide reliable estimates where observation points are sparse. In contrast, GNSS points are densely deployed along the southern coastline, facing Toyama Bay and Nanao Bay, resulting in smaller variances and more reliable interpolation.

For the EW component, the interpolated surface appears similar to that of the quasi-EW component derived from the 2.5D analysis. The variances are much smaller than those of the UD component because the interpolated surface exhibits less spatial variability than that of the UD component. The interpolated surface of the NS component has much smaller absolute values than that of the EW component. However, its directional and amplitude variations are large, resulting in variances greater than those of the EW component.

Spatial interpolation of GNSS observations is a simple and highly accurate method for estimating coseismic displacements. However, continuously operating reference stations (CORSs) are not yet sufficiently dense, even in a small but densely populated country such as Japan. Therefore, satellite SAR sensors remain essential for measuring crustal movements over wide areas associated with earthquakes and volcanic activities. Airborne LiDAR observations are expected to serve as a bridge between CORSs (high-accuracy point data) and satellite SAR (wide-area coverage).

5. Conclusions

In this study, three-dimensional (3D) ground-surface displacements were estimated using airborne LiDAR data acquired before and after the 2024 Noto Peninsula earthquake. Digital terrain (elevation) models (DTMs/DEMs) were generated from pre-earthquake point cloud data acquired by Ishikawa Prefecture and compared with post-earthquake DTMs developed by the Forestry Agency of Japan. Three-dimensional coseismic displacements were derived from the spatial correlation between pre- and post-event DTMs for 50 m × 50 m tiles using the Iterative Closest Point (ICP) algorithm.

The calculation results depend on tile size and are influenced by the amplitude and spatial distribution of coseismic ground movements. Therefore, moving-average windows of 250 m and 550 m were applied to the 50 m tiles to obtain continuous 3D displacement fields across the ground surface. A comparison between GNSS-measured displacements at 59 locations and the corresponding moving-average estimates for tiles containing triangulation points and GEONET CORS stations showed that the accuracy of the estimated displacements in all three components was within 0.2 m in terms of root mean square error (RMSE). The proposed method provides accurate and spatially continuous estimates of crustal deformation, effectively linking high-accuracy GNSS observations with wide-area satellite SAR-based measurements.