<figure><img src="image_1.png" /></figure>

HAL
open science

# Theoretical relationship between back-projection imaging and classical linear inverse solutions

Yukitoshi Fukahata, Yugi Yagi, Luis Rivera

► To cite this version:

Yukitoshi Fukahata, Yugi Yagi, Luis Rivera. Theoretical relationship between back-projection imaging and classical linear inverse solutions. Geophysical Journal International, 2013, 196 (1), pp.552-559. 10.1093/gji/ggt392. hal-00932059

HAL Id: hal-00932059

https://hal.science/hal-00932059

Submitted on 8 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Theoretical relationship between back-projection imaging and classical linear inverse solutions

Yukitoshi Fukahata,<fnref n="1" />,* Yuji Yagi<fnref n="2" /> and Luis Rivera<fnref n="3" />

Accepted 2013 September 24. Received 2013 August 19; in original form 2013 May 10

## SUMMARY

We clarified the theoretical relationship between the back-projection (BP) imaging and classical linear inverse solutions via the hybrid back-projection (HBP) imaging. In the HBP imaging, which is mathematically similar to the time-reversal source imaging, cross correlations of observed waveforms with the corresponding Green's functions are calculated. The key condition for BP to work well is that the Green's function is sufficiently close to the delta function. Then, the BP image represents the slip motion on a fault, and approximately equals to the least-squares solution (LSS). In HBP, instead of the Green's function in BP, the stacked autocorrelation function of the Green's function must be close to the delta function to obtain a fine image. Because the autocorrelation function is usually closer to the delta function than the original function, we can expect that HBP works better than BP, if we can reasonably estimate the Green's function. With additional condition that the stacked cross-correlation function of the Green's functions for different source locations is small enough, the HBP image is also approximately equal to the LSS. If these assumptions are not satisfied, however, the HBP image corresponds to a damped LSS with an extremely large damping parameter, which is clearly inferior to usual inverse solutions. On the other hand, the advantages of the BP method are much less computation and no necessity of Green's functions.

**Key words:** Inverse theory; Earthquake dynamics; Earthquake source observations; Theoretical seismology.

## 1 INTRODUCTION

In order to clarify the rupture process of earthquakes, inversion analyses of seismic waveform data have traditionally been carried out (e.g. Trifunac 1974; Olson & Apsel 1982; Hartzell & Heaton 1983; Ji et al. 2002; Piatanesi et al. 2007; Yagi & Fukahata 2011a). In the inversion analysis, observed data are inverted to estimate model parameters, using some theoretical relationship between the observed data and model parameters, commonly on the basis of least squares principle. In the inverse problem to estimate the rupture process, the model parameters are generally taken to represent slip motion on the fault.

In recent years, the back-projection (BP) analysis has become a popular tool to image the rupture process of large earthquakes since the success of Ishii et al. (2005), which clarified the northward rupture progression of the 2004 Sumatra–Andaman earthquake. In the BP analysis, observed waveform data are not inverted, but simply stacked with appropriate time-shifts. In these days, once a

large earthquake happened, the BP image has been presented quite soon, for example, for the 2005 Nias (Walker et al. 2005), the 2008 Wenchuan (Xu et al. 2009), the 2010 Chili (Kiser & Ishii 2011), the 2011 Tohoku (Koper et al. 2011; Meng et al. 2012; Yagi et al. 2012) and the 2012 far-off northern Sumatra (Wang et al. 2012). In some of these studies, improvements of the BP method were also carried out. Xu et al. (2009) investigated the effects of array distributions and stacking methods. Yagi et al. (2012) proposed a hybrid back-projection (HBP) method in which cross correlations of observed waveforms with the corresponding Green's functions were calculated. Kiser et al. (2011) and Wang et al. (2012) also added the information of Green's functions to obtain an image with better resolution.

We now see so many images obtained by the BP method. On the other hand, it has not been clear what the BP image represents physically. In addition to this, the relationship between the inverted solution and the BP image has also not been understood yet. Ishii et al. (2005) stated that the obtained BP image represents energy release. Of course, the BP image must be related to energy release in some way. However, no explicit relationship was shown. Recently, Yao et al. (2012) showed that the time-evolving BP image

<footnote n="1">¹Disaster Prevention Research Institute, Kyoto University, Kyoto, Japan. E-mail: fukahata@rcep.dpri.kyoto-u.ac.jp</footnote>

<footnote n="2">²Graduate School of Life and Environmental Sciences, University of Tsukuba, Japan</footnote>

<footnote n="3">³Institut de Physique du Globe de Strasbourg, Université de Strasbourg/CNRS, France</footnote>

<footnote>*Now at: Department of Earth Sciences, University of Oxford, UK.</footnote>
approximately corresponds to the source time function at the corresponding source location. However, the situation is a bit more complicated. In this paper, we more clearly show the conditions under which that approximation is valid and we will see the BP analysis intends to obtain the time-evolving image of slip motion on the fault.

Once what the BP image represents has become clear, it is now possible to compare the BP image with classical linear inverse solutions. However, the BP method is too simple to be directly compared to inverse solutions, although the simplicity is the superior point of the BP analysis. Therefore, in order to compare BP with inverse solutions, we pass through the HBP method. That is to say, we make comparisons of BP to HBP and HBP to inversion.

## 2 METHODS OF SEISMIC SOURCE IMAGING

In general, seismic waveform *d* observed at a station *k* is related to fault slip *a* through a Green's function *G*:

$$d_k(t) = \int_S [a(\xi) * G_k(\xi)](t) d\xi = \int_S \int_{-\infty}^{\infty} a(\xi, \tau) G_k(\xi, t-\tau) d\tau d\xi, \quad (1)$$

