Method for determining the accuracy of measured environmental data parameters

Abstract

A method of estimating expected errors of environmental data parameters based on radiance measurements obtained from visible infra red radiometric satellite sensors (VIIRS) orbiting the earth, comprising the steps of: obtaining N radiance measurements of a surface body I.sub.I, . . . ,I.sub.N defining matrix I depending on p unknown surface and atmospheric parameters T.sub.i, . . . T.sub.p, defining a matrix T; generating a forward model I=f(T) for obtaining the I radiance measurements from the p parameters; choosing an initial set of values for the p parameters and linearizing f(T) about the initial values to obtain a linearized forward I=s+H.theta. as f.sub.i (T.sub.01, T.sub.02, . . . T.sub.0p)+.SIGMA.H.sub.ij .theta..sub.j where I=1,2, . . . N and .theta..sub.j =T.sub.oj and H.sub.ij =.delta.f.sub.i /.delta.T.sub.j and where .theta. is a column matrix .theta..sub.l, . . . .theta..sub.p and H is a matrix of H.sub.ij values; adding measurement noise vector w of noise values to the forward model; determining the covariance of the measurement noise w to obtain a covariance matrix C; and manipulating the matrices H and C according to the equation C.sub.EDR =(H.sup.T C.sup.-1 H).sup.-1 to obtain matrix element C.sub.EDR indicative of the expected errors in the values of T.sub.1, . . . T.sub.p parameters.

Description

FIELD OF THE INVENTION

The present invention relates generally to satellite radiometry measurements for measuring environmental data parameters, and more particularly, to a method for determining the accuracy with which instruments providing said radiometry measurements yield the true values of the measured parameters at the earth's surface.

BACKGROUND OF THE INVENTION

There are certain phenomena on our earth which have great significance from military and/or commercial points of view. One such phenomenon is the water temperature of the oceans or other large bodies of water. Knowing the temperature of the water over large areas is of importance for military purposes, as well as from a commercial point of view, since ocean water temperature affects commercial fishing. Theoretically, one can take direct measurements of the water temperature at specific locations. However, due to the large surface area of the oceans, taking direct water temperature measurements over any large surface area of hundreds of thousands of square miles is clearly impractical.

Similarly, cloud measurements in general, and more particularly, ice particle size measurements of cirrus clouds, are also important in imaging certain bodies in order to remove or limit interfering clouds from any satellite images, as well as to obtain measurements of precipitation occurring within the atmosphere. Satellite sensor instruments have been used to measure specified radiance of certain ground-based bodies to within a specified error. That is, the satellite sensors themselves make known radiance measurements at the top of the atmosphere of bodies located at the earth's surface, such as water temperature.

However, the goal is to determine a given Environmental Data Record (EDR) or surface parameter to within a specified error. It is therefore crucial to relate the radiance measurement errors, which we know how to predict, to the expected errors in the EDR's. For many EDR's it is relatively easy to create forward models which, given specific values of the EDR parameters, predict the radiance leaving the top of the atmosphere. What is not so easy, however, is to exercise these models over a range of environmental parameters large enough to find the level of radiance error which produces adequate EDR estimates. Fortunately, as will be shown herein, signal processing formulas applied to linearized versions of an EDR forward model, permit one to proceed directly from the known radiance error (i.e. the satellite sensor noise of the point design) to the expected errors in the EDR. These formulas also specify algorithms, based on the linearized forward model, which can be used to estimate EDR's from their associated radiances.

As will be shown in the detailed description below, a basic statistical methodology for connecting EDR errors to radiance measurement errors will be presented in Section 1, followed by two examples in Sections 2 and 3 showing how to deal with VIIRS-based EDR's when the forward models can be written down algebraically. The Sea-Surface Temperature EDR analyzed in Section 2 shows that the predicted EDR error using the method disclosed herein matches NOAA's past experience. This provides support that the statistical methodology works as expected. The cirrus cloud EDR analyzed in Section 3 shows how the methodology can be used to analyze the performance different versions of an EDR algorithm. Section 4 shows how to handle EDR's where only a computerized version of the forward model (such as MODTRAN or FASCODE) exists. The technique described in Section 4 applies equally well, of course, to those VIIRS-based EDR's where only computerized versions of the forward model are easily available.

SUMMARY OF THE INVENTION

It is an object of the present invention to provide a method for estimating expected errors of environmental data parameters based on radiance measurements obtained from visible infra red radiometric satellite sensors (VIIRS) orbiting the earth, comprising the steps of obtaining N radiance measurements of a surface body I.sub.I, . . . I.sub.N defining matrix I depending on p unknown surface and atmospheric parameters T.sub.i, . . . T.sub.p, defining a matrix T; generating a forward model I=f(T) for obtaining the I radiance measurements from the p parameters, where f(T) comprises a matrix; choosing an initial set of values for the p parameters and linearizing f(T) about the initial values to obtain a linearized forward I=s+H.theta. as f.sub.i (T.sub.01, T.sub.02, . . . T.sub.0p)+.SIGMA.H.sub.ij .theta..sub.j where I=1,2, . . . N and .theta..sub.j =T.sub.j -T.sub.oj and H.sub.ij =.delta.f.sub.i /.delta.T.sub.j and where .theta. is a column matrix .theta..sub.1, . . . .theta..sub.p and H is a matrix of H.sub.ij values; adding measurement noise vector w of noise values to the forward model; determining the covariance of the measurement noise w to obtain a covariance matrix C; and manipulating the matrices H and C according to the equation C.sub.EDR =(H.sup.T C.sup.-1 H).sup.-1 to obtain matrix element C.sub.EDR indicative of the expected errors in the values of T.sub.1, . . . T.sub.p parameters.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a standard two-layer transmission-emission model for the ocean surface and atmosphere according to the present invention.

FIG. 2 illustrates the transmission curves as a function of wavenumber listed in Table 1 and shows how well they match the straight-line approximations generated by the values of M and A.sub.atm.

FIG. 3 represents a graph of functions P and .alpha. for T=290 deg K., 300 deg K., 310 deg K. and 800 cm.sup.-1 .ltoreq.v.ltoreq.1000 cm.sup.-1.

FIG. 4 shows the relative error involved in approximating P by .alpha. for three different temperatures.

FIG. 5 illustrates the physical model for measuring cirrus cloud temperature and ice-particle size using radiance measurements according to the present invention.

FIG. 6 provides an overall flow chart for estimating the expected errors of environmental data parameters based on radiance measurements from VIIRS sensors according to the present invention.

FIG. 7 is a flow chart depicting the steps for predicting errors in sea surface temperature measurements according to the present invention.

FIGS. 8A-B represent a flow chart depicting a method for estimating error involved in determining ice particle size of cirrus clouds using VIIRS radiance measurements according to the present invention.

FIG. 9 contains Table 1 which provides a set a of ratios for midlatitude winter to midlatitude summer and for tropical to midlatitude summer temperature values a function of wavenumber.

FIG. 10 contains Table 2 which provides a ration of the slope to the y-intercept from the value of Table 1.

FIG. 11 contains Table 3 which lists the expected errors in sea-surface temperature .DELTA.T, based on the VIIRS point design, for different values of the sea-surface temperature.

FIG. 12 contains Table 4 which provides a list of expected sea-surface errors based on AVHRR specifications for different values of sea-surface temperature.

FIG. 13 contains Table 5 which illustrates parameter values representing bands 3 and 4 of the VIIRS point design for the 3.7 .mu.m bands of the 10-12 .mu.m bands.

FIG. 14 contains Table 6 which represents ice particle size measurement values, and the absolute and relative errors associated with those measurements over a range of T.sub.C temperature values.

DETAILED DESCRIPTION OF THE INVENTION

Section 1: The Statistical Methodology

With reference now to FIG. 6, as part of the statistical methodology for VIIRS-based EDR parameters, we suppose there exist N measurements (e.g. N in-band sensor signals) which depend on p unknown atmospheric and surface parameters. In order to make the mathematical notation compact the N measurements I.sub.1, I.sub.2, . . . , I.sub.N are stacked in a column and treated as a column vector I. ##EQU1##

The p unknown surface and atmospheric parameters T.sub.1, T.sub.2, . . . , T.sub.p are treated similarly. ##EQU2##

In general the T are the EDR parameters which are to be determined from the satellite sensor measurements I. For example, there may be N=2, p=2 for a sea-surface temperature EDR, where I.sub.1,2 are two radiance measurements in the 10-12 .mu.m atmospheric window, and T.sub.1,2 are respectively the sea-surface temperature and an atmospheric transmission-emission parameter as shown in equations (2.8a-d) below. Note that for CrIS-based atmospheric profiles there may be N.congruent.1500 for the spectral channels and p.congruent.70 for a 35-layer atmospheric model with temperature and water-vapor concentrations to be determined for each of the layers.

The forward model going from T to I can be written as

I=f(T) where f(T) stands for a column vector of functions all depending on T. ##EQU3## The next step is to choose a standard set of values for the surface and atmospheric parameters. ##EQU4## and linearize f.sub.1, f.sub.2, . . . , f.sub.N about these values. ##EQU5## where i=1,2, . . . , N, .theta..sub.j =T.sub.j -T.sub.0j, (1.3c) and ##EQU6## If one defines ##EQU7## and H=matrix of H.sub.ij values (1.4c) then the linearized forward model in equation (1.2a) can be written as I=H.theta.+s (1.5)

