MODELING LAND DEFORMATIONS IN MOUNTAINS BY COMBINING TIME-SERIES L-BAND SAR IMAGES AND SPATIOTEMPORAL STATISTICAL MODELS

: Predicting land deformations using remote-sensing technology is extremely challenging because the simulation models that are involved have many parameters, including geological features, that must be calibrated against areas of interest, and no universal model is available. Instead, this paper proposes a method for modeling land deformations with time-series L-band synthetic aperture radar (SAR) images and spatiotemporal statistical models for prediction. First, time-series SAR analysis generates temporal land deformation. Then, water-related products are obtained from the Today’ Earth system of the Japan Aerospace Exploration Agency. The three-day consecutive maximum values are calculated and then interpolated spatially. Finally, a descriptive or dynamic spatiotemporal statistical model is applied. The low-rank Gaussian process model, a descriptive model, is found to be highly accurate for explaining land deformations, even for actual slid areas, as long as the areas of analysis are limited, e.g., 150 m × 50 m. It is concluded that the proposed method is effective for predicting land deformations spatially and temporally and so can be used to identify possible landslide locations.


INTRODUCTION
In Japan, be they large or small, several landslides occur near expressways every year, and some damage the neighboring expressway, thereby hampering traffic.For example, in August 2019, heavy rain triggered a large landslide with an area of approximately 170 m × 50 m near the Nagasaki Expressway operated by the West Nippon Expressway Company Limited (NEXCO-West), and it took almost a year to restore the expressway and operate the traffic service again (NEXCO-West, 2020).This area has a very gentle slope (~17°), and company staff mentioned that they did not imagine that a large landslide could happen on such a gentle slope.
For the past few decades, satellite synthetic aperture radar (SAR) has been a focus for measuring land deformation by using timeseries SAR images.PSInSAR (Ferretti, 2000) is one of the most popular algorithms for measuring land deformation.Among various types of SAR, L-band SAR is reported to be suitable for monitoring land surface conditions and therefore deformation in mountains because it penetrates into forests and can receive the backscatter from the land surface.One of the best known applications of L-band SAR to forests is the monitoring of deforestation in the Amazon rainforest by using the Phased Array type-L band SAR (PALSAR) onboard the Advanced Land Observing Satellite (ALOS) (Reiche et al., 2013;Pereira et al., 2016), while other frequencies (e.g., C, X, or P band) have potential applications to forest monitoring with deep-learning approaches (Bueso-Bello et al., 2021;Ortega et al., 2021).In a preliminary examination of the aforementioned landslide, we conducted time-series SAR analysis of the area by using ALOS2/PALSAR2, and the results revealed that some parts of the slid area had experienced considerable uplift and other parts had experienced considerable subsidence after heavy rain a year earlier in August 2018.This indicates that SAR images can be used to monitor a gradually ongoing landslide.
In the area of statistics, spatiotemporal analysis is very popular (Iaco et al., 2001;Castruccio and Genton, 2018), and for applying spatiotemporal models to environmental issues, a Bayesian hierarchical model is highly effective.Using the Markov-chain Monte Carlo (MCMC) method, it provides a data model based on observations and a process model that estimates the hidden condition of the phenomena (Guhaniyogi and Banerjee, 2019;Manago et al., 2019).Such a framework assists in flexibly expressing the relationship between an objective variable and explanatory variables.Our interest fits such a framework because land deformations are subject to many factors; for example, some of the internal factors are geological, terrain, and surface conditions, while the external factors can be meteorological data, mainly precipitation.One of the difficulties in predicting a landslide based on simulation is that we need many internal physical parameters related to the site.Moreover, calibrating those parameters is difficult, partly because their contribution to the landslide is site-specific.However, the spatiotemporal Bayesian hierarchical model has the potential to provide sitespecific spatiotemporal models for predicting land deformation.
Therefore, the aim of this research was to develop a spatiotemporal statistical model for predicting land deformations in mountains in order to take countermeasures against landslides.In the research framework, land-deformation data-the objective variable in the spatiotemporal statistical model-are derived by time-series SAR analysis, and the explanatory variables in the model are water-related physical quantities, e.g., precipitation, soil moisture, and surface runoff.As explained in Section 3.2, the hourly products are open to the public.The remainder of this paper is organized as follows.The location and data used for the research are described in Sections 2 and 3, respectively.Section 4 explains the existing spatiotemporal statistical models.Section 5 reports the experimental results including the time-series SAR analysis results and the preprocessing results of water-related products.The implications of these results and the validity of the algorithm are then discussed in Section 6.Finally, Section 7 concludes the paper.