where *S* is a fault surface and *\* denotes convolution. In a far field, *a* represents slip velocity on the fault for displacement data *d* or slip acceleration on the fault for velocity data *d*, if *G* is the Green's function for displacement waveform. The Green's function $G_k(\xi, t - \tau)$ means the deformation (e.g. displacement or velocity) at the observation station *k* at time *t* generated by a unit impulse slip on the fault plane *S* at the point $\xi$ at time $\tau$. For simplicity, we consider only one component for *d* and *a*, respectively, and so we omitted the suffixes which represent the directions of deformation at the observation station and slip motion at the source in eq. (1). The origin of time (*t* = 0) is set at the initiation of an earthquake. Then, we may take the range of the integral with respect to $\tau$ to be [0, *t*) due to causality, instead of the infinite integral in eq. (1). We may rewrite eq. (1) by taking a summation for the integration with respect to the fault surface *S* as

$$d_k(t) = \sum_l (a_l * G_{kl})(t), \quad (2)$$

where $a_l$ represents slip motion at the *l*th source grid on the fault plane *S* and $G_{kl}$ means the deformation at the observation station *k* generated by a unit impulse slip at the *l*th source grid on the fault plane *S*.

### 2.1 Back-projection

The basic equation of the BP imaging is expressed as follows (Ishii et al. 2005):

$$s_l^{\text{BP}}(t) = \sum_k c_k d_k(t + t_{kl}^p), \quad (3)$$

where $s_l^{\text{BP}}$ represents the stacked BP image at the *l*th potential source grid, $t_{kl}^p$ is the predicted traveltime between the grid *l* and the station *k* and $c_k$ is the normalizing (weighting) factor for each seismogram. Concerning observed data $d_k$, only the *P*-wave portion is usually used. Although the sign in front of $t_{kl}^p$ is taken to be minus in Ishii et al. (2005, 2007), the sign must be plus if we use a common time axis for the source and observation points, because the first seismic wave arrives at the station *k* at $t = t_{kl}^p$ and $d_k(t)$ is zero for $t < t_{kl}^p$.

By substituting eq. (2) into (3), we obtain

$$s_l^{\text{BP}}(t) = \sum_k c_k \sum_{\ell} (a_{\ell} * G_{k\ell})(t + t_{kl}^p). \quad (4)$$

In the BP image $s_l^{\text{BP}}(t)$, it is expected that the waveforms due to slips on gridpoints except for *l* are cancelled out each other (Ishii et al. 2005). This condition that holds true for arbitrary $a_{\ell}$ can be expressed as

$$\sum_k c_k G_{k\ell}(t + t_{kl}^p) \approx 0 \quad (l \neq \ell, \forall t). \quad (5)$$

Then, eq. (4) becomes

$$s_l^{\text{BP}}(t) \approx \sum_k c_k (a_l * G_{kl})(t + t_{kl}^p). \quad (6)$$

It should be noted that $s_l^{\text{BP}}(t)$ includes the contribution of the slip on gridpoints other than *l*, when the condition (eq. 5) does not hold true.

In the basic equation of BP (eq. 3), the theory of elastodynamics is used only for the predicted traveltime $t_{kl}^p$. Information concerning the later phases, inelastic attenuation and amplitude is neglected. Therefore, we can consider that in BP analyses the Green's function is implicitly assumed to be like the delta function with an amplitude that does not depend on the source location and mechanisms

$$G_{kl}(t + t_{kl}^p) \approx c\delta(t), \quad (7)$$

where *c* is a constant proportional to the amplitude of the Green's function, which may depend on the station location *k*. The condition (eq. 7) can be loosened to

$$\sum_k c_k G_{kl}(t + t_{kl}^p) \approx \delta(t), \quad (8)$$

where the proportional constant *c* is included in the weighting factor $c_k$. It has been reported that the images of BP are clearer for deep earthquakes (Suzuki & Yagi 2011). This is reasonable because the assumption of eq. (8) is more appropriate for such cases. If eq. (8) is not a good approximation, $s_l^{\text{BP}}(t)$ shows a smeared image in time.

By substituting eq. (8) into (6), we finally obtain

$$s_l^{\text{BP}}(t) \approx a_l(t). \quad (9)$$

Eq. (9) shows that the BP image $s_l^{\text{BP}}(t)$ is directly related to the motion on the fault. A similar result has already been obtained by Yao et al. (2012). In brief, BP intends to obtain the time-evolving image of slip motion for an earthquake, under the assumptions of eqs (5) and (8). In the BP analysis, however, the Green's function is not explicitly used. So, we never know the BP image represents slip velocity or acceleration. This point will be further discussed in a later section.

### 2.2 Hybrid back-projection

Yagi et al. (2012) proposed an improvement of BP by taking the cross correlation of observed waveforms and the Green's function

$$s_l^{\text{HBP}}(t) = \sum_k c'_k (d_k \widehat{\times} G_{kl})(t) = \sum_k c'_k \int_{-\infty}^{\infty} d_k(\tau) G_{kl}(\tau - t) d\tau, \quad (10)$$

where $\widehat{\times}$ denotes cross correlation and $c'_k$ is a normalizing (weighting) factor for each seismogram. Although an *N*th root stack is used in Yagi et al. (2012), we used the same stacking procedure as BP in eq. (10) to focus on the essential difference between them. Here, it
should be noted that we can rewrite the basic equation of BP (eq. 3) as

$$s_l^{BP}(t) = \sum_k c_k (d_k \widehat{\times} \delta)(t + t_{kl}^p). \quad (11)$$

This means that the BP method is equivalent to stacking the cross correlations of observed waveforms with the shifted delta function. This observation lead Yagi et al. (2012) to propose stacking the cross correlations of observed waveforms with the Green's function, because the synthetic waveform for an impulse slip is the Green's function. Yagi et al. (2012) called the imaging method expressed by eq. (10) the hybrid BP, because this method is relatively closer to inversion analyses. However, the method is essentially an improvement of BP.

From eqs (2) and (10), we obtain

$$s_l^{HBP}(t) = \sum_k c'_k \left( \sum_l (a_l * G_{kl}) \widehat{\times} G_{kl} \right) (t). \quad (12)$$

HBP also expects that the waveforms due to slips on gridpoints except for $l$ are cancelled out each other by stacking. This condition that holds true for arbitrary $a_l$ is expressed as

