The forum is locked.
The Ocean Color Forum has transitioned over to the Earthdata Forum (https://forum.earthdata.nasa.gov/). The information existing below will be retained for historical reference. Please sign into the Earthdata Forum for active user support.
All,
I need the power of the IOP Workshop Brain Trust.
I've been investigating the issue of satellite IOP retrieval sensitivity to inversion method (Levenburg-Marquart, SVD, etc.). Jeremy has also done some analysis using completely independent code applied to NOMAD data, which largely confirms the satellite processing results. I have posted a description of the analyses along with the results/conclusions here:
http://oceancolor.gsfc.nasa.gov/MEETINGS/OOXIX/IOP/inversion_methods.html
The bottom-line is this: the inversion method matters, and SVD matrix inversion appears to be much more sensitive to the input band suite.
For example, SVD results change dramatically (for the worse) in blue water when SeaWiFS Rrs(670) is included in a GSM inversion. Blue-water Rrs(670) retrievals have small absolute value (near zero) but large relative noise (digitization). The iterative methods (L-M, Amoeba) do not care if 670 is included in blue water. This is what I would expect, since we're minimizing sum ( retrieved_Rrs[lambda] - modeled_Rrs[lambda] )^2, and Rrs[670] is very small (so absolute noise is small and 670 noise does not contribute much to the merit function).
My hypothesis is that SVD is more sensitive because the linearization process (as in Hoge & Lyon 1996) inverts the spectral weighting. The overdetermined matrix inversion methods do not minimize Rrs differences; rather, they minimize (Ax - b), where b=(-aw[lambda]-bbw[lambda]*v), u=bb/a+bb, v=1-1/u. Here's a random case:
band Rrs[0-] u v b
412 0.030356 0.262304 -2.812377 0.004716
443 0.023846 0.213230 -3.689772 0.001896
490 0.015602 0.146454 -5.828084 -0.006041
510 0.009068 0.088935 -10.244158 -0.018931
555 0.004298 0.043689 -21.889184 -0.039295
670 0.000757 0.007923 -125.212256 -0.387961
So instead of minimizing Rrs differences where 670 is two orders of magnitude smaller than 412, we're minimizing Ax and b differences where 670 is two orders of magnitude greater than 412 (in an absolute sense). Presumably, measured-to-modeled differences and any noise contribution are equally magnified.
Can anybody confirm, deny, or provide alternative explanations to what we are seeing?
-- bryan
I need the power of the IOP Workshop Brain Trust.
I've been investigating the issue of satellite IOP retrieval sensitivity to inversion method (Levenburg-Marquart, SVD, etc.). Jeremy has also done some analysis using completely independent code applied to NOMAD data, which largely confirms the satellite processing results. I have posted a description of the analyses along with the results/conclusions here:
http://oceancolor.gsfc.nasa.gov/MEETINGS/OOXIX/IOP/inversion_methods.html
The bottom-line is this: the inversion method matters, and SVD matrix inversion appears to be much more sensitive to the input band suite.
For example, SVD results change dramatically (for the worse) in blue water when SeaWiFS Rrs(670) is included in a GSM inversion. Blue-water Rrs(670) retrievals have small absolute value (near zero) but large relative noise (digitization). The iterative methods (L-M, Amoeba) do not care if 670 is included in blue water. This is what I would expect, since we're minimizing sum ( retrieved_Rrs[lambda] - modeled_Rrs[lambda] )^2, and Rrs[670] is very small (so absolute noise is small and 670 noise does not contribute much to the merit function).
My hypothesis is that SVD is more sensitive because the linearization process (as in Hoge & Lyon 1996) inverts the spectral weighting. The overdetermined matrix inversion methods do not minimize Rrs differences; rather, they minimize (Ax - b), where b=(-aw[lambda]-bbw[lambda]*v), u=bb/a+bb, v=1-1/u. Here's a random case:
band Rrs[0-] u v b
412 0.030356 0.262304 -2.812377 0.004716
443 0.023846 0.213230 -3.689772 0.001896
490 0.015602 0.146454 -5.828084 -0.006041
510 0.009068 0.088935 -10.244158 -0.018931
555 0.004298 0.043689 -21.889184 -0.039295
670 0.000757 0.007923 -125.212256 -0.387961
So instead of minimizing Rrs differences where 670 is two orders of magnitude smaller than 412, we're minimizing Ax and b differences where 670 is two orders of magnitude greater than 412 (in an absolute sense). Presumably, measured-to-modeled differences and any noise contribution are equally magnified.
Can anybody confirm, deny, or provide alternative explanations to what we are seeing?
-- bryan
Bryan,
I see no flaw in your analysis.
What needs to be done in the matrix inversion approach is a weighting that is based on S/N ratios at each wavelength.
I will think on how best to incorporate it (I will not be surprise if it is addressed in Numerical Recipes).
Thanks!
Emmanuel
I see no flaw in your analysis.
What needs to be done in the matrix inversion approach is a weighting that is based on S/N ratios at each wavelength.
I will think on how best to incorporate it (I will not be surprise if it is addressed in Numerical Recipes).
Thanks!
Emmanuel
You can use LUD to do over determined direct inversions, you just must do a matrix multiplication trick to square off the the matrices. The over-determined direct inversions (not svd) that I use are least squares inversions that use LUD on the squared off matrices. I can send you the code. I don't have the math on the top of my head, but I can look it up and post that here later...
Paul
Paul
Paul,
I can try it (after the workshop), but do you have any reason to think that this would yield a different result than SVD? There are a multitude of alternatives available within the GSL library:
http://www.gnu.org/software/gsl/manual/html_node/Linear-Algebra.html
I don't have the knowledge to choose one method over another, but my suspicion is they would all yield nearly the same result (unless we change the weighting).
-- bryan
I can try it (after the workshop), but do you have any reason to think that this would yield a different result than SVD? There are a multitude of alternatives available within the GSL library:
http://www.gnu.org/software/gsl/manual/html_node/Linear-Algebra.html
I don't have the knowledge to choose one method over another, but my suspicion is they would all yield nearly the same result (unless we change the weighting).
-- bryan
Bryan:
You are correct that the over determined SVD will probably retrieve the same values as the linear over determined. However, it is critically important how the aph, adg and bbt spectral models are formulated. Most notably, I think I said this many times before, so forgive my repeating it, but including bands that don't co-vary without modifying the way the IOPs are formulated will tend to make the retrievals less accurate as they could have been. For example, the aph(443) value cannot be modeled with a fixed relationship between 412,490 and 555. If aph(443) is retrieved separately from aph(412,490,555), the additional benefit in retrieving the other IOPs from adding 443 can be realized. If a fixed model is used for aph(412,443,490,555,...) the error in matching the aph(443) will pile into something, in a least squares sense as the inversion attempts to best fit the data.
Weighting is one way to handle this, but some people are very interested in how aph(443) changes relative to the other channels, and with weighting, you'd get one aph(lambda_ref) value which you could use with your spectral model to get aph(X). If X is 443, it will be driven by the fixed spectral model and the least squares fit of the input data in the inversion and will not reflect the value that would be obtained in the function aph(443)=F(adg(443),bbt(443),aw(443),bbw(443),nLw(443). Know what I mean, or am I just babbling?
Paul
Man now I really wish I was able to attend the Italy meeting. It looks like it will be an excellent exchange of ideas.
You are correct that the over determined SVD will probably retrieve the same values as the linear over determined. However, it is critically important how the aph, adg and bbt spectral models are formulated. Most notably, I think I said this many times before, so forgive my repeating it, but including bands that don't co-vary without modifying the way the IOPs are formulated will tend to make the retrievals less accurate as they could have been. For example, the aph(443) value cannot be modeled with a fixed relationship between 412,490 and 555. If aph(443) is retrieved separately from aph(412,490,555), the additional benefit in retrieving the other IOPs from adding 443 can be realized. If a fixed model is used for aph(412,443,490,555,...) the error in matching the aph(443) will pile into something, in a least squares sense as the inversion attempts to best fit the data.
Weighting is one way to handle this, but some people are very interested in how aph(443) changes relative to the other channels, and with weighting, you'd get one aph(lambda_ref) value which you could use with your spectral model to get aph(X). If X is 443, it will be driven by the fixed spectral model and the least squares fit of the input data in the inversion and will not reflect the value that would be obtained in the function aph(443)=F(adg(443),bbt(443),aw(443),bbw(443),nLw(443). Know what I mean, or am I just babbling?
Paul
Man now I really wish I was able to attend the Italy meeting. It looks like it will be an excellent exchange of ideas.
Paul,
I do plan to implement an option for allowing one or two bands to float in aph. It makes sense to me. I don't think it will be difficult to implement, but I had to cut-off the development to run some of these analyses. If you have any results that you can share (e.g., fitting your gaussian with a floating aph(443)), that would be very useful. I also plan to implement a Boss and Roesler-type scheme as a wrapper around our generalized IOP model (giop). This would be another way to vary the aph spectrum on a per-pixel basis. No doubt there are other approaches. Perhaps it's time for somebody to start a thread on the topic of aph spectral models and selection methods.
I agree that down-weighting of non-covarying bands to minimize cross-spectral contamination due to a poorly chosen aph spectrum is a fool's errand.
- bryan
I do plan to implement an option for allowing one or two bands to float in aph. It makes sense to me. I don't think it will be difficult to implement, but I had to cut-off the development to run some of these analyses. If you have any results that you can share (e.g., fitting your gaussian with a floating aph(443)), that would be very useful. I also plan to implement a Boss and Roesler-type scheme as a wrapper around our generalized IOP model (giop). This would be another way to vary the aph spectrum on a per-pixel basis. No doubt there are other approaches. Perhaps it's time for somebody to start a thread on the topic of aph spectral models and selection methods.
I agree that down-weighting of non-covarying bands to minimize cross-spectral contamination due to a poorly chosen aph spectrum is a fool's errand.
- bryan
Hi Bryan,
Mutterings from a non-algorithm expert:
Your interpretation seems logical- fitting prescribed shapes to noisy data means adjusting the fit-parameters to the noise, whereas iteratively minimising a cost-function fits a smooth prescribed function through the noise. But going into detail with the implementation, the test may not be doing the matrix inversion any favours: Frank Hoge's 1996 paper showed aph and bbp retrieval errors to increase with increasing radiance error at increasing wavelength, and v.v. for adg. They also rejected 443nm and used [412, 490, 555] in the 3-band approach. It seems likely that the matrix inversion method will be more sensitive to the band-widths chosen because you're assuming distinct spectral differences in the IOPs. Which is not to say that you can't cajole either approach into giving good results..
Er, did that help?!
You're giving us plenty of food for thought!
jill.
Mutterings from a non-algorithm expert:
Your interpretation seems logical- fitting prescribed shapes to noisy data means adjusting the fit-parameters to the noise, whereas iteratively minimising a cost-function fits a smooth prescribed function through the noise. But going into detail with the implementation, the test may not be doing the matrix inversion any favours: Frank Hoge's 1996 paper showed aph and bbp retrieval errors to increase with increasing radiance error at increasing wavelength, and v.v. for adg. They also rejected 443nm and used [412, 490, 555] in the 3-band approach. It seems likely that the matrix inversion method will be more sensitive to the band-widths chosen because you're assuming distinct spectral differences in the IOPs. Which is not to say that you can't cajole either approach into giving good results..
Er, did that help?!
You're giving us plenty of food for thought!
jill.
All,
To be clear: I think the major problem is not with SVD or matrix inversion per se. I think it is with the way the problem has been formulated. The linearization is traditionally done by finding u = bb/(a+bb) and then defining v = 1 - 1/u = -a/bb such that -a = bb*v. I believe the problem is in the transformation to 1/u, which effectively means we're now fitting in 1/Rrs space. The small signal becomes the big signal and the small error becomes the big error.
I had an epiphany on the drive home. We just reformulate the problem as u*a = (1-u)*bb. Initial tests suggest this largely solves my 670 problem. 6-Band SVD looks much more like 6-Band LM (based on one scene histogram). More analysis to follow ...
-- bryan
To be clear: I think the major problem is not with SVD or matrix inversion per se. I think it is with the way the problem has been formulated. The linearization is traditionally done by finding u = bb/(a+bb) and then defining v = 1 - 1/u = -a/bb such that -a = bb*v. I believe the problem is in the transformation to 1/u, which effectively means we're now fitting in 1/Rrs space. The small signal becomes the big signal and the small error becomes the big error.
I had an epiphany on the drive home. We just reformulate the problem as u*a = (1-u)*bb. Initial tests suggest this largely solves my 670 problem. 6-Band SVD looks much more like 6-Band LM (based on one scene histogram). More analysis to follow ...
-- bryan
Bryan,
While the fix you find solves the 1/Rrs issue, yet the problem you found brings out an important point pertaining to all the methods used.
In the solution method used we should give different weights to wavelength based on signal to noise considerations (as one should do when fitting a line to data with variable S/N). The current methods we use weigh the Rrs (or (1/Rrs) at all wavelengths the same. There is no reason we could not include S/N considerations (e.g. based on specific sensor radiometric information and propagating error in the Rrs determinations, such as from atmospheric correction etc') in the inversion thus weighing more information we are more confident about.
Cheers,
Emmanuel
While the fix you find solves the 1/Rrs issue, yet the problem you found brings out an important point pertaining to all the methods used.
In the solution method used we should give different weights to wavelength based on signal to noise considerations (as one should do when fitting a line to data with variable S/N). The current methods we use weigh the Rrs (or (1/Rrs) at all wavelengths the same. There is no reason we could not include S/N considerations (e.g. based on specific sensor radiometric information and propagating error in the Rrs determinations, such as from atmospheric correction etc') in the inversion thus weighing more information we are more confident about.
Cheers,
Emmanuel
Emmanuel,
The GlobColour dataset is produced with an approach that accounts for S/N consideration:
In practice, the method we use in GSM/GlobColour weighs the Rrs according to an estimate of their uncertainty. The characterisation of Rrs uncertainties is based on the consideration of:
1. Statistical errors between satellite derived Rrs and in-situ measurements
2. Statistical errors of the semi-analytical modeling
We get an assessment of the quality of the error budget, through the following considerations:
1. Computation of chi-2 (should be close to 1 if error budget is correct).
2. Computation of residuals (residual is the difference between observed Rrs and modeled Rrs rebuilt with GSM direct model) (rms of residuals should be comparable with uncertainties estimates, assuming no bias).
3. Computation of covariance matrix of retrieved parameters (diagonal of the matrix should compare well with uncertainty estimates derived from statistical analysis with in situ ; off-diagonal terms give an indication of the correlation of retrieved parameters - the less the better)
Cheers
Odile
The GlobColour dataset is produced with an approach that accounts for S/N consideration:
In practice, the method we use in GSM/GlobColour weighs the Rrs according to an estimate of their uncertainty. The characterisation of Rrs uncertainties is based on the consideration of:
1. Statistical errors between satellite derived Rrs and in-situ measurements
2. Statistical errors of the semi-analytical modeling
We get an assessment of the quality of the error budget, through the following considerations:
1. Computation of chi-2 (should be close to 1 if error budget is correct).
2. Computation of residuals (residual is the difference between observed Rrs and modeled Rrs rebuilt with GSM direct model) (rms of residuals should be comparable with uncertainties estimates, assuming no bias).
3. Computation of covariance matrix of retrieved parameters (diagonal of the matrix should compare well with uncertainty estimates derived from statistical analysis with in situ ; off-diagonal terms give an indication of the correlation of retrieved parameters - the less the better)
Cheers
Odile
Dear Odile,
I think that Bryan pointed to a problem that is not answered by any of the current approaches.
When the number of wavelength=number of unknown all the algorithm will give the same exact solution (within convergence criteria error) for the same Rrs and same choice of eigenfunction. When we have more wavelengths than unknowns the problem is over determined and the solution with each method (linear or not) is found through some minimization. The issue is how to set this minimization to weigh the wavelengths based on the uncertainties we have (e.g. in Rrs due to S/N in space as well as error introduce through the atmospheric correction).
The cost function should be something of the form: chi(lambda) ~ |reconstructed Rrs (lambda)- space-based Rrs (lambda)|/(uncertainty in space-based Rrs(lambda)).
Cheers,
Emmanuel
I think that Bryan pointed to a problem that is not answered by any of the current approaches.
When the number of wavelength=number of unknown all the algorithm will give the same exact solution (within convergence criteria error) for the same Rrs and same choice of eigenfunction. When we have more wavelengths than unknowns the problem is over determined and the solution with each method (linear or not) is found through some minimization. The issue is how to set this minimization to weigh the wavelengths based on the uncertainties we have (e.g. in Rrs due to S/N in space as well as error introduce through the atmospheric correction).
The cost function should be something of the form: chi(lambda) ~ |reconstructed Rrs (lambda)- space-based Rrs (lambda)|/(uncertainty in space-based Rrs(lambda)).
Cheers,
Emmanuel
All,
I agree that we need to:
1) assign input uncertainties on Rrs
2) derive output uncertainties on model parameters
I understand that GlobColour has done this using in situ match-up statistics for (1) and L-M confidence intervals for (2). This is how I would do it. I believe the cost (merit) function is as Emmanuel describes, and if the input uncertainties are meaningful, and the model is sufficiently flexible, the chi-square will be close to 1. I hope that Odile or one of her colleagues can provide some small presentation on the efficacy of this approach.
Of course, assigning uncertainties on satellite Rrs can easily become an industry in itself. We don't really know the error on the atmospheric correction, and in any case it's likely to be spectrally systematic and geometry/geography dependent. An IOP model approach with sufficient degrees of freedom (e.g., Boss-Roesler) can probably get an equally good fit when faced with systematic error such as bad aerosol model selection. If we used the true error budget, our chi-square would likely be much less than 1, since it's a measure of goodness of fit, not true error.
Another source for satellite Rrs uncertainty is the Level-3 spatial/temporal variance (standard deviation on the binned mean). Does GlobColour make use of this information? If processing from Level-3 Rrs to Level-3 IOPs, standard deviation and n is available in the NASA products. Alternatively, it can be pre-tabulated by analysis of Level-3 climatologies, and temporally and geographically gridded look-up tables can be provided to the Level-2 or Level-3 IOP models.
I will look at implementing input/output uncertainies in the NASA codes.
-- bryan
I agree that we need to:
1) assign input uncertainties on Rrs
2) derive output uncertainties on model parameters
I understand that GlobColour has done this using in situ match-up statistics for (1) and L-M confidence intervals for (2). This is how I would do it. I believe the cost (merit) function is as Emmanuel describes, and if the input uncertainties are meaningful, and the model is sufficiently flexible, the chi-square will be close to 1. I hope that Odile or one of her colleagues can provide some small presentation on the efficacy of this approach.
Of course, assigning uncertainties on satellite Rrs can easily become an industry in itself. We don't really know the error on the atmospheric correction, and in any case it's likely to be spectrally systematic and geometry/geography dependent. An IOP model approach with sufficient degrees of freedom (e.g., Boss-Roesler) can probably get an equally good fit when faced with systematic error such as bad aerosol model selection. If we used the true error budget, our chi-square would likely be much less than 1, since it's a measure of goodness of fit, not true error.
Another source for satellite Rrs uncertainty is the Level-3 spatial/temporal variance (standard deviation on the binned mean). Does GlobColour make use of this information? If processing from Level-3 Rrs to Level-3 IOPs, standard deviation and n is available in the NASA products. Alternatively, it can be pre-tabulated by analysis of Level-3 climatologies, and temporally and geographically gridded look-up tables can be provided to the Level-2 or Level-3 IOP models.
I will look at implementing input/output uncertainies in the NASA codes.
-- bryan
Hi Bryan,
I completely agree with you about the satellite Rrs uncertainties and the fact that they are spectrally and geometry/geography (and time/season !) dependent. I've pointed that out in all my ocean color talks in the past few years (and that's the reason why I decided not to weigh Lwns in GSM after a while. But the GlobCOLOUR approach is actually pretty nice, in my opinion).
Anyway, as you know, we (at ICESS) are also producing confidence intervals for the GSM products. The method uses the projection of the sum of the Rrs residuals in solution space (Bates and Watts). Let me know if you are interested in implementing it.
Best,
Stéphane
I completely agree with you about the satellite Rrs uncertainties and the fact that they are spectrally and geometry/geography (and time/season !) dependent. I've pointed that out in all my ocean color talks in the past few years (and that's the reason why I decided not to weigh Lwns in GSM after a while. But the GlobCOLOUR approach is actually pretty nice, in my opinion).
Anyway, as you know, we (at ICESS) are also producing confidence intervals for the GSM products. The method uses the projection of the sum of the Rrs residuals in solution space (Bates and Watts). Let me know if you are interested in implementing it.
Best,
Stéphane
Out of curiosity, reformatted the problem in my IDL environment per Bryan's suggestion and re-ran the related NOMAD analyses.
In the parlance of this thread, the form to be fit was simplified from [ a + bb (1 - 1/u) = 0 ] to [ a u + bb (u - 1) = 0 ]. Assuming my implementation is correct, the results look encouraging.
Compare the original comparisons with these new comparisons.
No weighting applied, although I agree with its utility.
In the parlance of this thread, the form to be fit was simplified from [ a + bb (1 - 1/u) = 0 ] to [ a u + bb (u - 1) = 0 ]. Assuming my implementation is correct, the results look encouraging.
Compare the original comparisons with these new comparisons.
No weighting applied, although I agree with its utility.
It is interesting to note that the most dramatic improvement with increase in the number of wavelengths is in the decomposition of the absorption coefficient. Relatively little change is observed in the skill of the backscattering inversion.
Cheers,
Emmanuel
Cheers,
Emmanuel
Anybody care to speculate why the 4 and 5-band SVD underpredict the absorption highs. We see this in both the satellite and NOMAD results.
-- bryan
-- bryan
SeaWiFS scene processed with original and revised linearization:
http://oceancolor.gsfc.nasa.gov/MEETINGS/OOXIX/IOP/inversion/inver_lm_svd_scenes/
http://oceancolor.gsfc.nasa.gov/MEETINGS/OOXIX/IOP/inversion/inver_lm_svd_scenes_new/
I feel better now..
-- bryan
http://oceancolor.gsfc.nasa.gov/MEETINGS/OOXIX/IOP/inversion/inver_lm_svd_scenes/
http://oceancolor.gsfc.nasa.gov/MEETINGS/OOXIX/IOP/inversion/inver_lm_svd_scenes_new/
I feel better now..
-- bryan