Consistent estimation of strain-rate fields from GNSS velocity data using basis function expansion with ABIC

Present day crustal displacement rates can be accurately observed at stations of global navigation satellite system (GNSS), and crustal deformation has been investigated by estimating strain-rate fields from discrete GNSS data. For this purpose, a modified least-square inversion method was proposed by Shen et al. (J Geophys Res 101:27957–27980, 1996). This method offers a simple formulation for simultaneously estimating smooth velocity and strain-rate fields from GNSS data, and it has contributed to clarify crustal deformation fields in many regions all over the world. However, we notice three theoretical points to be examined when we apply the method: mathematical inconsistency between estimated velocity and strain-rate fields, difficulty in objectively determining the optimal value of a hyperparameter that controls smoothness, and inappropriate estimation of uncertainty. In this study, we propose a method of basis function expansion with Akaike’s Bayesian information criterion (ABIC), which overcomes the above difficulties. Application of the two methods to GNSS data in Japan reveals that the inconsistency in the method of Shen et al. is generally insignificant, but could be clear in regions with sparser observation stations such as in islet areas. The method of basis function expansion with ABIC shows a significantly better performance than the method of Shen et al. in terms of the trade-off curve between the residual of fitting and the roughness of velocity field. The estimated strain-rate field with the basis function expansion clearly exhibits a low strain-rate zone in the forearc from the southern Tohoku district to central Japan. We also find that the Ou Backbone Range has several contractive spots around active volcanoes and that these locations well correspond to the subsidence areas detected by InSAR after the 2011 Tohoku-oki earthquake. Thus, the method of basis function expansion with ABIC would serve as an effective tool for estimating strain-rate fields from GNSS data.