$$\sum_k c'_k (G_{kl} \widehat{\times} G_{kl})(t) \approx 0 \quad (l \neq l, \forall t). \quad (13)$$

With the condition (13), eq. (12) becomes

$$s_l^{HBP}(t) \approx \sum_k c'_k (a_l * Au_{G_{kl}})(t), \quad (14)$$

where

$$Au_{G_{kl}}(t) = (G_{kl} \widehat{\times} G_{kl})(t) = \int_{-\infty}^{\infty} G_{kl}(\tau) G_{kl}(\tau - t) d\tau, \quad (15)$$

where $Au_{G_{kl}}(t)$ is the autocorrelation function of $G_{kl}(t)$. If we can further assume that

$$Au_{G_{kl}}(t) \approx c' \delta(t) \quad (16)$$

or more precisely

$$\sum_k c'_k Au_{G_{kl}}(t) \approx \delta(t) \quad (17)$$

then we have the same expression for HBP as eq. (9) for BP

$$s_l^{HBP}(t) \approx a_l(t), \quad (18)$$

where $c'$ is a proportional constant. From eq. (18), we can consider that HBP also intends to estimate the slip motion on the fault. Because the Green's function is explicitly given in HBP, there is no ambiguity as to what the HBP image represents physically. The HBP image represents slip velocity on the fault for displacement data $d$ or slip acceleration on the fault for velocity data $d$, if we use the Green's function for displacement waveform.

Based on the comparison of eqs (5) and (13), and (8) and (17), we can judge the goodness of the images obtained by BP and HBP. The autocorrelation function always has its maximum at $t = 0$. Therefore, we can expect that the approximation of (17) is better than (8), especially for shallow earthquakes, for which reflection phases at the Earth's surface ($pP$ and $sP$) are more important. However, it should be noted that the two Green's functions included in the auto- and cross-correlation functions are not identical: one is theoretically or empirically estimated and the other is true and unknown. This means that the obtained image becomes worse if we use a Green's function that is not close to the true one.

## 2.3 Time reversal

In the time reversal (TR) imaging of seismic sources, the time reversed observed waveforms are back propagated to the source region to search the locations where the waves are constructively interfered. It can be written as (Larmat et al. 2006; Kawakatsu & Montagner 2008)

$$s_l^{TR}(t) = \sum_k c_k'' \int_{-\infty}^{\infty} d_k(t_0 - \tau) \tilde{G}_{lk}(t - \tau) d\tau, \quad (19)$$

where $t_0$ is an arbitrary reference time and $c_k''$ is a weighting factor. Here, it should be noted that the Green's function $\tilde{G}_{lk}$ is not for a displacement discontinuity across a fault surface as in eq. (1), but for a single force. $\tilde{G}_{lk}(t - \tau)$ means the displacement at the point $l$ at time $t$ generated by a single force applied at the point $k$ at time $\tau$, and $d_k$ is the displacement waveform at a station $k$. By changing the integral variable in eq. (19), we obtain

$$s_l^{TR}(t) = \sum_k c_k'' \int_{-\infty}^{\infty} d_k(\tau) \tilde{G}_{kl}(\tau - t_0 + t) d\tau \\ = \sum_k c_k'' (d_k \widehat{\times} \tilde{G}_{kl})(t_0 - t), \quad (20)$$

where we used the reciprocity of the Green's function. Eq. (20) shows TR also stacks the cross-correlations of observed waveforms with the Green's function as in eq. (10), although the direction of the time axis is reversed. The relation of the times between HBP (eq. 10) and TR (eq. 19) is expressed as

$$t_{HBP} = t_0 - t_{TR}. \quad (21)$$

As can be seen in eq. (20), HBP is closely related with TR, although the definition of the Green's function and original ideas are different.

# 3 RELATION TO INVERSE SOLUTIONS

## 3.1 Vector expression of seismic source imaging

In order to compare seismic source images of the BP and HBP analyses with inverse solutions, we introduce vector expressions. When we consider only one source gridpoint $l$ and one observed waveform at a station $k$, eq. (2) may be written in the following form with a discretization of time:

$$\mathbf{d}_k = \mathbf{G}_{kl} \mathbf{a}_l, \quad (22)$$

where the components of the vectors and matrix are as follows:

$$(\mathbf{d}_k)_i = d_k(t_i), \quad (\mathbf{G}_{kl})_{ij} = G_{kl}(t_i - t_j), \quad (\mathbf{a}_l)_j = a_l(\tau_j), \quad (23)$$

where $(\cdot)_i$ denotes the $i$th component of the vector and $(\cdot)_{ij}$ denotes the $ij$ component of the matrix. When the numbers of the sources and observation stations are not limited to one, eq. (22) is generalized as

$$\mathbf{d} = \mathbf{G} \mathbf{a}, \quad (24)$$

with

$$\mathbf{d} = \begin{pmatrix} \mathbf{d}_1 \\ \mathbf{d}_2 \\ \vdots \\ \mathbf{d}_K \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} \mathbf{G}_{11} & \mathbf{G}_{12} & \cdots & \mathbf{G}_{1L} \\ \mathbf{G}_{21} & \mathbf{G}_{22} & \cdots & \mathbf{G}_{2L} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{G}_{K1} & \mathbf{G}_{K2} & \cdots & \mathbf{G}_{KL} \end{pmatrix}, \quad \mathbf{a} = \begin{pmatrix} \mathbf{a}_1 \\ \mathbf{a}_2 \\ \vdots \\ \mathbf{a}_L \end{pmatrix}, \quad (25)$$
The vectors **d** and **a** and the matrix **G** are composed of the subvectors **d**<sub>k</sub> and **a**<sub>l</sub> and the submatrices **G**<sub>kl</sub>.

Similarly, we can express the basic equation of HBP (eq. 10) as

$$
s_l^{\text{HBP}} = \sum_k c_k' \mathbf{G}_{kl}^T \mathbf{d}_k \qquad (26)
$$

with

$$
(s_l^{\text{HBP}})_j = s_l^{\text{HBP}}(\tau_j), \quad (27)
$$

