Optimal interpolation of spatially discretized geodetic data

We present an algorithm to calculate horizontal strains (or strain rates) through interpolation of geodetically derived displacements (or velocities). This is an underdetermined inverse problem to derive smoothly distributed strains (or strain rates) using spatially discretized geodetic observations. A priori information, in the form of weighted smoothing, is critical to facilitate the solution. At a given site, the horizontal displacement (or velocity) field in its vicinity is approximated by a bilinear function and represented by rigid block translation, rotation, and strains (or their rates). The weighted displacement (or velocity) data in the neighborhood are used to estimate the field parameters through a least-squares inversion procedure. Optimal weightings are prescribed for the neighboring data, based on their distances to the interpolation site and their spatial coverage. Nonelastic strains resulted from surface fault rupture and creep may also be excluded from the solution. We apply this method to the Southern California Earthquake Center Crustal Motion Map version 4.0 velocity field and derive the strain-rate field in southern California. Our result shows that (1) distance-dependent weighting can be optimally achieved by employing either a Gaussian or quadratic decay function, with the former offering a slightly sharper result than the latter. (2) Spatially dependent weighting is important to improve the interpolation, and can be done by invoking either an azimuthal weighting or a Voronoi cell areal weighting function. (3) The strain-rate pattern in southern California is dominated by dextral shear of the San Andreas fault (SAF) system, and the secondary faults surrounding the Big Bend of the SAF strike at oblique angles with respect to the maximum shear direction, suggesting that tectonic deformation field on and off the SAF is dominated by mechanic processes of the SAF.