Measurement noise associated with the N radiance measurements is represented as a vector w of N random numbers w.sub.1, W.sub.2, . . . , W.sub.N having a jointly gaussian probability distribution.

The probability that w.sub.1 has a value between w.sub.1 and w.sub.1 +dw.sub.1, w.sub.2 has a value between w.sub.2 and w.sub.2 +dw.sub.2, . . . , w.sub.N and W.sub.N +dw.sub.N is p(w.sub.1, w.sub.2, . . . w.sub.N)dw.sub.1 dw.sub.2. . . dw.sub.N where ##EQU8## Note that this is a zero-mean gaussian distribution. In equation (1.7a) C is the covariance matrix of the w.sub.1, W.sub.2, . . . , w.sub.N random numbers, det C is the determinant of C, and w.sup.T =(w.sub.2 w,..sub.N w ) is the transpose of w (which is, of course, a row vector). The ij'th component of matrix C is (C).sub.i =C.sub.i =E((w.sub.i -E(w.sub.i)).multidot.(w.sub.j -E(w.sub.j)))=E(w.sub.i w.sub.j) (1.7b) where E stands for the usual expectation operator. In the last step of equation (1.7b) we use that E(w.sub.i)=0 for all values of I because the w.sub.i obey a zero-mean gaussian distribution. For any two random numbers r.sub.1, r.sub.2 (which do not necessarily obey a zero mean distribution) we have E(r.sub.1,2)=mean value of r.sub.1,2 (1.7c) E((r.sub.1,2 -E(r.sub.1,2)).sup.2)=variance of r.sub.1,2 (1.7d) E((r.sub.1 -E(r.sub.1)).(r.sub.2 -E(r.sub.2)))=covariance of r.sub.1 and r.sub.2 (1.7e)

One can then construct a complete linearized forward model, including measurement noise, by adding w, the vector or random measurement errors, to the right-hand side of equation (1.5) to obtain:

I=H.theta.+s+wtm (1.8) The error in the estimate of .theta..sub.1, .theta..sub.2, . . . , .theta..sub.p is then given by the covariance matrix C.sub.EDR =(H.sup.T I.sup.-1 H).sup.-1 (1.9a) where the "T" superscript denotes the matrix transpose and the "-1" superscript denotes the matrix inverse. The covariance matrix C.sub.EDR follows from the theorem established in Steven M. Kay's, Fundamentals of Statistical Signal Processing: Estimation Theory, page 97, PTR Prentiss-Hall, Inc., Englewood Cliff, N.J., 1993 and incorporated herein by reference. The matrix elements of C.sub.EDR specify the expected errors in the values of .theta..sub.1,.theta..sub.2, . . . .theta..sub.p, which is the same as the expected errors in the values of T.sub.1,T.sub.2, . . . T.sub.p because the .theta.'s are just the differences of the p surface and atmospheric parameters from their standard values (see equation (1.3c)). The minimum-variance, unbiased estimator of .theta. given noise-contaminated values of I0 (the noise amplitude is specified by the covariance matrix C) and known values of s is .theta..sub.estimate =(H.sup.T C.sup.-1 H).sup.-1 H.sup.T C.sup.-1 (I-s)(1.91b) Equation (1.9a) specifies the variance, that is, the error, in the estimate of .theta. given by equation (1.9b). If the forward model in equation (1.8) is strictly true rather than a linearized version of equation (1.2a), then equation (1.9b) is the best possible way to estimate .theta. and equation (1.9a) gives the lowest possible error with which .theta. can be estimated.

The constant s in equation (1.8) and (1.9b) is a known quantity so we can always define a new set of measurements D such that

D-I-s (1.10a) and D=H.theta.+w (1.10b) If .theta. has been measured often enough in the past that one can estimate the probability that the next value of .theta. will take on any particular value, and if it is known that this probability distribution for .theta. looks gaussian with a covariance matrix C.sub..theta. such that (C.sub.0).sub.ij =E((.theta..sub.i -E(.theta..sub.i))(.theta..sub.j -E(.theta..sub.j))). (1.11a) then it is also known that the expected error in our estimate of .theta. will now be C.sub.EDR =(C.sub..theta..sup.-1 +H.sup.T C.sup.-1 H).sup.-1(1.11b) when using the estimation formula (see page 391 of Fundamentals of Statistical Processing: Estimation Theory) .theta..sub.estimate =E(.theta.)+[(C.sub.0.sup.-1 +H.sup.T C.sup.-1 H).sup.-1 H.sup.T C.sup.-1 ].(D-H.E( .theta.)) (1.11c) If the D values are unknown so as to be unable to provide information regarding the values of .theta., one may proceed by obtaining a best estimate as to the true values of .theta. as E(.theta.)=(E(.theta..sub.1), E(.theta..sub.2), . . . , E(.theta..sub.3)).sup.T, the mean of the set of previously-determined .theta. values. However, we do in fact have the .theta. values generated by the true .theta. values, the right-hand side of equation (1.11c) uses the information provided by the D values to adjust the already known values of E(.theta.) to get a better estimate of .theta.. In equations (1.3a,b,c) above, we can choose to linearize the forward model about a T.sub.0 vector such that E(.theta.)=0 (indeed that would be the most natural way of linearizing the forward model). This simplifies equation (1.11c) to .theta..sub.estimate =[(C.sub..theta..sup.-1 +H.sup.T C.sup.-1 H) .sup.-1 H.sup.T C.sup.-1 ]. D (1.11d) Equation (1.11 d) represents a standard linear algorithm used to retrieve temperature and water vapor from a set of FTIR spectra. The algorithm of this type is disclosed in pages 12-14 of J. P. Kerekes, A Retrieval Error Analysis Technique for Passive Infrared Atmospheric Sounders, Lincoln Laboratory MIT, Technical Report 978, Jul. 8, 1993, ESC-TR-92-204, incorporated herein by reference. Note that when (1.11d) is used, equation (1.11b) automatically gives us the expected error in the retrieved temperature and water-vapor profiles. If no information is available about which temperatures and water-vapor concentrations are most likely to occur, then all the elements of C.sub..theta. are zero and equations (1.11b,d) reduce to (1.9a,b). Section 2: Sea-Surface Temperature

A forward model for determining the sea-surface temperature from satellite observations of a cloud-free region of the ocean is given as

I(v,atm)=B(v,T.sub.s).multidot..tau.(v,atm)+B(v,T.sub.a).multidot.[1-.tau.(v,atm)] (2. 1) where I is the radiance leaving the earth's atmosphere (in watt/cm.sup.2 /sr/cm.sup.-1), T.sub.s is the sea-surface temperature (in deg K), T.sub.a is an effective temperature (in deg K) describing the atmospheric self-radiance in the 10-12 .mu.m window region, .tau. is the dimensionless atmospheric transmission from the ocean surface to the top of the atmosphere, v is a wavenumber (in cm.sup.-1) in the 10-12 .mu.m window region, "atm" is shorthand for a long list of atmospheric parameters determining the value of the transmittance .tau. at wavenumber v, and B(v,T) is the Planck function B(v,T)=C.sub.1 v.sup.3 P(C.sub.2 vT.sup.-1) (2.2a) where ##EQU9##

Equation (2.1) represents a standard two-layer transmission-emission model for the ocean surface and atmosphere as shown in FIG. 1. Empirical data has shown that T.sub.a varies by less than 1 deg K in the 10-12 .mu.m window region, so the dependence of T.sub.a on v can be neglected in equation (2.1). Another good approximation is that in the 10-12 .mu.m window

.tau.(v,atm).congruent.1-A.sub.atm .multidot.K(v) (2.3) where K is a function only of wavenumber and the effect of all the atmospheric parameters can be lumped into a single constant A.sub.atm. As a check to equation (2.3), a MODTRAN simulation run for 35 wavenumbers between 820 cm.sup.-1 and 990 cm.sup.-1 has been performed using three different standard atmospheres: tropical, midlatitude summer, and midlatitude winter. Table 1 (FIG. 9) gives the ratios of [1-.tau.(v,atm)] for midlatitude winter to midlatitude summer and for tropical to midlatitude summer, and we see that while some correlation exists with the approximation given in (2.3), there are definitely wavenumbers where it does not work very well. We will, however, be interested in wide-bandwidth radiance signals, so we should also examine how well equation (2.3) works as a moving average relationship. Performing a linear regression of [1-.tau.(v,atm)] against v for the values of 1-.tau. and v given in Table 1 yields the ratio of the slope to the y-intercept as shown in Table 2 (FIG. 10). As one can ascertain, the ratios are very similar for all three atmospheres (in fact they are equal to within 4%). Section 5 shows that this is exactly the result to be expected if equation (2.3) held true, so equation (2.3) can be used with some confidence as long as it is applied only to wide-bandwidth radiance signals. FIG. 2 gives the transmission curves .tau.(v,atm) listed in Table 1 and shows how well they match the straight-line approximations generated by the values of M and A.sub.atm in equations (2.4a,b). .tau.(v,atm).congruent.1-A.sub.atm .multidot.K.sub.L (v) (2.4a) where K.sub.L (v)=1-M.multidot.v (2.4b) for M=8.0.times.10.sup.-4 cm,

1.8 for standard tropical atm

and

A.sub.atm =1.2 for standard midlatitude summer atm

0.27 for standard midlatitude winter atm