where the superscript *T* denotes transpose. When we consider more
than one gridpoints for the source, eq. (26) becomes

$$
\mathbf{s}^{\text{HBP}} = \mathbf{G}^{\text{T}} \mathbf{C}' \mathbf{d} \quad (28)
$$

with

$$
\mathbf{s}^{\text{HBP}} = \begin{pmatrix} \mathbf{s}_1^{\text{HBP}} \\ \mathbf{s}_2^{\text{HBP}} \\ \vdots \\ \mathbf{s}_L^{\text{HBP}} \end{pmatrix}, \quad \mathbf{C}' = \begin{pmatrix} c_1' \mathbf{I}_1 & & & \mathbf{0} \\ & c_2' \mathbf{I}_2 & & \\ & & \ddots & \\ & & & c_K' \mathbf{I}_K \end{pmatrix}, \quad (29)
$$

where $\mathbf{I}_k$ is the unit matrix with a suitable size.

In the same way, with an appropriate discretization of time, we
can express the basic equation of BP (eq. 3) for multiple source
gridpoints as

$$
\mathbf{s}^{\text{BP}} = \mathbf{B}^{\text{T}} \mathbf{C} \mathbf{d} \quad (30)
$$

with

$$
\begin{align}
\mathbf{s}^{\text{BP}} &= \begin{pmatrix} \mathbf{s}_1^{\text{BP}} \\ \mathbf{s}_2^{\text{BP}} \\ \vdots \\ \mathbf{s}_L^{\text{BP}} \end{pmatrix}, &
\mathbf{B} &= \begin{pmatrix} \mathbf{B}_{11} & \mathbf{B}_{12} & \cdots & \mathbf{B}_{1L} \\ \mathbf{B}_{21} & \mathbf{B}_{22} & \cdots & \mathbf{B}_{2L} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{B}_{K1} & \mathbf{B}_{K2} & \cdots & \mathbf{B}_{KL} \end{pmatrix}, \nonumber \\
\mathbf{C} &= \begin{pmatrix} c_1 \mathbf{I}_1 & & & \mathbf{0} \\ & c_2 \mathbf{I}_2 & & \\ & & \ddots & \\ \mathbf{0} & & & c_K \mathbf{I}_K \end{pmatrix} \tag{31}
\end{align}
$$

and

$$
(\mathbf{B}_{kl})_{ij} = \begin{cases} 1 & (t_i - \tau_j = t_{kl}^P) \\ 0 & (t_i - \tau_j \neq t_{kl}^P) \end{cases} . \quad (32)
$$

In the BP image $s_l^{\text{BP}}$, observed data are simply stacked with time-shifts (eq. 3). The matrix $\mathbf{B}_{kl}$ merely controls the time-shifts. Incidentally, the expression of $\mathbf{B}_{kl}$ in eq. (32) can be obtained from $\mathbf{G}_{kl}$ in eq. (23) with the condition that the Green's function is approximated to the delta function, as can be seen in eqs (10) and (11).

3.2 Inverse solutions

The observation eq. (24), **d** = **Ga**, can be solved by inversion anal-
yses, in which we estimate the model parameter **a** from observed
data **d**.

The simplest way to estimate **a** practically used would be the
following expression (e.g. Claerbout 2001):

$$
\hat{\mathbf{a}} = \mathbf{G}^{\mathrm{T}} \mathbf{d}, \qquad (33)
$$

where the hat denotes the estimate. As can be seen in eq. (28), the
HBP image is nothing different from this estimate, if we take the
weighting factor $c_k$ to be common for all observed data. Note that
the BP and HBP images have a degree of freedom in amplification.
Some readers may notice that eq. (33) represents an estimation
with the adjoint operator ($\mathbf{G}^T$), which is widely used to obtain the
gradient of an evaluation function in improving the inverse solution
of non-linear problems (e.g. Tromp *et al.* 2005).

The most famous inverse solution is the least-squares solution
(LSS)

$$
\hat{\mathbf{a}} = (\mathbf{G}^T \mathbf{E}^{-1} \mathbf{G})^{-1} \mathbf{G}^T \mathbf{E}^{-1} \mathbf{d}, \quad (34)
$$

where it is assumed that the misfit, **d** − **Ga**, follows the Gaussian
distribution with mean zero and covariance **E**. The covariance ma-
trix **E** controls the weight of each data, and so **E** has a similar
function to **C** and **C'**. If the variance of the misfit for each data is
common and independent of each other (E⁻¹ = σ²**I**), the LSS is
simplified as

$$
\hat{\mathbf{a}} = (\mathbf{G}^T \mathbf{G})^{-1} \mathbf{G}^T \mathbf{d}. \quad (35)
$$

When we do not have enough amount of data with good accuracy,
as it is well known, the LSS (eqs 34 or 35) does not give a good
estimate of the model parameters. We can never obtain a meaning-
ful solution only by enhancing the fit to data, when the data are
insufficient. In such a situation, some studies (e.g. Aki et al. 1977)
have used a damped LSS

$$
\hat{\mathbf{a}} = (\mathbf{G}^T \mathbf{E}^{-1} \mathbf{G} + \varepsilon^2 \mathbf{I})^{-1} \mathbf{G}^T \mathbf{E}^{-1} \mathbf{d}, \quad (36)
$$

where $\varepsilon^2$ is called a damping parameter.

