% B.I.Stepanov Institute of Physics % National Academy of Sciences of Belarus % Aerosol backscatter and depolarization retrieval algorithm utilizing % polarized lidar signals. function polar global k nu; global hStep Nbeg Nfin; global betaMref gamma Lmeas Ymeas Omega betaM; global beta2Coeff leakage; global tauAMatrix smoothOp2; global iter; global betaArr dArr RArr QArr Psi1Arr Psi2Arr global LmeasRArr LcalcArr YmeasQArr YcalcArr; %------------------------------------------------------------------------- disp('Reading input data...'); % Lidar signal weighting coefficient (by channel). k = hdf5read('data.h5', '/input/k'); % Backscatter smoothness weighting coefficient (by channel). nu = hdf5read('data.h5', '/input/nu'); % Height grid step, in meters. hStep = hdf5read('data.h5', '/input/hStep'); % First index of lidar signal data selected for retrieval. % Add 1 to obtain Matlab's 1-based indices. Nbeg = hdf5read('data.h5', '/input/Nbeg') + 1; % Reference point index. Nref = hdf5read('data.h5', '/input/Nref') + 1; % Last index of lidar signal data selected for retrieval. Nfin = hdf5read('data.h5', '/input/Nfin') + 1; % Molecular backscatter coefficient at the reference point. betaMref = hdf5read('data.h5', '/input/betaMref'); % Aerosol lidar ratio. gamma = hdf5read('data.h5', '/input/gamma'); % Measured lidar signal, divided by its average value at 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). tauM = hdf5read('data.h5', '/input/tauM'); % 1.0 if first lidar channel is unpolarized, and 0.0 if parallel-polarized. beta2Coeff = hdf5read('data.h5', '/input/beta2Coeff'); % Parallel leakage of the cross-polarized lidar channel. leakage = hdf5read('data.h5', '/input/leakage'); % Zenith angle of the lidar track, in degrees (by height). Z = hdf5read('data.h5', '/input/Z'); % Initial approximation to parallel backscatter profile. beta0 = hdf5read('data.h5', '/input/beta0'); % Initial approximation to depolarization ratio profile. d0 = hdf5read('data.h5', '/input/d0'); % Middle value for aerosol backscatter ratio at the reference point. R0 = hdf5read('data.h5', '/input/R0'); % Middle value for the cross-polarized lidar correction coefficient. Q0 = hdf5read('data.h5', '/input/Q0'); % Tolerance value for aerosol backscatter ratio at the reference point. Rtolerance = hdf5read('data.h5', '/input/Rtolerance'); % Tolerance value for the cross-polarized lidar correction coefficient. Qtolerance = hdf5read('data.h5', '/input/Qtolerance'); % Initial approximation to aerosol backscatter ratio at the reference % point. Rinit = hdf5read('data.h5', '/input/Rinit'); % Initial approximation to the cross-polarized lidar correction % coefficient. Qinit = hdf5read('data.h5', '/input/Qinit'); TolFun = hdf5read('data.h5', '/input/TolFun'); TolX = hdf5read('data.h5', '/input/TolX'); %------------------------------------------------------------------------- % Check consistency of input array dimensions. assert(all(size(k) == [2, 1])); assert(all(size(nu) == [2, 1])); assert(all(size(hStep) == [1, 1]) && hStep > 0); assert(all(size(Nbeg) == [1, 1])); assert(all(size(Nref) == [1, 1])); assert(all(size(Nfin) == [1, 1])); assert(all(size(betaMref) == [1, 1])); assert(all(size(gamma) == [1, 1])); assert(all(size(S) == [Nfin, 2])); assert(all(size(Omega) == [Nfin, 2])); assert(all(size(betaM) == [Nfin, 2])); assert(all(size(tauM) == [Nfin, 1])); assert(all(size(Z) == [Nfin, 1])); assert(all(size(beta2Coeff) == [1, 1]) && ... (beta2Coeff == 0.0 || beta2Coeff == 1.0)); assert(all(size(leakage) == [1, 1])); assert(all(size(beta0) == [Nfin, 1])); assert(all(size(d0) == [Nfin, 1])); assert(all(size(R0) == [1, 1]) && R0 > 0); assert(all(size(Q0) == [1, 1]) && Q0 > 0); assert(all(size(Rtolerance) == [1, 1]) && Rtolerance > 0); assert(all(size(Qtolerance) == [1, 1]) && Qtolerance > 0); assert(all(size(Rinit) == [1, 1]) && Rinit > 0); assert(all(size(Qinit) == [1, 1]) && Qinit > 0); 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(Nfin, 2); for j = 1:2 Lmeas(Nbeg:Nfin, j) = S(Nbeg:Nfin, j) .* exp(-2 * tauM(Nbeg:Nfin)); end % Extend arrays related to lidar signal equations below the lower % signal boundaries. for j = 1:2 assert(all(Omega(Nbeg:Nfin) > 0)) Lmeas(1:Nbeg-1, j) = Lmeas(Nbeg, j); Omega(1:Nbeg-1, j) = Omega(Nbeg, j); betaM(1:Nbeg-1, j) = betaM(Nbeg, j); end Z(1:Nbeg-1) = Z(Nbeg); % Measured cross-polarized lidar signal ratio profile (by height). Ymeas = Lmeas(:, 2) ./ Lmeas(:, 1); % Steps along the lidar track, with zenith angles taken into account. effStep = hStep ./ cos(degtorad(Z(1:Nfin))); % Matrix that transforms extinction coefficient vectors to optical % thickness vectors, relative to the reference point. tauAMatrixTmp = tauIntegrationMatrix(Nbeg, Nref, Nfin, effStep); tauAMatrix = tauAMatrixTmp(Nbeg:Nfin, Nbeg:Nfin); % Matrix that transforms concentration profiles into vectors evaluating % their smoothness. smoothOp2 = secondDifferenceOperator(Nfin-Nbeg+1, 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 parallel backscatter profile. betaArr = zeros(Nfin, 1, 2); % Aerosol depolarization ratio profile. dArr = zeros(Nfin, 1, 2); % Aerosol backscatter ratio at the reference point. RArr = zeros(1, 1, 2); % Cross-polarized lidar correction coefficient. QArr = zeros(1, 1, 2); % Residuals for lidar equations. Psi1Arr = zeros(1, 2, 2); % Residuals for smoothness criteria. Psi2Arr = zeros(1, 2, 2); % Corrected measured lidar signal, divided by its average value around the % reference point, with molecular component removed (by height). LmeasRArr = zeros(Nfin, 1, 2); % Calculated lidar signal, divided by its value at the reference point, % with molecular component removed (by height). LcalcArr = zeros(Nfin, 1, 2); % Corrected measured cross-polarized lidar signal ratio profile (by % height). YmeasQArr = zeros(Nfin, 1, 2); % Calculated cross-polarized lidar signal ratio profile (by height). YcalcArr = zeros(Nfin, 1, 2); %------------------------------------------------------------------------- disp('Starting optimization...'); x0 = [beta0(Nbeg:Nfin); d0(Nbeg:Nfin); Rinit; Qinit]; l0 = zeros(size(x0)); u0 = ones(size(x0)) * inf; % Set boundary values for the correction coefficients. l0((Nfin-Nbeg+1)*2 + 1) = R0 * (1 - Rtolerance); u0((Nfin-Nbeg+1)*2 + 1) = R0 * (1 + Rtolerance); l0((Nfin-Nbeg+1)*2 + 2) = Q0 * (1 - Qtolerance); u0((Nfin-Nbeg+1)*2 + 2) = Q0 * (1 + Qtolerance); optopt = optimset('Display', 'Iter', 'DerivativeCheck', 'off', ... 'TolFun', TolFun, 'TolX', TolX, 'Jacobian', 'on'); % Invoke optimizer. lsqnonlin(@optfun, x0, l0, u0, optopt); %------------------------------------------------------------------------- disp('Saving output data...') hdf5write('data.h5', '/output/iterCount', iter, 'WriteMode', 'append'); hdf5write('data.h5', '/output/beta', betaArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/d', dArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/R', RArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Q', QArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Psi1', Psi1Arr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Psi2', Psi2Arr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/LmeasR', LmeasRArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Lcalc', LcalcArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/YmeasQ', YmeasQArr, 'WriteMode', 'append'); hdf5write('data.h5', '/output/Ycalc', YcalcArr, '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); % Remove constraints for the first and the last nodes, as % approximations for the second order derivative cannot be obtained % there. diffOp(1, 1:N) = 0; diffOp(N, 1:N) = 0; diffOp = diffOp / step^2; end %------------------------------------------------------------------------- function [F, Jac] = optfun(x) global k nu; global hStep Nbeg Nfin; global betaMref gamma Lmeas Ymeas Omega betaM; global beta2Coeff leakage; global tauAMatrix smoothOp2; global iter; global betaArr dArr RArr QArr Psi1Arr Psi2Arr global LmeasRArr LcalcArr YmeasQArr YcalcArr; iter = iter + 1; % Current parallel backscatter profile. beta = reshape(x(1 : (Nfin-Nbeg+1)), Nfin-Nbeg+1, 1); % Current depolarization ratio profile. d = reshape(x((Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2), Nfin-Nbeg+1, 1); % Current aerosol backscatter ratio at the reference point. R = x((Nfin-Nbeg+1)*2 + 1); % Current cross-polarized lidar correction coefficient. Q = x((Nfin-Nbeg+1)*2 + 2); betaArr(Nbeg:Nfin, 1, iter) = beta; dArr(Nbeg:Nfin, 1, iter) = d; RArr(1, 1, iter) = R; QArr(1, 1, iter) = Q; % Lidar portion of the equations. F1 = zeros((Nfin-Nbeg+1)*2, 1); Jac1 = zeros((Nfin-Nbeg+1)*2, (Nfin-Nbeg+1)*2 + 2); sigmaA = beta .* (ones(Nfin-Nbeg+1, 1) + d) * gamma; betaA = beta .* (ones(Nfin-Nbeg+1, 1) + beta2Coeff * d); % Calculated aerosol optical thickness, relative to the reference % point. tauA = tauAMatrix * sigmaA; LmeasR = Lmeas(Nbeg:Nfin, 1) * R; % Calculated lidar signal, divided by its value at the reference % point, with molecular component removed. Lcalc = (betaA + betaM(Nbeg:Nfin, 1)) / betaMref .* exp(2 * tauA); % Weighting coefficients for unpolarized lidar equations. coeff1 = k(1) * sqrt(hStep ./ Omega(Nbeg:Nfin, 1)); F1(1:(Nfin-Nbeg+1)) = (LmeasR - Lcalc) .* coeff1; Psi1Arr(1, 1, iter) = norm(F1(1 : (Nfin-Nbeg+1))); LmeasRArr(Nbeg:Nfin, 1, iter) = LmeasR; LcalcArr(Nbeg:Nfin, 1, iter) = Lcalc; % Derivative with respect to aerosol parallel backscatter beta. % Use 'sparse' for diagonal matrices to improve performance. Jac1(1:(Nfin-Nbeg+1), 1:(Nfin-Nbeg+1)) = ... -sparse(1:(Nfin-Nbeg+1), 1:(Nfin-Nbeg+1), ... (ones(Nfin-Nbeg+1, 1) + beta2Coeff * d) .* ... exp(2 * tauA) / betaMref .* coeff1) - ... sparse(1:(Nfin-Nbeg+1), 1:(Nfin-Nbeg+1), ... 2 * gamma * (ones(Nfin-Nbeg+1, 1) + d) .* ... Lcalc .* coeff1) * tauAMatrix; % Derivative with respect to aerosol depolarization ratio d. % Use 'sparse' for diagonal matrices to improve performance. Jac1(1:(Nfin-Nbeg+1), (Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2) = ... -sparse(1:(Nfin-Nbeg+1), 1:(Nfin-Nbeg+1), ... beta2Coeff * beta .* exp(2 * tauA) / ... betaMref .* coeff1) - ... sparse(1:(Nfin-Nbeg+1), 1:(Nfin-Nbeg+1), ... 2 * gamma * beta .* Lcalc .* coeff1) * tauAMatrix; % Derivative with respect to R. Jac1(1:(Nfin-Nbeg+1), (Nfin-Nbeg+1)*2 + 1) = ... Lmeas(Nbeg:Nfin, 1) .* coeff1; % Equation for the ratio between cross-polarized and unpolarized lidar % channels. Ynumerator = beta .* (d + ones(Nfin-Nbeg+1, 1) * leakage) + ... betaM(Nbeg:Nfin, 2); Ydenominator = beta .* (ones(Nfin-Nbeg+1, 1) + beta2Coeff * d) + ... betaM(Nbeg:Nfin, 1); YmeasQ = Ymeas(Nbeg:Nfin) * Q; Ycalc = Ynumerator ./ Ydenominator; % Weighting coefficients for the lidar profile ratio equation. coeff2 = k(2) * sqrt(hStep ./ (Omega(Nbeg:Nfin, 1) + ... Omega(Nbeg:Nfin, 2))); F1((Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2) = (YmeasQ - Ycalc) .* coeff2; Psi1Arr(1, 2, iter) = norm(F1((Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2)); YmeasQArr(Nbeg:Nfin, 1, iter) = YmeasQ; YcalcArr(Nbeg:Nfin, 1, iter) = Ycalc; % Derivative with respect to beta. Jac1((Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2, 1:(Nfin-Nbeg+1)) = ... -sparse(1:(Nfin-Nbeg+1), 1:(Nfin-Nbeg+1), ... ((d + ones(Nfin-Nbeg+1, 1) * leakage) .* Ydenominator - ... (ones(Nfin-Nbeg+1, 1) + beta2Coeff * d) .* Ynumerator) ./ ... (Ydenominator .* Ydenominator) .* coeff2); % Derivative with respect to d. Jac1((Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2, ... (Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2) = ... -sparse(1:(Nfin-Nbeg+1), 1:(Nfin-Nbeg+1), ... (beta .* Ydenominator - beta2Coeff * beta .* Ynumerator) ./ ... (Ydenominator .* Ydenominator) .* coeff2); % Derivative with respect to Q. Jac1((Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2, (Nfin-Nbeg+1)*2 + 2) = ... Ymeas(Nbeg:Nfin) .* coeff2; % Smoothness portion of the equations. % The matrices are sparse here, to improve performance. F2 = zeros((Nfin-Nbeg+1)*2, 1); Jac2Matrix = cell(2, 1); F2(1:(Nfin-Nbeg+1)) = ... nu(1) * sqrt(hStep) * smoothOp2 * beta(:); F2((Nfin-Nbeg+1)+1 : (Nfin-Nbeg+1)*2) = ... nu(2) * sqrt(hStep) * smoothOp2 * d(:); for i = 1:2 Psi2Arr(1, i, iter) = ... norm(F2((Nfin-Nbeg+1)*(i-1)+1 : (Nfin-Nbeg+1)*i)); Jac2Matrix{i} = nu(i) * sqrt(hStep) * smoothOp2; end Jac2 = [blkdiag(Jac2Matrix{1:2}), zeros((Nfin-Nbeg+1)*2, 2)]; F = [F1; F2]; Jac = [Jac1; Jac2]; end