Introduction
Earth's crust is deformed due to relative motions of tectonic plates. In particular, we observe significant crustal deformation in plate convergence zones. Accumulation of tectonic stress causes earthquakes, and large-scale crustal deformation leads to tectonic landforms. Precise knowledge on crustal deformation plays a fundamental role in understanding the tectonics and dynamics of subject regions. The present day crustal motions can be accurately measured with space geodetic techniques. In particular, networks of permanent GNSS stations are operated worldwide and crustal velocity at the stations can be obtained very precisely, which contributes to monitor ongoing crustal deformation (e.g., Kreemer et al. 2014;Nishimura et al. 2014;Murray et al. 2020). Velocity data generally depend on a reference station or a reference frame and include rigid rotations. Therefore, strain rates, obtained by differentiating velocity fields with Okazaki et al. Earth, Planets and Space (2021) 73:153 respect to space, are more suitable to interpret internal deformation of Earth's crust.
Estimation of a strain-rate field from spatially discrete velocity data is a longstanding issue to quantify crustal deformation, and many methods have been proposed before the operation of space geodetic measurements. Classical methods divide a region into a triangulated network to estimate a mean strain rate within each cell using triangulation survey data (Frank 1966;Prescott 1976) as well as trilateration survey and GNSS data (Feigl et al. 1990). In this approach, however, the estimated strains are discontinuous at cell boundaries and depend on the way of partition. Later, continuous interpolation methods were developed (e.g., Haines and Holt 1993;Shen et al. 1996). These methods impose a certain degree of smoothness on strain-rate fields to stabilize the estimation from discrete velocity data without knowledge on major faults or block motions. They are still in progress using advanced mathematical tools. For example, Tape et al. (2009) applied spherical wavelets to estimate velocity and strain-rate fields, which enables to separate crustal deformation into different length scales. Sandwell and Wessel (2016) developed an interpolation method of discrete 2-D vector data on the basis of Green's functions of an elastic body, which ensures the coupling of the two components.
In this direction, Shen et al. (1996) proposed a modified least-square inversion method to investigate regional crustal deformation in the seismically active Los Angeles basin. Assuming local uniformity of a strain-rate field, they simultaneously estimated velocity and strain-rate fields through a bilinear fitting with weighted contributions of data according to the distance to an estimation point. The distance dependence of weights is controlled by a hyperparameter called the distance decaying constant (DDC). This method is easy to understand and implement, and has been widely applied to clarify characteristics of crustal deformation fields. For example, applying the method to GNSS data in Japan, Sagiya et al. (2000) found a high strain-rate region named the Niigata-Kobe tectonic zone (NKTZ) that passes through central Japan. Furthermore, through the analysis of GNSS data in this region before and after the 2011 Tohoku-oki earthquake, Meneses-Gutierrez and Sagiya (2016) succeeded in separating elastic strain and inelastic strain; Fukahata et al. (2020) further succeeded in separating the inelastic strain into plastic and viscous strains. Nishimura et al. (2018) clarified strain partitioning between inland active faults and the interplate coupling in southwest Japan. Other than Japan, their method has also been used to reveal crustal deformation related to regional tectonics and seismic activities, for example, in Taiwan (Lin et al. 2010), mainland China (Wang and Shen 2020), Spain (Stich et al. 2006), Italy (Devoti et al. 2011), and Greece (Chousianitis et al. 2015). The strain-rate field obtained by their method has also been used to long-term earthquake prediction (Shen et al. 2007). Some improvements have still been proposed for the method (Shen et al. 2015), as described in the next section. In brief, the method of Shen et al. (1996) has made great contribution to geophysics.
However, we notice three theoretical points to be examined when we apply the method. First, the simultaneously estimated velocity and strain-rate fields are inconsistent mathematically: the estimated strain-rate field cannot be obtained by differentiating the estimated velocity field, as shown in Appendix. Secondly, there is no criterion to objectively determine the optimal value of DDC, though the value of DDC greatly affects the estimation result. Thirdly, as pointed out by Shen et al. (2015), it is difficult to properly estimate uncertainties of velocity and strain-rate fields.
As an alternative, in this study, we propose a method of basis function expansion to estimate a velocity field from spatially discrete geodetic data, in which a velocity field is expressed by a linear combination of basis functions. Coefficients of the basis functions, which determine the weights of respective basis functions, are obtained through a procedure of an inversion analysis. Once we obtain a velocity field, the associated strain-rate field can be analytically calculated by spatially differentiating the velocity field.
Techniques of basis function expansion have been broadly used in the waveform inversion to estimate seismic source processes (e.g., Olson and Aspel 1982;Hartzell and Heaton 1983;Ide et al. 1996;Yagi et al. 2004). In these studies, however, boxcar functions have commonly been used as basis functions. Therefore, we cannot differentiate them at cell boundaries. Yabuki and Matsu'ura (1992) introduced bicubic B-spline functions as basis functions to estimate coseismic slip distribution from geodetic data. Cubic B-splines are continuous until the second derivative. Their method has been widely used to estimate not only coseismic slip distribution (e.g., Fukahata and Wright 2008;Funning et al. 2014), but also interseismic slip distribution (e.g., Yoshioka et al. 1993;Sagiya 1999;Fukahata et al. 2004). Fukahata et al. (1996) applied the method to reconstruct temporal variation of vertical crustal motions from levelling data. In this study, we essentially follow the formulation of Yabuki and Matsu'ura (1992) and Fukahata et al. (1996) to estimate a velocity field from GNSS data. In this approach, as mentioned above, the strain-rate field is obtained by spatially differentiating the estimated velocity field.
In Yabuki and Matsu'ura (1992), which used a framework of Bayesian inversion, smoothness constraint was used as a prior constraint and the relative importance of it to observed data was objectively determined from observed data based on Akaike's Bayesian information criterion (ABIC; Akaike 1980). A statistically rigorous framework of Yabuki and Matsu'ura (1992) also enables us to evaluate estimation errors appropriately. Therefore, the above-mentioned three points on the method of Shen et al. (1996Shen et al. ( , 2015 are all overcome in the method of this study. Luo et al. (2016) also used ABIC to reconstruct a surface deformation field from InSAR and GNSS data, but they used boxcars as basis functions, and so, their method was not suitable to estimate strain-rate fields.
In the following, we first explain the method of Shen et al. (1996Shen et al. ( , 2015 and basis function expansion with ABIC to estimate a strain-rate field. Next, we apply these methods to GNSS data in Japan, and compare the results with discussion about the characteristics of these methods. Finally, we examine the strain-rate field in Japan based on the results obtained by the method of basis function expansion.

Methods to estimate a strain-rate field from spatially discrete geodetic data
Method of Shen et al. Shen et al. (1996) proposed a method to estimate a velocity field and a strain-rate field simultaneously from spatially discrete velocity data. In their method, assuming local uniformity of a strain-rate field, horizontal velocity components (u, v) , strain rates e xx , e xy , e xy = ∂ x u, ∂ x v + ∂ y u /2, ∂ y v and a rotation rate ω = ∂ x v − ∂ y u 2 at an arbitrary point x, y are related with observed velocity data (u i , v i ) at i th station of a coordinate x i , y i by the following relation: where (�x i , �y i ) = x i − x, y i − y are the relative position of the observation station to the estimation point. δu i and δv i are called observation errors in Sagiya et al. (2000), but it would be more appropriate to regard them as fitting errors. Note that, following the convention (e.g., Segall 2010), anticlockwise rotation is taken to be positive for the rotation rate in Eq. (1), although clockwise rotation is taken to be positive in the original formulation of Shen et al. (1996). Using the definitions of the strain rate and rotation rate, which are given above, Eq. (1) can be decoupled to Equation (2) corresponds to the Taylor expansion of the velocity field at the first order. Equation (2) indicates that, as long as the errors (δu i , δv i ) are mutually independent, the estimations of u, ∂ x u, ∂ y u and v, ∂ x v, ∂ y v are independent of each other. Therefore, we explain the method for one velocity component in the following.
We have velocity data v i with variance (observation errors) σ 2 i at (x i , y i ) of N stations ( i = 1, . . . , N ). We consider a multiple regression model at each estimation point P x, y : Model parameters to be estimated are velocity v and its spatial derivatives ∂ x v and ∂ y v at a position P x, y . A key point in the estimation is that the fitting errors δv i are weighted according to the distance between the observation point (x i , y i ) and the estimation point P x, y as where the notation N m, σ 2 represents the Gaussian distribution with mean m and variance σ 2 . Equation (4) means that data at remote points from P do not contribute to the estimation at P. D is the distance decaying constant (DDC), which is a hyperparameter in this analysis and controls the spatial range of significance. Since the weight of each station varies continuously, this method leads to a smooth estimation of velocity and strain-rate fields.
The observation Eq.
(3) is expressed in a matrix form as where "diag" represents a diagonal matrix. We attach the subscript P to note that these quantities depend on an estimation point P. Since e P follows the Gaussian distribution, the relation of Eq. (5) can be expressed in the form of probability distribution: which can be regarded as the likelihood function for the model parameter a P . By maximizing the likelihood for given data d and the hyperparameter D , we obtain the optimal values of a P with its uncertainties as Shen et al. (2015) investigated the model structure further: they applied a quadratic decay besides the Gaussian decay (Eq. 4), made compensation for uneven azimuthal distribution of data, and determined the value of D at each point by assigning the total weighting of data W , which is defined as They selected an optimal scheme by mapping the difference of estimated strain-rate fields.
It is a good advantage of the method of Shen et al. (1996Shen et al. ( , 2015 that both velocity and strain rate are simultaneously computed at any point and these fields vary continuously in space. In addition, the implementation is simple, and so, this method has been used by many researchers and contributed to the advancement of geophysical research (e.g., Sagiya et al. 2000). As mentioned in "Introduction", however, there are three points to be examined in this method. First of all, because the velocity v and its derivatives ∂ x v and ∂ y v are estimated at each point independently, the velocity and strain-rate fields are mutually inconsistent in general. That is to say, the obtained strain-rate field is different from the spatial derivative of the obtained velocity field except for special cases (see Appendix). Secondly, it is difficult to objectively determine the optimal value of the hyperparameter D . This is because the regression model, Eq. (7), is independently constructed at each estimation point and we do not treat all (unweighted) data at once. This prevents us from determining the hyperparameter using a single objective function. Thirdly, as pointed out by Shen et al. (2015), estimation of uncertainties is inappropriate. This is because the uncertainty (Eq. 8) depends only on the (7) variance σ 2 i of velocity data at each station, so it does not reflect actual fitting to velocity data v i .

Basis function expansion with ABIC
We first formulate for one velocity component, and then describe how to deal with two components. We express the velocity field as a linear combination of a set of fixed basis functions Φ j x, y M j=1 : where a is a model parameter vector to be determined from observation data. As mentioned in "Introduction", we use bicubic B-splines for basis functions Φ j . Once we obtain the velocity field v , we can analytically calculate the strain-rate field using its spatial derivatives (e.g., In our problem to reconstruct a velocity field from spatially discrete velocity data, Eq. (10) corresponds to the observation equation, which is written in a matrix form as with We assume the errors to be isotropic Gaussian e ∼ N 0, σ 2 I , where σ 2 is an unknown scale factor (hyperparameter) of the variance, and I is the N × N unit matrix. Then, we obtain the data distribution as Following Yabuki and Matsu'ura (1992), we construct a Bayesian model incorporating prior information that the velocity field should change smoothly in space. Although velocity discontinuity could exist at boundaries of rigid motions, gradual variation of velocity is usually detected by dense observation networks in deformation zones. The method of Shen et al. (1996Shen et al. ( , 2015 also imposes local uniformity on strain rates. The roughness of a velocity field can be measured by the following quantity: In order to realize a smooth velocity field, the roughness r should be small. In terms of model parameters, it is rewritten in a positive-semidefinite quadratic form as where with With the roughness defined in Eq. (15), we can introduce prior constraints on the smoothness of a velocity field in the form of the degenerate Gaussian distribution (Fukahata 2012): where P is the rank of R , and | P | is the product of the non-zero eigenvalues of R . ρ 2 is a hyperparameter that controls the strength of smoothness constraint.
We combine the two distributions, Eqs. (13) and (18), through Bayes' theorem to derive the posterior distribution of model parameter a for given observation data d: where c is a normalization constant independent of model parameter a and hyperparameters σ 2 and ρ 2 . For fixed values of the hyperparameters, the model parameter a and its covariance are obtained by maximizing Eq. (19) as The optimal values of hyperparameters can be objectively selected by minimizing Akaike's Bayesian information criterion (ABIC) introduced by Akaike (1980). ABIC is defined by where N h is the number of hyperparameters ( N h = 2 for the current modeling). By minimizing ABIC, the optimal value of σ 2 can be analytically solved as with Then, ABIC is expressed in terms of α 2 as The search for the optimal value of α 2 is performed numerically. Substituting the optimal value α 2 into Eqs. (22) and (20) gives optimal values of σ 2 and the model parameters.
It is now straightforward to treat two velocity components, where we assume that the degree of smoothness is common to both components. Components of data, model parameters and errors are distinguished by subscripts as while matrices H and R and hyperparameters σ 2 and ρ 2 are common in both components. The observation equation and roughness are written as The optimal model parameters and their covariance matrices are given by which is identical to Eq. (20). ABIC is calculated as with Here, note that the number of data and model parameters are 2N and 2M , respectively, in this case. The optimal value of σ 2 is given by The two velocity components are mostly independent, but interacted with each other through the common hyperparameters.

Data and model setting
We apply the above two methods to GNSS velocity data in Japan to estimate velocity and strain-rate fields. The number of the used stations is 1336 in the region ranging 128°-146° E and 30°-46° N ( Fig. 1; Additional file 1). Raw data of these stations are archived at the Geospatial Information Authority of Japan (GSI), the Japan Coast Guard, Kyoto University, the International GNSS service (IGS), and UNAVCO. We estimate daily coordinates of continuous GNSS stations using precise point positioning with ambiguity resolution (Bertiger et al. 2010a, b) implemented in the GIPSY Ver. 6.4 software (https:// gipsy-oasis. jpl. nasa. gov/). The coordinates are transformed into the IGS14 reference frame (http:// www. igs. org/ artic le/ igs14-refer ence-frame-trans ition), which is GNSS realization of the International Terrestrial Reference Frame (ITRF) 2014 (Altamimi et al. 2016) using the transformation parameters provided by the Jet Propulsion Laboratory. The daily coordinates from January 2006 to December 2009 are used in the following analysis. Coordinate offsets associated with 14 large earthquakes (Table 1), 2 dyke intrusion events in eastern Izu volcanoes (Ueno et al. 2012), and maintenance of equipment referring to catalogues provided by GSI and Japan Meteorological Agency (JMA) are removed. The velocity vectors are estimated with a conventional procedure similar to that in Sagiya et al. (2000): linear trend, sinusoidal annual variation are fitted by the least-square method for each component separately, and then the coefficient of the linear term is used as the velocity component.
In the method of Shen et al. (1996Shen et al. ( , 2015, or shortly Shen's method, the form of the weighting function must be chosen. We adopt the original formulation of Shen's method (Shen et al. 1996), because the hyperparameter D is more intuitive and the trade-off curve indicates a better performance than the method proposed by Shen et al. (2015), which uses W instead of D as shown in Discussion (see Fig. 11). An appropriate value of the hyperparameter D generally depends on the observation density and the length-scale of crustal deformation, and it is actually determined through comparison of estimated results. In consequence, different values of D have been adopted in different studies: e.g., 25 km in Shen et al. (1996) and 35 km in Sagiya et al. (2000). Therefore, the results of several values of hyperparameters, specifically D = 15, 25 and 35 km, are compared in this study.
In the basis function expansion with ABIC, or shortly the ABIC method or ABIC, we set the Cartesian coordinates with the origin at (137° E, 38° N), and take a rectangle of 1800 km (north-south) by 1600 km (east-west) as

Comparison of the results of the two methods
The velocity fields estimated by both methods are shown in Fig. 2, in which the results are presented for the places where three or more observation stations exist within the radius of 50 km. All results are similar in most regions. Looking into the details, however, Shen's method with D = 15 km yields an unstable eastward velocity field along the Pacific coast in northeast Japan. This suggests an overfitting to observation data, and D = 15 km is considered to be too short to perform a stable estimation. Figure 3 is an enlarged view of eastward velocity around the Izu Islands. Observation data generally have westward motions, but Shen's method commonly estimates eastward motions in the east offshore of the islets. This is because the method extrapolates the large gradient of velocity data observed in Miyakejima Island, where adjacent stations have significantly different velocities. On the other hand, ABIC estimates a reasonable velocity field even in the offshore region. In brief, both methods usually give similar velocity fields, but Shen's method can be unstable on the periphery of observation networks.
Estimation bias, i.e., mean residual of estimations to observations at the observation sites Iwate-Miyagi inland earthquake, which is also shown in Fig. 1 Okazaki et al. Earth, Planets and Space (2021) Table 2. The bias of ABIC is very tiny (the order of 10 −15 , which may be numerical errors). This is an advantage of treating all data simultaneously, which leads to the cancellation of the bias. On the other hand, in Shen's method, since velocities are estimated independently at each data point, biases are not cancelled out and accumulate randomly.
Estimated strain-rate fields, specifically the dilatation rate = e xx + e yy and maximum shear-strain rate � = e 2 xy + e xx − e yy 2 /4 , are presented in Fig. 4. Here, � > 0 (warm colors) represents expansive deformation, while � < 0 (cold colors) represents contractive deformation. Differences between the estimation methods are clearer in the strain-rate field than in the velocity field. The scale of spatial variation differs by the value of DDC in Shen's method: the result of D = 15 km shows shorter length-scale variations, while that of D = 35 km shows a much smoother field. It is a difficult point that there is no objective criterion to determine the value of D. ABIC estimates the strain-rate field similar to Shen's method with D = 25 km. Differences between the two methods can be better recognized around the focal region of the 2008 M7.2 Iwate-Miyagi inland earthquake (140.9° E, 39.0° N), the epicenter of which is marked by a green star in Fig. 1, and an enlarged view is shown in Fig. 5. A pair of weak positive and strong negative dilatation rates is observed in the result of ABIC. In Shen's method, the positive dilation is unclear for D = 25 km, and disappears for D = 35 km; the result of D = 15 km shows the pair of positive and negative dilatation rates, similar to ABIC, but there exist many such pairs in other areas and it is difficult to distinguish meaningful signals.
To compare the estimated results more quantitatively, the velocity and dilatation-rate fields along three cross sections (see Fig. 2) are plotted in Fig. 6. Figure 6a, b shows results along north-south and east-west sections, respectively. Velocity fields are almost identical among all methods and different D values except for the outside of the observation network. There exists an outlier at 33.5° N on the 133° E section (Fig. 6a). The velocity profile of Shen's method is not affected by it, while that of ABIC is slightly, but not significantly, dragged. Both methods can be said to be sufficiently robust to a single outlier in estimating a velocity field. Strain-rate fields are more sensitive to estimation methods. In particular, along the 35.3° N section that passes through Mt. Fuji around 138.7° E, the magnitude of the dilatation rate is significantly different among the methods (Fig. 6b). Figure 6c presents the velocity and dilatation-rate profiles that pass through the focal region of the 2008 Iwate-Miyagi inland earthquake. Observed eastward velocities generally decrease from west to east, which indicates contractive motion in the east-west direction. However, there are a few data with faster eastward velocities in 140.5°-141.0° E; the motion causes larger east-west contraction around the epicenter and small expansion in the outside (Ohzono et al. 2012). ABIC captures this   . 4 Comparison of the estimated strain-rate fields. Results of the ABIC method and Shen's method with different D values are presented. Positive dilatation means expansion tendency smoothly. In Shen's method, the eastward velocity monotonically decreases for D = 35 km and, is barely stagnant for D = 25 km; the velocity profile for D = 15 km shows a better fit to the data in the focal region, but the strain-rate field oscillates in a small scale. ABIC exhibits sharp crustal deformation in the focal region while keeping smooth variation in the outside.

Setting of basis functions in the ABIC method
In the analysis using ABIC, the interval L of basis functions has been typically fixed beforehand (e.g., Fukahata et al. 1996), but it can have significant influence on the fitting performance. Therefore, we examine its effect in the following. As L gets smaller, the fitting generally improves and saturates when it reaches fine enough to resolve the variation of data, whereas the computational cost increases dramatically. The computation time of the inverse matrix (the most expensive portion in the analysis) is proportional to the cube of the number of basis functions M , which is inversely proportional to L 2 ; if the interval decreases by half, the computational cost becomes heavier by approximately 2 2×3 = 64 times. Therefore, we need to take a balance between finer fitting and computational cost.
The dependence of M and the value of ABIC (Eq. 29) on L is shown in Fig. 7. Figure 7a illustrates that the  number of basis functions, and hence computational cost, increases rapidly as L decreases. As shown in Fig. 7b, ABIC tends to increase with L , because of the lack of resolution to represent the variation of data. Fluctuation of the plot implies that the results could be affected by the relative position between observation data and nodes of basis functions. This fluctuation also suggests that the interval is too coarse to fit the data. The values of ABIC stay almost constant for L ≤ 30 km, which implies sufficient resolution in this range.
The estimated velocity and strain-rate fields for different L are presented in Figs. 8 and 9. The velocity field is almost identical for L ≤ 30 km. Although the estimated strain-rate changes around an outlier (Fig. 9a) and a focal region of the 2008 Iwate-Miyagi inland earthquake (Fig. 9c), the difference is much smaller than that of Shen's method with different D (Fig. 6). In particular, the results of L = 15 and 20 km are very similar to each other. In brief, results are stable for the change of L in the ABIC method, when L is taken to be sufficiently small. In this study, we basically present the result of L = 20 km, because it achieves a sufficient resolution and has moderate computational cost.
We next investigate the effect of the form of basis functions at the boundary of the analysis region (Fig. 10). Yabuki and Matsu'ura (1992) and Fukahata et al. (1996) used basis functions that take zero at the edge of the analysis region (Fig. 10a), which forces the field to be zero at the boundary. This model assumption is reasonable in expressing slip distribution of earthquakes, but inappropriate in expressing a horizontal velocity field. If we use these basis functions, the velocity field gradually approaches zero in the region outside the observation network (Fig. 10c) to mitigate large roughness near the boundary (Fig. 10g). This leads to the reversal of the dilatation rate from a coastal area to the sea area (Fig. 10e). In Fig. 10e, contraction is dominant in the land area, which is consistent with compressional tectonics in Japan, while expansion is dominant in the sea area. We observe the reversal of the dilatation rate more clearly, where observation sites are closer to the boundary (e.g., north and east of Hokkaido, west of Kyushu).
To remove such artificial deformation outside the observation network, in this study, we use basis functions truncated at the boundary shown in Fig. 10b; this strategy was partly used by Fukahata and Wright (2008) to express coseismic slip distribution at the Earth's surface. This changes the number of model parameters and components of the prior constraint R (Eq. 16) near the boundary. The obtained result shows natural extrapolation of the velocity field to the sea area (Fig. 10d), and the reversal of dilatation rates disappears (Fig. 10f ). It is critical for reasonable modeling to set appropriate basis functions.

Discussion on the problems in Shen's method
We mentioned three theoretical points to be examined in the method of Shen et al. (1996Shen et al. ( , 2015: velocity and strain-rate fields do not satisfy the relationship of differentiation, the value of the hyperparameter D (distance decaying constant) cannot be objectively determined, and the estimated uncertainty is unreliable. We discuss how much impact these factors have on the estimation results.
First, we consider the optimization of hyperparameters. It is crucial for fitting problems to take a balance between resolution and certainty (robustness), which are in a reciprocal relationship (Backus and Gilbert 1970;Menke 2012). In Shen's method, the hyperparameter D plays a role of regularization; smaller values yield a highresolution but unstable (unrobust) solutions, while larger values yield a stable (robust) but low-resolution solutions. In the basis function expansion, the smoothness of solutions are represented by prior information and its significance is controlled by ρ 2 ; the optimal value of it is objectively determined from observed data using ABIC.
To quantify the resolution and certainty, we use RMS of fitting errors (residual) and roughness (Eq. 14) of the velocity field, respectively. The lower these indices, the better the model is. The relation between the residual and roughness of the estimated results is plotted in Fig. 11. For Shen's method, we also show the case where W (Eq. 9) is used as a hyperparameter instead of D. In Shen's method, the residual and roughness show a clear reciprocal relationship. The original method (Shen et al. 1996), which changes D, gives slightly better performance than the revised method (Shen et al. 2015), which changes W, in this criterion. Roughness increases rapidly below a certain value of the hyperparameters ( D ∼ 25 km and W ∼ 6 ), which implies overfitting to observation data. On the other hand, the results with large D or W strongly flatten the estimated velocity field, and has too small roughness, which leads to increase of the residual. The indices of Shen's methods with D = 25 km are closest to those of ABIC with L ≤ 30 km, which is in agreement with a visual comparison of Fig. 4. A difficult problem is how to determine the optimal point on the trade-off curve in the absence of other criteria. In ABIC, the results of L ≥ 50 km are comparative or worse than those of Shen's method, which suggests that basis functions are too coarse to fit the data variation. On the other hand, the indices (residual and roughness) of L ≤ 30 km, for which the value of ABIC is sufficiently small (Fig. 7b), are located in the lower left part compared to the tradeoff curves of Shen's method. This indicates the superiority of the ABIC method in this criterion. Although the curve of the ABIC method does not converge for small L , the dependence of the indices on L is much milder than that on D and W of Shen's method. It is a manageable property of ABIC that a smaller L generally gives a better result and that a sufficiently small L gives a similar result. This is because a smaller L ensures a higher resolution, and the certainty is maintained by determining the optimal values of the hyperparameters with ABIC. This point is a clear contrast with Shen's method; results are sensitive to the value of D, and both too large or too small D yield poor results.
Second, as shown in Appendix, the velocity and strain-rate fields estimated by Shen's method do not satisfy the relationship of differentiation. The discrepancy between the dilatation rate directly estimated by Shen's method and that calculated from the estimated velocity field through differentiation is presented in Fig. 12. The maximum discrepancy is 159.13, 30.03 and 14.83 nanostrain/year for D = 15 , 25 and 35 km, respectively. The discrepancy is generally small except for an overfitting model D = 15 km. In practice, no serious problem would occur in most regions. However, the discrepancy is evident around the Izu Islands (the bottom panels in Fig. 12), where observation stations are sparse, even for D = 25 km. A large discrepancy exists not only in the extrapolated sea area, but also in the vicinity of the GNSS stations on the islets. We should keep in mind the possibility of inconsistency in strain-rate fields. The discrepancy does not occur in the ABIC method because the Finally, we consider the problem of uncertainties of the estimated velocity and strain-rate fields. Figure 13 shows uncertainties of the absolute value of velocity vectors estimated by both methods. In Shen's method, the uncertainty is small in the whole land area and increases exponentially to the outside of the observation network. We also observe a smaller D results in a larger uncertainty despite a smaller residual (Fig. 11). This is because the uncertainty, defined by Eq. (8), depends only on input uncertainties σ i estimated at each observation site, and it does not reflect the actual fitting to velocity data v i . Equations (4) and (6) indicate that a smaller D yields larger errors at each data point (i.e., larger δv i and components of E P ), which leads to a larger estimated uncertainty (Eq. 8). As pointed out by Shen et al. (2015), this value itself cannot be used to measure the error of estimations. On the other hand, in the ABIC method, the covariance of the model parameters is given by Eq. (28) in which the value of σ is determined through the balance between data fit and smoothness of a model, which are related to the observed data d as well as the model parameter a (Eqs. 30 and 31). Therefore, the estimated uncertainty reflects the fitting accuracy of velocity fields. Figure 14 shows uncertainties of the dilatation-rate field. Characteristics are similar to those of the velocity field; the uncertainty monotonically decreases with D in Shen's method, while the ABIC method estimates larger uncertainty than Shen's method in most land areas. The change of uncertainty with D in Shen's method is even clearer in the strain-rate field. In Hokkaido, the result of the ABIC method shows a slightly larger uncertainty owing to sparser distribution of observation stations, while this characteristic is not so clear in the results of Shen's method of D = 25 km and 35 km.

Strain-rate field in Japan
In this section, we overview characteristics of the strain-rate fields in Japan revealed by the ABIC method (Figs. 15, 16 and 17). Although very high strain rates are observed along the Pacific coast particularly in southwestern Japan, we neglect them in the following discussion, because this crustal deformation is chiefly caused by interplate coupling and considered to be mostly cyclic.
The results of the principal strain rates (Fig. 15) and dilatation rates (Fig. 16) indicate that EW contraction dominates the Japanese Islands, which is consistent with previous studies of GNSS data analysis (e.g., Sagiya et al. 2000). This EW contraction is also in harmony with the stress field of Japan estimated from seismological data (e.g., Terakawa and Matsu'ura 2010) and geologically detected deformation in recent a few million years (e.g., Hujita 1980). Although the EW contraction is dominant in the Japanese Islands, the obtained strain-rate fields show large spatial variation. It is also observed that areas of high dilatation rates usually correspond to those of high shear-strain rates (Fig. 17), though there are exceptions.
A conspicuous high strain-rate zone passes through from the eastern margin of Japan Sea in the Tohoku district, via the Niigata Plain to Kobe, which is known as NKTZ (Sagiya et al. 2000). The high shear-strain rate zone further extends from Kobe to middle Kyushu along the Median Tectonic Line active fault system in Shikoku, which is consistent with the most active boundary in southwestern Japan revealed by a block modeling of GNSS data (Nishimura et al. 2018). We also observe a branch of a high strain-rate zone from around Hakusan volcano (136.8° E, 36.2° N) to Ise Bay along 137° E. This region has already been suggested as the contractive boundary between northeastern Japan (the North American plate) and southwestern Japan (the Amurian plate) from a GNSS data analysis (Heki and Miyazaki 2001).
In contrast, we have found that a low strain-rate zone extends in the forearc, from the Pacific coast of the southern Tohoku district (Fukushima Prefecture) to central Japan (Aichi Prefecture or possibly Nara Prefecture). This low strain-rate zone was briefly mentioned in Sagiya (2004), but the result of this study elucidates it more sharply due to denser GNSS stations and the improvement of the analysis method (Figs. 4 and 11). This forearc low strain-rate zone, which basically well corresponds to high seismic-velocity zone in the shallow crust (Nishida et al. 2008), would be attributed to a low geothermal gradient in the forearc (Tanaka et al. 2004) and less effects of interplate coupling between continental and oceanic  plates, that is, a low coupling ratio for the Pacific plate (Noda et al. 2013) and a low plate convergence rate due to strain partitioning of both continental and oceanic plates in central Japan (Heki and Miyazaki 2001;Nishimura et al. 2018).
The Ou Backbone Range in the Tohoku district has several contractive spots around active volcanoes, although the postseismic deformation of the 2008 Iwate-Miyagi inland earthquake is also included. It is interesting that the locations of these contractive spots well correspond Fig. 15 The principal axes of the strain-rate field estimated using the ABIC method with L = 20 km. Red and blue arrows represent expansive and contractive strain rates, respectively. Green stars represent the epicenters of earthquakes whose coordinate offsets are removed from GNSS data. Yellow lines trace MTL and OBR. Black lines represent surface traces of major active faults (Headquarters for Earthquake Research Promotion 2017). Red triangles represent Mt. Fuji, Hakusan and Sakurajima, and black triangles represent other active volcanos (JMA 2013). The following acronyms are used: MTL Median Tectonic Line active fault system, OBR Ou Backbone Range, IB Ise Bay, A Aichi, F Fukushima, I Ibaraki, K Kobe, Na Nara, Ni Niigata to the subsidence areas detected by InSAR after the 2011 Tohoku-oki earthquake , although an earlier period (1997)(1998)(1999)(2000) of GNSS data did not show such a good coincidence (Miura et al. 2002). This coincidence suggests that strain-rate variations with such short wavelengths catch meaningful signals. This contraction may be attributed to a high geothermal gradient under the stress condition of EW compression (Shibazaki et al. 2008).
Although contraction is dominant in the Japanese Islands, there are several regions with areal expansion, which are better captured owing to improved resolution than in previous studies. The most distinctive expansion occurs in the Izu Islands, where backarc spreading is ongoing (Nishimura 2011). The dilatation is also apparent  The maximum shear-strain rate field estimated from the ABIC method with L = 20 km. Green stars represent the epicenters of earthquakes whose coordinate offsets are removed from GNSS data. Black lines represent surface traces of major active faults (Headquarters for Earthquake Research Promotion 2017). White triangles represent Mt. Fuji, Hakusan and Sakurajima, and black triangles represent other active volcanos (JMA 2013) earthquakes is likely to be the cause of these weak dilatation and high shear-strain rates. Weak positive dilatation along the Pacific coast in Fukushima and Ibaraki prefectures might also be related to the 2008 Fukushima-oki (M6.9) and the 2008 Ibaraki-oki (M7.0) earthquakes as well as decadal transient deformation before the 2011 Tohoku-oki earthquake (Mavrommatis et al. 2014), though shear-strain rates are low in these regions. It is interesting that positive dilatation occurred along the Pacific coast of Fukushima prefecture before the 2011 Tohoku-oki earthquake, where a large normal fault event (M7.0), which is very rare in Honshu, Japan, occurred after the 2011 Tohoku-oki earthquake ).

Conclusions
Among various methods for estimating strain-rate fields, we investigated theoretical properties of a widely applied method of Shen et al. (1996Shen et al. ( , 2015. In this method, the estimated strain-rate field is mathematically not consistent with the estimated velocity field (Appendix). Moreover, the value of a hyperparameter D (Eq. 4) that controls the degree of smoothness must be manually selected, and estimated uncertainty of the velocity and strain-rate fields is unreliable. We then proposed a method using basis function expansion with ABIC to resolve these difficulties.
We applied the two methods to GNSS data (2006)(2007)(2008)(2009) in Japan, and examined the characteristics of these methods. The basis function expansion with ABIC yielded a better result than Shen's method in terms of the trade-off curve between the residual and roughness (Fig. 11). If we use an appropriate value of D = 25 km, the results of the ABIC method and Shen's method are quite similar, but significant differences exist in islet and coastal regions, where the density of observation stations drastically changes, and in local large deformation areas, where postseismic deformation and volcanic inflations occur (Figs. 2, 3, 4, 5 and 6). In Shen's method, the self-inconsistency between the estimated velocity and strain-rate fields does not have serious effects on most regions (Fig. 12), but the estimated uncertainty is certainly unreliable (Figs. 13 and 14).
A main practical concern of Shen's method is the determination of the value of the hyperparameter D, which significantly changes the estimation result and interpretation of crustal deformation. Other mathematical interpolation methods also require adjustment of hyperparameters to specify the smoothness of solutions. The optimal values are usually determined by checking whether the obtained result is in harmony with expected forms of crustal deformation. That is to say, subjective judgement by analysts affects the selection of a final result. This may allow us to show goodlooking results, but could lead to a biased solution. On the other hand, ABIC objectively determines the optimal values of hyperparameters in the method of basis function expansion. Although the interval L of basis functions must be chosen, the dependence is simple: estimation results are reasonable and hardly change for sufficiently small L (Figs. 7 and 8). Therefore, once we check that the value of ABIC and estimation result do not significantly change below a certain L 0 (approximately 30 km in the current dataset), any value of L ≤ L 0 can be chosen as long as computational cost allows.
The strain-rate field in Japan (Figs. 15, 16 and 17) estimated from the basis function expansion with ABIC detected a forearc low strain-rate zone that passes from the southern Tohoku district to central Japan. The result also shows several contractive spots around active volcanoes in the Ou Backbone Range, the locations of which well correspond to the subsidence areas detected by InSAR after the 2011 Tohoku-oki earthquake . Thus, the method of basis function expansion with ABIC would serve as an effective tool for estimating strain-rate fields from GNSS data.
The difference between w(x) and the derivative of v(x) is calculated from Eq. (34)  We calculate F (x) by specifying the number of data N . For N = 2 , computation leads to F (x) = 0 . The strainrate field is spatially uniform for N = 2 , so the consistency is maintained. For N = 3 , we obtain with This indicates that the relationship of differentiation is not satisfied except for the two cases (i) (x 2 − x 3 )(x 3 − x 1 )(x 1 − x 2 ) = 0 , and (ii) (x 2 − x 3 )y 1 + (x 3 − x 1 )y 2 + (x 1 − x 2 )y 3 = 0 . The (36) (38) F (x) = f 1 f 3 − f 2 2 g 0 + f 1 f 2 − f 0 f 3 g 1 + f 0 f 2 − f 2 1 g 2 . (42) condition (i) indicates that the positions of two data points (e.g., x 1 and x 2 ) are identical. Data at the same position should be treated as one data, so this virtually corresponds to the case of N = 2 . The condition (ii) indicates the collinearity equation of the three points. In this case, all the three data points align on a straight line, and so the strain-rate field is spatially uniform. In other words, when a strain-rate field changes spatially, the relationship of differentiation between the estimated velocity and strain-rate fields is not satisfied. It is expected that the relationship does not generally hold for N ≥ 4 . Formally, if x 3 = · · · = x N and the three points x 1 , x 2 and x 3 are not collinear, then F (x) = 0 as in the case of N = 3 . Since F defined by Eq. (41) is continuous with respect to (x, x 1 , . . . , x n ) , sufficiently small perturbation of the positions of the data points ( x 1 , . . . , x n ) maintains the inequality F (x) = 0 . These data points constitute examples of F (x) = 0 for truly N different data points. In fact, as shown in Fig. 12, the relationship of differentiation is not satisfied for the GNSS data analyzed in this study.