Transcript Document

Variational methods for
retrieving cloud, rain and
hail properties combining
radar, lidar and radiometers
Robin Hogan
Julien Delanoe
Department of Meteorology, University of Reading, UK
Outline
• Increasingly in active remote sensing, many instruments are
being deployed together, and individual instruments may
measure many variables
– We want to retrieve an “optimum” estimate of the state of the
atmosphere that is consistent with all the measurements
– But most algorithms use at most only two instruments/variables and
don’t take proper account of instrumental errors
• The “variational” approach (a.k.a. optimal estimation theory) is
standard in data assimilation and passive sounding, but has
only recently been applied to radar retrieval problems
– It is mathematically rigorous and takes full account of errors
– Straightforward to add extra constraints and extra instruments
• In this talk, two applications will be demonstrated
– Polarization radar retrieval of rain rate and hail intensity
– Retrieving cloud microphysical profiles from the A-train of satellites
(the CloudSat radar, the Calipso lidar and the MODIS radiometer)
Polarization
radar: Z
• We need to retrieve rain rate for accurate flood forecasts
• Conventional radar estimates rain-rate R from radar reflectivity
factor Z using Z=aRb
– Around a factor of 2 error in retrievals due to variations in raindrop
size and number concentration
– Attenuation through heavy rain must be corrected for, but gate-bygate methods are intrinsically unstable
– Hail contamination can lead to large overestimates in rain rate
Polarization
radar: Zdr
• Differential reflectivity Zdr is a measure of drop shape, and
hence drop size: Zdr = 10 log10 (ZH /ZV)
– In principle allows rain rate to be retrieved to 25%
– Can assist in correction for attenuation
• But
– Too noisy to use at each range-gate
– Needs to be accurately calibrated
– Degraded by hail
ZH
1 mm
ZV
3 mm
4.5 mm
Polarization
radar: fdp
• Differential phase shift fdp is a propagation effect caused by the
difference in speed of the H and V waves through oblate drops
– Can use to estimate attenuation
– Calibration not required
– Low sensitivity to hail
• But
– Need high rain rate
– Low resolution information:
need to take derivative but far too noisy
to use at each gate: derivative can be negative!
phase shift
• How can we make the best use of the Zdr and fdp information?
Variational method
• Start with a first guess of coefficient a in Z=aR1.5
• Z/a implies a drop size: use this in a forward model to predict
the observations of Zdr and fdp
– Include all the relevant physics, such as attenuation etc.
• Compare observations with forward-model values, and refine a
by minimizing a cost function:
n
J 
i 1
Z
dr ,i
Z

2
Z dr
  f
fwd 2
dr ,i
dp ,i
f
f
  a  a 
fwd 2
dp ,i
2
dp
Observational errors are
explicitly included, and the
solution is weighted accordingly
ap 2
i

2
a ap
+ Smoothness
constraints
For a sensible solution at
low rainrate, add an a priori
constraint on coefficient a
How do we solve this?
• Generalize: suppose I have two estimates of variable x:
– m with error m (from measurements)
– b with error b (“background” or “a priori” knowledge of the PDF of x)

m  x
J
2
• The best estimate of x minimizes a cost
function:
• At minimum of J, dJ/dx=0, which leads to:
– The least-squares solution is simply a
weighted average of m and b, weighting each
by the inverse of its error variance
2m

b  x

2
b2
m
b

 2m b2
x
1
1

 2m b2
m  xi b  xi
 2
2
• Can also be written in terms of difference
m
b
x

x

of m and b from initial guess xi:
i 1
i
1
1

 2m b2
The Gauss-Newton method
• We often don’t directly observe the variable we want to
retrieve, but instead some related quantity y (e.g. we observe
Zdr and fdp but not a) so the cost function becomes
2
2

y  H ( x) b  x 
J

2
y
b2
– H(x) is the forward model predicting the observations y from state x
and may be complex and non-analytic: difficult to minimize J
• Solution: linearize forward model about a first guess xi
y
y
H ( x)  H ( xi ) 
 x  xi 
x xi
– The x corresponding to y=H(x), is
equivalent to a direct measurement m:
m  xi 
y  H ( xi )
y / x
…with error:  m 
Observation
y
y / x
y
x
xi+2 xi+1 xi
(or m)
• Substitute into prev. equation:
– If it is straightforward to
calculate y/x then iterate this
formula to find the optimum x
xi 1  xi 
y / x y  H ( xi ) b  xi
 2
2
y
b
y / x 2 
 2y
• If we have many observations
and many variables to retrieve
then write this in matrix form:


xi 1  xi  A1 HT R1y  H (xi ) B1b  xi 
Where the Hessian matrix is
A  H T R 1H  B 1
– The matrices and vectors are
defined by:
 x1 
 b1 
 y1 
 
 
 
