Is it correct for converting relative perturbation to absolute perturbation?

[YoushanLiu]

On lines 211-213 of workflow/migration.py:
211 # Merge to vector and convert to absolute perturbations:
212 # log dm --> dm (see Eq.13 Tromp et al 2005)
213 gradient.update(vector=gradient.vector * model.vector)

Maybe there is a problem involved in the conversion on line 213. In SPECFEM, after calling finizalize_simulation, the outputs are relative kernels.

In Eq. 13 Tromp et al., 2005, K_relative * delta_(ln rho) = K_relative * delta_rho / rho.
Thus, K_absolute = K_relative / rho, i.e., the relative perturbations can be converted to absolute perturbations divided by the material parameters (rho) instead of multiplied by corresponding parameters.

Please check line 213 in migration.py.

Alternatively, we can update models using the relative perturbations directly by:
m_(k+1) = m_k * exp(alpha * p) because ln(m_(k+1)) = ln(m_k) + alpha * p

[bch0w]

Hi @YoushanLiu, sorry for taking so long to get to this. Looking at the code, and the equation from the paper, it looks to me like multiplication is correct here.

• SPECFEM outputs relative perturbation, e.g., for density: $\delta\ln\rho = \delta \rho / \rho$
• But to update the model we want the absolute perturbation $\delta \rho = (\delta \ln\rho) \rho$
• So we need to multiply the relative perturbation from SPECFEM ($\delta\ln\rho$) by the model vector( $\rho$) to get the absolute perturbation.

Let me know if that makes sense or if I might be missing something here.

---

Below I just reiterate some of your points with references so that the background is captured appropriately as I was trying to convince myself of this. The relevant section in Tromp et al. (2005) on Pg 197 calculates the Frechet derviatives $\delta\chi$

In the code the conversion from relative to absolute kernels occurs here:

seisflows/seisflows/workflow/migration.py

Lines 319 to 323 in 9ddaedc

	# Merge to vector and convert to absolute perturbations:
	# log dm --> dm (see Eq.13 Tromp et al 2005)
	logger.info("scaling gradient to absolute model perturbations")
	gradient.update(vector=gradient.vector * model.vector)
	gradient.write(path=os.path.join(self.path.eval_grad, "gradient"))

Going back to the legacy version of the code this is how things have always been done:

seisflows/seisflows/postprocess/base.py

Lines 66 to 68 in 5c596ed

	# convert to absolute perturbations, log dm --> dm
	# see Eq.13 Tromp et al 2005
	gradient *= solver.merge(solver.load(path +'/'+ 'model'))

In SPECFEM save_adjoint_kernels we are told that the output kernels are given as the relative kernel $\delta \ln \rho$ and not the absolute perturbation ($\delta \rho$)
https://github.com/SPECFEM/specfem2d/blob/98741db1a0c8082ca57364f5b17ea95df6cbf1c2/src/specfem2D/save_adjoint_kernels.f90#L124-L142

          ! for parameterization (rho,mu,kappa): "primary" kernels
          !
          ! kernel expressions below are for relative perturbations
          ! for example using dln(rho) = drho/rho (and not for absolute perturbations drho)
          !
          ! dChi = \int_V [ K_rho dln(rho) + K_mu dln(mu) + K_kappa dln(kappa) ] d^3x
          !
          ! and K_rho   = - rho   int_0^T [ s^adj * \partial_t^2 s ] dt
          !     K_mu    = - 2 mu  int_0^T [ D^adj : D ] dt               (where D = 1/2[grad(s) + grad(s)^T] - 1/3 div(s) I )
          !     K_kappa = - kappa int_0^T [ div(s^adj) div(s) ] dt
          !
          ! (see e.g. Peter et al. 2011, eq. 15, 18, 19)

          ! density kernel
          rho_kl(i,j,ispec) =  - rhol * (deltat * NTSTEP_BETWEEN_COMPUTE_KERNELS) * rho_kl(i,j,ispec)
          ! shear modulus kernel
          mu_kl(i,j,ispec) =  - TWO * mul * (deltat * NTSTEP_BETWEEN_COMPUTE_KERNELS) * mu_kl(i,j,ispec)
          ! bulk modulus kernel
          kappa_kl(i,j,ispec) = - kappal * (deltat * NTSTEP_BETWEEN_COMPUTE_KERNELS) * kappa_kl(i,j,ispec)