Approximating the function K(v) in equation (2.3) by the straight line K.sub.L (v) given in (2.4b) and substitute (2.4a,b) into equations (2.1) yields.

I(v,atm)=B(v,T.sub.s)+K.sub.L (v).multidot.A.sub.atm .multidot.[B(v,T.sub.a)-B(v,T.sub.s)] (2.5) From page 5042 of C. Prabhakara, G. Dalu, and V. G. Kunde's article on "Estimation of Sea Surface Temperature From Remote Sensing in the 11- to 13 .mu.m Window Region", J. Of Geophysical Research, vol. 79, No. 33, 1974, pp. 5039-5044, incorporated herein by reference, we know that T.sub.s and T.sub.a are almost always within 10 deg K of each other, so for a given patch of ocean at a specified time of year, the Planck function B(v,T) will only need to be evaluated over a limited range of V and T values. A reasonable range for midlatitude summer is 285 deg K.ltoreq.T.ltoreq.315 deg K (2.6a) 800 cm.sup.-1 .ltoreq.v.ltoreq.1000 cm.sup.-1 (2.6b) For this range of v and T function P(C.sub.2 vT.sup.-1) can be approximated as P(C.sub.2 vT.sup.-1).congruent..alpha.(C.sub.2 vT.sup.-1)=.OMEGA..multidot.(C.sub.2 vT.sup.-1).sup.-n (2.7a) where .OMEGA.=8.2527 (2.7b) and n=4.4 (2.7c) FIG. 3 shows functions P and .alpha. for T=290 deg K, 300 deg K, 310 deg K and 800 cm.sup.-1 .ltoreq.v.ltoreq.1000 cm.sup.-1 ; one can see that a does a reasonable job of approximating P. FIG. 4 shows the relative error involved in approximating P by .alpha., and for the three temperatures graphed in FIG. 3, the relative error is 3% or less. We substitute (2.2a) and (2.7a) into equation (2.5) to get I(v,atm)=C.sub.b v.sup.-3-n .theta..sub.s +C.sub.B K.sub.L (v)v.sup.3-n .theta..sub.a (2.8a) where C.sub.B =.OMEGA.C.sub.1 C.sub.2.sup.-n .congruent.1.98.times.10.sup.-12 watt/cm.sup.2 /sr/(cm.sup.-1).sup.4 /(cm-degK).sup.n, (2.8b) .theta..sub.s =T.sub.s.sup.n (2.8c) and .theta..sub.a =A.sub.atm .multidot.(T.sub.a.sup.n -T.sub.s.sup.n)(2.8d) Equations (2.8a-d) reduce the nonlinear problem of equation (2.5) into a linear problem in two unknowns: .theta..sub.s the n'th power of the sea-surface temperature, and .theta., an atmospheric variable describing the combined effect of T.sub.s, T.sub.a, and A.sub.atm. Therefore, for this particular VIIRS-based EDR we do not have to linearize the forward model to use equations (1.9a,b), although that could have been done had it turned out to be necessary.

The VIIRS point design specifies two bandwidths in the 10-12 .mu.m atmospheric window which are labeled band 1 and band 2. Band 1 ranges from wavenumbers V.sub.1a to v.sub.1b and band 2 ranges from v.sub.2a to v.sub.2b where

v.sub.1a =800 cm.sup.-1, v.sub.1b =870 cm.sup.-1 (2.9a) V.sub.2a =885 cm.sup.-1, v.sub.2b =971 cm.sup.-1 (2.9b) The noise amplitudes in band 1 and 2, based on the point design specification of NedT=0.1 deg K at T=293 deg K, are .sigma..sub.1 and .sigma..sub.2 respectively where .sigma..sub.1 =1.211.times.10.sup.-6 watt/cm.sup.2 /sr (2.10a) .sigma..sub.2 =1.284.times.10.sup.-6 watt/cm.sup.2 /sr (2.10b) One obtains .sigma..sub.1 and .sigma..sub.2 from the formula (I=1,2) ##EQU10## where the second step in equation (2.10c) comes from the approximation for B(v,T) given by equations (2.2a,b) and (2.7a-c).

Integrating equation (2.8a) over bands 1 and 2 gives (I=1,2) ##EQU11## where equation (2.4b) has been used to go from (2.11a) to (2.11b). Equation (2.1 lb) is written in matrix form as

I=H.theta. (2.12a) where ##EQU12## The forward model with additive noise w=(w.sub.1, w.sub.2).sup.T, where w.sub.1 is the random error in band 1 and w.sub.2 is the random error in band 2, can be written as I=H.theta.+w (2.13a) From equations (2.10a,b) we see that, assuming gaussian signal noise with independent noise pulses in bands 1 and 2, the noise covariance matrix representing the measurement errors specified by w is ##EQU13## Equation (2.13a) is identical to equation (1.8) in Section 1 above (with s=0). Therefore, the C.sub.EDR covariance matrix may be written as ##EQU14## The part of C.sub.EDR of interest to the present application is the "1,1" component specifying the expected variance in the sea-surface variable .theta..sub.s =T.sub.s.sup.n. The error in T.sub.s, called .increment.t.sub.s, comes from ##EQU15## Equations (2.7c) and (2.15 b) are used in Table 3 (FIG. 11) to list .increment.T, for different values of the sea-surface temperature. The values in Table 3 run from 0.47 deg K to 0.59 deg K and are a good match to NOAA's past experience with AVHRR data (typical errors are 0.5 deg K to 0.7 deg K, see reference 14). This is not surprising because the VIIRS point design closely resembles the AVHRR specifications in the 10-12 .mu.m atmospheric window. Repeating our calculations for the AVHRR bands in the 10-12 .mu.m atmospheric window using AVHRR noise specifications, .sigma..sub.1,AVHRR =1.470.times.10.sup.-6 watt/cm.sup.2 /sr .sigma..sub.2,AVHRR =1.740.times.10.sup.-6 watt/cm.sup.2 /sr the results show that the expected temperature errors range from 0.61 deg K to 0.78 deg K (see Table 4, FIG. 12). Again, there exists a reasonable match to NOAA's past experience in retrieving sea-surface temperatures. As shown in section 6, an algorithm may be implemented to calculate the temperature errors as shown in Tables 3 and 4. A flow chart depicting the above operation is illustrated in FIG. 7. Section 3: Ice Particle Size in Cirrus Clouds

Here there is shown that the statistical methodology can be used to investigate the accuracy of different versions of an EDR algorithm. In this analysis, the processing steps involve directly linearizing the forward model to obtain an accurate error estimation rather than re-arranging into linear format as was done in the Section 2 example above. Note that the steps involved in obtaining an accurate error estimate may be implemented by means of a computer program resident on a computer host such as SPARC station or other commercially available computer systems. A flow chart depicting this method is illustrated in FIG. 8.

As is well known, the phenomenological relation between cirrus cloud temperature and ice-particle size is ##EQU16## where

L=average ice-particle size (in .mu.m)

T.sub.c =cirrus=cloud temperature (in deg K)

c.sub.0 =326.3 .mu.m

c.sub.1 =12.42 .mu.m/deg K

c.sub.2 =0.197 .mu.m/(deg K).sup.2

c.sub.3 =0.0012 .mu.m/(deg K).sup.3

Equation (3.1a) determines the ice-particle size once the cirrus clouds temperature is known. The forward model for the cirrus cloud temperature T.sub.c is I.sub.3 (v)=(1-.epsilon..sub.3 (v)).multidot.I.sub.G3 (v)+.epsilon..sub.3 (V).multidot.B(v,T.sub.c) (3.2a) for v a wavenumber (in cm.sup.-1) in the 3.7 .mu.m region of the earthshine radiance spectrum, I.sub.3 the radiance (in watt/sm.sup.2 /sr/cm).sup.1 in the 3.7 .mu.m region leaving the cloud, IG.sub.G3 the radiance (in watt/cm.sup.2 /sr/cm.sup.-1) in the 3.7 .mu.m region coming up from the earth beneath the cloud, .epsilon..sub.3, the dimensionless cloud emissivity in the 3.7 .mu.m region, and B the Planck function defined above in equations (2.2a,b). A similar equation can be written for the radiance in the 10-12 .mu.m region of the earthshine spectrum. I.sub.4 (v)=(1-.epsilon..sub.4)).multidot.I.sub.G4 (v)+.epsilon..sub.4 (v).multidot.B(v,T.sub.c) (3.2b) with I.sub.4, I.sub.G4, and .epsilon..sub.4 the same as I.sub.3, I.sub.G3, and E.sub.3 except that they apply to the 10-12 .mu.m instead of the 3.7 .mu.m spectral region. FIG. 5 shows that the conceptual basis of equations (3.2a,b) is even simpler than that used for the sea-surface temperature analysis. That is, the present model is a one layer rather than a two-layer transmission-emission model.

Equations (3.2a,b) as written show the two measured quantities I.sub.3 and I.sub.4 as depending on six unknowns: I.sub.G3, I.sub.G4, .epsilon..sub.3, .epsilon..sub.4, and T.sub.c. The values of I.sub.G3 and I.sub.G4 however, can be determined from nearby image pixels where there are no cirrus clouds present; and if .epsilon..sub.3 (v) and .epsilon..sub.4 (v) are replaced by band-averaged constant emissivities .epsilon..sub.3 and .epsilon..sub.4, there is a known functional relationship between the two emissivities as evidenced in the article by S. C. Ou, K. N. Liou, W. M. Gooch, and Y. Takano, "Remote Sensing of Cirrus Cloud Parameters Using Advanced Very-High-Resolution Radiometer 3.7- and 10-9 .mu.m Channels", Applied Optics, Vol. 32, No. 12, 1993, pp. 2171-2180, incorporated herein by reference.