x   b    y   
y 
x 
b 
 m
 n
 n
State vector, a priori vector
and observation vector
 y1

 x1
H 
 ym
 x
 1
y1 


xn 

 
ym 

xn 

The Jacobian
1
b2
  2y1



R 



2 

ym 

 b21



B



2 

bn 

Error covariance
matrices of observations
and background
Finding the solution
New ray of data
First guess of x
– In this problem, the
observation vector y and
state vector x are:
Forward model
Predict measurements y and
Jacobian H from state vector x
using forward model H(x)
Compare measurements to
forward model
Has the solution converged?
2 convergence test
 ln a1 


x  
 ln x 
n

 Z dr 1 


  
Z 
dr m

y 
 f dp1 
  


 f dp 
 m
No
Yes
Calculate error in retrieval
The solution error covariance
matrix is S=A-1
Proceed to next ray
Gauss-Newton iteration step
Predict new state vector:
xi+1= xi+A-1{HTR-1[y-H(xi)]
+B-1(b-xi)}
where the Hessian is
A=HTR-1H+B-1
Chilbolton example
3-GHz radar
25-m dish
• Observations
• Retrieval
Forward-model values at final iteration are essentially leastsquares fits to the observations, but without instrument noise
A ray of data
•
Zdr and fdp are well fitted by the
forward model at the final
iteration of the minimization of
the cost function
• The scheme also reports the
error in the retrieved values
• Retrieved coefficient a is forced
to vary smoothly
– Represented by cubic spline basis
functions
Enforcing smoothness 1
• Cubic-spline basis functions
– Let state vector x contain the
amplitudes of a set of basis functions
– Cubic splines ensure that the solution
is continuous in itself and its first and
second derivatives
– Fewer elements in x  more efficient!
Representing
a 50-point
function by
10 control
points
Forward model
Convert state vector to high
resolution: xhr=Wx
Predict measurements y and
high-resolution Jacobian Hhr
from xhr using forward model
H(xhr)
Convert Jacobian to low
resolution: H=HhrW
W
The
weighting
matrix
Enforcing smoothness 2
• Background error covariance matrix
– To smooth beyond the range of individual basis
functions, recognise that errors in the a priori
B
estimate are correlated
– Add off-diagonal elements to B assuming an
exponential decay of the correlations with range
– The retrieved a now doesn’t return immediately
to the a priori value in low rain rates
• Kalman smoother in azimuth
– Each ray is retrieved separately, so how do we
ensure smoothness in azimuth as well?
– Two-pass solution:
• First pass: use one ray as a constraint on the
retrieval at the next (a bit like an a priori)
• Second pass: repeat in the reverse direction,
constraining each ray both by the retrieval at
the previous ray, and by the first-pass
retrieval from the ray on the other side
Response to observational errors
Nominal Zdr error of ±0.2 dB
Additional random error of ±1 dB
Retrieved a
Zdr
and
fdp
Retrieval error
What if we use
only Zdr or fdp ?
Very similar retrievals: in
moderate rain rates, much
Zdr
only
more useful information
obtained from Zdr than fdp
fdp
Where observations
provide no information,
retrieval tends to a priori
value (and its error)
only
fdp only useful where
there is appreciable
gradient with range
Difficult case: differential attenuation of
1 dB and differential phase shift of 80º
Heavy rain
and
hail
• Observations
• Retrieval
How is hail retrieved?
• Hail is nearly spherical
– High Z but much lower Zdr than
would get for rain
– Forward model cannot match both
Zdr and fdp
• First pass of the algorithm
– Increase error on Zdr so that rain
information comes from fdp
– Hail is where Zdrfwd-Zdr > 1.5 dB and
Z > 35 dBZ
• Second pass of algorithm
– Use original Zdr error
– At each hail gate, retrieve the
fraction of the measured Z that is
due to hail, as well as a.
– Now the retrieval can match both
Zdr and fdp
Distribution of hail
Retrieved a
Retrieval error
Retrieved hail fraction
– Retrieved rain rate much lower
in hail regions: high Z no longer
attributed to rain
– Can avoid false-alarm flood
warnings
– Use Twomey method for
smoothness of hail retrieval
Enforcing smoothness 3
• Twomey matrix, for when we have no useful a priori information
– Add a term to the cost function to penalize curvature in the solution:
ld2x/dr2 (where r is range and l is a smoothing coefficient)
– Implemented by adding “Twomey” matrix T to the matrix equations

xi 1  xi  A1 HT R1y  H (xi ) B1b  xi  lTx i
A  H T R 1H  B 1  lT
 1 2 1



 2 5  4 1

 1 4 6 4 1


T
1  4 6  4 



1

4
6




  


