% B.I.Stepanov Institute of Physics % National Academy of Sciences of Belarus % Aerosol profile retrieval algorithm utilizing combined lidar and CIMEL % Sun photometer data. function liric global J M; global k f d; global b a V; global hStep Nfin N; global betaMref Omega betaM Lmeas; global tauAMatrix smoothOp; global iter; global CArr RArr Psi1Arr Psi2Arr Psi3Arr LmeasRArr LcalcArr; %------------------------------------------------------------------------- disp('Reading input data...'); % Number of lidar channels. J = hdf5read('data.h5', '/input/J'); % Number of aerosol modes. M = hdf5read('data.h5', '/input/M'); % Lidar signal weighting coefficient (by channel). k = hdf5read('data.h5', '/input/k'); % Aerosol column volume weighting coefficient (by aerosol mode). f = hdf5read('data.h5', '/input/f'); % Profile smoothness weighting coefficient (by aerosol mode). d = hdf5read('data.h5', '/input/d'); % Aerosol backscatter coefficient (by channel, aerosol mode). b = hdf5read('data.h5', '/input/b'); % Aerosol extinction coefficient (by channel, aerosol mode). a = hdf5read('data.h5', '/input/a'); % Aerosol column volume, in micrometers (by aerosol mode). V = hdf5read('data.h5', '/input/V'); % Height step of the lidar data grid, in meters. hStep = hdf5read('data.h5', '/input/hStep'); % First index of lidar signal data selected for retrieval (by channel). % Add 1 to obtain Matlab's 1-based indices. Nbeg = hdf5read('data.h5', '/input/Nbeg') + 1; % Reference point index (by channel). Nref = hdf5read('data.h5', '/input/Nref') + 1; % Last index of lidar signal data selected for retrieval (by channel). Nfin = hdf5read('data.h5', '/input/Nfin') + 1; % Molecular backscatter coefficient at the reference point (by channel). betaMref = hdf5read('data.h5', '/input/betaMref'); % Measured lidar signal, divided by its average value around the reference % point (by height, channel). S = hdf5read('data.h5', '/input/S'); % Weighting coefficients for lidar signal (by height, channel). Omega = hdf5read('data.h5', '/input/Omega'); % Molecular backscatter coefficient (by height, channel). betaM = hdf5read('data.h5', '/input/betaM'); % Molecular optical thickness, relative to the reference point, with % positive values below it (by height, channel). tauM = hdf5read('data.h5', '/input/tauM'); % Zenith angle of the lidar track, in degrees (by height, channel). Z = hdf5read('data.h5', '/input/Z'); % Initial approximation to aerosol concentration, in ppb (by height, % aerosol mode). C0 = hdf5read('data.h5', '/input/C0'); % Initial approximation to backscatter ratio at reference point (by channel). R0 = hdf5read('data.h5', '/input/R0'); TolFun = hdf5read('data.h5', '/input/TolFun'); TolX = hdf5read('data.h5', '/input/TolX'); % Size of the aerosol concentration profile to be retrieved. N = max(Nfin); %------------------------------------------------------------------------- % Check consistency of input array dimensions. assert(all(size(J) == [1, 1]) && J > 0); assert(all(size(M) == [1, 1]) && M > 0); assert(all(size(k) == [J, 1])); assert(all(size(f) == [M, 1])); assert(all(size(d) == [M, 1])); assert(all(size(b) == [J, M])); assert(all(size(a) == [J, M])); assert(all(size(V) == [M, 1])); assert(all(size(hStep) == [1, 1]) && hStep > 0); assert(all(size(Nbeg) == [J, 1])); assert(all(size(Nref) == [J, 1])); assert(all(size(Nfin) == [J, 1])); assert(all(size(betaMref) == [J, 1])); assert(all(size(S) == [N, J])); assert(all(size(Omega) == [N, J])); assert(all(size(betaM) == [N, J])); assert(all(size(tauM) == [N, J])); assert(all(size(Z) == [N, J])); assert(all(size(C0) == [N, M])); assert(all(size(R0) == [J, 1])); assert(all(size(TolFun) == [1, 1]) && TolFun > 0); assert(all(size(TolX) == [1, 1]) && TolX > 0); %------------------------------------------------------------------------- % Measured lidar signal, divided by its average value around the reference % point, with molecular component removed, without the backscatter ratio % correction factor (by height, channel). Lmeas = zeros(N, J); for j = 1:J Lmeas(Nbeg(j):Nfin(j), j) = S(Nbeg(j):Nfin(j), j) .* ... exp(-2 * tauM(Nbeg(j):Nfin(j), j)); assert(all(Omega(Nbeg(j):Nfin(j)) > 0)) % Extend arrays related to lidar signal equations below the lower % signal boundaries. Thus signal equations in the near-to-the-ground % region will be the same, and retrieved concentrations in that region % will tend to be near-constant. Lmeas(1:Nbeg(j)-1, j) = Lmeas(Nbeg(j), j); Omega(1:Nbeg(j)-1, j) = Omega(Nbeg(j), j); betaM(1:Nbeg(j)-1, j) = betaM(Nbeg(j), j); Z(1:Nbeg(j)-1, j) = Z(Nbeg(j), j); end % Matrices that transform extinction coefficient vectors to optical % thickness vectors, relative to the reference points (by channel). tauAMatrix = cell(J, 1); for j = 1:J % Steps along the lidar track, with zenith angles taken into account. effStep = hStep ./ cos(degtorad(Z(1:Nfin(j), j))); tauAMatrix{j} = tauIntegrationMatrix(... Nbeg(j), Nref(j), Nfin(j), effStep); end % Matrix that transforms a concentration profile to a vector evaluating % its somoothness. smoothOp = secondDifferenceOperator(N, hStep); %------------------------------------------------------------------------- % Current iteration number + 1 (1 means no iterations have completed yet). iter = 0; % The following arrays will be used to store snapshots of the optimization % process at every iteration as well as the retrieval results. They will % be 3-dimensional; new layers will be added to these arrays in 'optfun'. % Preallocate space for 2 iterations to make sure that the arrays will remain % 3-dimensional even in the non-realistic case where initial approximation % itself triggers termination of the retrieval process. % Aerosol concentrations. CArr = zeros(N, M, 2); % Backscatter ratios. RArr = zeros(1, J, 2); % Residuals for lidar equations. Psi1Arr = zeros(1, J, 2); % Residuals for photometer equations. Psi2Arr = zeros(1, M, 2); % Residuals for smoothness criteria. Psi3Arr = zeros(1, M, 2); % Measured lidar signal, divided by its average value around the reference % point, with molecular component removed (by height, channel). LmeasRArr = zeros(N, J, 2); % Calculated lidar signal, divided by its value at the reference point, % with molecular component removed (by height, channel). LcalcArr = zeros(N, J, 2); %------------------------------------------------------------------------- disp('Starting optimization...'); x0 = [reshape(C0, N*M, 1); R0]; l0 = zeros(size(x0)); optopt = optimset('Display', 'Iter', 'DerivativeCheck', 'off', ... 'TolFun', TolFun, 'TolX', TolX, 'Jacobian', 'on'); % Invoke optimizer. lsqnonlin(@optfun, x0, l0, [], optopt); %------------------------------------------------------------------------- disp('Saving output data...') hdf5write('data.h5', '/output/iterCount', iter, 'WriteMode', 'append'); hdf5write('data.h5', '/output/C', CArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/R', RArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Psi1', Psi1Arr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Psi2', Psi2Arr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Psi3', Psi3Arr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/LmeasR', LmeasRArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Lcalc', LcalcArr, 'WriteMode', 'append'); disp('Done.') end %------------------------------------------------------------------------- % Optical thickness integration operator. % % Returns a matrix that transforms a vector of extinction coefficients % defined on a grid with the given variable step to the vector of optical % thicknesses, relative to the reference point 'Nref' (with positive % values below it). 'Nfin' defines the size of the matrix. Optical % thicknesses below the 'Nbeg' point will be the same and equal to that at % 'Nbeg'. function tauMatrix = tauIntegrationMatrix(Nbeg, Nref, Nfin, stepVector) % Assume that 'stepVector(i)' is the distance between nodes i and i+1. % Rectangular rule integration matrix, using left corner function % values. This matrix will produce optical thicknesses relative to the % ground. Using a sparse matrix speeds up the calculation a bit. direct = tril(ones(Nfin), -1) * sparse(1:Nfin, 1:Nfin, stepVector); % Adjust optical thicknesses below 'Nbeg'. direct(1:Nbeg-1, :) = repmat(direct(Nbeg, :), Nbeg-1, 1); % Vector to produce optical thickness at the reference point relative % to the ground, using the same integration rule as in 'direct'. reverseVector = stepVector'; if Nref <= Nfin reverseVector(Nref:Nfin) = 0.0; end % To obtain the relative optical thickness, subtract the absolute one % from that in the reference point. tauMatrix = repmat(reverseVector, Nfin, 1) - direct; end % ------------------------------------------------------------------------ % Second order finite difference operator for concentration profiles. % % Returns a matrix that transforms a function defined on a grid with the % given step to an approximation to that function's second derivative % defined on the same grid. function diffOp = secondDifferenceOperator(N, step) % A matrix with '2's on its main diagonal and '-1's on the diagonals % just below and just above the main one. % % Using sparse matrices speeds up the calculation a bit. diffOp = sparse(1:N, 1:N, 2) + ... sparse(2:N, 1:N-1, -1, N, N) + sparse(1:N-1, 2:N, -1, N, N); % Assume that the function's first order derivatives at the segment % boundaries are zeros. diffOp(1,1) = 1; diffOp(N,N) = 1; diffOp = diffOp / step^2; end %------------------------------------------------------------------------- function [F, Jac] = optfun(x) global J M; global k f d; global b a V; global hStep Nfin N; global betaMref Omega betaM Lmeas; global tauAMatrix smoothOp; global iter; global CArr RArr Psi1Arr Psi2Arr Psi3Arr LmeasRArr LcalcArr; iter = iter + 1; C = reshape(x(1:N*M), N, M); % Current aerosol concentrations. R = x(N*M+1 : N*M+J); % Current backscatter ratios. CArr(1:N, 1:M, iter) = C; RArr(1, 1:J, iter) = R; % Lidar portion of the equations. F1 = zeros(sum(Nfin), 1); Jac1 = zeros(sum(Nfin), N*M + J); for j = 1:J % Equation indices for the current channel. eqnRange = sum(Nfin(1:j-1))+1 : sum(Nfin(1:j)); % Calculated aerosol extinction and backscatter coefficients. sigmaA = 0.0; betaA = 0.0; for i = 1:M sigmaA = sigmaA + a(j, i) * C(1:Nfin(j), i); betaA = betaA + b(j, i) * C(1:Nfin(j), i); end % Calculated aerosol optical thickness, relative to the reference % point. tauA = tauAMatrix{j} * sigmaA; % Measured lidar signal, divided by its average value around the % reference point, with molecular component removed. LmeasR = Lmeas(1:Nfin(j), j) * R(j); % Calculated lidar signal, divided by its value at the reference % point, with molecular component removed. Lcalc = (betaA + betaM(1:Nfin(j), j)) ./ betaMref(j) .* exp(2 * tauA); % Weighting coefficients for lidar equations. coeff = k(j) * sqrt(hStep ./ Omega(1:Nfin(j), j)); F1(eqnRange) = (LmeasR - Lcalc) .* coeff; Psi1Arr(1, j, iter) = norm(F1(eqnRange)); LmeasRArr(1:Nfin(j), j, iter) = LmeasR; LcalcArr(1:Nfin(j), j, iter) = Lcalc; for i = 1:M % Derivatives with respect to aerosol concentrations C(i). % Use 'sparse' for diagonal matrices to improve performance. Jac1(eqnRange, N*(i-1)+1 : N*(i-1)+Nfin(j)) = ... -sparse(1:Nfin(j), 1:Nfin(j), ... b(j, i) * exp(2 * tauA) ./ betaMref(j) .* coeff) - ... sparse(1:Nfin(j), 1:Nfin(j), ... 2 * a(j, i) * Lcalc .* coeff) * tauAMatrix{j}; end % Derivatives with respect to backscatter ratio R(j). Jac1(eqnRange, N*M+j) = Lmeas(1:Nfin(j), j) .* coeff; end % Photometer portion of the equations. F2 = zeros(M, 1); Jac2 = zeros(M, N*M + J); for i = 1:M % Use 0.001 to convert meters * ppb to micrometers. F2(i) = (V(i) - sum(C(:, i)) * 0.001 * hStep) * f(i); Psi2Arr(1, i, iter) = abs(F2(i)); % Derivatives with respect to aerosol concentrations C(i). Jac2(i, N*(i-1)+1 : N*i) = ... -0.001 * hStep * f(i) * ones(1, N); end % Smoothness portion of the equations. % The matrices are sparse here, to improve performance. F3 = zeros(N*M, 1); Jac3Matrix = cell(3, 1); for i = 1:M F3(N*(i-1)+1 : N*i) = d(i) * sqrt(hStep) * smoothOp * C(:, i); Psi3Arr(1, i, iter) = norm(F3(N*(i-1)+1: N*i)); Jac3Matrix{i} = d(i) * sqrt(hStep) * smoothOp; end Jac3 = [blkdiag(Jac3Matrix{1:M}), zeros(N*M, J)]; F = [F1; F2; F3]; Jac = [Jac1; Jac2; Jac3]; end