(1-.epsilon..sub.3)=(1-.sub.4).sup.m(T.sbsp.c.sup.) (3.3a) where ##EQU17## Function L(T.sub.c) in equation (3.3b) is derived equation (3.1a) above. Now the forward model in equation (3.2a,b) can be approximated by I.sub.3 (v)=(1-.epsilon..sub.4).sup.m(T.sbsp.c.sup.) .multidot.I.sub.G3 (V)+[1-(1-.epsilon..sub.4).sup.m(T.sbsp.c.sup.) ].multidot.B(v,T.sub.c)(3.4a) I.sub.4 (v)=(1-.epsilon..sub.4 .multidot.I.sub.G4 (v)+.epsilon..sub.4 .multidot.B(v,T.sub.c) (3.4b)

Such approximations are validated within the following articles incorporated herein by reference. Robert P. d'Entremont, Michael K. Griffin, and James T. Bunting, "Retrieval of Cirrus Radiative Properties and Altitudes Using Multichannel Infrared Data", Fifth Conference on Satellite Meteorology and Oceanography, Boston, Mass., American Meteorology Society, Sept. 3-7, 1990, is pp. 4-9; Robert P. D'ENTREMONT, Donald P. Wylie, J. William Snow, Michael K. Griffin, and James T. Bunting, "Retrieval of Cirrus Radiative and Spatial Properties Using Independent Satellite Data Analysis Techniques:, Proceedings Sixth Conference on Satellite Meteorology and Oceanography, Atlanta, Ga., American Meteorology Society, Jan. 5-10, pp. 17-20; Robert P. d'Entremont, Donald P. Wylie, Daniel C. Peduzzi, and Joseph Doherty, "Retrieval of Cirrus Radiative and Spatial Properties Using Coincident AVHRR and HIRS Satellite Data", Passive Infrared Remote Sensing of Clouds and the Atmosphere. SPIE-The International Society for Optical Engineering, Vol, 1934, 1993, pp. 180-196; S. C. Ou, K. N. Liou, W. M. Gooch, and Y. Takano, "Remote Sensing of Cirrus Cloud Parameters Using Advanced Very-High-Resolution Radiometer 3.7- and 10.9-.mu.m Channels", Applied Optics, Vol. 32, No. 12, 1993, pp. 2171-2180; Robert P. d'Entremont and Gary B. Bustafson, "Support of Environmental Requirements for Cloud Analysis and Archive: Detection and Analysis of Cirrus Clouds Using Passive Infrared Satellite Data", Memo from Remote Sensing Group, Atmospheric and Environmental Research, Inc., 840 Memorial Drive, Cambridge, Mass., March 1996; and Algorithm Summaries, Robert d'Entremont, Jan. 22, 1997, VIIRS-40.4.6, pp. 35-37.

When equations (3.4a,b) are integrated over two bands, one in the 3.7 mm region and one in the 10-12 .mu.m region, the four measured values of .intg.I.sub.3,dv, .intg.I.sub.4,dv, .intg.I.sub.G3 dv, and .intg.I.sub.G4 dv give two equations in two unknowns, providing enough information to use the statistical methodology of Section 1. For the 3.7 .mu.m band, which we call band 3, equation (3.4a) is integrated from v.sub.3a to v.sub.3b, and for the 10-.mu.m band, which we call band 4, equation (3.4b) is integrated from v.sub.4a to v.sub.4b. ##EQU18##

We next define two unknowns .theta..sub.e and .theta..sub.c using

.epsilon..sub.4 =.epsilon..sub.04 +.theta..sub.c (3.6a) T.sub.c =T.sub.0 +.theta.c (3.6b) where .epsilon..sub.04 and T.sub.0 are typical (or average) values for, respectively, the cloud emissivity in band 4 and the cloud temperature. They are chosen so that the true values of the cloud temperature and emissivity lie reasonably close to T.sub.0 and .epsilon..sub.04, making .theta..sub.c and .theta..sub.e small quantities.

The integrals of B(v,T.sub.c) between v.sub.ia and v.sub.ib (for I=3,4) are then linearized with respect to .theta..sub.c. ##EQU19## Equations (3.7a) come from a Taylor series expansion of B(v,T.sub.c) about T.sub.c =T.sub.0, and equations (2.2a,b) are used to write the integrals in equations (3.7b,c). The next expression which has to be linearized is f(.epsilon..sub.4,T.sub.c). We calculate

f.sub.00 =f(.epsilon..sub.04, T.sub.0) (3.8a) ##EQU20## with function L(T.sub.c) defined in equation (3.1a) above. Now f(.epsilon..sub.4, T.sub.c) can be approximated as f(.epsilon..sub.4,T.sub.c).congruent.f.sub.00 +f.sub..epsilon.0 .theta..sub.e +f.sub.TO .theta..sub.c (3.8e) using .theta..sub.c and .theta..sub.e as defined in equations (3.6a,b).

The last step in the linearization process requires the definition of two more variables .theta..sub.G3 and .theta..sub.G4.

I.sub.G3 =I.sub.0G3 +.theta..sub.G3 (3.9a) I.sub.G4 =I.sub.0G4 +.theta..sub.G4 (3.9b) In effect, the determination of I.sub.G3 and I.sub.G4 from nearby cloudless pixels is treated as two more measurements. The values of I.sub.0G3 and I.sub.0G4 are chosen to make .theta..sub.G3 and .theta..sub.G4 small, so they must be typical (or average) values of the I.sub.G3, I.sub.G4 radiances.

Equations (3.5a,b) and (3.9a,b) are combined to obtain the complete set of forward equations using (3.6b), (3.7a), and (3.8e) to represent .epsilon..sub.4, f, and the integrals over B.

I.sub.3 =(f.sub.00 f.sub..epsilon.0 .theta..sub.e +f.sub.TO .theta..sub.e)(I.sub.0G3 +.theta..sub.G3)+[1-f.sub.00 -f.sub..epsilon.0 .theta..sub.e -f.sub.TO .theta..sub.c ](b.sub.3 +.beta..sub.3 .theta..sub.c) (3.10a) I.sub.4 =(1-.epsilon..sub.04 -.theta..sub.e)(I.sub.0G4 +.theta..sub.G4)+(.epsilon..sub.04 +.theta..sub.e)(b.sub.4 +.beta..sub.4 .theta..sub.c) (3.10b) I.sub.G3 =I.sub.0G3 +.theta..sub.G3 (3.10c) I.sub.G4 =I.sub.0G4 +.theta..sub.G4 (3.10d)

Equations (3.10a-d) are linearized by keeping terms of 0(.theta..sub.c), 0(.theta..sub.e), 0(.theta..sub.G3), 0(.theta..sub.G4); terms of 0(.theta..sub.e .theta..sub.c), 0(.theta..sub.e,.theta..sub.G3), 0(.theta..sub.e .theta..sub.G4), 0(.theta..sub.c .theta..sub.G3), 0(.theta..sub.c .theta..sub.G4), 0(.theta..sub.G3 .theta..sub.G4) are neglected, as are terms of 0(.theta..sub.c.sup.n), 0(.theta..sub.e.sup.n), 0(.theta..sub.G3.sup.n), 0(.theta..sub.G4.sup.n) for n>1. The method followed in the linearization process is that terms which are constant or linear in the .theta.'s are kept, and terms which depend on higher powers of the .theta.'s are dropped.

I.sub.3 =.theta..sub.c [.epsilon.(I.sub.0G3 -b.sub.3)+.beta..sub.3 (1-f.sub.00)]+.theta..sub.e [f.sub..epsilon.0 (I.sub.0G3 -b.sub.3)]+.theta..sub.G3 [f.sub.00 ]+[f.sub.00 (I.sub.0G3 -b.sub.3)+b.sub.3 ] (3.11a) I.sub.4 =.theta..sub.c [.beta..sub.4 .epsilon..sub.04 ]+.theta..sub.e [b.sub.4 -I.sub.0G4 ]+.theta..sub.G4 [1-.epsilon..sub.04 ]+[I.sub.0G4 (1-.epsilon..sub.04)+b.sub.4 .epsilon..sub.04 ] (3.11b) I.sub.G3 =I.sub.0G3 .theta..sub.G3 (3.11c) I.sub.G4 =I.sub.0G4 .theta..sub.G4 (3.11d) Equations (3.11a-d) can be written in matrix form as ##EQU21## To check the accuracy of equations (3.12a-e) we represent bands 3 and 4 by the VIIRS point design for the 3.7 .mu.m band and one of the 10-12 .mu.m bands, as shown in Table 5. A noise vector w=(w.sub.1,W.sub.2,W.sub.3,W.sub.4).sup.T of zero-mean gaussian random numbers is added to the right-hand side of (3.12a). I=H.theta.+s+w (3.13a) The noise amplitudes in bands 3 and 4 are .sigma..sub.3 and .sigma..sub.4 respectively (see Table 5 FIG. 13), so the covariance matrix for the w vector is ##EQU22## Equation (3.13b) assumes the I.sub.3, I.sub.G3 measurements have the same .sigma..sub.3 noise amplitude; it assumes that the I.sub.4, I.sub.G4 measurements have the same .sigma.4 noise amplitude; and it assumes the errors in all four measurements I.sub.3, I.sub.4, I.sub.G3, I.sub.G4 are uncorrelated. Equation (3.13a) has the form given in equation (1.8) above, so equation (1.9a) can be used to calculate C=(H.sup.T C.sup.-1 H).sup.-1 to get the expected error in our estimate of .theta..sub.c, the cirrus cloud temperature. As will be shown below, and algorithm (see section 7 may be implemented which uses (1.9a) with the constants and matrices specified in Table 5 and equations (3.12e), (3.13b) to calculate .increment.T.sub.c =0.836 deg K when T.sub.0 =230 deg K, .epsilon..sub.04 =0.4, and the values of I.sub.0G3, I.sub.0G4 are the band 3 and 4 radiances of a black body at 290 deg K. The corresponding error in the ice-particle size of the cirrus cloud is given by ##EQU23## Table 6 (FIG. 14) is derived from the last page of Section 7 and shows that for a range of T.sub.c values running between 220 deg K and 240 deg K the ice-particle size L runs between 41 .mu.m and 88 .mu.m, the absolute error in L is between 1.4 .mu.m and 2.8 .mu.m, and the relative error is between 2.9% and 3.3%. These values are a reasonable match to the threshold requirements of the ice-particle size EDR (5% relative error and 2 .mu.m absolute error for ice particles between 0 and 50 .mu.m in size, see reference 15). Again, as in Section 2, the calculations have produced believable results.