STUDY AREA
As the study area, we selected the mountains and flat areas near the Takeo Junction (JCT) of the Nagasaki Expressway operated by NEXCO-West; the geolocation of Takeo JCT is 33.1479°N, 129.9956°E.Figures 1(a The slid area has a slope of ~17°, and the reported mechanism for the landslide is that because of the extreme precipitation, rainfall penetrated into the weathered tuff breccia, the groundwater level rose rapidly, and so the ground was weakened (NEXCO-West, 2019).

ALOS2/PALSAR2
In this research, we used ALOS2/PALSAR2 launched and operated by the Japan Aerospace Exploration Agency (JAXA).Table 1 gives the details of the 15 PALSAR2 images used for this research.The observation mode was UBS, and the observation direction was right.Period Ti is defined as that between observations i and i+1; with 15 PALSAR2 images, we defined the 14 periods of T1 to T14.

Today's Earth Japan
JAXA developed and operates the Today's Earth land-surface and river simulation system, which implements a land-surface model and a river-routing model.The input to the system are meteorological parameters derived from satellite observations, e.g., the Himawari meteorological satellites of the Japan Meteorological Agency, and the Global Precipitation Measurement Mission.The water and energy interactions between land surface and atmosphere are calculated, and then the surface runoff and baseflow are estimated independently (JAXA, 2022); see Takata et al. (2003), Yamazaki et al. (2011), Beven et al. (1984), and Bates et al. (2010) for details about the models used in the system.Today's Earth is divided into two sub-systems, i.e., Today's Earth Global (TE-Global) and Today's Earth Japan (TE-Japan), the spatial resolutions of which are 1/4 and 1/60 of a degree, respectively.Because our study area is within Japan, we used the products of TE-Japan, the temporal resolutions for which are hourly, daily, and monthly.
TE-Japan provides approximately 70 products related to heat and water including the meteorological input derived from satellite data.From these we selected and examined several water-related products, and herein we report the results for precipitation (GPRCT), total soil moisture (GLWtot), and surface runoff (SRUNOF).
We now introduce a covariance function that represents spatial autocorrelation and temporal correlation and can be expressed as a decay function with respect to distance and time difference (Iaco et al., 2001;Christopher et al., 2019): The covariance function must be positive-definite, which means that for any (, ) ∈  ×  , any  , ∈ ℝ,  = 1, ⋯ ,   ,  = 1, ⋯ ,   , and any positive integers   and   ,  , must satisfy When both   (ℎ  ) and   (ℎ  ) are positive-definite, so are their sum and product.In this research, we used the following two covariance functions: Eq. ( 4) is called the separable covariance with the product of   (ℎ  ) and   (ℎ  ) , and Eq. ( 5) is called the product-sum covariance, where k is a coefficient representing the crosscorrelation between   (ℎ  ) and   (ℎ  ) .The prediction of spatiotemporal kriging with these covariances are expressed as   =    +   ,   ~(0, c� , ,  , �). (6)

Basis-function representation
The R package mgcv can implement spatiotemporal prediction by using the low-rank Gaussian process (GP).The latter approximates full-rank GP variables η whose number is N with the regression coefficients whose number is L less than N, expressed by (; )  + (; ), In the low-rank GP, the coefficients in the basis function are assumed to satisfy a Gaussian distribution, i.e.,  ∼ (,    ′ +   ), Here, we focus on (; ) in Eq. ( 7) and approximate it as (; )  ,   ~(0,  2 ).

Time-series SAR analysis
We implemented time-series SAR analysis on the 15 PALSAR2 images described in Table 1, using PSInSAR (Ferretti, 2000) to estimate the deformation rate.the analysis, we used the publically available AW3D digital elevation model, whose pixel spacing is 30 m (EORC, 2022).We referred to Global Navigation Satellite System (GNSS) data obtained at the Nakahara, Saga, Yamauchi, Sagakashima, and Kawatana GNSS-based control stations of the Geospatial Information Authority of Japan (GSI) [Ref.9] around the study area, as shown in Figure 2(a).The Yamauchi control station is the nearest one to the study area.
Figures 2(b)-(d) show the original and 15-day moving mean deformation of Yamauchi, which was calculated as the difference against the geolocation on January 1, 2017. Figure 3(a) shows the geolocations of regions of interest (ROIs) 1, 2, and 3: ROI1 corresponds to the landslide area; ROI2 is a mountainous area, but no landslide occurred there; the main landcover in ROI3 is agricultural, i.e., paddy fields.
In PSInSAR analysis, we set the threshold for F-value in Eq. ( 17) [Anahara and Shimada, 2017] as -0.95 in order to select stable permanent scatterer (PS), where m and n are scene number, t is acquisition time,  is observed interferometric phase,  is wavelength, v is deformation rate,  is incidence angle, r is slant range, Δ is DEM error height,  ⊥ is perpendicular baseline.
It was found that the results of PSInSAR analysis by using the 15 PALSAR2 images did not have sufficient number of PSs.Therefore, instead of lowering the threshold for F-value, we divided the 15 PALSAR2 images into 3 datasets; 6 images from February 25, 2017 (No. 1 in Table 1

Water-related products of TE-Japan
In the case of a landslide, the 72-h precipitation before it can be one of the most effective variables for explaining land deformation (Yamamoto et al., 2019).Therefore, we calculated the maximum value of the consecutive three-day precipitation within a period of two satellite image observations.We used daily rather than hourly products because of the much smaller file volume.From the precipitation, we calculated the maximum values of total soil moisture and surface runoff, and we examined the periods of the three maximum values obtained.In most of the cases, the starting and ending days were either the same or differed by a day.We judged this processing to be reasonable.
The study area has an area of approximately 1.8 km × 1.8 km.
The spatial resolution of the TE-Japan products is 1/60 of a degree, which is equivalent to ~1.6 km at a latitude of 33°.We used the interpolation technique of ordinary kriging to generate data with a spatial resolution of 0.0001°, which is ~9.3 m and accords with the spatial resolution of the land-deformation estimate.

Descriptive modeling
We used mgcv in R for the descriptive modeling.The geolocations of the PSs were originally in degrees, and we converted them into plane orthogonal coordinates under the world geodetic system with meters as the unit of distance.The water-related products were normalized by subtracting the mean and dividing by the standard deviation.In R, we set the function as where the function te forms the product of X, Y, and Time, the "50" and "10" in (50,10) are the numbers of basis functions in space and time, respectively, the "2" and "1" in (2,1) are the dimensions in space and time, respectively, and the function s represents the spline smoothing function.As shown in Figure 4(a), we have 3 ROIs.Each ROI is divided into sub-ROIs (e.g., ROIs 1_1, 1_2, …, 3_4), each of which has an area of approximately 150 m × 50 m.The width of 50 m originated from our finding that the slid areas show different features in terms of deformation in 50-m intervals along the cross section against the expressway (Kusunose et al., 2022).The PSs in each sub-ROI are separated into one class for modeling and another for validation.For ROI1, Table 2 gives the statistics of the validation results, Figure 4(b) shows the geolocations of ROIs 1_1 to 1_4, and Figure 5 shows the temporal deformations predicted by three descriptive models, i.e., a low-rank GP model, a product-sum covariance model, and a separable variance model.

Low-rank Gauss
Product Table 2 Results of applying spatiotemporal models to ROI1.

Time-series SAR analysis
The estimated land deformations were found to be reasonable and supported by experts at the expressway companies, even though the results were not validated because of a lack of field survey data.It was found that parts of the landslide area have no PSs because of the low spatiotemporal coherence, and most of the PSs show uplift.Figure 3 (b) shows the uplift in the grid close to the expressway of the slid area after the heavy rain in 2018.This finding is consistent with in situ observations.We visited the slid areas, and the surrounding mountains have mixed forests and trees higher than 10 m.It was confirmed that L-band SAR has the potential to detect land deformations in mountains.We used the one-dimensional land-deformation results obtained from PALSAR2 images obtained on a descending part of the orbit.However, PALSAR2 images are available from both ascending and descending parts of the orbit, so we examined the results of time-series deformation analysis for each set of images; we found that the descending-orbit images are more suitable for the analysis because the results indicate the uplift and subsidence in the slid areas before the large landslide (Kusunose et al., 2022).We also examined the three-dimensional analysis by combining the two deformation vectors obtained from two sets of images and another deformation vector from the GNSS-based control stations operated by the GSI.However, the periods of both orbit observations do not overlap well, and so the obtained results did not show consistent deformation information.This issue of the gap between the periods of two-orbit observations should be investigated in the near future.

Modeling of land deformation
In this research, as the objective variable we used the relative land deformation from one observed SAR image to the next.We examined the cumulative land deformation as the objective variable, and the prediction by using the cumulative land deformation was much worse than that by using the relative land deformation.This may be because the water-related products do not reflect the cumulative land deformation even though the term to express the temporal correlation is included.
We now discuss the period for producing the water-related variables from the TE-Japan products.Some previous work has reported that 72-h integration of precipitation is suitable for explaining the degree of landslide damage (Yamamoto, 2019).We also examined the results obtained by setting two-day and four-day integration of water-related variables, and the statistics including R 2 and root-mean-square error (RMSE) show that the accuracies were lower than those in the case of three-day integration.For example, R 2 in two-day and four-day integration cases were approximately 0.01 to 0.02 smaller than those in three-day integration cases.This could be because two days (i.e., 48 h) are insufficient and four days (i.e., 96 h) are redundant for explaining the heavy rain and the accompanying surface runoff,.Therefore, we recommend using the 72-h period for preprocessing the water-related products.
We used other water-related products and terrain-related parameters (e.g., slope angle and elevation) as explanatory variables in the spatiotemporal statistical models, but none were found to be statistically significant.As for the three variables of precipitation, soil moisture, and surface runoff, their levels of statistical significance depend on the sub-ROI, and in some cases only one of them is regarded as statistically insignificant.These findings indicate that it may be enough to use three variables in the spatiotemporal models for predicting land deformation, and the optimal model depends on the area of interest.This dependency of the model on the area may reflect the geological features and underground structure.Because such information is quite difficult to measure or obtain, the spatiotemporal statistical models are very effective for predicting land deformation.

CONCLUSIONS
In this paper, we presented a method for predicting land deformations by combining time-series L-band SAR images and spatiotemporal statistical models.Temporal land deformations were estimated by applying PSInSAR, one of the most popular time-series SAR analysis algorithms, and the land deformations were used as observed data.Then, the water-related products of precipitation, soil moisture, and surface runoff obtained from TE-Japan of JAXA were processed into three-day consecutive maximum values within a period between two observations of ALOS/PALSAR2 images.After they were interpolated spatially, spatiotemporal statistical models were applied.As a descriptive model, the low-rank GP model showed high accuracies of explaining land deformation (even for actual slid areas) as long as the areas of analysis were limited, e.g., 150 m × 50 m.The proposed method is effective for modeling land deformation both spatially and temporally and so can be used to identify possible landslide locations.In future work, the prediction of land deformation should be investigated more to improve its predictive accuracy.
) and (b) show the study area.Heavy rain on August 27, 2019 triggered a landslide with an area of approximately 170 m × 60 m near Takeo JCT, and the slid soil damaged and uplifted the road surface of the expressway.It took approximately a year to restore the expressway completely, and the surface of the slid area was covered with concrete (NEXCO-West, 2020).The red rectangle in Figure 1(b) shows the landslide area, Figure 1(c) shows the damaged road just after the landslide, and Figure 1(d) shows the current surface of the slid area covered with concrete.
Study area: (a) Takeo Junction (JCT) in Japan; (b) Takeo JCT and landslide area shown by red rectangle; (c) expressway damaged by landslide on August 27, 2019 (NEXCO-West, 2020); (d) slid area covered by concrete.The photograph in (d) was taken by one of the authors on July 21, 2021.The images of (a) and background of (b) are from Google Earth.

Figure 2 .Figure 3 .
Figure 2. Global Navigation Satellite System (GNSS)-based control stations around study area: (a) Relative east-west movement and moving mean movement against geolocation of Yamauchi GNSS-based control station on January 1, 2017; (b) relative north-south movement and moving mean movement and (c) relative vertical movement and moving mean movement.The background image of (a) is from Google Earth.

Figure 4 .
Figure 4. Geolocations of regions of interest (ROIs).(a) 3 ROIs, and (b) Geolocation of ROI1_1 to 1_4.Yellow points indicate the geolocation of the examples shown in Figure 5.The background images of (a) and (b) are from Google Earth.

Table 1 .
Details of PALSAR2 images used around study area.The landslide near Takeo JCT occurred on August 27, 2019 in observations 10 and 11.