New hyperspectral calibration algorithm documentation

[czender]

Description

This is the draft improved algorithm to retrieve spectral reflectance based on existing and soon-to-be available inputs mentioned in #281. Suggestions and correction welcome (the more specific, the better). This draft is based on the algorithm proposed #88 from 2016.
For simplicity, the algorithm description currently omits the time dimension. It is implicit in all "real time" quantities. Zenith angle (theta) refers to the discrete set of angles at which spectralon imagery was taken.

Inputs and Outputs

Variable are denoted with the syntax type variable_name(dimensions) [units] where type, variable_name, and dimensions refer to their expected on-disk and in-memory storage structure in the hyperspectral workflow. The rank of some variables, notably the white and dark exposures, will be reduced between the Spectralon calibration and the HS workflow. The Spectralon exposures will be taken at a variety of zenith angles. Whether they are eventually stored one-zenith angle per file or aggregated into a new dimension (theta = zenith angle) is TBD, and does not change the mathematics of the procedure proposed below.

Required known or measured inputs:**

1. uint16 xps_img(wavelength,y,x) [cnt] =Exposure from experiment image (i.e., plants) (known, this is the "real time" imagery caputured by VNIR, SWIR cameras)
2. uint16 xps_wht(wavelength) [cnt] ="White" exposure from Spectralon panel (measured by VNIR, SWIR at multiple daytime zenith angles during dedicated calibration exercises in clear and overcast skies in Hyperspectral Calibration #208). The calibration reduction procedure currently averages over the spatial (x,y) dimensions of the Spectralon exposures prior to use in image calibration. This assumption can be relaxed if it is demonstrated that the Spectralon panel has systematically non-uniform spatial behavior.
3. uint16 xps_drk(wavelength) [cnt] ="Dark" exposure from Spectralon panel (measured once by VNIR, SWIR at night (preferably), or "with the lens cap on" in dedicated calibration exercise). Same comments on (x,y) as for xps_wht.
4. float32 rfl_wht(wavelength) [fraction] =Reflectance of "white" Spectralon panel (factory calibration) (assumed to be isotropic and time-constant, though should be periodically factory re-calibrated, and isotropicicity should be tested)
5. float32 flx_dwn(wavelength) [W m-2 m-1] =Downwelling spectral irradiance (measured by sensor on environmental logger)

Intermediate derived quantities:

1. float32 cst_cnv(wavelength) [cnt/(W m-2 m-1)] =Conversion constant expressing the proportionality between reflected spectral flux and Exposure (derived)
2. float32 flx_upw(wavelength) [W m-2 m-1] =Upwelling spectral flux (derived. and possibly measured for closure?)

Outputs

1. float32 rfl_img(wavelength,y,x) [fraction] =Reflectance of image (i.e., plants)

Proposed Algorithm to retrieve reflectance from measurements.

NB: The quantities below are conforming arrays of different ranks, so broadcasting (rank expansion) from lower to higher ranks is required to perform the arithmetic.

1. Assume incremental Spectralon exposure beyond dark counts is linear with incident (and reflected) spectral flux:
   • xps_wht-xps_drk=flx_dwn*rfl_wht*cst_cnv
2. Derive proportionality constant from calibrated Spectralon exposures:
   • cst_cnv=(xps_wht-xps_drk)/(flx_dwn*rfl_wht)
     Applying this procedure to Spectralon exposure data (xps_wht) at multiple zenith angles will determine whether cst_cnv is f(theta). Note that cst_cnv will implicitly account for any f(theta) dependence neglected in the available factory calibrated Spectralon reflectance (rfl_wht).
3. Assume proportionality constant for calibration sheet and plant image are identical:
   • xps_img-xps_drk=flx_dwn*rfl_img*cst_cnv
4. Derive plant reflectance from exposure
   • rfl_img=(xps_img-xps_drk)/(flx_dwn*cst_cnv)
     Here we use the cosine-zenith-angle-interpolated cst_cnv and time-interpolated flx_dwn
5. (Optional) Derive upwelling spectral flux from reflectance and compare it to measured upwelling spectral flux (if available) for closure/validation
   • flx_upw=flx_dwn*rfl_img
6. (Optional) Apply PAR-sensor SRF to downwelling irradiance, integrate, compare to measured PAR for closure (integration would require detailed information or assumptions about bandpass of each spectral channel).
7. More?

Alternative Calibration Procedure drafted and tested beginning 20170921:

Nomenclature:

cnt = Counts
cnv = Conversion
cst = Constant
dgr = Degrees
drk = Dark (not White)
drn = Duration
dwn = Downwelling
flx = Flux
frc = Fraction
img = Image (a "picture" of plants and soil, not a Target)
spc = Spectral (not broadband)
sza = Solar zenith angle (often used interchangably with Time)
trg = Target (i.e., a calibration target like the Spectralon 95% target, not an Image)
upw = Upwelling
wht = White (not Dark)
wvl = Wavelength

Coordinates:
float64 time [s] = Time
float64 x [m] = X distance
float64 y [m] = Y distance
float32 wvl [m] = Wavelength (300-1200 nm)
float32 drn [s] = Exposure duration (e.g., 20 ms, 30 ms, 40 ms)
float32 sza [dgr] = Solar zenith angle (sort of like a time)
float32 trg [frc] = Target (each targets has a nominal reflectance)

Required known or measured inputs:**

float32 flx_dwn_spc_trg(wvl) [W m-2 m-1] = Downwelling spectral irradiance towards target
float32 flx_dwn_spc_img(wvl) [W m-2 m-1] = Downwelling spectral irradiance towards image
float32 rfl_trg(trg,wvl) [frc] = Reflectance of target (measured at factory, time-invariant)
uint16 xps_trg_wht(sza,drn,trg,wvl) [cnt] = Exposure of calibration target during daylight
uint16 xps_trg_drk(trg,wvl) [cnt] = Exposure of calibration target at night (xps_trg_drk << xps_trg_wht)
uint16 xps_img(wvl,y,x) [cnt] = Exposure of realtime image containing plants, soil

Intermediate derived quantities:

float32 cst_cnv_trg(sza,trg,wvl) [cnt/(W m-2 m-1)] = Conversion constant expressing proportionality between incident spectral flux and exposure of target
float32 cst_cnv_img(sza,trg,wvl) [cnt/(W m-2 m-1)] = Conversion constant expressing proportionality between incident spectral flux and exposure of image
float32 flx_upw_spc_trg(sza,wvl) [W m-2 m-1] = Upwelling spectral flux (derived. and possibly measured for closure?)

Outputs:

float32 rfl_img(wvl,y,x) [frc] = Reflectance of image (i.e., plants)

Off-line Procedure:
Do this once off-line after every calibration campaign.
Store results in files or look-up table for use by on-line production procedure.

1. Clip each xps_trg_wht to rectangle and average over x,y to create scalar exposure.
   This removes inhomogeneities and best characterizes target.
   Throw away any calibration exposures that are saturated (artificially flat spectral response).
   Throw away any calibration exposures where the target is partially in direct light and partially shaded.

2. Find flx_dwn_spc_trg nearest in time/sza to xps_trg_wht. Should be within ~5 minutes.

3. Interpolate flx_dwn_spc_trg to VNIR/SWIR wavelength grid

4. Compute target conversion factor
   cst_cnv_trg=xps_trg_wht/flx_dwn_spc_trg
   cst_cnv_trg(drn,sza,trg,wvl) [cnt/(W m-2 m-1)]=(xps_trg_wht-xps_trg_drk)/flx_dwn_spc_trg
   cst_cnv_trg(drn,sza,trg,wvl) tells us, for a given duration and solar zenith angle, how many counts per unit incident flux to exptect from a surface with the known reflectance of that target. We subtract xps_trg_drk from xps_trg_wht to correct for dark counts. This is a small correction that ensures that black surfaces would have a retrieved reflectance of zero.

On-line Production Procedure:
Do this for each plant/soil image that needs calibration.

1. Find flx_dwn_spc_img nearest in time/sza to xps_img. Should be within ~5 minutes.

2. Interpolate flx_dwn_spc_img to VNIR/SWIR wavelength grid

3. Find cst_cnv nearest in sza to xps_img. Ideally within 5 degrees.

4. Compute image conversion factor
   cst_cnv_img=xps_img/flx_dwn_spc_img
   cst_cnv_img(drn,sza,trg,wvl) [cnt/(W m-2 m-1)] = (xps_img-xps_trg_drk)/flx_dwn_spc_img
   cst_cnv_img(drn,sza,trg,wvl) tells us, for the actual duration and solar zenith angle, how many counts per unit incident flux were measured from the plants with unknown reflectance. The small correction for dark counts must again be made.

5. Derive reflectance
   If the ratio of the image conversion factor to the target conversion factor is one, then the image must have the same reflectance as the target, so multiply this ratio by the target reflectance to get the unknown image reflectance:

   rfl_img=rfl_trg(trg,wvl)*cst_cnv_img/cst_cnv_trg(drn,sza,trg,wvl)
   rfl_img(trg,wvl,y,x)=rfl_trg(trg,wvl)*cst_cnv_img(drn,sza,wvl)/cst_cnv_trg(drn,sza,trg,wvl)

   This relies on reasonable yet sometimes subtle assumptions. A primary concern is that image and target conversion factors must be normalized to apply for the same duration. The spectral fluxes employed in each (i.e., flx_dwn_spc_trg and flx_dwn_spc_img) are already normalized per unit time. However, the measured count in each exposure (xps_trg_wht, xps_trg_drk, and xps_img) is sensitive to exposure duration. Each exposure may have been made with its own duration, and they must all be converted to the same equivalent duration. It is reasonable to assume that counts are linearly proportional to exposure duration for small changes in duration. Thus if the image exposure is 40 ms and the nearest valid target exposure was for 20 ms then the former should be divided by two or the latter multiplied by two when the conversion factors are applied.

6. Each production image takes place at a given solar zenith angle and for a given duration. The procedure retrieves a set of possible reflectances. Each member of the set uses the calibration curve of its respective target. Detector non-linearity (and drift or noise in factory calibrated reflectance) can cause retrieved reflectances to disagree. Some targets will have better characterized factory reflectances than others. Some targets will behave more uniformly with zenith angle than others. The optimal reflectance to report to the user is a weighted average of the retrieved reflectances. Ideally the weighting is characterized by long-term analysis of the results.

Storage is expensive and we can only store rfl_img(wvl,y,x) not rfl_img(trg,wvl,y,x). Hence, at least to begin, we will use only the most reliable target instead of the entire vector. This determination must still be made by experimenting with retrievals with each target, one at a time.

Next steps

More assumptions, input measurements, and/or more sophisticated algorithms would incorporate these additional sources of information:

1. BRDF (angular reflectance properties) of plant/leaves
2. 3D orientation of reflective surfaces (plant/leaves)
3. BRDF (angular reflectance) of calibration sheet
4. Direct/diffuse partitioning of downwelling flux

Notes

• Units of exposure are vague. Exposure is similar to a photon counter modulated by the spectral response function (SRF) of the sensor. Hence we call them "counts", denoted [cnt]. Counts range from [0..2^16-1] = [0..65535].
• add QAQC tests for sensitivity and saturation in each band
• cross validate with other sensors with known spectral response functions.

Completion Criteria

• Original algorithm modified to account for dark counts
• Routines to determine and locate Spectralon calibration from two zenith angles bounding a given HS image zenith angle
• Test in HS workflow, refine if necessary