This statistical methodology can also be used to examine the effect of modeling error. If the forward model for predicting the values of I.sub.3 and I.sub.4 in equations (3.10a,b) does not predict exactly the I.sub.3, I.sub.4 atmospheric radiances, this can be regarded as another source of random error and include the amount by which the model is incorrect in the noise vector w. The modeling error can be treated as measurement "noise" because what is desired of the instruments to measure is the I.sub.3, I.sub.4 values predicted by the inaccurate model of (3.10a,b). Instead, what is obtained is the radiance coming from the unknown but accurate model, i.e. physical reality. There is no way of knowing whether the modeling error in any given measurement acts to compensate for the measurement noise, or to make the discrepancy worse, so .sigma..sub.mod 3 and .sigma..sub.mod 4 is taken as the modeling error in bands 3 and 4 respectively, to be independent sources of error obeying a zero-mean gaussian distribution. The covariance matrix C is now written as ##EQU24## The I.sub.G3, I.sub.G4 errors have not been changed because these are the radiances reaching the bottom of the cirrus cloud layer and are measured directly from nearby cloudless pixels. They cannot be affected by modeling error. They might, however, have somewhat different values in the nearby pixels from what they in fact are in the cirrus cloud pixels. We can, therefore, specify that .sigma..sub.G3 and .sigma..sub.G4 are the errors in I.sub.G3, I.sub.G4 respectively coming from the use of nearby cloudless pixels to estimate the values of I.sub.G3, I.sub.G4 in the cirrus cloud pixels. Including them in the formula for the covariance matrix C gives ##EQU25## Estimating reasonable values for a .sigma..sub.mod 3, .sigma..sub.mod 4, .sigma..sub.G3, .sigma..sub.G4 and then using (3.15) to specify matrix C in equation (1.9a) gives us a quick and easy way of including the modeling and procedure errors in the error budget for the ice-particle and cirrus temperature EDR's. Many of the other VIIRS-based EDR's can have their modeling and procedure errors including in their C matrices in a similar way.

Note also how easy it is to find a performance improvement of an EDR algorithm caused by the addition of an extra measurement to the forward model. Suppose the suggestion is made to add a 6.7 .mu.m measurement, which shall be called band 5, to the cirrus EDR algorithm. Note that the relationship between .epsilon..sub.5, the cirrus emissivity in band 5, and .epsilon..sub.4, the cirrus emissivity in band 4 is given.

(1-.epsilon..sub.5)=(1-.epsilon..sub.4).sup.n(T.sub.c) (3.16) where n(T.sub.c) is a specified function of the cloud temperature T.sub.c. We define, following the procedure used on f(.epsilon..sub.4,T.sub.c), that (see equations (3.8a-e)) g(.epsilon..sub.4,T.sub.c)=(1-.epsilon..sub.4).sup.m(T.sbsp.c.sup.)(3.17a) g.sub.00 =g(.epsilon..sub.04, T.sub.0) (3.17b) ##EQU26## Following the pattern of equation (3.11a), the two equations shown below are combined, equations (3.18a,b), with equations (3.11a-d) to get the new forward model. I.sub.5 =.theta..sub.c [g.sub.TO (I.sub.0G5 -b.sub.5)+.beta..sub.5 (1-g.sub.00)]+.theta..sub.e [g.sub..epsilon.0 (I.sub.0G5 -b.sub.5)]+.theta..sub.G5 [g.sub.00 ]+[g.sub.00 (I.sub.0G5 -b.sub.5)+b.sub.5 ] (3.18a) I.sub.G5 =I.sub.0G5 .theta..sub.G5 (3.18b) We next define b.sub.5, .beta..sub.5 in terms of the wavenumber limits v.sub.5a, v.sub.5b of band 5. ##EQU27## The in-band signals for the band 5 radiances coming from pixels with and without the cirrus clouds are, respectively, I.sub.5 and I.sub.G5. The value of I.sub.0G5 is chosen to make .theta..sub.G5 small. Augmenting H, .theta., and s to include equations (3.18a,b) gives ##EQU28## The matrix C corresponding to equation (3.13b) is ##EQU29## where .sigma..sub.5 is the noise amplitude of the band 5 signal. Modeling and procedure errors can be added to (3.19d) the same way they were added to equation (3.15). Equation (1.9a) can now be used to find .increment.t.sub.c, the expected error in the cirrus cloud temperature, and equation (3.13c) then gives .increment.L, the expected error in the size of the cirrus ice particles. Section 4: Computer-Based Forward Models

Some forward models are too complex to be easily written down algebraically. The models connecting atmospheric temperature and water-vapor profiles with the CrIs infrared spectra are one obvious example, and the VIIRS-based ocean-color and aerosol EDR's may be another (see reference 16). Typically, what does exist in these cases is a complex, computer-based forward model going from the EDR inputs to a collection of output radiances. Fortunately, this is all the information needed to find matrix H. We outline below how this is down, defining the variables in terms of the CrIS-based EDR's because the best-known forward models for these EDR's, such as MODTRAN or FASCODE, are entirely computerized.

The computer-based forward model is represented by the vector ##EQU30## where C is the total number of radiance outputs of the forward model. For the CrIS-based EDR's, it is the total number of FTIR channels in the measured infrared spectra (C is approximately 1500 for the current point design). The forward model inputs are represented by the vector ##EQU31## where S is the total number of these inputs. For the CrIS-based EDR's it is about twice the total number of atmospheric layers because each layer is characterized by both a temperature and a water-vapor concentration (S is typically about 70 for this type of forward model). An exact representation of the forward model can be written as

f=f(T) (4.3a) or f.sub.i =f.sub.i (T)=f.sub.i (T.sub.1, T.sub.2, . . . ,T.sub.5)(4.3b) for I=1,2, . . . ,C.

The linearized forward equations are ##EQU32## where T.sub.0 is the collection of input values for some standard atmosphere which is a good approximation for the true temperature and water-vapor profiles. ##EQU33## We calculate f.sub.0, the vector of forward model outputs for the standard atmosphere T.sub.0.

f.sub.0 =f(T.sub.0) (4.5a) or ##EQU34## We define a new vector .delta.T.sup.[k] to be an S-dimensional column vector which is zero everywhere except for the d'th component which is .delta.T.sub.k. ##EQU35## We run the forward model for T.sub.0 +.delta.T.sup.[k], giving us a vector of outputs which we call f.sup.[k]. f.sup.[k] =f(T.sub.0 +.sigma.T.sup.[k]) (4.7a) or ##EQU36## Vector f.sub.0 represents the output of the computer-based forward model for the standard atmosphere inputs T.sub.0, and f.sup.[k] is the output of the forward model with all the standard atmosphere inputs left unchanged except for the k'th input, which has been changed by an amount .delta.T.sub.k. To create a complete collection of f.sup.[k] vectors the forward model must be run S times.

Returning to equation (4.4a) above, we recognize that it can be written as

f=H.theta.+f.sub.0 (4.8a) where .theta.=T-T.sub.0, (4.8b) and H.sub.ij, the ij'th component of matrix H, is given by ##EQU37## Equations (4.8a-c) represent a restatement of equations (1.3b-d) in Section 1 above. The forward model predictions f in equation (1.8a) have no errors; but the actual CrIS spectra, of course, will. The measurement errors are represented in the usual way as an additive random vector w on the right-hand side of equation (4.8a), giving us f.sub.measured =H.theta.+s+w (4.8d) where s=f.sub.0. The random vector w is characterized by a covariance matrix C coming from the point design of the CrIS. Equation (4.8d) is the same as equation (1.8), so as soon as matrix H is known, equation (1.9a) can be used to calculate the expected errors in the retrieved values of T.sub.0 +.theta.. From (4.8b) we see that ##EQU38## where (f.sup.[k] -f.sub.0).sub.I is the i'th component of the vector difference (f.sup.[k] -f.sub.0). Equation (4.9) shows that the k'th column of H can be calculated by subtracting f.sub.0, the standard atmosphere outputs of the forward model, from f[k], the forward-model outputs when only the k'th input of the standard atmosphere has been changed by .delta.T.sub.k. This procedure requires the forward model to be run S+q times, once to calculate f.sub.0 and S times to calculate f.sup.[k] for k=1,2, . . . ,S.

