Inverse Models

Size Distribution Inversion

Theory

The size distribution can be found from the inverse of the response function using Tikhonov regularization.

$L_1$ and $L_2$ are the Euclidian norms of the solution and the intial guess.

$L_1 = \left\lVert \bf{A}\rm{𝕟} - \rm{𝕣}\right\rVert_2$

$L_2 = \left\lVert\bf{L}(\rm{𝕟} - \rm{𝕟_i})\right\rVert_2$

where $\rm{𝕟}$ is the true size distribution, $\bf{A}$ is the convolution matrix, 𝕣 is the measured response vector, $\bf{L}$ is a weights matrix, and $\rm{𝕟_i}$ is an initial guess. To solve for the inverse, a balance is sought between the $L_1$ and $L_2$ norms.

$𝕟_{inv} = \arg \min\{L_1^2 + \lambda^2 L_2^2\}$

where $\lambda$ is the regularization parameter. Taking $\bf{L} = \bf{I}$ (weight matrix equals the identity Matrix 𝐈) and $\rm{𝕟_i} = \bf{S}^{-1}\rm{𝕣}$ as initial guess (Matrix 𝐒), the regularized inverse is computed via:

$𝕟_{inv} = (\bf{A}^{\rm{T}}\bf{A} + \lambda^{\rm{2}} \bf{I})^{\rm{-1}}(\bf{A}^{\rm{T}} \rm{𝕣} + \lambda^2\bf{S}^{-1} \rm{𝕣})$

The regularization parameter $\lambda$ "interpolates" between the noisy least-square inverse ($\lambda = 0$) and the initial guess ($\lambda >> 1$). The optimal regularization parameter is found using the L-curve method. The optimal $\lambda_{opt}$ is found using the L-curve method. The L-curve is defined as a plot of $\log(L_1)$ vs. $\log(L_2)$, where $L_1$ and $L_2$ are obtained for a series of discrete $\lambda$ values. The $\lambda_{opt}$ is found at the corner of the L-curve, which is mathematically defined as the point where where the curvature of the L-curve is maximum. Here, the corner of the L-curve is found using the iterative algorithm described in Talukdar and Swihart (2003). The curvature is calculated using Eq. (14) in Hansen (2000), which requires the first and second derivatives of $d\ln(L_i)^2/d\lambda$. These derivatives of $d\ln(L_i)^2/d\lambda$ are estimated numerically.

See session 3 of the Tutorial for an interactive explanation on how Tikhonov regularization works.

Example

The L-curve search method is implemented in the rinv function. Also check out the new rinv2 function, which has much faster performance and is based on generalized cross validation. It also allows to specify higher order inversions. The default implementation gives near identical results to rinv.

𝕟ⁱⁿᵛ = rinv(𝕣.N, δ, λ₁=0.1, λ₂=1.0);

The function takes the response vector to be inverted, the DMA grid to read the matrices, and an upper and lower search bound for the optimal regularization parameter.

df = CSV.read("example_data.csv", DataFrame)

t, p, lpm = 293.15, 940e2, 1.666e-5
r₁, r₂, l = 9.37e-3,1.961e-2,0.44369
Λ = DMAconfig(t,p,1lpm,4lpm,r₁,r₂,l,0.0,:+,6,:cylindrical)
δ = setupDMA(Λ, vtoz(Λ,10000), vtoz(Λ,10), 120)

𝕣 = (df,:Dp,:Rcn,δ) |> interpolateDataFrameOntoδ
𝕟ⁱⁿᵛ = rinv(𝕣.N, δ, λ₁=0.1, λ₂=1.0)
𝕟ⁱⁿᵛ² = rinv2(𝕣.N, δ, λ₁=0.1, λ₂=1.0)

The file example_data.csv is a comma delimited text files that contains diameter and a resonse vector. The function interpolateDataFrameOntoδ takes the columns :Dp and :Rcn from the DataFrame df and creates a response distribution of the type SizeDistribution that matches the DMA grid δ. Finally, 𝕟ⁱⁿᵛ is computed using the rinv function.

The method was compared to the inverted size distribution given by the manufacturer Aerosol Instrument Manager (AIM) software, assuming no diffusion correction in either method. Figure 3 in the manuscript shows excellent correspondence between this package and the commercial software.

CCN Inversion

Size resolved CCN measurements have a long history (e.g. Cruz and Pandis, 1997, Snider et al., 2006). Particles are typically dried, charge neutralized, and passed through the DMA.At the exit, the flow is split between a CPC that measures all particles and a CCN counter that measures particles that form cloud droplets at a water supersaturation set by the instrument. In this configuration the DMA voltage can either be stepped or continuously scanned. The ratio of CCN-to-CPC response function is used to determine the fraction of particles that activate at a given diameter. The diameter where 50% of the particles activate is taken to be the activation diameter. The activation of particles into cloud droplets is proportional to the volume of solute in the particle. Therefore, larger multiply charged particles activate first. This leads to a bias in the inferred D50 diameter if the activated fraction is computed from the ratio of the raw response functions (Petters et al., 2007, 2009).

The CCN transmission function is modeled using a cumulative Gauss integral

$T_{af} = \frac{1}{2}\left[1 + \mathrm{erf}\left(\frac{d-\mu}{\sigma}\right) \right]$

This function is applied to the mobility size distribution. Then the response function is computed. Empirically, activated fraction can be computed using from the ratio of size distributions and response functions.