Summary
• New scheme achieves a seamless transition between the
following separate algorithms:
– Drizzle. Zdr and fdp are both zero: use a-priori a coefficient
– Light rain. Useful information in Zdr only: retrieve a smoothly
varying a field (Illingworth and Thompson 2005)
– Heavy rain. Use fdp as well (e.g. Testud et al. 2000), but weight
the Zdr and fdp information according to their errors
– Weak attenuation. Use fdp to estimate attenuation (Holt 1988)
– Strong attenuation. Use differential attenuation, measured by
negative Zdr at far end of ray (Smyth and Illingworth 1998)
– Hail occurrence. Identify by inconsistency between Zdr and fdp
measurements (Smyth et al. 1999)
– Rain coexisting with hail. Estimate rain-rate in hail regions
from fdp alone (Sachidananda and Zrnic 1987)
Hogan (2006, submitted to J. Appl. Meteorol.)
The
A-train
• The CloudSat radar and the
Calipso lidar were launched
on 28th April 2006
• They join Aqua, hosting the
MODIS, CERES, AIRS and
AMSU radiometers
• An opportunity to tackle
questions concerning role of
clouds in climate
• Need to combine all these
observations to get an
optimum estimate of global
cloud properties
MODIS RGB composite
13.10 UTC
June 18th
Scotland
Lake
district
England
Isle of
Wight
France
MODIS Infrared window
13.10 UTC
June 18th
Scotland
Lake
district
England
Isle of
Wight
France
Met Office rain radar network
13.10 UTC
June 18th
Scotland
Lake
district
England
Isle of
Wight
France
• Calipso lidar
7 June 2006
Aerosol from
China?
Mixed-phase
altocumulus
Drizzling
stratocumulus
• CloudSat radar
Cirrus
Molecular
scattering
Non-drizzling
stratocumulus
Rain
East China Sea
Japan
Sea of Japan
Eastern Russia
5500
km
Motivation
• Why combine radar, lidar and radiometers?
– Radar ZD6, lidar b’D2 so the combination provides particle size
– Radiances ensure that the retrieved profiles can be used for radiative
transfer studies
• Some limitations of existing radar/lidar ice retrieval schemes
(Donovan et al. 2000, Tinel et al. 2005, Mitrescu et al. 2005)
– They only work in regions of cloud detected by both radar and lidar
– Noise in measurements results in noise in the retrieved variables
– Eloranta’s lidar multiple-scattering model is too slow to take to greater
than 3rd or 4th order scattering
– Other clouds in the profile are not included, e.g. liquid water clouds
– Difficult to make use of other measurements, e.g. passive radiances
– Difficult to also make use of lidar molecular scattering beyond the
cloud as an optical depth constraint
– Some methods need the unknown lidar ratio to be specified
• A “unified” variational scheme can solve all of these problems
Why not to invert the lidar separately
• “Standard method”: assume a value for the extinction-tobackscatter ratio, S, and use a gate-by-gate correction
– Problem: for optical depth d>2 is excessively sensitive to choice of S
– Exactly the same instability identified for radar in 1954
Implied
optical
depth is
infinite
• Better method (e.g. Donovan et al. 2000): retrieve the S that
is most consistent with the radar and other constraints
– For example, when combined with radar, it should produce a profile of
particle size or number concentration that varies least with range
Formulation of variational scheme
• Observation vector
– Elements may be missing
– Logarithms prevent
unphysical negative values
 ln b1 





 ln b 
n


Z1





y 

 Zm



d


 I 8.7 m 


 I 8.7 12.0m 
Attenuated lidar
backscatter profile
Radar reflectivity
factor profile (on
different grid)
Visible optical depth
Infrared radiance
Radiance difference
• State vector
 ln 1ice 


  
 ln  
n 

 ln N1ice 





 ln N mice 
x

S


 LW P 
1


liq
 ln N1 
  


aer
 ln 1 