For the CrIS-based EDR's, such as the temperature and water-vapor profiles, there usually exists a matrix C.sub..theta. describing the covariance of the expected EDR values between different levels of the atmosphere. Therefore, equation (1.11b) should be used to estimate the errors in these EDR's. In general, a covariance matrix C.sub..theta. describes what is known about the statistics of the desired EDR's before any measurements are made. By making conservative or aggressive assumptions about C.sub..theta., we can determine conservative or aggressive values for the expected errors in these Cr1S-based EDR's. When no assumptions are made about the EDR values, the C.sub..theta. matrix can be set to zero and equation (1.9a) used to find the expected errors in the retrieved EDR's. Equation (4.9) can be used to find H for any computer-based forward model; accordingly, the above detailed steps are also applicable to VIIRS-based EDR's.

Section 5

Suppose we have two straight lines written in functional form as f.sub.1 (v) and f.sub.2 (v) where

f.sub.1 (v)=m.sub.1 v+b.sub.1, f.sub.2 (v)=m.sub.2 v+b.sub.2, and we want to see whether f.sub.1 (v).congruent.A.multidot.f.sub.2 (v) for some constant A. One quick way of doing this is to see whether ##EQU39## because m.sub.1 v+b.sub.1 =A.multidot.(m.sub.2 v+b.sub.2) can only be true (for all values of v) when m.sub.1 =Am.sub.2 and b.sub.1 =Ab.sub.2, or ##EQU40##

Section 6

This algorithm is used to calculate the expected error in the sea surface temperature measured by satellite IR radiometers operating in the 800 to 1000 cm.sup.-1 atmospheric spectral window and utilizes formulas on pages 97 and 391 of Steven M. Kay's textbook Fundamentals of Statistical Signal Processing: Estimation Theory, PTR Prentice Hall, Englewood Cliffs, N.J. 1993, ISBN 0-13-345711-7.

Start by defining constants used in the calculations. Variable n is the power of v (wavenumber) over temperature used to match (exp(x)-1) (-1) in the Planck function for wavenumbers between 800 and 100 cm.sup.-1 and temperatures between 310 and 290 deg K. Variable CB is a constant used to normalize the measured inband radiances when applying the model equations. It has units of watt-cm2/sr/(cm-deg K) n. The atmospheric transmittance is modeled as 1-K(v)*A where K(v) is a straight line of slope M in units of cm. (This transmittance model assumes that the radiometer measurements are made over a wide enough band that we can neglect the rapid up-and-down variation of the atmospheric transmittance with respect to v.)

N:=4.4 CB:=1.98.multidot.10.sup.-12 M:=8.0.multidot.10.sup.-4 Define the F4 and F5 functions which will be used to calculate the elements of matrix H as defined on pages 97 and 391 of Steven Kay's text. Variables va and vb are, respectively, the beginning and ending wavenumbers of the spectral band being measured by the satellite radiometer. ##EQU41## Specify two bands for the radiometer measurement, band 1 between 800 and 870 cm.sup.-1 and band 2 between 885 and 971 cm.sup.-1. Variable .sigma.1 is the standard deviation of the measurement error in band 1 in units of watt/cm2/sr, and .sigma.2 is the standard deviation of the measurement error in band 2 in units of watt/sm2/sr. We set up matrices C and H as defined on pages 97 and 391 of Steven Kay's book. ##EQU42## Utilize the formula on page 97 of Steven Kay's book to calculate the expected error in the estimate of TS n, the n'th power of the sea-surface temperature in units of degK n. Matrix Cu is the covariance matrix for the estimate of Ts n and the unknown atmospheric transmittance parameter. ##EQU43## The corresponding error in the measured sea-surface temperature is .sigma.Ts, its standard deviation in degK. It is expected that the sea-surface temperature lies between 290 and 310 degK, so we calculate a table of corresponding .sigma.Ts values. ______________________________________Ts.sub.0 : = 290. Ts.sub.1 : = 295. Ts.sub.2 : = 300. Ts.sub.3 : = 305.Ts.sub.4 : = 310. i: = 0..4 ##STR1## ##STR2## Ts.sub.i .sigma.Ts.sub.i______________________________________ 290 0.7672964 295 0.7239716 300 0.6837606 305 0.6463934 310 0.6116269______________________________________

Section 7

This algorithm is used to estimate the error involved in finding the ice-particle size of cirrus clouds using VIIRS measurements (based on the 1997 NPOES proposal) in the 3.7 micron (band 3) and 10.9 micron (band 4) radiance bands. The analysis presented below is based on information about the dependence of the bands 3 and 4 radiances in a paper "Remote Sensing of Cirrus Cloud Parameters Using Advanced Very-High-Resolution Radiometer 3.7 and 10.9 .mu.m Channels" in Applied Optics, Vol. 32, No. 12, April 1993 by S. C. Ou, K. N. Liou, W. M. Gooch, and Y. Takanco, pp. 2171-2180, (see especially eqns. (1a,b) on p. 2172 and eqn. (4) on p. 2173), and incorporated herein by reference. The statistical analysis is based on eqns. (4.30) and (4.32) on page 97 of Fundamentals of Statistical Signal Processing: Estimation Theory by Steven M. Kay, PTR Prentice Hall, Englewood Cliffs, N.J. 1993 (ISBN 0-13-345711-7).