According to Jackson (1979) and Tarantola (2005), the damping term can be interpreted to express prior information about the problem. In this case, the damping term requires the solution length to be minimized. The relative weight to minimizing the misfit is controlled by the damping parameter $\varepsilon^2$. In inversion analyses for seismic source processes, smoothing constraint has been commonly used (e.g. Yagi & Fukahata 2011b). The smoothing constraint is also considered to be a kind of prior information (Matsu'ura et al. 2007).

3.3 Model resolution matrix

As shown in the previous section, there can be several inverse so-
lutions for the same problem. As a measure of the goodness of the
inverse solutions, we introduce a model resolution matrix (Menke
2012). The solution of the linear inverse problem (eq. 24) can be
symbolically expressed as

$$
\hat{\mathbf{a}} = \mathbf{G}^{-\mathbf{g}} \mathbf{d}, \qquad (37)
$$

where $\mathbf{G}^{-\mathbf{g}}$ is called the generalized inverse. By substituting eq. (24)
into eq. (37), we obtain

$$
\hat{\mathbf{a}} = \mathbf{G}^{-\mathbf{g}} \mathbf{G} \mathbf{a} = \mathbf{R} \mathbf{a}, \quad (38)
$$

where **R** is called the model resolution matrix. We can see that **â**
is a good estimate of **a**, when **R** is close to the identity matrix. It
should be noted that the model resolution matrix is independent of
the actual values of the data.

For the case of the LSS, the model resolution matrix is

$$
\mathbf{R}^{\text{LSS}} = \mathbf{I}. \tag{39}
$$

This means that $\hat{\mathbf{a}} = \mathbf{a}$ and each model parameter is perfectly re-
solved.

We have seen in Section 2 that $\mathbf{s}^{\text{HBP}}$ and $\mathbf{s}^{\text{BP}}$ are intended to
represent slip motion **a** on the fault. In other words, $\mathbf{s}^{\text{HBP}}$ and $\mathbf{s}^{\text{BP}}$
are a kind of estimates of the model parameter a. Therefore, by substituting eq. (24) into eqs (28) and (30), we can obtain the model resolution matrices for the HBP and BP images as

$$ \mathbf{R}^{\text{HBP}} = \mathbf{G}^{\text{T}} \mathbf{C}' \mathbf{G} \quad (40) $$

and

$$ \mathbf{R}^{\text{BP}} = \mathbf{B}^{\text{T}} \mathbf{C} \mathbf{G} \quad (41) $$

Because $\mathbf{R}^{\text{HBP}}$ and $\mathbf{R}^{\text{BP}}$ are not the identity matrix, the HBP and BP images are clearly inferior to the LSS in terms of the model resolution matrix.

However, if the model resolution matrices, $\mathbf{R}^{\text{BP}}$ and $\mathbf{R}^{\text{HBP}}$, are close to the identity matrix, it does not cause a serious problem. From simple algebra, we obtain

$$ \mathbf{R}^{\text{HBP}} = \begin{pmatrix} \sum_k c'_k \mathbf{G}_{k1}^T \mathbf{G}_{k1} & \sum_k c'_k \mathbf{G}_{k1}^T \mathbf{G}_{k2} & \cdots & \sum_k c'_k \mathbf{G}_{k1}^T \mathbf{G}_{kL} \\ \sum_k c'_k \mathbf{G}_{k2}^T \mathbf{G}_{k1} & \sum_k c'_k \mathbf{G}_{k2}^T \mathbf{G}_{k2} & \cdots & \sum_k c'_k \mathbf{G}_{k2}^T \mathbf{G}_{kL} \\ \vdots & \vdots & \ddots & \vdots \\ \sum_k c'_k \mathbf{G}_{kL}^T \mathbf{G}_{k1} & \sum_k c'_k \mathbf{G}_{kL}^T \mathbf{G}_{kL} \end{pmatrix} \quad (42) $$

The conditions that $\mathbf{R}^{\text{HBP}}$ is equal to the identity matrix are expressed as follows:

$$ \sum_k c'_k \mathbf{G}_{k\ell}^T \mathbf{G}_{k\ell} = \begin{cases} \mathbf{0} & (l \neq \ell) \\ \mathbf{I} & (l = \ell) \end{cases}, \quad (43) $$

The conditions of eq. (43) exactly correspond to eqs (13) and (17), which are the assumptions made in the HBP analysis in order for the HBP image to represent slip motion on the fault.

In the same way, the conditions that $\mathbf{R}^{\text{BP}}$ is equal to the identity matrix are

$$ \sum_k c_k \mathbf{B}_{k\ell}^T \mathbf{G}_{k\ell} = \begin{cases} \mathbf{0} & (l \neq \ell) \\ \mathbf{I} & (l = \ell) \end{cases}. \quad (44) $$

These conditions also exactly correspond to eqs (5) and (8). Incidentally, if the Green's function is like the delta function, the conditions of eq. (44) are approximately satisfied.

As shown in eqs (43) and (44), if the assumptions of (13) and (17) or (5) and (8) are satisfied, the HBP and BP images represent slip motion on the fault. However, these conditions are usually not satisfied. Then, what will happen? In the damped LSS (eq. 36), by taking $\epsilon^2$ to be extremely large, we obtain

$$ \hat{\mathbf{a}} \approx c \mathbf{G}^{\text{T}} \mathbf{E}^{-1} \mathbf{d}, \quad (45) $$

where c is a proportional constant. This solution is nearly equal to zero practically, but not exactly zero. Considering that $\mathbf{E}$ has a similar function to $\mathbf{C}'$ and that the HBP image has a degree of freedom in amplification, the estimate of eq. (45) is equivalent to the HBP image (eq. 28). That is to say, the HBP image corresponds to a damped LSS with an extremely large damping parameter. In the BP image, the Green's function in $\mathbf{G}$ in eq. (45) is further approximated by the delta function.

## 4 SIMPLE NUMERICAL EXAMPLES AND DISCUSSION

In this section, through simple numerical computations, we investigate the validity of the approximations used in the BP and HBP methods, and compare the images obtained by these methods with the damped LSS.

<figure><img src="image_3.png" /><figcaption>Figure 1. Time profiles of the stacked Green's function (top) and the stacked autocorrelation function of the Green's function (bottom) for a strike-slip point source on a vertical fault (solid lines). Those profiles for dip-slip point sources on a vertical fault (dotted lines) and on a thrust fault with a dip angle of 30° (broken lines) are also shown. All point sources are put at a depth of 15 km. Observation stations are distributed every 10° with the radial distance from 30° to 90° and all azimuthal direction. These time profiles represent the left-hand side of eq. (8) for top and that of eq. (17) for bottom. In order to obtain a sharp image by the BP and HBP analyses, these profiles must be close to the delta function. When an impulse slip was given at the point source, these profiles correspond to the time-evolving image of BP (top) and HBP (bottom) at the location of the point source.</figcaption></figure>

We put an impulse strike-slip point source on a vertical fault plane at a depth of 15 km and $t = 0$. With the method of Kikuchi & Kanamori (1991), we computed the synthetic velocity waveform at observation stations that are distributed every 10° with the radial distance from 30° to 90° and all azimuthal direction; the total number of the stations is 252. The synthetic waveforms were calculated every 0.05 s and no errors were added. Because an impulse slip was given, the synthetic velocity data $d_k$ is equal to the Green's function $G_{k\ell}$ for velocity waveform, where $l$ represents the point-source location. The source grid interval in the BP and HBP images were taken to be 1 km × 1 km. We gave an equal weight for each observed data ($c_k$ and $c'_k$ are constant).

