Inversion procedure

2.3.1 Linearisation of the discretised forward problem

To solve the forward problem efficiently, the leapfrog formulas (Eqs 2.57 and 2.58) are linearised with respect to the permittivity entryc_j at the permittivity distribu-tionϵas:

p_ℓ+1

∂c_j =∂p_ℓ

∂c_j +∆tC⁻¹

‚

h^diff_ℓ +∂b_ℓ+¹

2

∂c_j

Œ

(2.70)

∂q^(k)

ℓ+¹₂

∂c_j =

∂q^(k)

ℓ−¹₂

∂c_j +∆t

∂a^(k)

ℓ−¹₂

∂c_j , (2.71)

whereh^diff_ℓ is an auxiliary source function originating from the scattering point and defined by

h^diff_ℓ =−∂C

∂c_jC⁻¹(f_ℓ+b_ℓ+¹

2). (2.72)

The more detailed derivation of this linearised form is provided in [64]and in Publication III.

Based on the definition of the permittivity distribution (Eq. 2.69) and the entries of the matrixC(Eq. 2.50),

‚∂C

∂c_j

Œ

k,l

=

∫︂

T_j

ϕ_kϕ_ldΩ. (2.73)

This means that the differential entry(∂C/∂c_j)_k,l is non-zero only whereϕ_k and ϕ_l are supported onT_j, the elements of the fine meshT. The differential of the equation 2.73 can therefore be interpreted as a special case of the first-order BA (Eq.

2.67). Also, the multigrid formulation of the permittivity distribution is valid be-cause the differential update(∂C/∂c_j)to the permittivity distribution differs from zero only at the coarse mesh elementT_j^′∈ T^′.

In the reconstruction procedure, the relation between the permittivity vector x= (c₁,c₂, . . . ,c_M) and the discretised electric potential fieldpis approximated by

the linearised forward model

p[x] =p[x₀] +J[x₀](x−x₀), (2.74) in whichx₀is ana prioriestimate for the permittivity, andJ[x₀]the Jacobian matrix consisting of partial derivatives∂p_ℓ/∂c_j evaluated atx₀and formulated as

∂p_ℓ

∂c_j = ∑︂

⃗r_i∈T_j,i≤n

d^(i,_ℓ ^j). (2.75)

The entriesd_ℓ^i,j can be computed by solving an auxiliary system which follows from the leapfrog time integration system (Eqs 2.57 and 2.58) by placingh^diff_ℓ (Eq.

2.72) as the source. This procedure is formulated in[75]. The number of the terms in the sum of the equation 2.75 depends on the density of the fine meshT, relating to the number of nodes⃗r_i∈ T. To reduce the number of terms, and the computa-tional work requirement in simulating the model, the multigrid approach is used to redefine the sum 2.75 with respect to the coarse meshT^′as

∂p_ℓ

∂c_j ≈ ∑︂

⃗ri∈T_j^′,i≤N

d^′(_ℓ^i,j)=∑︂^d+1

k=1

d^′(_ℓ^ik,j), (2.76) where the number of termsd+1 is the number of nodesr⃗_ibelonging to the coarse mesh elementT_j^′

2.3.2 Deconvolution regularisation parameter

The selection of the contrast factorc(ρ)in the equations 2.42 and 2.46 is related to the permittivity perturbationρat the point⃗. As the source termr h˜^BA,n is linear with respect to the perturbation, and the BA is found through the Tikhonov-regularised deconvolution process, the selection of the magnitude of the permittivity perturba-tionρcan be associated with choosing an appropriate Tikhonov deconvolution reg-ularisation parameterν(Publication III). Hence, an update of the formρ→γρwith γ >0 corresponds toh^BA,n→γh^BA,n. If the estimate top_ℓis updated asp_ℓ→γp_ℓ, the corresponding deconvolution regularisation parameter is updated asν→γ⁻¹ν.

This means that the same effect which follows from a decrease in the perturbation

can be achieved by an increase in the regularisation parameterν. Therefore, in the numerical evaluation of the BA, any perturbation can be assumed to beρ=1, and the deconvolution regularisation parameter can be selected with respect to that.

2.3.3 Total variation

In the inversion stage, the linearised forward model (Eq. 2.74) can be written as the linear system