Create a function called PLnck(x) which calculates 1(exp(x)-1). When x is close to zero it replaces exp(x)-1 and returns 1/x exactly. ##EQU44## PLnck(x):=if(.vertline.x.vertline.>0.01, PLnck(x), PLnck2(x)) Generate C1 and C2 black body constants used to specify Planck's curve. Variable C1 has units of watt/sm2/sr/(cm-1)4 and variable C2 has units of degK/cm-1. C1:=1.191062.multidot.10.sup.12 C2:=1.multidot.438786 Specify Ibb.sub.-- Winvcminvsr as a function of wavenumber a (in cm-1) and temperature T (in degK). The units of Ibb.sub.-- Winvminvsr are Watts/cm/sr. ##EQU45## Specify IdT.sub.-- WinvcminvsrinvdegK, the derivative of Ibb.sub.-- Winvcminvsr, as a function of wavenumber .sigma. (in cm-1) and temperature T (in degK). The units of IdT.sub.-- WinvcminvsrinvdegK are W/cm/sr/degK. ##EQU46## Specify the beginning and ending wavenumbers (in cm-1) of band 3, called v3a and v3b respectively, and the beginning and ending wavenumbers of band 4, called v4a and v4b respectively. v3a.sub.-- invcm:=2545 v3b.sub.-- invcm:=2817 v4a.sub.-- invcm:=885 v4b.sub.-- invcm:=971 Variable nedt3.sub.-- degK is the NEDT (in degK) expected for the band 3 point design of VIIRS in the 1997 NPOES proposal. Variable nedt4.sub.-- degK is the NEDT (in degk) for the band 4 point design of VIIRS in the 1997 NPOES proposal. Tnedt3.sub.-- degK and Tnedt4.sub.-- degK are the band 3 and band 4 scene temperatures at which the NEDT values are evaluated. Nedt3.sub.-- degK:=0.2 Tnedt3.sub.-- degK:=271.0 nedt4.sub.-- degK:=0.1 Tnedt4.sub.-- degK:=293.0 Variable .sigma.3.sub.-- Winvcm2invsr is the expected measurement error (standard deviation of the additive noise) in band 3 in units of Watt/sm2/sr, and .sigma.4.sub.-- Winvcm2invsr is the expected measurement error (standard deviation) in band 4 in units of Watt/sm2/sr. ##EQU47## Calculate the constants a3.sub.-- Winvcm2invsr, .beta.3.sub.-- Winvcm2invsrinvdegK and a4.sub.-- Winvcm2invsr, .beta.4.sub.-- Wincm2invsrinvdegK which are used to linearize with respect to temperature the integrals over bands 3 and 4 of the Planck black-body function. The integrals are linearized about temperature TO.sub.-- degK (in units of degK). TO.sub.-- degK:=230.0 ##EQU48## Define function L.sub.-- .mu.m m, the cirrus cloud ice-particle size (in microns) as a function of cloud temperature Tc (in deg K). The best fit constants c0.sub.-- .mu.m, c1.sub.-- .mu.minvdegK, c2.sub.-- minvdegK2, and c3.sub.-- minvdegK3 come from curve fitting of cloud data. C0.sub.-- .mu.m:=326.3 c1.sub.-- .mu.minvdegK:=12.42 c2.sub.-- .mu.minvdegK2:=0.197 c3 .mu.minvdegK3:=0.0012 L.sub.-- .mu.m(tc):=c0.mu.m+c1.sub.-- .mu.minvdegK.multidot.(Tc-273)+c2.sub.-- .mu.minvdegK2 .multidot.(Tc-273).sup.2 +c3.sub.-- .mu.minvdegK3 (Tc-273).sup.3 Define function m (dimensionless), the exponent used in the forward radiance model for bands 3 and 4. It depends on L.sub.-- .mu.m, the cirrus cloud ice-particle size in microns. The constants a0s, a2s.sub.-- .mu.m, and A2s.sub.-- .mu.m2 define the dependence of m on L.sub.-- .mu.m. a0s:=0.722 a1s.sub.--.mu.m:= 55.08 a2s.sub.-- .mu.m2:=-174.12 ##EQU49## Define the derivative of function m with respect to its argument L.sub.-- .mu.m, mdL.sub.-- inv.mu.m (in units of .mu.m-1) ##EQU50## Define the derivative of L.sub.-- .mu.m with respect to its argument Tc, LdT.sub.-- .mu.minvdegK (in units of .mu.m/degK) LdT.sub.-- .mu.minvdegK(Tc):=c1.sub.-- .mu.minvdegK+2.multidot.c2.sub.-- .mu.minvdegK2.multidot.(Tc-273)+3.multidot.c3.sub.-- .mu.minvdegK3.multidot.(Tc-273).sup.2 Define the derivative of m with respect to the cirrus cloud temperature Tc, mdT.sub.-- invdegK (in degK-1). MdT.sub.-- invdegK(Tc):=mdL.sub.-- inv.mu.m(L.sub.-- .mu.m(Tc)).multidot.LdT.sub.-- .mu.minvdegK(Tc) Specify .epsilon.04, the dimensionless average emissivity of the cirrus clouds. .epsilon.04:=0.4 Define the values of the constants f00, f.sub..epsilon., 0, and fT0.sub.-- invdegK. F00:=(1-.epsilon.04).sup.m (T0.sub.-- degK) F00=0.586773 f.epsilon.0:=-m(T0.sub.-- degK).multidot.(1-.epsilon.04).sup.m(T0.sbsp.--.sup.degK)-1 f.epsilon.0=1.020631 fT0.sub.-- invdegK:=(1-.epsilon.04).sup.m(T0.sbsp.--.sup.degK) .multidot.In(1-.epsilon.044).multidot.mdT.sub.-- invdegK(T0.sub.-- degK) fT0.sub.-- invdegK=0.003404 Create typical values for the amount of in-band radiance in bands 2 and 3 reaching the cirrus clouds from the ground by estimating the surface temperature as Ts.sub.-- degK (in degK) and integrating the Planck function at this temperature over bands 3 and 4. Variable IOG3.sub.-- Winvcm2invsr and IOG4.sub.-- Winvcm2invsr are these in-band radiances for bands 3 and 4 respectively (in units of watts/cm2/sr). ##EQU51## Define the H matrix used on page 97 of Steven Kay's text (referenced at start of this spreadsheet). H11:=fT0.sub.-- invdegK-(I0G3.sub.-- Winvcm2invsr-a3.sub.-- Winvcm2invsr)+.beta.3.sub.-- Winvcm2invsrinvdegK.multidot.(1-f00) H12:=f.epsilon.0.multidot.(I0G3.sub.-- Winvcm2invsr-a3.sub.-- Winvcm2invsr) H13:=f00 H21:=.beta.4.sub.-- Winvcm2invsrinvdegK.multidot..epsilon.04 H22:=a4.sub.-- Winvcm2invsr-I0G4.sub.-- Winvcm2invsr H24:=I-.epsilon.04 ##EQU52## Specify the noise covariance matrix C specified in page 97 of Steven Kay's text (referenced at the beginning of the spreadsheet) ##EQU53## Calculate the covariance matrix of the errors in the determination of unknown parameters: the cirrus cloud temperatures, the emissivity of the cirrus in band 4, the deviation of the radiance reaching the cloud from the ground from its nominal values in band 3, and the same radiance deviation in band 4. Here we are only interest in the error in the estimate of the cirrus cloud temperature, whose variance will be given by the "1,1" value of the error covariance matrix C.theta.. C.theta.:=(H.sup.T .multidot.C.sup.-1 .multidot.H).sup.-1 ##EQU54## The standard deviation of the error in the cirrus temperature (in degK) is .sigma.Tc.sub.-- degK. ##EQU55## Use the previously-defined function L.sub.-- .mu.m, which depends on the cirrus cloud temperature in degK, to create a table of the error in L.sub.-- .mu.m expected given the above error .sigma. Tc.sub.-- degK in the estimated cloud temperatures. ______________________________________i: = 0,,21 Tc.sub.i : = 220 + (i - 1) .multidot. 1 Lmicron.sub.i : =L.sub.-- .mu.m(Tc.sub.i)Lmicron.sub.-- err.sub.i : = LdT.sub.-- .mu.minvdegK(Tcp.sub.i).multidot. .sigma.Tc.sub.-- degK ##STR3##Lmicron.sub.i Lmicron.sub.-- err.sub.i Rel.sub.-- errL.sub.i______________________________________41.1152 1.373158 0.03339842.7606 1.380519 0.03228544.4184 1.393903 0.03138146.0958 1.413309 0.0306647.8 1.438738 0.03009949.5382 1.47019 0.02967851.3176 1.507664 0.02937953.1454 1.55116 0.02918755.0288 1.60068 0.02908856.975 1.656222 0.02906958.9912 1.717786 0.02911961.0846 1.785374 0.02922863.2624 1.858984 0.02938565.5318 1.938616 0.02958367.9 2.024271 0.02981370.3742 2.115949 0.03006772.9616 2.213649 0.0303475.6694 2.317372 0.03062578.5048 2.427118 0.03091781.475 2.542886 0.03121184.5872 2.664677 0.03150287.8486 2.79249 0.031788______________________________________ Conclusion

The statistical methodology outlined in Section 1, and applied to an algebraic forward model for the sea-surface temperature in Section 2 and an algebraic forward model for the size of cirrus ice particles in Section 3 provides a measure of the accuracy in which a satellite may yield parameters of the earth's surface. In Section 2 the predicted EDR performance matched NOAA's past experience for two-band sea-surface temperature retrieval, and in Section 3 the predicted performance matched NOAA's threshold requirements. Thus, in both cases, the statistical methodology led to reasonable results. The performance of an augmented EDR algorithm can also be determined when an extra band is used to provide more spectral information, as well as how the effect of modeling error could be easily included in the covariance matrix describing the measurement uncertainties. Section 4 showed how to obtain the matrix elements of the linearized forward model (which the statistical methodology requires) when the forward model is only available as a computer program. The forward model for the CrIS-based EDR's, as well as some of the more complicated VIIRS-based EDR's, may in fact, only exist in computerized form.

Current techniques for deriving VIIRS-based EDR's seem to rely on creating forward models consisting of N equations in N unknowns and then solving for one or more of the unknowns to create the desired algorithm. Clearly, such a procedure cannot be applied to forward models consisting of N equations in M unknowns when N>M; yet the information provided by the extra equations might well lead to more accurate EDR algorithms. This sort of situation, where the forward model has more equations than unknowns, is exactly the sort of situation which is easily handled by the statistical methodology according to the present invention. The current VIIRS-based EDR algorithms also have no way of using correlations between the M unknowns of the forward mode, and this is additional extra information which can be easily handled by this methodology. While EDR algorithms based on the methodology will have come from linearized forward models and so may be less accurate than EDR algorithms based on the full non-linear forward models, their ability to use extra information may provide a performance advantage over their non-linear counterparts, especially when the desired EDR's exhibit only moderate variations about a climatic mean. For these reasons, algorithms derived from equations (1.9b) and (1.11 c,d) in Section 1 may well provide a significant advance over the current state of the art in VIIRS-based EDR retrieval.

Moreover, the statistical method according to the present invention provides a quick and easy means of estimating the expected performance of an EDR algorithm. Current methods for determining performance envisions a series of Monte-Carlo, end-to-end simulations based on combinations of EDR forward models with the simulated sensors and the proposed EDR algorithms. The inputs of the forward models would be changed randomly and compared again and again to the outputs coming from the EDR algorithms. This is a truly ambitious program of number crunching and may, in the end, provide only imprecise knowledge of the EDR algorithms' performance, especially if only a relatively few Monte-Carlo runs are performed. Although some type of end-to-end Monte-Carlo simulation will be required to verify performance predictions achieved by other means (such as the statistical methodology discussed in Sections 1-4 above), this sort of simulation is at best a clumsy way of examining different algorithms for the same EDR. As we have shown in Sections 1-4, the present statistical methodology only needs an EDR's linearized forward model to predict the EDR algorithm's behavior; compared to an end-to-end simulation it well more quickly and easily suggest which version of an EDR algorithm has the best performance. When it becomes time to perform an end-to-end simulation, this methodology will, by having already examined the algorithm's behavior, provide an approximate answer to the simulation's results, helping to avoid gross errors in the computer calculations.

Claims

1. A method of estimating expected errors of environmental data parameters based on radiance measurements obtained from visible infra red radiometric satellite sensors (VIIRS) orbiting the earth, comprising the steps of:

obtaining N radiance measurements of a surface body I.sub.I, . . . I.sub.N defining a matrix I depending on p unknown surface and atmospheric parameters T.sub.i, . . . T.sub.p, defining a matrix T;

generating a forward model I=f(T) for obtaining said I radiance measurements from said p parameters, where f(T) comprises an N element matrix;

choosing an initial set of values for said p parameters and linearizing f(T) about said initial values to obtain a linearized forward model I=s+H.theta. as f.sub.i (T.sub.01, T.sub.02, . . . T.sub.0p)+.SIGMA.H.sub.ij .theta..sub.j where I=1,2, . . . N and .theta..sub.j =T.sub.j -T.sub.oj and H.sub.ij =.delta..sub.i f/.delta.T.sub.j and where .theta. is a column matrix .sub.1 .theta., . . . .sub.p .theta. and H is a matrix of H.sub.ij values;