In Fig. 1, we show the time profiles of the stacked Green's function (top) and the stacked autocorrelation function of the Green's function (bottom) by solid lines. These profiles correspond to eq. (8) for top and eq. (17) for bottom. As shown in the previous sections, in order to obtain a sharp image by the BP and HBP methods, these profiles must be close to the delta function. The Green's function, however, commonly has later phases, such as $pP$ and $sP$, which are not easily cancelled out by stacking. In the case of a strike-slip on a vertical fault, the amplitude of a later phase after stacking is larger than the direct $P$ phase, which results in ghosts in the BP image. For the case of HBP, the situation is better. The highest peak is always taken at $t = 0$ owing to the nature of autocorrelation functions. However, the effect of later phases still remains and ghosts emerge even before $t = 0$. In Fig. 1, we also show dip-slip cases on a vertical fault (dotted lines) and on a thrust fault with a dip angle of 30° (broken lines). Even if we consider a dip slip, as shown in Fig. 1, the situation is similar; the stacked Green's function for BP and the stacked autocorrelation function of the Green's function for HBP are significantly different from the delta function, although the direct $P$ phase takes the largest positive value.
In Fig. 2, we show the images of BP (left-hand side) and HBP (middle) on the vertical fault plane for the case of the strike-slip point source. These images and the time profiles in Fig. 1 are normalized by the maximum value of each image. Since the observed waveform is equal to the Green's function, these images correspond to the model resolution matrices on the fault plane for the given point source. The obtained BP and HBP images are not confined to the black dot at $t = 0$, where the source was put, but widely smeared in space and time. This means the stacking procedure to cancel out the waveforms does not work efficiently, because the Green's functions in the smeared area are quite similar to that at the point source. Apparent movement of the source is also observed for both BP and HBP images. The source moves from top to bottom in BP, because the traveltime is longer for shallower grids. On the other hand, the source moves from bottom to top in HBP, because the traveltime of the sP phase, which is dominant in the computation of the autocorrelation functions in the present case, is shorter for shallower grids. Such artefacts may result in incorrect estimation of the rupture velocity. In the BP image, the largest peak emerges at a different time and location from the true ones. Even in an ideal situation, such as equally distributed many stations, no errors and a point source, the images show such artefacts and have limited resolution. The results show that the resolution matrices for BP (eq. 41) and HBP (eq. 40) are not well approximated by the identity matrix.

We also show a damped LSS in the right-hand column of Fig. 2. The image also corresponds to the model resolution matrix for the given point source. In obtaining the solution, we added a very small damping parameter to avoid numerical instability due to computational errors. As shown in Fig. 2, the damped LSS with a small damping parameter is completely different from the images obtained by the BP and HBP methods. The difference simply comes from the equations used in the estimations: $\mathbf{B}^T \mathbf{d}$ (eq. 30) for BP, $\mathbf{G}^T \mathbf{d}$ (eq. 28) for HBP and $(\mathbf{G}^T \mathbf{G} + \varepsilon^2 \mathbf{I})^{-1} \mathbf{G}^T \mathbf{d}$ (eq. 36) for the damped LSS. This example vividly shows the importance of equations used in estimations. If some errors are added in data, we should use a larger damping parameter, which results in a more smeared image of the damped LSS. However, the solution is still better than the image of HBP. Only in the exceptional condition that the damping parameter is extremely large, the damped LSS converges to the HBP image.

In order to enhance the correlation of waveforms, the BP analysis often uses a local network. When we confine the observation stations from 0° to 90° in the azimuth, the obtained images by both BP and HBP becomes horizontally asymmetrical (Fig. 3), because the traveltimes are different for source grids of the same depth (the traveltime is longer for grids in the positive direction of the horizontal axis.). Here, we kept the number of stations the same, that is, the observation stations are four times denser than the previous case (Fig. 2). Because an impulse slip was given at the point source, these images also correspond to the model resolution matrix for the given point source. Note that the model resolution matrix is independent of observed data, but depends on the network geometry. As shown in Fig. 3, if we use a local network, the resolution of the BP and HBP images becomes worse. Larger false signals can be seen. On the other hand, the resolution is still good for the case of a damped LSS with a small damping parameter.

In terms of the model resolution matrix, as shown in Figs 2 and 3, the damped LSS is clearly superior to the BP and HBP images. On the other hand, in order to obtain the damped LSS, we need significantly higher computational power than in the BP and HBP analyses. It is not easy to compute the inverse matrix of $\mathbf{G}^T \mathbf{G}$, when the number of model parameters are very large. Because of that, we had to restrict the model fault area as shown by the dotted line in

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time-evolving images of BP (left-hand side), HBP (middle) and a damped least-squares solution (right-hand side) on the vertical fault plane due to an impulse strike-slip point source. The setting of the computation is the same as in Fig. 1. Because an impulse slip was given at the point source, these images correspond to the model resolution matrix (eqs 39–41) for the given point source. The amplitude of each image is normalized so that the maximum is 1. The black dots in the left-hand and middle columns at $t = 0$ represent the source location. In the diagrams of the right-hand column, the dotted square represents the model fault area, only for which the model parameters were estimated in order to reduce the computation task. In the damped least-squares solution, a very small damping parameter was added to avoid numerical instability due to computational errors.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3. Time-evolving images of BP (left-hand side), HBP (middle) and a damped least-squares solution (right-hand side) on the vertical fault plane due to an impulse strike-slip point source. The distribution of observation stations is confined from 0° to 90° in the azimuthal direction, although the number of the stations is kept the same. Except the distribution of observation stations, the figure is calculated in the same way as in Fig. 2. Because an impulse slip was given at the point source, these images also correspond to the model resolution matrix for the given point source.</figcaption></figure>