y−y₀=L(x−x₀) +n, (2.77)

whereyis the actual data vector relating to the perturbed permittivity distribution, y₀relates to the data on the background permittivity distributionϵ_{b g},Lis the Jaco-bian matrix for the constanta prioriguessx₀, andnis a noise vector including both modelling and measurement errors. The regularised solution forxis then obtained by the iteration

x_ℓ+1=x₀+ (L^TL+αDΓ_ℓD)⁻¹L^T(y−y₀), (2.78) in which Γ_ℓ is a weighting matrix satisfyingΓ₀ = I andΓ_ℓ = diag(|D[x_ℓ−x₀]|)⁻¹ forℓ≥1. The constantαis a regularisation parameter which value is commonly adjusted to the same order of magnitude as||y−y₀||². The matrixDis a regularisation matrix of the form

D_i,j =βδ_i,j+ ℓ^(i,j)

max_(i,j)ℓ^(i,j)(2δ_i,j−1), (2.79) δ_i,j =

⎧

⎨

⎩

1 if j=i

0 otherwise. (2.80)

The first term of the equation 2.78 is a weighted identity operator which limits the total magnitude ofx, and the second term penalises the jumps ofxover the edges of the coarse meshT^′multiplied by the edge lengthℓ^(i,j)=∫︁

T_i^′∩T_j^′d s. If the regulari-sation constantβ=0, the regularisation procedure yields the total variation ofx. If β >0 then also the norm ofxis regularised at each iteration step.

Following the detailed derivation in[64], the inversion process minimises the regularising function

F(x) =||L(x−x₀)−(y−y₀)||²₂+2⎷

α||D(x−x₀)||₁, (2.81) which is the sum of the total variation of the estimate and the L2-norm penalty term.

The former of these corresponds to the Euclidean norm of the permittivity gradient integrated overD, while the latter penalises the total magnitude of the distribution.

The regularisation parameterαaffects the reconstruction quality, andβensures the numerical stability of the inversion process bounding the reconstruction values. A characteristic of the total variation (TV) regularisation procedure is to yield a low-noise reconstruction with large connected areas close to constant, because the length of the boundary curve between jumps is regularised. This is especially useful in the context where clear-contrast inclusions are sought from a constant background. In the present context, such inclusions would be internal voids, cracks, other higher-contrast details, and the surface layer of the asteroid target.

The multigrid approach[15, 75] can also be used in the inversion stage to re-construct coarse details before the finer ones. The coarse-to-fine inversion routine alternates between finding the coarse resolution inverse estimate to x^coarse_ℓ+1 and the fine resolution inverse estimatex^fine_ℓ+1, in which the variables in the equation 2.78 are computed for both the resolutions separately. The inverse estimate is the sum of the coarse and the fine solutionsx_ℓ+1 =x^coarse_ℓ+1 +x^fine_ℓ+1. The validity and details of the approach are explained in[75].

2.3.4 Tomographic backprojection

In tomographic reconstruction, a signal measured in a projection of the target is re-distributed to the image space. Backprojection is a central concept in this procedure.

In the presented problem of full-wave tomography, the Born approximation defined in the section 2.2.2 can be interpreted as a differential of the full wave signal with respect to the perturbed permittivity distribution ϵ(Eq. 2.41). The resulting BA matrix can be associated with an array of the formL= [d₁,d₂, . . . ,d_N], where each columnd_j denotes the BA evaluated at a perturbation pointr_j, as also specified by the Eq. 2.75. Consequently, by defining the perturbed permittivity distribution as

∆c= [∆c₁,∆c₂, . . . ,∆c_N]^T, the signal perturbation ∆s corresponding to the

per-mittivity perturbation can be approximated numerically by the product

∆s=L∆c. (2.82)

A rough backpropagated reconstruction can be obtained by interpreting the ad-joint operation of BA as the tomographic backprojection. The multiplicationL^T∆s is obtained from the equation

∆s^T(L∆c) =∆c^T(L^T∆s). (2.83) The nonzero entries of the resulting backpropagation matrix correspond to the lo-cations at which a permittivity perturbation contributes to the signal∆sat the spec-ified time interval.

In document Full-Wave Radar Tomography of Complex, High-Contrast Targets: Imaging the Interior Structure of a Small Solar System Body (sivua 42-46)