adding measurement noise vector w of noise values to said forward model;

determining the covariance of said measurement noise w to obtain a covariance matrix C; and

manipulating the matrices H and C according to the equation C.sub.EDR=(H.sup.T C.sup.-1 H).sup.-1 to obtain a matrix element of C.sub.EDR indicative of the expected errors in said values of T.sub.1. . . T.sub.p parameters.

2. The method according to claim 1, wherein said measurement noise w comprises N random numbers w.sub.1, . . . ,w.sub.N.

3. The method according to claim 2, wherein said N random numbers have a jointly gaussian distribution probability.

4. The method according to claim 3, wherein said jointly gaussian distribution probability is a zero mean gaussian distribution.

5. The method according to claim 1, wherein said matrix element C.sub.EDR specify expected errors in the values .theta..sub.1, . . . .theta..sub.p.

6. The method according to claim 5, wherein the expected errors in the values .theta..sub.1, .theta..sub.2, . . . .theta..sub.p represent the differences in said p surface and atmospheric parameters from said initial set of values.

7. A method of predicting errors in sea surface temperature measurements based on radiance measurements I obtained from visible infra red radiometric satellite sensor (VIIRS) orbiting the earth having at least a first bandwidth and a second bandwidth corresponding to a particular window region, comprising the steps of:

applying a forward model for determining sea-surface temperature from said satellite radiance measurements I over said particular wavelength window region as I(v,atm)=B(v,T.sub.s).multidot.t(V,atm)+B(v,T.sub.a).multidot.[1-t(v,atm)] where T.sub.s is measured sea-surface temperature, Ta represents a temperature indicative of atmospheric self-radiance in said wavelength window region, t represents a parameter indicative of transmissivity from said sea surface to the top of the atmosphere, v is a wavenumber in said window region, and atm is representative of atmospheric parameters determining the value of said transmittance at wavenumber v, and where B(V,T) is a Planck function;

approximating said transmissivity parameter t(v,atm) as 1-A.sub.atm .multidot.K(v) in said window region where K is a function only of wavenumber and where A.sub.atm is indicative of atmospheric parameters;

performing a straight line approximation of K(v) to obtain K.sub.L (V);

evaluating B(V,T) over a limited range of V and T parameter values to obtain I(V,atm)=C.sub.B is a constant and .theta..sub.s =(T.sub.s).sup.n, and .theta..sub.a =A.sub.atm (T.sub.a.sup.n -T.sub.s.sup.n);

measuring noise amplitudes over said first and second bandwidths to obtain a first noise amplitude .sigma..sub.1, and a second noise amplitude .sigma..sub.2 ;

integrating said radiance function measurements I(v,atm) over said first and second bandwidths to obtain said radiance I as a function of a matrix H and a matrix .theta. such that I=H.theta., where I is a matrix function of radiances I1 and I2 corresponding to said first and second bandwidths, and .theta. is a matrix representation of .theta..sub.s and .theta..sub.a, and where H is a function of V;

adding to said radiance I a noise matrix w corresponding to said noise amplitudes to obtain a radiance I=H.theta.+w;

determing the covariance of said noise matrix w to obtain a covariance matrix C;

and manipulating the matrices H and C according to the equation (H.sup.T C.sup.-1 H).sup.-1 to obtain a matrix C.sub.EDR having an element indicative of the expected variance in said sea surface temperature variable T.sub.s.sup.n.

8. The method according to claim 7, wherein said wavelength window region is 10-12 .mu.m.

9. The method according to claim 8, wherein said T.sub.a varies by less than 1 Kelvin within said 10-12 .mu.m window.

10. The method according to claim 9, wherein said limited range over which B(V,T) is evaluated further comprises a temperature range of between 285 K and 315 K.

11. The method according to claim 10, wherein said limited range over which B(V,T) is evaluated further comprises a wavenumber range of between 800 cm.sup.-1 and 1000 cm.sup.-1.

12. The method according to claim 11, wherein said Planck function B(V,T) is C.sub.1 V.sup.3 P(C.sub.2 VT.sup.-1) where C.sub.1 is 1.191.times.10.sup.-12 watt/cm.sup.2 /sr/(cm.sup.-1).sup.4, and where P(C.sub.2 VT.sup.-1)=1/(e.sup.c2VT-1 -1).

13. The method according to claim 7, wherein the step of approximating t(v,atm) as 1-A.sub.atm .multidot.K(v) is applied only to wide bandwidth radiance signals.

14. A method estimating an error involved in determining ice particle size of cirrus clouds using radiance measurements obtained from visible infrared radiometric satellite sensors (VIIRS) in a first and second radiance bands, said VIIRS sensors orbiting the earth and obtaining said radiance measurements from a body located on a surface of said earth, comprising the steps of:

creating a Planck function P(x) such that the value of said P(x) is 1/x when x approaches 0, and is 1/(e.sup.x -1) otherwise;

determining first and second black body constants C1 and C2 associated with specifying a curve of said Planck function P(x);

determining a first parameter as a function of wavenumber .sigma. and temperature .tau.;

determining a second parameter as a function of wavenumber .sigma. and temperature .tau.;

determining beginning and ending wavenumbers associated with said respective first and second radiance bands;

storing an expected temperature error NEDT expected for said first band associated with the design of said VIIRS sensor and a second expected temperature error NEDT4 expected for said second band associated with the design of said VIIRS sensor;

storing first and second temperature measurement values TNEDT and TNEDT4, respectively, with each of said first and second respective radiation bands over which said expected temperature errors NEDT and NEDT4 are evaluated;

obtaining an expected measurement error .sigma..sub.3 in said first band by integrating said first parameter as a function of said first expected temperature error and first temperature measurement value over said first radiance band and obtaining a second expected measurement error .sigma..sub.4 associated with said second radiance band by integrating said first parameter as a function of said second expected temperature error and said second temperature measurement value over said second radiance band;

determining a set of constants A3, B3, A4, and B4 for linearizing with respect to temperature integrals over said first and second bands of said Planck black body function P(x), wherein said constant A3 comprises integrating said first parameter as a function of temperature T.sub.0 over said first radiance band;

wherein said constant B3 is evaluated by integrating said second parameter at the temperature T.sub.0 over said first radiance band;

wherein said A4 constant evaluated by integrating said first parameter at said temperature T.sub.0 over said second radiance band;

and wherein said constant B4 is evaluated by integrating said second parameter at said temperature T.sub.0 over said second radiance band;

defining a cirrus cloud ice particle size parameter L as a function of cloud temperature T.sub.c based on a plurality of constants c0, c1, c2 and c3 derived from a curve fitting of cloud data;

defining a function M as a function of said L cirrus cloud ice particle size and constant values A0, A1, and A2 for associating said dependence of M on L;

defining a derivative dM of function M with respect to L and a derivative dL of L with respect to temperature T.sub.c ;

defining a derivative dTM of M with respect to cirrus cloud temperature T.sub.c ;

determining an average emissivity of said cirrus clouds as .epsilon..sub.04 ;

defining values associated with constants f.sub.00, f.sub..epsilon.0, and fT0 where f.sub.00 is characterized by the function f.sub.00 =(1-.epsilon..sub.04).sup.m(To)-1 and f.sub..epsilon.0 is characterized by the function -m(T.sub.0) (1-.epsilon..sub.04).sup.m(To)-1 and fT0 is characterized by the function FT0=(1-.epsilon..sub.04).sup.m To .multidot.1n(1.epsilon..sub.04).multidot.dTM(T.sub.0);

estimating the surface temperature T.sub.s and integrating said Planck function P(x) at said surface temperature Ts over each of said first and second radiance bands to obtain a first estimated radiance I.sub.G3 associated with said first band, and a second estimated radiance I.sub.G4 associated with said second band;

defining a matrix H having corresponding matrix elements indicative of said values fT0, f.sub..epsilon.0, and f.sub.00 ;

determining a noise covariance matrix C based on said first and second expected measurement errors .sigma..sub.3 and .sigma..sub.4 ; determining a covariance matrix of errors in the estimation of said cirrus cloud temperature C.sub..theta. where C.sub..theta. is characterized by the equation C.sub..theta. =(H.sup.T .multidot.C.sup.-1 .multidot.H).sup.-1 ;

determining a standard deviation .sigma.T of said estimated errors in C.theta.; and applying said .sigma.T to said L(Tc) particle size to obtain said estimate of error in said particle size.

15. The method according to claim 14, further comprising the step of creating a table of errors L given said above error estimate .sigma.T in said estimated cloud temperatures using said previously defined function L which depends on said cirrus cloud temperature.

16. The method according to claim 14, wherein said temperature measurements are measured in Kelvins.

17. The method according to claim 14, wherein said first band comprises a 3.7 .mu.m region and where said second band comprises a 10.9 .mu.m region radiance band.

18. The method according to claim 17, wherein said beginning and ending wave numbers in said first band are 2545 cm.sup.-1 and 2817 cm.sup.-1 respectively.

19. The method according to claim 18, wherein said beginning and ending wavenumbers for said second band comprise 885 cm.sup.-1 to 971 cm.sup.-1 respectively.

20. The method according to claim 19, wherein said first and second temperature measurement values are respectively 271 K and 293 K.