To predict the ratio of 𝕣ᶜᶜⁿ/𝕣ᶜⁿ in terms of the activated fraction model Taf:

$𝕟^{\mathrm{cn}} = (\bf{A}^{\rm{T}}\bf{A} + \lambda^{\rm{2}} \bf{I})^{\rm{-1}}(\bf{A}^{\rm{T}} \rm{𝕣^{\mathrm{cn}}} + \lambda^2 \bf{S}^{-1} \rm{𝕣^{\mathrm{cn}}})$

$\left( \frac{𝕣^{\mathrm{ccn}}}{𝕣^{\mathrm{cn}}} \right)_{\mathrm{model}} = \frac{\mathbf{A}[T_{af}\; .*\; 𝕟^{\mathrm{cn}}]}{\mathbf{A}𝕟^{\mathrm{cn}}}$

This mirrors the approach in Petters et al. (2007). The model + fit can be solved using this package as follows.

𝕟ᶜⁿ = (𝐀'𝐀 + λ^2𝐈)^(-1) * (𝐀'𝕣ᶜⁿ  + λ^2 * 𝐒^(-1)*𝕣ᶜⁿ)
model(x,p) = (δ.𝐀*(𝕟ᶜⁿ.N.*Taf(𝕟ᶜⁿ, p[1], p[2])))./(δ.𝐀*𝕟ᶜⁿ.N)
fit = curve_fit(model, δ.Dp, 𝕣ᶜᶜⁿ.N./𝕣ᶜⁿ.N, [65.0, 3.0])

The code loads the CCN (𝕣ᶜᶜⁿ) and CN (𝕣ᶜⁿ) response functions from file. The activated fraction is computed through the ratio of the distributions using one of the SizeDistribution Number Operators.

𝕒𝕗 = 𝕣ᶜᶜⁿ/𝕣ᶜⁿ

And here is the complete example showing how to fit the 𝕒𝕗 response curve.

t, p, lpm = 293.15, 940e2, 1.666e-5
r₁, r₂, l = 9.37e-3,1.961e-2,0.44369
Λ = DMAconfig(t,p,1lpm,4lpm,r₁,r₂,l,0.0,:+,6,:cylindrical)
δ = setupDMA(Λ, vtoz(Λ,10000), vtoz(Λ,10), 120)

# Load a simple comma delimited text file - file contains :Dp, :Rcn, :Rccn
df = CSV.read("example_data.csv", DataFrame)
𝕣ᶜⁿ = (df,:Dp,:Rcn,δ) |> interpolateDataFrameOntoδ        # CN response distribution
𝕣ᶜᶜⁿ = (df,:Dp,:Rccn,δ) |> interpolateDataFrameOntoδ;     # CCN response distribution
𝕒𝕗 = 𝕣ᶜᶜⁿ/𝕣ᶜⁿ
Taf(𝕟,μ,σ) = @. 0.5 * (1.0 + erf((𝕟.Dp - μ)./(sqrt(2.0σ))))
𝐈, 𝐒, 𝐀, λ =  δ.𝐈, δ.𝐒, δ.𝐀, 0.5
𝕟ᶜⁿ = (𝐀'𝐀 + λ^2𝐈)^(-1) * (𝐀'𝕣ᶜⁿ  + λ^2 * 𝐒^(-1)*𝕣ᶜⁿ)
model(x,p) = (𝐀 * (𝕟ᶜⁿ.N .* Taf(𝕟ᶜⁿ, p[1], p[2])))./( 𝐀 * 𝕟ᶜⁿ.N)
fit = curve_fit(model, 𝕒𝕗.Dp, 𝕒𝕗.N, [65.0, 3.0])
Ax = fit.param
afmodel = model(δ.Dp, Ax)

Tandem DMA Inversion

A complete example for tandem DMA inversion is provided in inversion3.jl in the examples folder, which reproduces Figure 3 in Petters (2021).

Briefly, the code for inversion contains these key features:

Λ₁, Λ₂, δ₁, δ₂ = initializeDMAs(Dd, k)
Ax = [[1300.0, 60.0, 1.4], [2000.0, 200.0, 1.6]]
𝕟ᶜⁿ = DMALognormalDistribution(Ax, δ₁)
gf, ge, 𝐀 = TDMAmatrix(𝕟ᶜⁿ, Dd, Λ₁, Λ₂, δ₂, k)
model = TDMA1Dpdf(𝕟ᶜⁿ, Λ₁, Λ₂, (Dd, 0.8, 5.0, k));

The initializeDMAs sets up the grids for DMA 1 and 2, Ax are the parameters for the assumed size distribution, 𝕟ᶜⁿ is the instantiation of that size distribution. gf, ge, and 𝐀 are the growth factor bin midpoints, growth factor bin edges and the convolution matrix, and model is the forward model function.

The code below is to invert the data. x₀ is the initial guess, lb,ub the lower and ope bounds, and xλ1 and xλ2 are the inverted solutions for the two main methods stated in the manuscript.

x₀ = N1./sum(N1)
lb, ub = zeros(k), ones(k)
xλ1 = invert(𝐀, N1, Lₖx₀B(2, x₀, lb, ub))
e1 = @> sqrt.(sum((xλ1 .- f).^2.0)./k) round(digits = 3)
xλ2 = invert(𝐀, N1, LₖDₓB(0, 0.001, lb, ub))
e2 = @> sqrt.(sum((xλ2 .- f).^2.0)./k) round(digits = 3)