Figs 2 and 3, in order to avoid the overburden of computation. If we take the same model fault area as in the BP and HBP analyses, we have to use a larger fault patch and a time step, which results in low resolution of the solution. However, the resolution of the solution would be still better than those of the BP and HBP images for most cases. In the BP analysis, we can take a very small grid size, as in the computation of Fig. 2. Because of that, it has been sometimes considered that a higher resolution image can be obtained by BP. However, the resolution of BP is not restricted by the grid size in computation, but by the conditions of eqs (5) and (8) or (44). In other words, the BP analysis is useful for a gigantic earthquake, such as the 2004 Sumatra-Andaman earthquake.

Another important point of BP is that we do not need information on Green's functions. Without Green's functions, we can obtain a rough rupture image of earthquakes, which is an incomparable advantage of the BP analysis. This simplicity, however, can be the source of ambiguity in physically interpreting the BP image. As described in Section 2, seismic waveform d is related to fault slip a through a Green's function G

$$d_k^d(t) = \sum_l (a_l^v * G_{kl}^d)(t), \quad (46)$$

where the superscripts d and v denote displacement and velocity, respectively (i.e. $a_l^v$ represents slip velocity on the fault patch l). The BP analysis usually uses velocity data for seismogram (e.g. Yao et al. 2012). When we take time differential of eq. (46), there are at least two candidates for the right-hand side:

$$d_k^v(t) = \sum_l (a_l^a * G_{kl}^d)(t) = \sum_l (a_l^v * G_{kl}^v)(t), \quad (47)$$

where the superscript a denotes acceleration. If we consider a homogeneous elastic whole space with no attenuation, the Green's function for displacement waveform $G_{kl}^d$ is like the delta function. Hence, from a theoretical point of view, it may be consistent to interpret the BP image as the slip acceleration on the fault. In practice, however, the Earth has a free surface and is neither homogeneous nor completely elastic. As shown in eqs (28) and (30), the Green's function in HBP is approximated to the delta function in BP. However, it is difficult to say which Green's functions, $G_{kl}^d$ or $G_{kl}^v$, is approximated to the delta function. There can be other candidates. The BP image seems to have physically intrinsic ambiguity: we never know the BP image approximately represents slip velocity or slip acceleration on the fault.

## 5 CONCLUSIONS

Although physically intrinsic ambiguity consists in the BP image, which may represent slip velocity or slip acceleration, both the BP and HBP images intend to represent slip motion on the fault under the following two assumptions: (1) waveforms due to fault slips other than the target source grid l are cancelled out each other and (2) the stacked Green's function (for BP) or the stacked autocorrelation function of the Green's function (for HBP) is like the delta function. These conditions are mathematically expressed as follows:

$$\sum_k c_k G_{kl}(t + t_{kl}^p) \approx \delta(t)\delta_{\ell l} \quad (\text{BP}) \quad (48)$$

$$\sum_k c_k' (G_{kl} \hat{\times} G_{kl})(t) \approx \delta(t)\delta_{\ell l} \quad (\text{HBP}) \quad (49)$$

where $\delta_{\ell l}$ is the Kronecker delta. Because the autocorrelation function always has its maximum at $t = 0$, we can expect that HBP
works better than BP, if we can reasonably estimate the Green's function.

In estimating the slip motion on the fault, the following expression is used for each method:

$$s^{\text{BP}} = \mathbf{B}^{\text{T}} \mathbf{C} \mathbf{d} \quad (\text{BP}), \qquad (50)$$

$$s^{\text{HBP}} = \mathbf{G}^{\text{T}} \mathbf{C}' \mathbf{d} \quad (\text{HBP}), \qquad (51)$$

$$\hat{\mathbf{a}} = (\mathbf{G}^{\text{T}} \mathbf{E}^{-1} \mathbf{G} + \varepsilon^2 \mathbf{I})^{-1} \mathbf{G}^{\text{T}} \mathbf{E}^{-1} \mathbf{d} \quad (\text{damped LSS}). \qquad (52)$$

Because of that, the model resolution matrices for BP ($\mathbf{B}^{\text{T}}\mathbf{C}\mathbf{G}$) and HBP ($\mathbf{G}^{\text{T}}\mathbf{C}\mathbf{G}$) are not close to the identity matrix, which results in low resolution of the images (Figs 2 and 3). The conditions that the model resolution matrix is approximated to the identity matrix precisely echo the assumptions made in the analyses of BP (eq. 48) and HBP (eq. 49).

The HBP image corresponds to the damped LSS with an extremely large damping parameter. The Green's function in HBP is further approximated by the delta function in BP (eqs 50 and 51). In terms of the model resolution matrix, the damped least-squares method is clearly superior to the BP and HBP methods.

On the other hand, the damped least-squares method requires much larger computation, which gives the limitation of resolution. In addition to this, information on Green's functions and fault plane is needed. In other words, these points are advantages of the BP analysis.

## ACKNOWLEDGEMENTS

We thank David Al-Attar for useful information. We also thank Toshiro Tanimoto and an anonymous reviewer for their useful comments.

## REFERENCES

Aki, K., Christoffersson, A. & Husebye, E.S., 1977. Determination of the three-dimensional seismic structure of the lithosphere, *J. geophys. Res.*, **82**, 277–296.

Claerbout, J.F., 2001. *Basic Earth Imaging (Version 2.4)*, Stanford University.

