260 likes | 463 Views
Synergistic cloud, aerosol and precipitation products Progress so far in RATEC. Robin Hogan Julien Delano ë Nicola Pounder University of Reading. Synergy products as defined in CASPER. Other products follow from this. Develop “best estimate” first. Overview of our role in RATEC.
E N D
Synergistic cloud, aerosol and precipitation productsProgress so far in RATEC Robin Hogan JulienDelanoë Nicola Pounder University of Reading
Synergy products as defined in CASPER Other products follow from this Develop “best estimate” first
Overview of our role in RATEC • “Radiative transfer for EarthCARE” (RATEC) • Lead contractor: Howard Barker • Started August 23, 2009 • 15 month project • Our role is to start the development of synergy algorithms • Also define the required imager and target classification products • WP120: Review - first 3 months • Products and Algorithms Requirements Document (PARD) • Algorithm Development Plan (ADP) • WP 220: Design - next 5 months • Product Definition Document (PDD) • WP 320: Prototype and test - final 7 months • Algorithm Theoretical Basis Document (ATBD)
Overview of talk Laying the foundations for the algorithm... Summary of retrieval framework State variables (describing ice, liquid, rain and aerosol) Forward models (for radar, HSRL lidar, infrared & solar radiometers) Cost function Minimization techniques Gauss-Newton & Levenberg-Marquardt Gradient Descent (adjoint method) Ensemble methods Coding progress C++ library Scattering code Calculation and storage of error descriptors Solution error covariance matrix Averaging kernel Progress for specific atmospheric constituents Liquid cloud: gradient constraint
1. New ray of data: define state vector Use classification to specify variables describing each species at each gate Ice: extinction coefficient , N0’, lidar extinction-to-backscatter ratio Liquid: extinction coefficient and number concentration Rain: rain rate and mean drop diameter Aerosol: extinction coefficient, particle size and lidar ratio 6. Iteration method Derive a new state vector Retrieval framework 2. Convert state vector to radar-lidar resolution Often the state vector will contain a low resolution description of the profile 3. Forward model 3a. Radar model Including surface return and multiple scattering 3c. Radiance model Solar and IR channels 3b. Lidar model Including HSRL channels and multiple scattering Ingredients developed before Not yet developed Not converged 4. Compare to observations Check for convergence 5. Convert Jacobian to state-vector resolution Jacobian initially will be at the radar-lidar resolution Converged 7. Calculate retrieval error Error covariances and averaging kernel Proceed to next ray of data
Proposed state variables Ice clouds largely done in CASPER; Snow & riming in convective clouds needs to be added (VARSY?) Liquid clouds to be tackled in RATEC Basic rain to be added in RATEC; Full representation in VARSY? Basic aerosols added in RATEC; Full representation via collaboration (IRMA?)
Forward model components • From state vector x to forward modelled observations H(x)... Jacobian matrix H=y/x x Ice & snow Liquid cloud Rain Aerosol Lookup tables to obtain profiles of extinction, scattering & backscatter coefficients, asymmetry factor Ice/radar Ice/lidar Ice/radiometer Liquid/radar Liquid/lidar Liquid/radiometer Lots of expensive matrix multiplications: likely to be the most expensive part of the entire algorithm Rain/radar Rain/lidar Rain/radiometer Aerosol/radiometer Aerosol/lidar Sum the contributions from each constituent Lidar scattering profile Radar scattering profile Radiometer scattering profile Radiative transfer models Gradient of radar measurements with respect to radar inputs Gradient of lidar measurements with respect to lidar inputs Gradient of radiometer measurements with respect to radiometer inputs H(x) Lidar forward modelled obs Radar forward modelled obs Radiometer forward modelled obs
Solution methods (1/3) • We want to minimize this cost function: • Possibly using the gradient (a vector): • …and the second derivative (a matrix): 1. Gauss-Newton method (Rodgers p85): • Advantage: rapid convergence (instant convergence for linear problems) • Another advantage: get the error covariance of the solution “for free” • Disadvantage: need the Jacobian of every forward model: can be expensive
Solution methods (2/3) 2. Gradient descent method(e.g. ECMWF data assimilation system): • Advantage: we don’t need to calculate the Jacobian so forward model is cheaper! • Quasi-Newton method to get search direction (conjugate gradient, BGFS etc) • Disadvantage: more iterations needed since we don’t know curvature of J(x) • Disadvantage: don’t get the error “for free” at the end • Why don’t we need the Jacobian H? • The “adjoint” of a forward model takes as input the vector {.} and outputs the vector Jobs without needing to calculate the matrix H on the way • Adjoint can be coded to be only ~3 times slower than original forward model • Tricky coding for newcomers, although some automatic code generators exist Simple steepest descent requires many steps Conjugate gradient method more intelligent…
Solution methods (3/3) 3. Levenberg-Marquardt (hybrid of Gauss-Newton & steepest descent): • If cost function reduces, keep h small: rapid convergence • If cost function increases, increase h so that always go downhill • More robust than Gauss-Newton, but still need the Jacobian 4. Ensemble method (e.g. Ensemble Kalman Filter): • Create ensemble of “trial solutions” or “particles” in state space • Spread of results allows Jacobian to be calculated “automatically” to create better set of particles: repeat until convergence • Applied to GPM retrievals by Mircea Grecu (NASA) • Advantage: no Jacobian required • Estimate of retrieval error from final spread • Disadvantage: may need many ensemble members so probably slower than adjoint method 5. Simulated annealing (e.g. Donovan in ECSIM) • Suitable for very non-linear problems • Too slow for this application? . . . . . .
Computational cost can scale with number of points describing vertical profile N; we can cope with an N2dependencebut not N3 Radiative transfer forward models • Lidar will use PVC+TDTS, radar will use single scattering+TDTS • Jacobian of TDTS is too expensive: need to develop reduced-resolution Jacobian or alternatively adjoint code • Also need depolarization forward model with multiple scattering • Infrared will probably use RTTOV, solar radiances will use RADIANT
Coding progress • Coding has started on the underpinning libraries • Use C++ because of key language features • Matrix and Vector Expression Library • Operator overloading and expression templates allows complex operations on matrices to be written on one line, with similar efficiency as Fortran-90 • Important feature is to allow a matrix object to behave as part of another matrix – important when different parts of the Jacobian refer to different instruments and need to be passed to functions in bits • Optimal Synergy Library • Object orientation and polymorphism: Radar and Lidar classes inherit from generic Observation class; IceCloud and Aerosol inherit from generic Constituent class • Enables code to be written to be completely flexible: observations can be added and removed without needing to keep track of indices to matrices, so same code can be applied to many platforms • Also written code to provide scattering libraries in NetCDF files
Storage of error descriptors If ice is present at N points and M variables are retrieved: MN retrieved variables to report with MN errors Error covariance and averaging kernel matrices also needed A standard output from a variational algorithm But ~(MN)2/2 values in each matrix, much larger than the retrieved variables and their errors Matrix is a different size each ray! Need a suitable compression strategy... “Measurements” “Measurements” Retrievals Retrievals
Retrieval error covariance matrix (Extinction, lidar ratio, N’) Top -> Height -> bottom N’ S a Top -> Height -> bottom a N’ Top -> Height -> bottom Top -> Height -> bottom
Retrieval error correlation matrix (Extinction, lidar ratio, N’) Top -> Height -> bottom N’ S a Top -> Height -> bottom a N’ Top -> Height -> bottom Top -> Height -> bottom
For autocorrelation of extinction... • 3 Gaussians requires 8N or fewer values to approximate the error covariance matrix • Need to work on correlations between different variables
Averaging kernel • The averaging kernel matrix expresses how the retrieval at a point depends on the truth at every other point • The sum of each column expresses how much the retrieval at that point derives from the measurements (rather than the prior): 1 is good! • The shape of each column expresses the effective vertical resolution of the retrieval: more peaked is good! • Easy to parameterize: • Store integral and peak • Or store integral and the width for each column Height (m) Information content
Gradient constraint • We have a good constraint on the gradient of the state variables with height for: • LWC in stratocu (adiabatic profile, particularly near cloud base) • Rain rate (fast falling so little variation with height expected) • Not suitable for the usual “a priori” constraint • Solution: add an extra term to the cost function to penalize deviations from gradient c: A. Slingo, S. Nichols and J. Schmetz, Q. J. R. Met. Soc. 1982
Example in liquid clouds • Using simulated observations: • Triangular cloud observed by a 1- or 2-field-of-view lidar • Retrieval uses Levenberg-Marquardt minimization with Hogan and Battaglia (2008) model for lidar multiple scattering Two FOVs: very good performance One 10-m footprint: saturates at optical depth=5 One 100-m footprint: multiple scattering helps! One 10-m footprint with gradient constraint: can extrapolate downwards successfully Optical depth=30 Footprint=100m Footprint=600m
Ongoing work • Implement various instrument forward models • Not too difficult: interfacing different codes and languages • Write adjoint for slowest models • Solver • Test gradient descent algorithms (conjugate gradient, BGFS etc) • Implementation of different atmospheric constituents • Ice cloud: we know how to do this • Liquid cloud: to be done in RATEC • Aerosol and precipitation: basic implementation in RATEC; longer-term development needs collaboration with DAME, IRMA etc... • Test datasets • Several airborne comparisons for validation (TC4, POLARCAT etc) • A-Train • ECSIM (after RATEC) • Two-instrument synergy products • To be tested using same code removing one instrument (after RATEC)
Aerosol refractive index 355-nm lidar Solar MSIchannels • Values from Ellie Highwood • Strong dependence on aerosol type • Varies with hydration • Can we use one continuous variable to represent different aerosol properties? • Need to make simplifications if information from multiple wavelengths is to be combined Thermal MSI channels
Matrix and Vector Expression Library • Intuitive interface and ability to do complex expressions Matrix M(3,3); Vector x(3) = 1.0, 2.0, 3.0; // Simple initialization M = N * R + exp(P); // ‘*’ invokes element multiplication y = transpose(x) ()* inverse(R); // ‘()*’ invokes matrix multiplication x = solve(A, b); // Linear algebra • Important aspects not in any other library (to my knowledge) S.subset(M, 3, 30, 2, 20); // Matrix subsetting M(find(M < 0.0)) = 0.0; // Matlab-style indexing BitfieldMatrix b; b.bit(4) = 1; // Handling of classification data • Efficient and robust implementation • “Expression templates” to avoid temporary objects in long expressions • Reference counting for allocated objects protects against memory leaks
Key classes • State • Methods: • solve() • add_obs(Observation) • Data: • List of Observations • List ofConstituents • Vector x, y • Matrix Jacobian etc. BackgroundConstituents • Constituent • Methods: • calc_scat_props(wavelength) • Data: • Vector x, x_prior • Observation • Methods: • calc_scat_props() • virtual forward_model() • Data: • Real wavelength • Vector y, y_error • Matrix Jacobian IceCloud LiquidCloud Rain Aerosol • ActiveInstrument • Methods: • forward_model() • Radiometer • Methods: • forward_model() Radar Lidar
Forward model part 1 State vector x xice . . . . . xliq . . . To radiative transfer xradar xradar,gas xradar,ice xradar,liq bradar . . . aradar . . . 0 0 0 0 aradar . . . bradar . . . aradar . . . bradar . aradar . = + + xhighres,ice W1 xice lnavis . . . lnN0* . . . xhighres,ice/ xice Basis functions for variables represented at reduced resolution lnavis . . . lnN0’ . W3 Sum the contributions from each atmospheric constituent = xradar / xradar,ice xradar,ice W2 xradar,ice / xhighres,ice bradar . . . aradar . . . Scattering lookup table
Forward model part 2 Radiative transfer y'radar Hradar y’radar / xradar Z . . . Hradar,ice Hradar W2 W1 W3 y’radar /xice y’radar / xradar xradar / xradar,ice xhighres,ice / xice xradar,ice / xhighres,ice = y' H Z . . . b' . . . I y’/ x To solver Forward modelled observations and Jacobian