So it would make sense that we want to convert from relative to absolute we need to multiply by the model vector

[YoushanLiu]

Hi, @bch0w, Thank you for your reply! Yes, I know that the SPECFEM outputs the relative kernels after calling finalize_simulation.

But, I noticed that your derivation constructed the relationship between relative- and absolute-parameter perturbations. However, we use the sensitivity kernels (actually some forms of their preconditioning, scaling versions) instead of parameter perturbations to update the models. Therefore, we should derive the relationship between relative- and absolute-sensitivity kernels.

For the absolute kernels,
$\delta\chi = \int{[K_{\rho}\delta\rho + K_{\kappa}\delta\kappa + K_{\mu}\delta\mu] d^{3} \textbf{x}}$
$= \int{[\rho K_{\rho}\delta\rho/\rho + \kappa K_{\kappa}\delta\kappa/\kappa + \mu K_{\mu}\delta\mu/\mu] d^{3} \textbf{x}}$
$= \int{[\rho K_{\rho}\delta\ln\rho + \kappa K_{\kappa}\delta\ln\kappa + \mu K_{\mu}\delta\ln\mu] d^{3} \textbf{x}}$
$= \int{[K_{\ln\rho}\delta\ln\rho + K_{\ln\kappa}\delta\ln\kappa + K_{\ln\mu}\delta\ln\mu] d^{3} \textbf{x}}$

We can obtain,
$K_{\ln\rho} = \rho K_{\rho}$
where, $K_{\ln\rho}$ is the relative kernel, while $K_{\rho}$ is the absolute kernel.

Thus,
$K_{\rho} = K_{\ln\rho}/\rho$

As a result, the material parameters should be divided instead of multiplied when converting relative kernels to absolute kernels.

The following references are two relative papers. The first one presented the relative kernels (Tromp et al., 2005), while the second one presented the absolute kernels (Kohn et al., 2012).

References:
Tromp J., Tape C., Liu Q. 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels. Geophys. J. Int., 160, 195-216.
Kohn D., De Nil D., Kurzmann A., Przebindowska A., and Bohlen T. 2012. On the influence of model parameterization in elastic full waveform tomography. Geophys. J. Int., 191, 325-345.

Therefore, the model update and sensitivity kernels are not inconsistent in your code!

As I said before, the division is not stable. You can change your code to update models by using the relative kernels, i.e., $\textbf{m}_{k+1} = \textbf{m}_{k}\cdot\exp(\alpha_{k}\textbf{d}_{k})$

[YoushanLiu]

A simple test can check the correctness of the kernels.
For relative kernels, the magnitudes of kernels for $\rho$, $\kappa$ and $\mu$ are comparable.
However, the magnitudes of kernel of $\rho$ is larger than those of $\kappa$ and $\mu$.
You can check them.

[bch0w]

Hi @YoushanLiu, sorry for the long delay in reply here. We have been looking into this in our group to double check the derivations and math but have not yet come to a conclusion.

[YoushanLiu]

@bch0w Thank you for yor reply! The derivation should be invariable for $\delta\chi$ just as I show above.
I think the current implementation is possibly problematic. For relative kernels, the amplitudes of kernels for $\rho$, $\kappa$ and $\mu$ should be comparable, while the amplitudes of kernels for $\rho$ is larger than those of $\kappa$ and $\mu$ for absolute kernels. You can make a check. Thanks a lot!