Hartzell, S. & Heaton, T.H., 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California earthquake, *Bull. seism. Soc. Am.*, **73**, 1553–1583.

Ishii, M., Shearer, P.M., Houston, H. & Vidale, J.E., 2005. Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the Hi-Net array, *Nature*, **435**, 933–936.

Ishii, M., Shearer, P.M., Houston, H. & Vidale, J.E., 2007. Teleseismic P wave imaging of the 26 December 2004 Sumatra-Andaman and 28 March 2005 Sumatra earthquake ruptures using the Hi-net array, *J. geophys. Res.*, **112**, B11307, doi:10.1029/2006JB004700.

Jackson, D.D., 1979. The use of a priori data to resolve non-uniqueness in linear inversion, *Geophys. J. R. astr. Soc.*, **57**, 137–157.

Ji, C., Wald, D.J. & Helmberger, D.V., 2002. Source description of the 1999 Hector Mine, California, Earthquake, part II: complexity of slip history, *Bull. seism. Soc. Am.*, **92**, 1208–1226.

Kawakatsu, H. & Montagner, J.P., 2008. Time-reversal seismic-source imaging and moment-tensor inversion, *Geophys. J. Int.*, **175**, 686–688.

Kikuchi, M. & Kanamori, H., 1991. Inversion of complex body wave-III, *Bull. seism. Soc. Am.*, **81**, 2335–2350.

Kiser, E. & Ishii, M., 2011. The 2010 Mw 8.8 Chile earthquake: triggering on multiple segments and frequency-dependent rupture behavior, *Geophys. Res. Lett.*, **38**, L07301, doi:10.1029/2011GL047140.

Kiser, E., Ishii, M., Langmuir, C.H., Shearer, P.M. & Hirose, H., 2011. Insights into the mechanism of intermediate-depth earthquakes from source properties as imaged by back projection of multiple seismic phases, *J. geophys. Res.*, **116**, B06310, doi:10.1029/2010JB007831.

Koper, K.D., Hutko, A.R., Lay, T., Ammon, C.J. & Kanamori, H., 2011. Frequency-dependent rupture process of the 2011 Mw 9.0 Tohoku Earthquake: comparison of short-period P wave backprojection images and broadband seismic rupture models, *Earth Planets Space*, **63**, 599–602.

Larmat, C., Montagner, J.P., Fink, M., Capdeville, Y., Tourin, A. & Clévédé, E., 2006. Time-reversal imaging of seismic sources and application to the great Sumatra earthquake, *Geophys. Res. Lett.*, **33**, L19312, doi:10.1029/2006GL026336.

Matsu'ura, M., Noda, A. & Fukahata, Y., 2007. Geodetic data inversion based on Bayesian formulation with direct and indirect prior information, *Geophys. J. Int.*, **171**, 1342–1351.

Meng, L., Ampuero, J.P., Luo, Y., Wu, W. & Ni, S., 2012. Mitigating artifacts in back-projection source imaging with implications for frequency-dependent properties of the Tohoku-Oki earthquake, *Earth Planets Space*, **64**, 1101–1109.

Menke, W., 2012. *Geophysical Data Analysis: Discrete Inverse Theory*, 3rd edn, Elsevier.

Olson, A.H. & Apsel, R.J., 1982. Finite faults and inverse theory with applications to the 1979 Imperial Valley earthquake, *Bull. seism. Soc. Am.*, **72**, 1969–2001.

Piatanesi, A., Cirella, A., Spudich, P. & Cocco, M., 2007. A global search inversion for earthquake kinematic rupture history: application to the 2000 Western Tottori, Japan earthquake, *J. geophys. Res.*, **112**, B07314, doi:10.1029/2006JB004821.

Suzuki, M. & Yagi, Y., 2011. Depth dependence of rupture velocity in deep earthquakes, *Geophys. Res. Lett.*, **38**, L05308, doi:10.1029/2011GL046807.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, SIAM.

Trifunac, M.D., 1974. A three-dimensional dislocation model for the San Fernando, California, earthquake of February 9, 1971, *Bull. seism. Soc. Am.*, **64**, 149–172.

Tromp, J., Tape, C. & Liu, Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**, 195–216.

Walker, K.T., Ishii, M. & Shearer, P.M., 2005. Rupture details of the 28 March 2005 Sumatra Mw 8.6 earthquake imaged with teleseismic P waves, *Geophys. Res. Lett.*, **32**, L24303, doi:10.1029/2005GL024395.

Wang, D., Mori, J. & Uchide, T., 2012. Supershear rupture on multiple faults for the Mw 8.6 Off Northern Sumatra, Indonesia earthquake of April 11, 2012, *Geophys. Res. Lett.*, **39**, L21307, doi:10.1029/2012GL053622.

Xu, Y., Koper, K.D., Sufri, O., Zhu, L. & Hutko, A., 2009. Rupture imaging of the Mw 7.9 12 May 2008 Wenchuan earthquake from back projection of teleseismic P waves, *Geochem. Geophys. Geosyst.*, **10**, Q04006, doi:10.1029/2008GC002335.

Yagi, Y. & Fukahata, Y., 2011a. Introduction of uncertainty of Green's function into waveform inversion for seismic source processes, *Geophys. J. Int.*, **186**, 711–720.

Yagi, Y. & Fukahata, Y., 2011b. Rupture process of the 2011 Tohoku-oki earthquake and absolute elastic strain release, *Geophys. Res. Lett.*, **38**, L19307, doi:10.1029/2011GL048701.

Yagi, Y., Nakao, A. & Kasahara, A., 2012. Smooth and rapid slip near the Japan Trench during the 2011 Tohoku-oki earthquake revealed by a hybrid back-projection method, *Earth planet. Sci. Lett.*, **355–356**, 94–101.

Yao, H., Shearer, P.M. & Gerstoft, P., 2012. Subevent location and rupture imaging using iterative backprojection for the 2011 Tohoku Mw 9.0 earthquake, *Geophys. J. Int.*, **190**, 1152–1168.