Ice visible extinction
coefficient profile
Ice normalized
number conc. profile
Extinction/backscatter
ratio for ice
Liquid water path and
number conc. for each
liquid layer
Aerosol visible
extinction coefficient
profile
Radar forward model and a priori
• Create lookup tables
– Gamma size distributions
– Choose mass-area-size relationships
– Mie theory for 94-GHz reflectivity
• Define normalized number
concentration parameter
– “The N0 that an exponential
distribution would have with same
IWC and D0 as actual distribution”
– Forward model predicts Z from
extinction and N0
– Effective radius from lookup table
• N0 has strong T dependence
– Use Field et al. power-law as a-priori
– When no lidar signal, retrieval
relaxes to one based on Z and T
(Liu and Illingworth 2000, Hogan et
al. 2006)
Field et al. (2005)
Lidar forward model: multiple scattering
• 90-m footprint of Calipso
means that multiple
scattering is a problem
• Eloranta’s (1998) model
forward
scattered
photons escape
– O (N m/m !) efficient for N
points in profile and m-order
scattering
– Too expensive to take to more
than 3rd or 4th order in
retrieval (not enough)
• New method: treats third
and higher orders together
– O (N 2) efficient
– As accurate as Eloranta when
taken to ~6th order
– 3-4 orders of magnitude
faster for N =50 (~ 0.1 ms)
Narrow
field-of-view:
Wide field-ofview:
forward
scattered
photons may be
returned
Ice cloud
Molecules
Liquid cloud
Aerosol
Hogan (2006, Applied Optics, in press). Code: www.met.rdg.ac.uk/clouds
Radiance forward model
• MODIS solar channels provide an estimate of optical depth
– Only very weakly dependent on vertical location of cloud so we simply
use the MODIS optical depth product as a constraint
– Only available in daylight
• MODIS, Calipso and SEVIRI each have 3 thermal infrared
channels in atmospheric window region
– Radiance depends on vertical distribution of microphysical properties
– Single channel: information on extinction near cloud top
– Pair of channels: ice particle size information near cloud top
• Radiance model uses the 2-stream source function method
–
–
–
–
Efficient yet sufficiently accurate method that includes scattering
Provides important constraint for ice clouds detected only by lidar
Ice single-scatter properties from Anthony Baran’s aggregate model
Correlated-k-distribution for gaseous absorption (from David Donovan
and Seiji Kato)
Ice cloud: non-variational retrieval
Observations
State
variables
Aircraftsimulated
profiles with
noise (from
Hogan et al.
2006)
Retrieval is
accurate
but not
perfectly
stable
where lidar
loses signal
Derived
variables
• Donovan et al. (2000) algorithm can only be applied where
both lidar and radar have signal
Variational radar/lidar retrieval
Observations
Lidar noise
matched by
retrieval
State
variables
Derived
variables
Noise
feeds
through to
other
variables
• Noise in lidar backscatter feeds through to retrieved extinction
…add smoothness constraint
Observations
State
variables
Derived
variables
Retrieval
reverts to
a-priori N0
Extinction
and IWC
too low in
radar-only
region
• Smoothness constraint: add a term to cost function to penalize
curvature in the solution (J’ = l Si d2i/dz2)
…add a-priori error correlation
Observations
State
variables
Derived
variables
Vertical
correlation
of error in
N0
Extinction
and IWC
now more
accurate
• Use B (the a priori error covariance matrix) to smooth the N0
information in the vertical
…add visible optical depth constraint
Observations
State
variables
Derived
variables
Slight
refinement
to
extinction
and IWC
• Integrated extinction now constrained by the MODIS-derived
visible optical depth
…add infrared radiances
Observations
State
variables
Derived
variables
• Better fit to IWC and re at cloud top
Poorer fit
to Z at
cloud top:
information
here now
from
radiances
Convergence
• The solution generally
converges after two or
three iterations
– When formulated in terms of
ln(), ln(b’) rather than , b’,
the forward model is much
more linear so the minimum of
the cost function is reached
rapidly
Radar-only retrieval
Observations
State
variables
Derived
variables
• Retrieval is poorer if the lidar is not used
Use a priori
as no other
information
on N0
Profile poor
near cloud
top: no
lidar for
the small
crystals
Radar plus optical depth
Observations
State
variables
Derived
variables
Optical
depth
constraint
distributed
evenly
through
the cloud
profile
• Note that often radar will not see all the way to cloud top
Observed
94-GHz
radar
reflectivity
Observed
905-nm lidar
backscatter
Forward
model radar
reflectivity
Forward
model lidar
backscatter
Ground-based example
Lidar fails to
penetrate
deep ice cloud
Retrieved
extinction
coefficient 
Retrieved
effective
radius re
Retrieved
normalized
number conc.
parameter
N0
Error in
retrieved
extinction

Radar only:
retrieval tends
towards a-priori
Lower error in
regions with both
radar and lidar
Conclusions and ongoing work
• Variational methods have been described for retrieving cloud,
rain and hail, from combined active and passive sensors
– Appropriate choice of state vector and smoothness constraints ensures
the retrievals are accurate and efficient
– Could provide the basis for cloud/rain data assimilation
• Ongoing work: cloud
– Test radiance part of cloud retrieval using geostationary-satellite
radiances from Meteosat/SEVIRI above ground-based radar & lidar
– Retrieve properties of liquid-water layers, drizzle and aerosol
– Apply to A-train data and validate using in-situ underflights
– Use to evaluate forecast/climate models
– Quantify radiative errors in representation of different sorts of cloud
• Ongoing work: rain
– Validate the retrieved drop-size information, e.g. using a distrometer
– Apply to operational C-band (5.6 GHz) radars: more attenuation!
– Apply to other problems, e.g. the radar refractivity method
Sd
Banda Sea
An island of
Indonesia
Antarctic ice sheet
Southern Ocean