- Email: [email protected]
A multispectrum nonlinear least squares fitting technique
J. Quanr. Specwosr. Rodialinr.Tramfer Vol. 53, No. 6. pp. 705-721. 1995 Copyright f; 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0022-4073(95)00015-1
0022-4073/95
A MULTISPECTRUM NONLINEAR LEAST FITTING TECHNIQUE D. CHRIS
59.50 + 0.00
SQUARES
BENNER,1_.$ CURTIS P. RINSLAND,$ V. MALATHY MARY ANN H. SMITH,8 and DAVID ATKINS7
DEVIJ
$Department of Physics, College of William and Mary, Williamsburg, VA 23187-8795, §Atmospheric Sciences Division, NASA Langley Research Center, Mail Stop 40lA, Hampton, VA 23681-0001, and r509 Constitution Avenue, NE, 4, Washington, DC 20002, U.S.A. (Received I June 1994; received for publication 23 Februury 1995)
Abstract-An extension of the nonlinear least squares spectrum fitting technique has been developed for the simultaneous fitting of multiple spectra. This procedure’s smaller number of fitted parameters as compared to fitting one spectrum at a time improves the determination of spectroscopic parameters. A more reliable evaluation of the errors associated with the solution is possible. Correlations among fitted parameters may preclude their determination from a single spectrum fit. If the correlations differ from spectrum to spectrum, however, separation of these parameters is often possible when including the spectra in one solution. Overfitting is also avoided when combining spectra with sufficiently different experimental physical conditions.
INTRODUCTION Several groups’-’ have employed a nonlinear least squares spectrum fitting technique to determine spectroscopic parameters from high resolution infrared laboratory molecular spectra. Based upon known experimental conditions, these programs fit the position, intensity, width and perhaps other parameters for each spectral line. An iterative procedure is used to minimize the sum of the squares of the differences between the observed spectrum and the spectrum calculated with the parameters. To account for the spectral distortions produced by the spectrometer, additional parameters such as the continuum (100% transmittance level) and instrumental line shape must be either specified or determined as part of the least squares solution. Each spectrum is fitted separately, and the results are later combined or averaged to derive the final values. We have employed this technique in many studies.2.8-‘9 Atmospheric spectra are commonly fitted with a similar technique.30-32 In such cases, the atmospheric technique assumes knowledge of spectral line parameters and derives atmospheric parameters such as mixing ratios of molecular species. Again, knowledge or derivation of the spectrometer characteristics is required. Often in a series of spectra the layers of the atmosphere are weighted differently in each spectrum. Physical conditions such as temperature, pressure and the mixing ratios of various gas molecules present in an atmospheric layer are derived by fitting a single spectrum. Physical conditions of the remaining layers are fixed to values determined previously from fits to other spectra in the series. A disadvantage of this technique is that small errors in the derived parameters in one layer of the atmosphere propagate to other layers. The least-squares algorithm compensates for these errors by adjusting the values of the parameters in the other layers during subsequent fits. This error propagation is not rigorously quantifiable. To avoid this problem, Carlotti33 developed a “global fit” least squares technique that determines atmospheric parameters by simultaneously fitting a series of spectra. The ratio of the number of observed points to derived parameters is not improved tTo
whom all correspondence
should
be addressed. 705
D. Chris Benner et al
706
(unless instrumental parameters are assumed equal or related on different spectra), but there is less error propagation. We have adopted the “global fitting” approach to analyze laboratory spectra. Series of spectra recorded over a range of experimental conditions are fitted simultaneously to derive the spectroscopic line parameters. Since the number of derived spectroscopic line parameters does not increase with the number of spectra fitted, the advantage of this technique is even greater for fitting laboratory spectra than for atmospheric spectra. However, the number of spectrometer-related parameters such as continuum level or channel spectrum parameters may increase proportionally to the number of spectra. INVERSION
TECHNIQUE
A modified Levenberg-Marquardt34”5 algorithm minimizes the sum of the squares of the residuals between the experimental spectra and the spectra calculated using the fixed and fitted parameters. Each fitted parameter is initialized with a starting value that is updated after each iteration. Calculated spectra are computed in two steps, namely, the generation of monochromatic spectra and then correction for spectrometer-induced distortions. The Appendix contains more details. The least squares technique requires the calculation of the derivative of each spectrum point with respect to each fitted parameter to form the Jacobian. One approach to the calculation of these derivatives is to recompute each spectrum with a small change in the fitted parameter [see Eq. (2) of the Appendix]. This technique is generally reliable, though it may fail sometimes if “small change” is not suitably defined or the calculation of the spectrum is not sufficiently precise. The computations by this method are, however, slow. The Appendix presents methods that are faster and/or more accurate. These methods often reduce computing time by orders of magnitude. A spectrum with more sample points will, of course, tend to have a larger weight in the solution simply because it has more points to be fitted. Oversampling of some spectra may especially pose a problem. Different spectra may also have different signal-to-noise ratios. These difficulties may be overcome by assigning an appropriate weight to each spectrum. ADVANTAGES
OF
THE
TECHNIQUE
There are several advantages to analyzing multiple spectra simultaneously. The technique is very versatile. For example, it is possible to fit spectra recorded with a Fourier transform spectrometer and a tunable diode laser spectrometer simultaneously, or to fit spectra obtained at different spectral resolutions. Thus, a high resolution spectrum may be appropriate at low pressure to resolve closely spaced lines, but a lower resolution spectrum may be adequate for high pressure since Lorentz broadening washes out fine detail. Analysis of absorption features blended into one profile because their line widths are greater than their line separation is often ambiguous. This is due to the strong correlation among the various retrieved parameters. In such instances, inclusion of a low pressure spectrum with unsaturated, partially or fully resolved absorption features will help determine the positions and intensities of these lines.
Fig. 1. Opposite.
Fig. 1. Nonlinear least squares spectrum fits to the P(6) manifold of the Y, band of “CH,. The observed spectra are normalized to the highest point in the spectral interval and all four panels have this same scaling for the ordinate. The color for each spectrum is the same from panel to panel. (a) The five observed spectra. All have a 5 cm path of methane and are obtained at room temperature. The orange line is approx. 0.4 torr of pure methane. The other spectra have an approx. I % mixing ratio of methane in dry air. The gray spectrum was obtained with a pressure of 75 torr, the light blue with a pressure of I50 torr, the dark blue with a pressure of 300 torr and the red with a pressure of 500 torr. (b) Amplified residuals of the spectra when fit individually with Voigt profiles. All five spectra are fit to the noise level, but the spectral line parameters are not forced to be consistent from fit to fit. (c) Amplified residuals of a simultaneous fit of all five spectra to Voigt profiles. The spectral line parameters are constrained to be consistent. (d) Amplified residuals of a simultaneous fit of all five using line profiles which include collisional narrowing and line coupling in the Rosenkrantz approximation.
0.01 --
(d)
Simultaneous - With Coupling
CA
z
5
0.00
2 -0.01
0.01
Single spectrum - Voigt
2957
2959
2958
Wavenumber (cm-‘) Fig. 1. Legend opposite. 707
2960
708
D. Chris
Benner
et al
0.01
0 3 4
0.00
2
”
II
-0.01
I
,I
I
iI
I
I
I
1.0
3747.0
3741.2
3147.4 Wavenumber
3147.6 (cm-t)
Fig. 2. Nonlinear least squares fit to a water vapor line at 3747.4 cm-‘. The line profile is a combination of residual water in the evacuated interferometer and the residual water vapor in a nitrogen purge path between the source and the entrance slit to the interferometer. The amount of water is somewhat different in each spectrum. The lower panel contains the observed spectra and the upper panel the residuals in the same units and with corresponding colors.
Torr
-3 -3 -3
-3
-3
-3
6a244i415
4.26(l) -12 -12
_a’
68
97 (51’
24
of solutions
64912 i6i’
Torr
.64725(15!' 6 5.21(O) _p?.O613(4j
--O.lS(lj
_?
(Dark Blue)
300
are
12
12
6.3
.6518(11)”
given in this of the last
table digit
(5)
are accurate displayed.
_;~.0587 (31
-14
.65:‘E3 -0.0075121 6.6Oiij O.O63i(ij
Voigt
to
;4j
the
levei
4.53(23 $.0623!?j \-0.084(143
68252 -0:00778(19j
.6507 3) -0.0067 141 6.77(5 0.065c; -0.0@81 i
-0.17:1 -C-l.17[! -C.2Dil -0.17il
-O.lEtl
nonVoigt
Multispectrum Solution
Air induced pressure shift coefficient of spectr-al line at 296K (~10.~cm l/arm! Parameter not obtainable by the single spectrum technique and only one path Spectral line intensity at 296K [x~O~~~cm~'/(molea~les/cm'j] 'Arithmetic mean and standard deviation of single spectrum solutions 'Air broadened Lorentz halfwidth coefficient at 296K of spectral line (cm '/atm) .
[uoted by the formal error estimates. The formal error estimates continuum line (x10-*/cm-‘) Slope of the fitted Not a parameter of the solution Zero pressure position of spectral line in cm-l plus 2955 cm m1 the effect of the pressure shift Includes By the technrque of ref. 14
(2j
4.1(2) 0.0535(111 12
_/6402W'
Results Combi 11 ed
0.01Ei3)6 6.O(llj'" 3.056(?J6 Ii
in Fig. 1.’
-0.17!1!
-3
_?
in units
6@094!6i' ."7550(14]' -8' -6 4.56(2! 4.72(5j -12 -12 -12 -12
5.47(3] 0.0615(4~ -12
-8’
150 ToI7 (Light Blue)
displayed
Due to modeling problems, none of the values
Ii
Ii
-6. 4.39 1)
Ii
Ii
parameters
Single Spectrum Solution
1. Some derived
6501?(4)' -8' 6.52(l)
_3
-0.16(l)
-3
65073(3)" _a. 6.40(1 ) -12 -12
Torr
(Gray)
75
-0.15(l)
(Grange)
3.4
This case is for illustration.
n solution”
lumber of jarameters fit
'osition' :hift' :ntensity¶ Ialfwidth" ,ine Mixing"
'osition‘ :hift7 ntensity' lalfwidth" ,ine Mixing"
j.4 Torr '5 Torr !50 Tort!3O Torr iO0 Tcrr
;lope2
Parameter
Table
710
D. Chris Benner et al
Figure l(a) displays five spectra of the P(6) manifold of the vj band of “CH, broadened by air and Table 1 lists some of the numerical results derived from fits to these spectra. The 0.01~cm-’ resolution spectra were obtained with the Fourier transform spectrometer at the National Solar Observatory’s McMath-Pierce solar telescope facility on Kitt Peak and a 5 cm path of methane. These spectra were used by Benner et alz9 to derive the widths and shifts of the methane spectral lines. Figure l(b) shows the expanded residuals (observed minus calculated) when each spectrum is fitted individually to between 20 and 26 coefficients (spectral line positions, intensities and/or halfwidths for up to 10 different spectral lines; continuum level and slope; spectrometer resolution; spectrometer phase error and/or zero level). The retrieved continuum slope is presented in Table 1 as an example of a spectrum parameter. This slope is caused primarily by the relative spectral sensitivity of the spectrometer itself and by the experimental setup, but is not very sensitive to the gas amount or other physical conditions such as temperature or pressure. The single spectrum analyses determine this parameter well because the fitted spectral interval is much wider than the manifold containing the strong methane lines and the wings of the lines are sufficiently weak that even at the higher pressures they have little effect upon the ends of the spectral interval. The change in the continuum slope in the multispectrum fit is within the estimated uncertainty of the single spectrum solution. Changing models to use a nonvoigt line profile (as discussed below) also had little effect upon the continuum slope. If a shorter spectral interval had been chosen, significant correlations between this parameter and parameters such as the halfwidths of the spectral lines would have made the slope more uncertain. In Fig. 1, the spectral feature at 2958.67 cm-’ consists of two lines only partially resolved even in the low pressure spectra. These two spectral lines are further detailed in Table 1. From the low pressure spectra the Lorentz halfwidth and pressure shift coefficients cannot be determined well because of their small effects upon the spectra. The higher pressure spectra, however, contain considerable information about these two parameters. The line parameters (position, intensity, Lorentz halfwidth and pressure shift) are not determinable, however, due to their strong correlations with each other. A nonunique fit is displayed in the table. The halfwidth of the higher wavenumber line was assumed to be 0.0625 cm-‘/atm29 and the halfwidth of the lower wavenumber line was similarly constrained to 0.0543 cm-‘/atm on the two lowest pressure spectra. The intensities derived for the four lowest pressure spectra in the single spectrum analyses are all within 5% of the mean. The highest pressure spectrum, however, gives intensities which are in error by a factor of about 1.5. The sum of the intensities of these two lines is close to that obtained in the other spectra, but the partitioning of the intensity is inaccurate. An attempt to match the wings with the halfwidth of one of the lines fixed to an inaccurate value moved the line positions within the blend which, in turn caused the repartitioning of the intensity. The halfwidth of the lower wavenumber line is also affected. The position of that line is not consistent with values derived from the fits to the other spectra in that the relationship between the position and the pressure is not very linear. As a result, the derived zero pressure positions and pressure shifts are poorly determined. The estimated uncertainties (see the following section) are also overly optimistic since the assumed model includes the uncertain value of at least one of the halfwidths of these two lines. By simultaneously fitting the five spectra, the halfwidth and pressure shift coefficients can be determined since the zero pressure line position and the line intensity are well constrained by the low pressure spectra. Furthermore, the five spectra are fit to a total of 47 coefficients rather than the total of 120 coefficients for the five individual spectrum fits (each spectral line is required to have only one value of each parameter which is valid for all five spectra). This reduces the overfitting that caused the residuals to look so small on the single spectrum fits [Fig. l(b)] even when the line parameters were poorly determined. The zero pressure position and intensity are now determined to an accuracy consistent with that of the low pressure spectrum fits. The pressure shifts and halfwidths are determined well for both lines. The drawback, however, is the much poorer fit to the spectrum [Fig. l(c)]. It is clear that there is something which is modeled incorrectly. An additional multispectrum fit was tried that included approximate corrections for collisional narrowing36 and line coupling. 37The residuals, displayed in Fig. I(d), from this fit show a marked improvement, but are still above the noise level. We believe that the remaining residuals are due to the inadequacy of the Rosenkrantz approximation of the line coupling and to small uncertainties
Nonlinear least squares
711
in physical parameters which are masked in a single spectrum fit by adjustment of the spectral line parameters. Notice, though, that the spectral line parameters have changed from those derived with Voigt profiles. In some cases (such as the halfwidths) the changes are far greater than the estimated uncertainties. This demonstrates the principle that the internal error estimate does not include uncertainties in the modeling. An additional advantage to fitting several spectra simultaneously is the ability to solve directly for parameters that are not determinable from a single spectrum fitting. Examples of these parameters are pressure shifts and temperature dependent parameters such as lower state energies of the transitions, temperature dependence of the Lorentz halfwidth and temperature dependence of the pressure shift coefficients. For example, the expression for the Lorentz halfwidth of a line is given in Eq. (45) of the Appendix. For a single spectrum with a single path there is a single value of the halfwidth that applies to the spectrum. It is, therefore, only possible to derive a single constant in that equation (generally at, the Lorentz halfwidth at a standard pressure and temperature). If, however, there are two or more spectra in the solution with different pressures and temperatures, there is a halfwidth appropriate for each spectrum and it is possible to derive as many halfwidth parameters as there are spectra. Since the pressures and temperatures are known there are only two parameters to be derived, namely at and n, the temperature dependence exponent. We have demonstrated the technique of deriving both parameters simultaneously (as well as position, pressure shift and temperature dependence of the pressure shift simultaneously) with the multispectrum technique in Ref. 38. Determination of the temperature dependence exponent of the halfwidth using a spectrum-by-spectrum analysis method requires a fitting of the results obtained from each of several spectra to Eq. (45) of the Appendix.14 This method can give inaccurate results unless the proper uncertainties in the retrieved parameters obtained from each spectrum are applied. Even so, different correlations among parameters on different spectra may constrain a multispectrum solution in ways not easily obtainable otherwise. ERROR
ANALYSIS
Error analysis is a difficult problem in least squares fitting of spectroscopic data.39-4’ A formal error estimate for each parameter can be obtained by the method of Whittaker and Robinson.42 This assumes that all second (and higher) derivatives of each fitted point with respect to each fitted parameter are zero and the data points are independent of each other. It also assumes no error in the modeling of the spectrum. These errors are reasonable estimates of the random error if none of the parameters is strongly correlated with any of the others and there are no significant modeling uncertainties. Off-diagonal elements of the covariance matrix can be used to calculate a correlation coefficient between each pair of parameters. If two parameters are strongly correlated (or anticorrelated), one or both may need to be constrained to a reasonable value known from other considerations such as theoretical constraints. This method for calculating errors is better than error calculations deduced from parameters obtained from fits to individual spectra (e.g., deriving broadening coefficients by least squares fitting measurements of Lorentz widths at a series of pressures). However, these derived uncertainties include only the uncertainties due to the noise level of the spectrum and the parameters derived in the solution. Not included are uncertainties due to errors in the measurement of the experimental conditions, the calculation of the partition function, incomplete modeling of the spectrometer and spectral line parameters not fitted in the solution. These unaccounted uncertainties often dominate a multispectrum solution since a large number of fitted points can decrease the formal error estimates to small values. CALIBRATION In a single spectrum fit the wavelength scale calibration is not required before fitting the data. The line positions may be corrected afterwards when an improved calibration is available. If absolute line positions are not needed, only a correct relative calibration is required to derive pressure shift information.
712
D. Chris Benner et al
The multispectrum technique, however, requires at least a relative calibration of the wavelength scale before the spectra are fitted. The relative wavelength scales must be precisely aligned or else the spectral lines will not be positioned properly on each spectrum. An absolute calibration adjustment can be applied after the fit to the spectral line positions if a better absolute calibration becomes available. Experimental parameters such as temperature, pressure, path length and mixing ratio must also be known to sufficient relative accuracy before performing multispectrum fits. Usually this means that careful laboratory measurement of these parameters is required. However, it is sometimes possible to use the spectra themselves to improve the relative calibration. For example, in a group of spectra obtained for measuring Lorentz foreign broadening coefficients, but not intensities of spectral lines, small mixing ratios are commonly used. If the mixing ratios are not known accurately enough, the derived intensity of the same spectral line will be different from spectrum to spectrum. Since the mixing ratios are small, the widths may not be significantly affected by self broadening. Thus, small adjustments to the nominal mixing ratios can make the mass paths of the spectra consistent. This is achieved when fits including all spectra no longer include intensity related residuals, Then the derived intensities of the spectral lines will be in error by a constant factor (due to uncertainty in the absolute mixing ratios). The derived widths, however, will be correct except for the usually insignificant error due to the self broadening correction. IMPLEMENTATION
The multispectrum nonlinear least squares fitting technique has been implemented into a FORTRAN program that runs on an IBM 80386 or 80486 compatible computer. It may be run in a batch mode, a completely interactive mode or a mode in which the initial input is completely specified as in the batch mode, but further modifications are allowed as in the interactive mode. In the interactive mode the user can modify the solution after a specified number of iterations or after meeting one or more convergence criteria, whichever occurs first. Changes entered through a series of menus allow the user to modify fixed parameters, change parameters for which a solution is found, add or delete parameters from the solution (such as line position or continuum level), add spectra to the solution and/or add or delete spectral lines. A backup option allows discarding the last set of iterations if the results were unsatisfactory. Continuation of the solution proceeds with a specified number of additional iterations if requested. Results of the error analysis, graphic output of the fits and output of derivatives are all supplied to help the user in evaluating the quality of the solution. Computation of the monochromatic spectrum is performed by standard line-by-line methods. Voigt profiles (including self and foreign pressure broadened Lorentz halfwidths), pressure shifts (self and foreign pressure induced), temperature dependence of Lorentz halfwidths and pressure shifts, temperature dependence of the line intensity (lower state energy of the transition), Dicke narrowing and line coupling (Rosenkrantz approximation) have been implemented. Calculation of the total internal partition sums is done by the method of Gamache et a1.43These are calculated separately for each gas and molecular isotope included in the solution. Flexible modeling of spectrometer-induced distortions is also included. Zero level, channel spectra and a polynomial fit to the continuum may be fitted to each experimental spectrum. For a Fourier transform spectrum, the effective apodization and phase error may be fitted. The maximum path difference of the interferometer and the solid angle of the entrance aperture are supplied by the user. For a tunable diode laser spectrum, the time constant of the electronics and the laser spectral width may be fitted. This multispectrum fitting program allows more than one absorbing path to be included in the solution. The physical conditions of these paths may be quite different and, therefore, the internal partition sums are calculated separately for each path, gas and isotope combination. An example of this situation is displayed in Fig. 2. These spectra were also obtained using the same Fourier transform spectrometer of the National Solar Observatory. The cell contained methane which does not have any appreciable absorption in this spectral region. A residual water vapor line arising from both within the evacuated interferometer tank and a dry nitrogen purged atmospheric path between the source and the interferometer entrance slit (in which the methane cell is placed) is present.
Nonlinear least squares
713
Table 2. Derived spectral line parameters from Fig. 2.’ Zero
Pressure
Line Position
Halfwidth
Coefficient
N,-Induced Pressure Shift Coefficient _.___
The water absorption caused by the residual gas in the interferometer tank shows as a narrow line. Not quite centered on the narrow line is the nitrogen-broadened water vapor line arising from the purged atmospheric path, As a result, these spectra provide N,-induced pressure shifts without the uncertainty of wavelength calibration. Several of these spectra are included in the solution (relative calibration of wavelength scale is necessary in this case) with different amounts of water in each spectrum. The zero pressure spectral line position, the N,-broadened Lorentz halfwidth coefficient and the N,-induced pressure shift coefficient are determined and given in Table 2. The amount of water vapor in each path is not well known and there is some unmodeled modification of the low pressure line shape due to the changing column of water in the path of the movable mirror of the FTS. These modeling problems cause the error analysis to produce values somewhat smaller than the true accuracy of the measurement. Indeed, the estimated uncertainties of the position and halfwidth are near the calibration accuracy of the wavenumber scale and the knowledge of the pressure in the purged atmospheric path. This illustrates the inadequacy of the error analysis to deal with systematic sources of errors when there are enough data included in the solution that the random sources of error are reduced too far. FUTURE
DIRECTIONS
Experimentally derived spectral line parameters are often used to determine more fundamental molecular parameters. For example, zero pressure positions of many lines in a single band may be combined to obtain vibrational energy and rotational constants of a molecule and the intensities of the same lines fit to a model with only a small number of free parameters. Ideally, residuals (observed minus calculated spectra) should be minimized not for the best derived spectral positions and intensities, but for the best vibrational energies, rotational constants, band intensity parameters and width parameters. This has been done to a greater or lesser extent for a single spectrum analysis by the whole band technique of Hoke and Shaw4 and by Pine and Looney.4s Carrying this idea even farther, if the spectra contain a sufficient number of vibration-rotation bands, the residuals should be minimized for the best spectroscopic constants and other theoretical parameters ultimately required. This involves fitting larger spectral intervals than those we have implemented so far and would require large amounts of computer processing time. On the other hand, theoretical constraints can also be imposed upon the spectrum on a much smaller scale. This might include constraints like equal intensity for doublets or the same spectrum parameters (such as channeling) for two spectra recorded close in time. Acknowledgements-Many people have contributed to this effort at one time or another. Donald J. Richardson of Hughes STX did some of the original programming of the single spectrum technique. Robert C. Reisweber of Computer Sciences Corporation assisted with information on previously existing software for the least squares technique, fast Fourier transforms and error analysis. C. Plymate and J. Wagner of the National Solar Observatory provided assistance in obtaining the data used for the examples and G. Ladd of the National Solar Observatory and C. H. Sutton of ST Systems Corporation helped in the initial processing of the data. Research at the College of William and Mary is supported by cooperative agreements with the National Aeronautics and Space Administration. The National Solar Observatory is operated by the Association of Universities for Research in Astronomy under contract with the National Science Foundation.
D. Chris Benner et al
714
REFERENCES I. Y. S. Chang and J. H. Shaw, Appl. Spectrosc. 31, 213 (1977). C. P. Rinsland, D. C. Benner, D. J. Richardson, and R. A. Toth, Appl. Opt. 22, 3805 (1983). L. R. Brown, J. S. Margolis, R. H. Norton, and B. D. Stedry, Appl. Spectrosc. 37, 287 (1983). J. W. C. Johns, J. Mol. Spectrosc. 125, 442 (1987). V. Dana and A. Valentin, Appt. Opt. 27, 4450 (1988). V. Dana and J.-Y. Mandin, JQSRT 48, 725 (1992). V. Dana, J.-Y. Mandin, C. Camy-Peyret, J.-M. Flaud, J.-P. Chevillard, R. L. Hawkins, and J.-L. Delfau,
2. 3. 4. 5. 6. 7.
Appl. Opt. 31, 1928 (1992). 8. V. Malathy Devi, C. P. Rinsland, and D. C. Benner, Appl. Opt. 23, 4067 (1984). 9. C. P. Rinsland and D. C. Benner, Appl. Opt. 23, 4523 (1984). IO. C. P. Rinsland, D. C. Benner, and V. Malathy Devi, Appl. Opt. 24, 1644 (1985). 11. D. C. Benner and C. P. Rinsland, J. Mol. Spectrosc. 112, 18 (1985). 12. V. Malathy Devi, C. P. Rinsland, M. A. H. Smith, and D. C. Benner, Appl. Opt. 24, 2788 (1985). 13. C. P. Rinsland, D. C. Benner, and V. Malathy Devi, Appl. Opt. 25, 1204 (1986). 14. C. P. Rinsland, V. Malathy Devi, M. A. H. Smith, and D. C. Benner, Appl. Opt. 27, 631 (1988). 15. V. Malathy Devi, C. P. Rinsland, M. A. H. Smith, and D. C. Benner, Appl. Opt. 27, 2296 (1988). 16. M. A. H. Smith, C. P. Rinsland, V. Malathy Devi, D. C. Benner, and K. B. Thakur, J. Opt. Sot. Am. B 5, 585 (1988). 17. D. C. Benner, V. Malathy Devi, C. P. Rinsland, and P. S. Ferry-Leeper, Appl. Opt. 27, 1588 (1988). 18. C. P. Rinsland, V. Malathy Devi, M. A. H. Smith, and D. C. Benner, Appl. Opt. 28, 211 (1989). 19. V. Malathy Devi, D. C. Benner, M. A. H. Smith, and C. P. Rinsland, Appl. Opt. 30,287 (1991). Erratum Appl. Opt. 30, 2928 (1991). 20. C. P. Rinsland, M. A. H. Smith, V. Malathy Devi, and D. C. Benner, J. Mol. Spectrosc. 150, 640
(1991). 21. C. P. Rinsland, M. A. H. Smith, V. Malathy Devi, and D. C. Benner, J. Mol. Spectrosc. 150, 173 (1991). 22. M. A. H. Smith, C. P. Rinsland, and V. Malathy Devi, J. Mol. Spectrosc. 147, 142 (1991). 23. C. P. Rinsland, A. Goldman, M. A. H. Smith, and V. Malathy Devi, Appl. Opt. 30, 1427 (1991). 24. V. Malathy Devi, D. C. Benner, C. P. Rinsland, and M. A. H. Smith, JQSRT 48, 581 (1992). 25. M. A. H. Smith, C. P. Rinsland, V. Malathy Devi, and D. C. Benner, Spectrochim. Acta 48A, 1257 (1992). 26. V. Malathy Devi, D. Chris Benner, M. A. H. Smith, and C. P. Rinsland, J. Mol. Spectrosc. 155, 333 (1992). C. P. Rinsland, M. A. H. Smith, V. Malathy Devi, and D. C. Benner, J. Mol. Spectrosc. 156, 507 (1992). V. Malathy Devi, D. C. Benner, M. A. H. Smith, and C. P. Rinsland, J. Mol. Spectrosc. 157, 95 (1993). D. C. Benner, V. Malathy Devi, M. A. H. Smith, and C. P. Rinsland, JQSRT 50, 65 (1993). A. Goldman and R. S. Saunders, JQSRT 21, 155 (1979). E. Niple, W. G. Mankin, A. Goldman, D. G. Murcray, and F. J. Murcray, Geophys. Res. Lett. 7, 489 (1980). 32. C. P. Rinsland, M. R. Gunson, J. C. Foster, R. A. Toth, C. B. Farmer, and R. Zander, J. Geophys. Res.
27. 28. 29. 30. 31.
96, 33. M. 34. K. 35. D.
1057 (1991). Carlotti, Appl. Opt. 27, 3250 (1988). Levenberg, Quart. Appl. Math. 2, 164 (1944). W. Marquardt, J. Sot. Indust. Appl. Math. 11, 431 (1963).
36. A. S. Pine and J. P. Looney, J. MoI. Spectrosc. 122, 41 (1987). 37. 38. 39. 40. 41. 42.
43. 44. 45. 46.
B. Gentry and L. L. Strow, J. Chem. Phys. 86, 5722 (1987).
V. Malathy Devi, D. C. Benner, M. A. H. Smith, and C. P. Rinsland, JQSRT
51, 439 (1994).
J. H. Shaw, N. Tu, and D. L. Agresta, Appl. Opt. 24, 2437 (1985). V. Dana, J.-Y. Mandin, and A. Hamdouni, Appl. Opt. 31, 1937 (1992). C. L. Brummel and L. A. Philips, J. Mol. Spectrosc. 159, 287 (1993). E. T. Whittaker and G. Robinson, The Calculus of Observations, Chap. 9, Blackie, Glasgow (1926). R. R. Gamache, R. L. Hawkins, and L. S. Rothman, J. Mol. Spectrosc. 142, 205 (1990). M. L. Hoke and J. H. Shaw, Appl. Opt. 21, 929 (1982). A. S. Pine and J. P. Looney, J. Chem. Phys. 93, 6942 (1990). F. Schreier, JQSRT 48, 743 (1992). APPENDIX Calculation
of Derivatives
Computationally, the most time consuming section of the nonlinear least squares technique is the calculation of derivatives. Often, these derivatives are calculated by the method of finite differences (defined below). In this Appendix some techniques that improve the accuracy of the derivatives and/or speed of the calculation significantly over this procedure are presented.
Nonlinear least squares
715
The spectrum S as a function of wavenumber, v, is calculated from the continuum B(v), channel spectrum C(v), gas spectrum convolved with the instrumental line shape b(v) and zero level offset z as: S(v) = B(v)C(v)b(v)
+ z.
(1)
The derivatives may be obtained numerically by using finite differences.
WV, 4) _ S(v,4 a4 -
+
6) - S(v,4) .!I
Here q is the parameter with respect to which the derivative is to be taken and c. is small enough that the derivative of S with respect to q is not significantly different at q and q + L, but large enough to avoid precision problems in the subtraction. In our approach, S(v) is first calculated with the current value of each parameter to be derived. Storing several intermediate steps of this computation allows many steps to be skipped in some derivative calculations. One step stored from the calculation of S(v) is the sum of the cross sections [defined in Eq. (24)] as a function of v for the spectral lines whose line parameters are all fixed in the analysis. This summation need not be repeated unless there is a change in these line parameters. Spectrometer Efects The derivative of the calculated signal with respect to each spectrometer effect is always possible without the complete recalculation of b(v). The recalculation of J(v) as defined in Eqs. (21) and (24) is never necessary since J(v) is always the same as it is in the nominal case. Also, in every case the derivative is only nonzero for the one spectrum to which the parameter applies. A discussion of each derivative in our implementation is given here. Zero level
The derivative of the calculated signal S(v) with respect to the zero level offset is trivial.
Continuum level The derivatives with respect to the continuum are expressed as:
as(v) = a4
-
C(v)b(v)
WV) at7 .
-
C(v) and b(v) are independent of the continuum parameters. Most of the computation is involved in the calculation of b(v). When b(v) is computed for the nominal case, it may be stored and retrieved at this point. C(v) could be treated in the same way, but recomputation of C(v) is not too time consuming. The continuum may be modeled as a polynomial of order N in terms of the starting wavenumber of the fitted interval, vStartrand a series of coefficients. B(v) = 2 Uv II,=0 The parameters to be derived [corresponding
(5)
to q in Eq. (4)] are B,. = (v -
gp
- v,,,,)“‘.
vstarty.
m
If there is more than one coefficient in the solution, it can be shown that the following rapidly calculable relationship holds.
D. Chris Eknner et al
716
Channel spectra
For the derivatives with respect to the channel spectrum parameters, a relationship similar to Eq. (4) may be obtained. y
= B(V)b(“) y
.
(8)
Again, b(v) is by far the most time-consuming part of the equation for a computer, but has already been computed for the nominal case using Eq. (1). B(v) has also been calculated for the nominal case, but can be quickly recomputed. The channel spectra are computed using an equation equivalent to that of Niple.3’ C(V) = P(V)=/{1
~(v){Re[b(v)lcos[e(v)l- Im[Wv)lWVv)l),
-Tca,[
1 -COS$(V
-&)]~+{~~~aiSin$(V
(9) -&))‘y
(lo)
(11)
khan
z, aisin
Nchiln 1 ai i= I is the number of channel spectra and ai, Pi and di are the amplitude, period In these equations NChan and phase of channel i. The derivatives with respect to the amplitude of channel i can be readily calculated. 1-
WV) -=~[j?(cosd,aai
l)-sind,]+C(v) 3
{Im[b(v)]sin d, + Re[b(v)][cos di - l]},
(12)
(13) p = tan
8,
(14)
p(v){Re[b(v)]sin 8 + Im[b(v)]cos 0} II=
9
P(1 + P’) p =
(15)
1 - C a,(1 - cos di).
(16)
i=l
The derivatives with respect to the period of channel i can also be determined.
at(v) diai
-=F api
, 1
rl(cosdi+/Isindi)+
sin di - 5 cos di)
5=%Z4” C %(P
ai sin di.
1,
(17)
(18)
i=l
The derivative with respect to the phase of channel i is particularly simple if the derivative with respect to the period of the same channel was calculated.
acw pi aqv) -= V+aP,' ah
(19)
Instrumental line shape
A derivative with respect to a parameter such as Fourier transform spectrometer phase error involved in the convolution with the instrument profile or for a spectral line parameter is again analogous to Eq. (4).
as(v) = a4
-
B(v)C(v)
WV)
-.
a4
(20)
117
Nonlinear least squares
B(v) and C(v) are simple enough that they can be easily recomputed in case they are not saved from the evaluation of the nominal case. The function b(v) can be expressed as a convolution of the line profile, w(v), and the monochromatic spectrum, J(v). Jc w(v)J(v) dv. s --3cI
b(v)=
(21)
For parameters contributing to w(v) there are two approaches possible for determining the derivative. First, note that the partial derivative of b(v) may be expressed as a modified convolution with the partial derivative of w(v) replacing w(v).
s
WV) a9
= +
-=
J(v)-
wv) a4
dv
(22)
.
If the partial derivative of w(v) is known, this convolution is accurate and, because the values of J(v) can be stored from the nominal case, faster than Eq. (2). If the partial derivative is not known, Eq. (2) may be used. Considerable computation time may be saved even here since J(v) need not be recomputed. Monochromatic Spectrum For parameters derived for b(v).
contributing
to the monochromatic
spectrum,
a similar integral can be
dJ(v>dv, w(v)---co 84 s
db(v) -= a4
=
(23)
This leaves only the determination of aJ(v)/aq. In most cases this can be handled without resorting to the method of finite differences [Eq. (4)]. Except in complex cases (for example, line mixing cases not simplified by too many approximations), the monochromatic spectrum can be expressed in terms of a line shape for each individual spectral line using Eq. (24). (24) A,, = mass path of gas i in path j, ZGk= intenstiy of line k of gas i in path j, rcVk = line shape factor for line k of gas i in path j. Here the Cpntity in the solution.
A,,z#
Kuk(v)
is the cross section. For laboratory spectra only ZVkand icOkare varied
Intensity The intensity Z,,,,can be expressed as a function of the line intensity Zoa at temperature T,,, the lower state energy Ei, the temperature q, the second radiation constant C, and several other factors (partition function, Honl-London factors, etc.), R,, which do not depend upon either Zoik or E;;.
The derivative with respect to the intensity of a line at temperature
aJ(v> -J(v)~~~A,Q~(v)$, -= azoik
aJ(v)
-=
azOik
,rA,hjkK,k(“).
J(v) N ‘11, Olk
(26) Olk
j=l
-r
TO can now be calculated.
(27)
, - I
Since J(v) and Key are calculated for the nominal case, the right side of this equation may be readily evaluated. The lower state energy Ei in the evaluation of S(v) only contributes to Zijk. The partial derivatives of J(v) with respect to E;; may then be calculated similarly.
718
D. Chris Benner et al
aJ(v) -= aE;;.
- J(v)C, ‘f” AijIUk
(28)
j=l
In the special case that there is only one path or that all paths have the same temperature T, these two derivatives are related closely. (29) Placing this relation into Eqs. (20) and (23), the following computationally found.
rapid relationship is
(30) Line shape
Some spectral line parameters (such as position and width parameters) contribute to only the line shape r$(v). aJw _=
-J(v)JkNfh/l,y. j=
a%k
I
(31)
I
A combination of the real and imaginary parts of the complex error function is used for most line shapes. For example, the Voigt function is simply the real part of the complex error function. The complex error function depends on only two parameters, x and y. The functions K(x, y) and L(x, y) are then defined in terms of the complex error function E(x,y). E(x, y) = KG, y) + iL(x, y).
(32)
Schreier46 reviews many of the techniques that have been proposed for calculation of E(x, y) and shows that the partial derivatives of K(x, y) and L(x, y) with respect to x and y can be expressed in terms of K(x, y), L(x, y), x and y only. The following relations can be derived from his equations.
W-T Y) ax
aw, Y 1 =2 ay
(33)
=~[YW,Y)-~W,Y)I,
yK(x, y) + xL(x, y) - -!-
1
Jn '
(34)
(35) aw,~)
ay
aK(x,~)
=ax’
(36)
Line position and pressure shift
The line position Vline is expressed in one of the models of Smith et alZSin terms of the zero pressure position, vo, pressure shift at temperature To, 6, temperature dependence of the pressure shift, 6’, gas sample pressure P and temperature T. Vline= ~0 + [6 + 6 ‘(T
- To)]P.
(37)
For the Voigt profile, r$(V) is just K(x, y) for the appropriate gas, path and spectral line. The parameter y is independent of vline9 vo, 6, and 6’ while x depends upon v, Vline and the Doppler width of the spectral line, an. x =
Jln 2(VaD
Vline)
(38)
719
Nonlinear least squares
The derivatives of +.(v) with respect to v,,, 6 and 6’ are then expressible in terms of K(x,, , yijk) and L(x,, , yijk). Consider, for example, the form for vO.
aKuk (V> aK(Xvkhk 1
---=
aVO,L
)
axok av+
axuk
m%k
3 Yijk) aY,ik
04
ayr/k avO,k
(39) ’
The second term on the right-hand side of Eq. (39) is zero. This equation can be expressed in a more convenient form by using Eqs. (33) (37) and (38). -2JiZ aKvk(v) -----= (4) b,k L (xijk 7 Yyk > $k K(-%,k , yi/k )I. a%, %,k Similar derivatives may be found for 6 and 6 ‘. General relations among these three derivatives hold. aKvk(v)
=
aslk
p
aKvk(V> ’
dvo,i
(41) ’
Under certain conditions (such as a spectrum obtained with only one path containing the absorbing gas), even more general relationships hold. (43) (44) If pressure shifts induced by different gases are separated in the solution, there are similar relationships that can be formed with these three derivatives that will involve the mixing ratios of the corresponding gases. Lorentz
halfwidth
A similar procedure may be followed with the Lorentz halfwidth, CI~,and the variables contained in it, namely, the Lorentz halfwidth at temperature To, the pressure P and the temperature dependence of the halfwidth n.
For the Voigt profile an equation analogous to Eq. (39) may be obtained. Since x does not depend upon CX~or n, considerable simplification is possible. The relationship between y and these two parameters is given by: (46)
The derivative of the line shape factor again derived.
Key
with respect to tx[ for the Voigt profile case can be
Jn
The derivative with respect to the temperature derivative.
1 .
I’,,kK(x,k,,y,,)+x,L(x,,y,)-l
(47)
dependence of the halfwidth is related to this
(48)
If only one path contains the absorbing gas or all paths have the same temperature fundamental relation can be derived.
T, a more
(49)
D. Chris Benner et al
720
Relations similar to Eqs. (47)-(49) may be derived for cases in which broadening due to different gases are handled separately. These equations involve the mixing ratios of the gases. Doppler
width
The derivative with respect to Doppler width a: at temperature TO, is somewhat more complicated since both x and y depend upon a:, but this can be simplified by relating it to previously derived derivatives. (50) If all paths are at the same temperature T and pressure P (or the pressure shift is zero at T) or if there is only one path, a more general relation holds. as(v) I = -c aa “D,, 1[
In our implementation,
(V-
Vline,k > y
+
a& 91.
O,k
(51) f-d
we solve for the molecular mass, M, rather than for the Doppler width. a0
D
=v”
(52)
c
The derivatives with respect to a: and A4 can be related.
a~,,,(4 ah(v) aa& aM,=x abf,'
aa; _=-_ aM
(53)
4
(54)
2~’
An analogous expression is derivable for the partial derivative of S(v). Nonvoigt
spectral lines
Nonvoigt line shapes can be treated in the same way if Key can still be expressed in terms of the complex error function. For example, Pine and Looney36 define a line profile that includes pressure narrowing. K(V) = Re
E(x, Y + e)
1
(55) [ 1 - JneE(x, y + e) ’ E was defined in terms of K and L in Eq. (32). The parameter e defines the pressure narrowing. This function can be differentiated with respect to x and y to obtain the same derivatives previously defined for the Voigt function. In addition, the derivative with respect to e may also be determined. All these derivatives are expressible in terms of K, L, x, y, e and v only. Similarly, the line shape for the Rosenkrantz approximation to line coupling is given in terms of E(x,y) by Gentry and Strow.37 rc(v) = Re[E(x, v)] + PYIm[E(x, r)].
(56)
Again this can be differentiated with respect to x and y to obtain the derivatives derived above for the Voigt case or be differentiated with respect to Yik,,the line coupling coefficient. To combine collisional narrowing and line coupling, a first approximation is to use the complex form of Key from Eq. (55) for E(x, y) in Eq. (56). All derivatives can still be expressed in terms of K, L, x, y, e, Y and v. NOMENCLATURE
A-Mass path of gas a-Amplitude of the channel spectrum B-Continuum &-Continuum polynomial coefficient b--Gas spectrum convolved with instrument line shape C-Channel spectrum
Nonlinear least squares
C,-Second radiation constant c--Speed of light d-Auxillary channel spectrum variable E-Complex error function E “-Lower state energy e-Collisional narrowing coefficient of spectral line Z-Spectral line intensity Z&Spectral line intensity at standard temperature i-Square root of - 1 [Eq. (32) only] i-Index variable J-Monochromatic gas spectrum j-Index variable (over physical paths) K-Real part of the complex error function k-Index variable (over spectral lines) k-Boltzman’s constant [Eq. (52) only] L--Imaginary part of the complex error function M-Mass of the absorbing molecule m-Continuum polynomial order N-Polynomial order of continuum Nchan-Number of channel spectra Npath-Number of physical paths in the solution n-Temperature dependence of the halfwidth exponent P-Period of channel spectrum [Eqs. (10-19) only] P-Gas pressure q-Parameter to be determined by the least squares fit R-Collection of miscellaneous factors S-Spectral amplitude T-Gas temperature TO-Standard temperature w-Instrumental line profile x-Voigt x parameter Y-Spectral line coupling coefficient (Rosenkrantz approximation) y-Voigt y parameter z-Zero level offset cm-Doppler halfwidth of spectral line cc,-Lorentz halfwidth of spectral line c$-Lorentz halfwidth coefficient of spectral line B-Auxillary channel spectrum variable b-Pressure shift coefficient 6’-Temperature dependence of the pressure shift coefficient c-Small finite number q-Auxillary channel spectrum variable 8-Auxillary channel spectrum variable k--Spectral line shape function p-Auxillary channel spectrum variable v-Wavenumber v,i,,-Wavenumber position of spectral line V ,,,,,-Wavenumber of start of fitted interval v,-Zero pressure wavenumber position of spectral line <-Channel spectrum auxillary variable X-Ratio of circumference to diameter of a circle p-Channel spectrum auxillary variable b--Wavenumber of a channel spectrum peak
721
0022-4073/95
A MULTISPECTRUM NONLINEAR LEAST FITTING TECHNIQUE D. CHRIS
59.50 + 0.00
SQUARES
BENNER,1_.$ CURTIS P. RINSLAND,$ V. MALATHY MARY ANN H. SMITH,8 and DAVID ATKINS7
DEVIJ
$Department of Physics, College of William and Mary, Williamsburg, VA 23187-8795, §Atmospheric Sciences Division, NASA Langley Research Center, Mail Stop 40lA, Hampton, VA 23681-0001, and r509 Constitution Avenue, NE, 4, Washington, DC 20002, U.S.A. (Received I June 1994; received for publication 23 Februury 1995)
Abstract-An extension of the nonlinear least squares spectrum fitting technique has been developed for the simultaneous fitting of multiple spectra. This procedure’s smaller number of fitted parameters as compared to fitting one spectrum at a time improves the determination of spectroscopic parameters. A more reliable evaluation of the errors associated with the solution is possible. Correlations among fitted parameters may preclude their determination from a single spectrum fit. If the correlations differ from spectrum to spectrum, however, separation of these parameters is often possible when including the spectra in one solution. Overfitting is also avoided when combining spectra with sufficiently different experimental physical conditions.
INTRODUCTION Several groups’-’ have employed a nonlinear least squares spectrum fitting technique to determine spectroscopic parameters from high resolution infrared laboratory molecular spectra. Based upon known experimental conditions, these programs fit the position, intensity, width and perhaps other parameters for each spectral line. An iterative procedure is used to minimize the sum of the squares of the differences between the observed spectrum and the spectrum calculated with the parameters. To account for the spectral distortions produced by the spectrometer, additional parameters such as the continuum (100% transmittance level) and instrumental line shape must be either specified or determined as part of the least squares solution. Each spectrum is fitted separately, and the results are later combined or averaged to derive the final values. We have employed this technique in many studies.2.8-‘9 Atmospheric spectra are commonly fitted with a similar technique.30-32 In such cases, the atmospheric technique assumes knowledge of spectral line parameters and derives atmospheric parameters such as mixing ratios of molecular species. Again, knowledge or derivation of the spectrometer characteristics is required. Often in a series of spectra the layers of the atmosphere are weighted differently in each spectrum. Physical conditions such as temperature, pressure and the mixing ratios of various gas molecules present in an atmospheric layer are derived by fitting a single spectrum. Physical conditions of the remaining layers are fixed to values determined previously from fits to other spectra in the series. A disadvantage of this technique is that small errors in the derived parameters in one layer of the atmosphere propagate to other layers. The least-squares algorithm compensates for these errors by adjusting the values of the parameters in the other layers during subsequent fits. This error propagation is not rigorously quantifiable. To avoid this problem, Carlotti33 developed a “global fit” least squares technique that determines atmospheric parameters by simultaneously fitting a series of spectra. The ratio of the number of observed points to derived parameters is not improved tTo
whom all correspondence
should
be addressed. 705
D. Chris Benner et al
706
(unless instrumental parameters are assumed equal or related on different spectra), but there is less error propagation. We have adopted the “global fitting” approach to analyze laboratory spectra. Series of spectra recorded over a range of experimental conditions are fitted simultaneously to derive the spectroscopic line parameters. Since the number of derived spectroscopic line parameters does not increase with the number of spectra fitted, the advantage of this technique is even greater for fitting laboratory spectra than for atmospheric spectra. However, the number of spectrometer-related parameters such as continuum level or channel spectrum parameters may increase proportionally to the number of spectra. INVERSION
TECHNIQUE
A modified Levenberg-Marquardt34”5 algorithm minimizes the sum of the squares of the residuals between the experimental spectra and the spectra calculated using the fixed and fitted parameters. Each fitted parameter is initialized with a starting value that is updated after each iteration. Calculated spectra are computed in two steps, namely, the generation of monochromatic spectra and then correction for spectrometer-induced distortions. The Appendix contains more details. The least squares technique requires the calculation of the derivative of each spectrum point with respect to each fitted parameter to form the Jacobian. One approach to the calculation of these derivatives is to recompute each spectrum with a small change in the fitted parameter [see Eq. (2) of the Appendix]. This technique is generally reliable, though it may fail sometimes if “small change” is not suitably defined or the calculation of the spectrum is not sufficiently precise. The computations by this method are, however, slow. The Appendix presents methods that are faster and/or more accurate. These methods often reduce computing time by orders of magnitude. A spectrum with more sample points will, of course, tend to have a larger weight in the solution simply because it has more points to be fitted. Oversampling of some spectra may especially pose a problem. Different spectra may also have different signal-to-noise ratios. These difficulties may be overcome by assigning an appropriate weight to each spectrum. ADVANTAGES
OF
THE
TECHNIQUE
There are several advantages to analyzing multiple spectra simultaneously. The technique is very versatile. For example, it is possible to fit spectra recorded with a Fourier transform spectrometer and a tunable diode laser spectrometer simultaneously, or to fit spectra obtained at different spectral resolutions. Thus, a high resolution spectrum may be appropriate at low pressure to resolve closely spaced lines, but a lower resolution spectrum may be adequate for high pressure since Lorentz broadening washes out fine detail. Analysis of absorption features blended into one profile because their line widths are greater than their line separation is often ambiguous. This is due to the strong correlation among the various retrieved parameters. In such instances, inclusion of a low pressure spectrum with unsaturated, partially or fully resolved absorption features will help determine the positions and intensities of these lines.
Fig. 1. Opposite.
Fig. 1. Nonlinear least squares spectrum fits to the P(6) manifold of the Y, band of “CH,. The observed spectra are normalized to the highest point in the spectral interval and all four panels have this same scaling for the ordinate. The color for each spectrum is the same from panel to panel. (a) The five observed spectra. All have a 5 cm path of methane and are obtained at room temperature. The orange line is approx. 0.4 torr of pure methane. The other spectra have an approx. I % mixing ratio of methane in dry air. The gray spectrum was obtained with a pressure of 75 torr, the light blue with a pressure of I50 torr, the dark blue with a pressure of 300 torr and the red with a pressure of 500 torr. (b) Amplified residuals of the spectra when fit individually with Voigt profiles. All five spectra are fit to the noise level, but the spectral line parameters are not forced to be consistent from fit to fit. (c) Amplified residuals of a simultaneous fit of all five spectra to Voigt profiles. The spectral line parameters are constrained to be consistent. (d) Amplified residuals of a simultaneous fit of all five using line profiles which include collisional narrowing and line coupling in the Rosenkrantz approximation.
0.01 --
(d)
Simultaneous - With Coupling
CA
z
5
0.00
2 -0.01
0.01
Single spectrum - Voigt
2957
2959
2958
Wavenumber (cm-‘) Fig. 1. Legend opposite. 707
2960
708
D. Chris
Benner
et al
0.01
0 3 4
0.00
2
”
II
-0.01
I
,I
I
iI
I
I
I
1.0
3747.0
3741.2
3147.4 Wavenumber
3147.6 (cm-t)
Fig. 2. Nonlinear least squares fit to a water vapor line at 3747.4 cm-‘. The line profile is a combination of residual water in the evacuated interferometer and the residual water vapor in a nitrogen purge path between the source and the entrance slit to the interferometer. The amount of water is somewhat different in each spectrum. The lower panel contains the observed spectra and the upper panel the residuals in the same units and with corresponding colors.
Torr
-3 -3 -3
-3
-3
-3
6a244i415
4.26(l) -12 -12
_a’
68
97 (51’
24
of solutions
64912 i6i’
Torr
.64725(15!' 6 5.21(O) _p?.O613(4j
--O.lS(lj
_?
(Dark Blue)
300
are
12
12
6.3
.6518(11)”
given in this of the last
table digit
(5)
are accurate displayed.
_;~.0587 (31
-14
.65:‘E3 -0.0075121 6.6Oiij O.O63i(ij
Voigt
to
;4j
the
levei
4.53(23 $.0623!?j \-0.084(143
68252 -0:00778(19j
.6507 3) -0.0067 141 6.77(5 0.065c; -0.0@81 i
-0.17:1 -C-l.17[! -C.2Dil -0.17il
-O.lEtl
nonVoigt
Multispectrum Solution
Air induced pressure shift coefficient of spectr-al line at 296K (~10.~cm l/arm! Parameter not obtainable by the single spectrum technique and only one path Spectral line intensity at 296K [x~O~~~cm~'/(molea~les/cm'j] 'Arithmetic mean and standard deviation of single spectrum solutions 'Air broadened Lorentz halfwidth coefficient at 296K of spectral line (cm '/atm) .
[uoted by the formal error estimates. The formal error estimates continuum line (x10-*/cm-‘) Slope of the fitted Not a parameter of the solution Zero pressure position of spectral line in cm-l plus 2955 cm m1 the effect of the pressure shift Includes By the technrque of ref. 14
(2j
4.1(2) 0.0535(111 12
_/6402W'
Results Combi 11 ed
0.01Ei3)6 6.O(llj'" 3.056(?J6 Ii
in Fig. 1.’
-0.17!1!
-3
_?
in units
6@094!6i' ."7550(14]' -8' -6 4.56(2! 4.72(5j -12 -12 -12 -12
5.47(3] 0.0615(4~ -12
-8’
150 ToI7 (Light Blue)
displayed
Due to modeling problems, none of the values
Ii
Ii
-6. 4.39 1)
Ii
Ii
parameters
Single Spectrum Solution
1. Some derived
6501?(4)' -8' 6.52(l)
_3
-0.16(l)
-3
65073(3)" _a. 6.40(1 ) -12 -12
Torr
(Gray)
75
-0.15(l)
(Grange)
3.4
This case is for illustration.
n solution”
lumber of jarameters fit
'osition' :hift' :ntensity¶ Ialfwidth" ,ine Mixing"
'osition‘ :hift7 ntensity' lalfwidth" ,ine Mixing"
j.4 Torr '5 Torr !50 Tort!3O Torr iO0 Tcrr
;lope2
Parameter
Table
710
D. Chris Benner et al
Figure l(a) displays five spectra of the P(6) manifold of the vj band of “CH, broadened by air and Table 1 lists some of the numerical results derived from fits to these spectra. The 0.01~cm-’ resolution spectra were obtained with the Fourier transform spectrometer at the National Solar Observatory’s McMath-Pierce solar telescope facility on Kitt Peak and a 5 cm path of methane. These spectra were used by Benner et alz9 to derive the widths and shifts of the methane spectral lines. Figure l(b) shows the expanded residuals (observed minus calculated) when each spectrum is fitted individually to between 20 and 26 coefficients (spectral line positions, intensities and/or halfwidths for up to 10 different spectral lines; continuum level and slope; spectrometer resolution; spectrometer phase error and/or zero level). The retrieved continuum slope is presented in Table 1 as an example of a spectrum parameter. This slope is caused primarily by the relative spectral sensitivity of the spectrometer itself and by the experimental setup, but is not very sensitive to the gas amount or other physical conditions such as temperature or pressure. The single spectrum analyses determine this parameter well because the fitted spectral interval is much wider than the manifold containing the strong methane lines and the wings of the lines are sufficiently weak that even at the higher pressures they have little effect upon the ends of the spectral interval. The change in the continuum slope in the multispectrum fit is within the estimated uncertainty of the single spectrum solution. Changing models to use a nonvoigt line profile (as discussed below) also had little effect upon the continuum slope. If a shorter spectral interval had been chosen, significant correlations between this parameter and parameters such as the halfwidths of the spectral lines would have made the slope more uncertain. In Fig. 1, the spectral feature at 2958.67 cm-’ consists of two lines only partially resolved even in the low pressure spectra. These two spectral lines are further detailed in Table 1. From the low pressure spectra the Lorentz halfwidth and pressure shift coefficients cannot be determined well because of their small effects upon the spectra. The higher pressure spectra, however, contain considerable information about these two parameters. The line parameters (position, intensity, Lorentz halfwidth and pressure shift) are not determinable, however, due to their strong correlations with each other. A nonunique fit is displayed in the table. The halfwidth of the higher wavenumber line was assumed to be 0.0625 cm-‘/atm29 and the halfwidth of the lower wavenumber line was similarly constrained to 0.0543 cm-‘/atm on the two lowest pressure spectra. The intensities derived for the four lowest pressure spectra in the single spectrum analyses are all within 5% of the mean. The highest pressure spectrum, however, gives intensities which are in error by a factor of about 1.5. The sum of the intensities of these two lines is close to that obtained in the other spectra, but the partitioning of the intensity is inaccurate. An attempt to match the wings with the halfwidth of one of the lines fixed to an inaccurate value moved the line positions within the blend which, in turn caused the repartitioning of the intensity. The halfwidth of the lower wavenumber line is also affected. The position of that line is not consistent with values derived from the fits to the other spectra in that the relationship between the position and the pressure is not very linear. As a result, the derived zero pressure positions and pressure shifts are poorly determined. The estimated uncertainties (see the following section) are also overly optimistic since the assumed model includes the uncertain value of at least one of the halfwidths of these two lines. By simultaneously fitting the five spectra, the halfwidth and pressure shift coefficients can be determined since the zero pressure line position and the line intensity are well constrained by the low pressure spectra. Furthermore, the five spectra are fit to a total of 47 coefficients rather than the total of 120 coefficients for the five individual spectrum fits (each spectral line is required to have only one value of each parameter which is valid for all five spectra). This reduces the overfitting that caused the residuals to look so small on the single spectrum fits [Fig. l(b)] even when the line parameters were poorly determined. The zero pressure position and intensity are now determined to an accuracy consistent with that of the low pressure spectrum fits. The pressure shifts and halfwidths are determined well for both lines. The drawback, however, is the much poorer fit to the spectrum [Fig. l(c)]. It is clear that there is something which is modeled incorrectly. An additional multispectrum fit was tried that included approximate corrections for collisional narrowing36 and line coupling. 37The residuals, displayed in Fig. I(d), from this fit show a marked improvement, but are still above the noise level. We believe that the remaining residuals are due to the inadequacy of the Rosenkrantz approximation of the line coupling and to small uncertainties
Nonlinear least squares
711
in physical parameters which are masked in a single spectrum fit by adjustment of the spectral line parameters. Notice, though, that the spectral line parameters have changed from those derived with Voigt profiles. In some cases (such as the halfwidths) the changes are far greater than the estimated uncertainties. This demonstrates the principle that the internal error estimate does not include uncertainties in the modeling. An additional advantage to fitting several spectra simultaneously is the ability to solve directly for parameters that are not determinable from a single spectrum fitting. Examples of these parameters are pressure shifts and temperature dependent parameters such as lower state energies of the transitions, temperature dependence of the Lorentz halfwidth and temperature dependence of the pressure shift coefficients. For example, the expression for the Lorentz halfwidth of a line is given in Eq. (45) of the Appendix. For a single spectrum with a single path there is a single value of the halfwidth that applies to the spectrum. It is, therefore, only possible to derive a single constant in that equation (generally at, the Lorentz halfwidth at a standard pressure and temperature). If, however, there are two or more spectra in the solution with different pressures and temperatures, there is a halfwidth appropriate for each spectrum and it is possible to derive as many halfwidth parameters as there are spectra. Since the pressures and temperatures are known there are only two parameters to be derived, namely at and n, the temperature dependence exponent. We have demonstrated the technique of deriving both parameters simultaneously (as well as position, pressure shift and temperature dependence of the pressure shift simultaneously) with the multispectrum technique in Ref. 38. Determination of the temperature dependence exponent of the halfwidth using a spectrum-by-spectrum analysis method requires a fitting of the results obtained from each of several spectra to Eq. (45) of the Appendix.14 This method can give inaccurate results unless the proper uncertainties in the retrieved parameters obtained from each spectrum are applied. Even so, different correlations among parameters on different spectra may constrain a multispectrum solution in ways not easily obtainable otherwise. ERROR
ANALYSIS
Error analysis is a difficult problem in least squares fitting of spectroscopic data.39-4’ A formal error estimate for each parameter can be obtained by the method of Whittaker and Robinson.42 This assumes that all second (and higher) derivatives of each fitted point with respect to each fitted parameter are zero and the data points are independent of each other. It also assumes no error in the modeling of the spectrum. These errors are reasonable estimates of the random error if none of the parameters is strongly correlated with any of the others and there are no significant modeling uncertainties. Off-diagonal elements of the covariance matrix can be used to calculate a correlation coefficient between each pair of parameters. If two parameters are strongly correlated (or anticorrelated), one or both may need to be constrained to a reasonable value known from other considerations such as theoretical constraints. This method for calculating errors is better than error calculations deduced from parameters obtained from fits to individual spectra (e.g., deriving broadening coefficients by least squares fitting measurements of Lorentz widths at a series of pressures). However, these derived uncertainties include only the uncertainties due to the noise level of the spectrum and the parameters derived in the solution. Not included are uncertainties due to errors in the measurement of the experimental conditions, the calculation of the partition function, incomplete modeling of the spectrometer and spectral line parameters not fitted in the solution. These unaccounted uncertainties often dominate a multispectrum solution since a large number of fitted points can decrease the formal error estimates to small values. CALIBRATION In a single spectrum fit the wavelength scale calibration is not required before fitting the data. The line positions may be corrected afterwards when an improved calibration is available. If absolute line positions are not needed, only a correct relative calibration is required to derive pressure shift information.
712
D. Chris Benner et al
The multispectrum technique, however, requires at least a relative calibration of the wavelength scale before the spectra are fitted. The relative wavelength scales must be precisely aligned or else the spectral lines will not be positioned properly on each spectrum. An absolute calibration adjustment can be applied after the fit to the spectral line positions if a better absolute calibration becomes available. Experimental parameters such as temperature, pressure, path length and mixing ratio must also be known to sufficient relative accuracy before performing multispectrum fits. Usually this means that careful laboratory measurement of these parameters is required. However, it is sometimes possible to use the spectra themselves to improve the relative calibration. For example, in a group of spectra obtained for measuring Lorentz foreign broadening coefficients, but not intensities of spectral lines, small mixing ratios are commonly used. If the mixing ratios are not known accurately enough, the derived intensity of the same spectral line will be different from spectrum to spectrum. Since the mixing ratios are small, the widths may not be significantly affected by self broadening. Thus, small adjustments to the nominal mixing ratios can make the mass paths of the spectra consistent. This is achieved when fits including all spectra no longer include intensity related residuals, Then the derived intensities of the spectral lines will be in error by a constant factor (due to uncertainty in the absolute mixing ratios). The derived widths, however, will be correct except for the usually insignificant error due to the self broadening correction. IMPLEMENTATION
The multispectrum nonlinear least squares fitting technique has been implemented into a FORTRAN program that runs on an IBM 80386 or 80486 compatible computer. It may be run in a batch mode, a completely interactive mode or a mode in which the initial input is completely specified as in the batch mode, but further modifications are allowed as in the interactive mode. In the interactive mode the user can modify the solution after a specified number of iterations or after meeting one or more convergence criteria, whichever occurs first. Changes entered through a series of menus allow the user to modify fixed parameters, change parameters for which a solution is found, add or delete parameters from the solution (such as line position or continuum level), add spectra to the solution and/or add or delete spectral lines. A backup option allows discarding the last set of iterations if the results were unsatisfactory. Continuation of the solution proceeds with a specified number of additional iterations if requested. Results of the error analysis, graphic output of the fits and output of derivatives are all supplied to help the user in evaluating the quality of the solution. Computation of the monochromatic spectrum is performed by standard line-by-line methods. Voigt profiles (including self and foreign pressure broadened Lorentz halfwidths), pressure shifts (self and foreign pressure induced), temperature dependence of Lorentz halfwidths and pressure shifts, temperature dependence of the line intensity (lower state energy of the transition), Dicke narrowing and line coupling (Rosenkrantz approximation) have been implemented. Calculation of the total internal partition sums is done by the method of Gamache et a1.43These are calculated separately for each gas and molecular isotope included in the solution. Flexible modeling of spectrometer-induced distortions is also included. Zero level, channel spectra and a polynomial fit to the continuum may be fitted to each experimental spectrum. For a Fourier transform spectrum, the effective apodization and phase error may be fitted. The maximum path difference of the interferometer and the solid angle of the entrance aperture are supplied by the user. For a tunable diode laser spectrum, the time constant of the electronics and the laser spectral width may be fitted. This multispectrum fitting program allows more than one absorbing path to be included in the solution. The physical conditions of these paths may be quite different and, therefore, the internal partition sums are calculated separately for each path, gas and isotope combination. An example of this situation is displayed in Fig. 2. These spectra were also obtained using the same Fourier transform spectrometer of the National Solar Observatory. The cell contained methane which does not have any appreciable absorption in this spectral region. A residual water vapor line arising from both within the evacuated interferometer tank and a dry nitrogen purged atmospheric path between the source and the interferometer entrance slit (in which the methane cell is placed) is present.
Nonlinear least squares
713
Table 2. Derived spectral line parameters from Fig. 2.’ Zero
Pressure
Line Position
Halfwidth
Coefficient
N,-Induced Pressure Shift Coefficient _.___
The water absorption caused by the residual gas in the interferometer tank shows as a narrow line. Not quite centered on the narrow line is the nitrogen-broadened water vapor line arising from the purged atmospheric path, As a result, these spectra provide N,-induced pressure shifts without the uncertainty of wavelength calibration. Several of these spectra are included in the solution (relative calibration of wavelength scale is necessary in this case) with different amounts of water in each spectrum. The zero pressure spectral line position, the N,-broadened Lorentz halfwidth coefficient and the N,-induced pressure shift coefficient are determined and given in Table 2. The amount of water vapor in each path is not well known and there is some unmodeled modification of the low pressure line shape due to the changing column of water in the path of the movable mirror of the FTS. These modeling problems cause the error analysis to produce values somewhat smaller than the true accuracy of the measurement. Indeed, the estimated uncertainties of the position and halfwidth are near the calibration accuracy of the wavenumber scale and the knowledge of the pressure in the purged atmospheric path. This illustrates the inadequacy of the error analysis to deal with systematic sources of errors when there are enough data included in the solution that the random sources of error are reduced too far. FUTURE
DIRECTIONS
Experimentally derived spectral line parameters are often used to determine more fundamental molecular parameters. For example, zero pressure positions of many lines in a single band may be combined to obtain vibrational energy and rotational constants of a molecule and the intensities of the same lines fit to a model with only a small number of free parameters. Ideally, residuals (observed minus calculated spectra) should be minimized not for the best derived spectral positions and intensities, but for the best vibrational energies, rotational constants, band intensity parameters and width parameters. This has been done to a greater or lesser extent for a single spectrum analysis by the whole band technique of Hoke and Shaw4 and by Pine and Looney.4s Carrying this idea even farther, if the spectra contain a sufficient number of vibration-rotation bands, the residuals should be minimized for the best spectroscopic constants and other theoretical parameters ultimately required. This involves fitting larger spectral intervals than those we have implemented so far and would require large amounts of computer processing time. On the other hand, theoretical constraints can also be imposed upon the spectrum on a much smaller scale. This might include constraints like equal intensity for doublets or the same spectrum parameters (such as channeling) for two spectra recorded close in time. Acknowledgements-Many people have contributed to this effort at one time or another. Donald J. Richardson of Hughes STX did some of the original programming of the single spectrum technique. Robert C. Reisweber of Computer Sciences Corporation assisted with information on previously existing software for the least squares technique, fast Fourier transforms and error analysis. C. Plymate and J. Wagner of the National Solar Observatory provided assistance in obtaining the data used for the examples and G. Ladd of the National Solar Observatory and C. H. Sutton of ST Systems Corporation helped in the initial processing of the data. Research at the College of William and Mary is supported by cooperative agreements with the National Aeronautics and Space Administration. The National Solar Observatory is operated by the Association of Universities for Research in Astronomy under contract with the National Science Foundation.
D. Chris Benner et al
714
REFERENCES I. Y. S. Chang and J. H. Shaw, Appl. Spectrosc. 31, 213 (1977). C. P. Rinsland, D. C. Benner, D. J. Richardson, and R. A. Toth, Appl. Opt. 22, 3805 (1983). L. R. Brown, J. S. Margolis, R. H. Norton, and B. D. Stedry, Appl. Spectrosc. 37, 287 (1983). J. W. C. Johns, J. Mol. Spectrosc. 125, 442 (1987). V. Dana and A. Valentin, Appt. Opt. 27, 4450 (1988). V. Dana and J.-Y. Mandin, JQSRT 48, 725 (1992). V. Dana, J.-Y. Mandin, C. Camy-Peyret, J.-M. Flaud, J.-P. Chevillard, R. L. Hawkins, and J.-L. Delfau,
2. 3. 4. 5. 6. 7.
Appl. Opt. 31, 1928 (1992). 8. V. Malathy Devi, C. P. Rinsland, and D. C. Benner, Appl. Opt. 23, 4067 (1984). 9. C. P. Rinsland and D. C. Benner, Appl. Opt. 23, 4523 (1984). IO. C. P. Rinsland, D. C. Benner, and V. Malathy Devi, Appl. Opt. 24, 1644 (1985). 11. D. C. Benner and C. P. Rinsland, J. Mol. Spectrosc. 112, 18 (1985). 12. V. Malathy Devi, C. P. Rinsland, M. A. H. Smith, and D. C. Benner, Appl. Opt. 24, 2788 (1985). 13. C. P. Rinsland, D. C. Benner, and V. Malathy Devi, Appl. Opt. 25, 1204 (1986). 14. C. P. Rinsland, V. Malathy Devi, M. A. H. Smith, and D. C. Benner, Appl. Opt. 27, 631 (1988). 15. V. Malathy Devi, C. P. Rinsland, M. A. H. Smith, and D. C. Benner, Appl. Opt. 27, 2296 (1988). 16. M. A. H. Smith, C. P. Rinsland, V. Malathy Devi, D. C. Benner, and K. B. Thakur, J. Opt. Sot. Am. B 5, 585 (1988). 17. D. C. Benner, V. Malathy Devi, C. P. Rinsland, and P. S. Ferry-Leeper, Appl. Opt. 27, 1588 (1988). 18. C. P. Rinsland, V. Malathy Devi, M. A. H. Smith, and D. C. Benner, Appl. Opt. 28, 211 (1989). 19. V. Malathy Devi, D. C. Benner, M. A. H. Smith, and C. P. Rinsland, Appl. Opt. 30,287 (1991). Erratum Appl. Opt. 30, 2928 (1991). 20. C. P. Rinsland, M. A. H. Smith, V. Malathy Devi, and D. C. Benner, J. Mol. Spectrosc. 150, 640
(1991). 21. C. P. Rinsland, M. A. H. Smith, V. Malathy Devi, and D. C. Benner, J. Mol. Spectrosc. 150, 173 (1991). 22. M. A. H. Smith, C. P. Rinsland, and V. Malathy Devi, J. Mol. Spectrosc. 147, 142 (1991). 23. C. P. Rinsland, A. Goldman, M. A. H. Smith, and V. Malathy Devi, Appl. Opt. 30, 1427 (1991). 24. V. Malathy Devi, D. C. Benner, C. P. Rinsland, and M. A. H. Smith, JQSRT 48, 581 (1992). 25. M. A. H. Smith, C. P. Rinsland, V. Malathy Devi, and D. C. Benner, Spectrochim. Acta 48A, 1257 (1992). 26. V. Malathy Devi, D. Chris Benner, M. A. H. Smith, and C. P. Rinsland, J. Mol. Spectrosc. 155, 333 (1992). C. P. Rinsland, M. A. H. Smith, V. Malathy Devi, and D. C. Benner, J. Mol. Spectrosc. 156, 507 (1992). V. Malathy Devi, D. C. Benner, M. A. H. Smith, and C. P. Rinsland, J. Mol. Spectrosc. 157, 95 (1993). D. C. Benner, V. Malathy Devi, M. A. H. Smith, and C. P. Rinsland, JQSRT 50, 65 (1993). A. Goldman and R. S. Saunders, JQSRT 21, 155 (1979). E. Niple, W. G. Mankin, A. Goldman, D. G. Murcray, and F. J. Murcray, Geophys. Res. Lett. 7, 489 (1980). 32. C. P. Rinsland, M. R. Gunson, J. C. Foster, R. A. Toth, C. B. Farmer, and R. Zander, J. Geophys. Res.
27. 28. 29. 30. 31.
96, 33. M. 34. K. 35. D.
1057 (1991). Carlotti, Appl. Opt. 27, 3250 (1988). Levenberg, Quart. Appl. Math. 2, 164 (1944). W. Marquardt, J. Sot. Indust. Appl. Math. 11, 431 (1963).
36. A. S. Pine and J. P. Looney, J. MoI. Spectrosc. 122, 41 (1987). 37. 38. 39. 40. 41. 42.
43. 44. 45. 46.
B. Gentry and L. L. Strow, J. Chem. Phys. 86, 5722 (1987).
V. Malathy Devi, D. C. Benner, M. A. H. Smith, and C. P. Rinsland, JQSRT
51, 439 (1994).
J. H. Shaw, N. Tu, and D. L. Agresta, Appl. Opt. 24, 2437 (1985). V. Dana, J.-Y. Mandin, and A. Hamdouni, Appl. Opt. 31, 1937 (1992). C. L. Brummel and L. A. Philips, J. Mol. Spectrosc. 159, 287 (1993). E. T. Whittaker and G. Robinson, The Calculus of Observations, Chap. 9, Blackie, Glasgow (1926). R. R. Gamache, R. L. Hawkins, and L. S. Rothman, J. Mol. Spectrosc. 142, 205 (1990). M. L. Hoke and J. H. Shaw, Appl. Opt. 21, 929 (1982). A. S. Pine and J. P. Looney, J. Chem. Phys. 93, 6942 (1990). F. Schreier, JQSRT 48, 743 (1992). APPENDIX Calculation
of Derivatives
Computationally, the most time consuming section of the nonlinear least squares technique is the calculation of derivatives. Often, these derivatives are calculated by the method of finite differences (defined below). In this Appendix some techniques that improve the accuracy of the derivatives and/or speed of the calculation significantly over this procedure are presented.
Nonlinear least squares
715
The spectrum S as a function of wavenumber, v, is calculated from the continuum B(v), channel spectrum C(v), gas spectrum convolved with the instrumental line shape b(v) and zero level offset z as: S(v) = B(v)C(v)b(v)
+ z.
(1)
The derivatives may be obtained numerically by using finite differences.
WV, 4) _ S(v,4 a4 -
+
6) - S(v,4) .!I
Here q is the parameter with respect to which the derivative is to be taken and c. is small enough that the derivative of S with respect to q is not significantly different at q and q + L, but large enough to avoid precision problems in the subtraction. In our approach, S(v) is first calculated with the current value of each parameter to be derived. Storing several intermediate steps of this computation allows many steps to be skipped in some derivative calculations. One step stored from the calculation of S(v) is the sum of the cross sections [defined in Eq. (24)] as a function of v for the spectral lines whose line parameters are all fixed in the analysis. This summation need not be repeated unless there is a change in these line parameters. Spectrometer Efects The derivative of the calculated signal with respect to each spectrometer effect is always possible without the complete recalculation of b(v). The recalculation of J(v) as defined in Eqs. (21) and (24) is never necessary since J(v) is always the same as it is in the nominal case. Also, in every case the derivative is only nonzero for the one spectrum to which the parameter applies. A discussion of each derivative in our implementation is given here. Zero level
The derivative of the calculated signal S(v) with respect to the zero level offset is trivial.
Continuum level The derivatives with respect to the continuum are expressed as:
as(v) = a4
-
C(v)b(v)
WV) at7 .
-
C(v) and b(v) are independent of the continuum parameters. Most of the computation is involved in the calculation of b(v). When b(v) is computed for the nominal case, it may be stored and retrieved at this point. C(v) could be treated in the same way, but recomputation of C(v) is not too time consuming. The continuum may be modeled as a polynomial of order N in terms of the starting wavenumber of the fitted interval, vStartrand a series of coefficients. B(v) = 2 Uv II,=0 The parameters to be derived [corresponding
(5)
to q in Eq. (4)] are B,. = (v -
gp
- v,,,,)“‘.
vstarty.
m
If there is more than one coefficient in the solution, it can be shown that the following rapidly calculable relationship holds.
D. Chris Eknner et al
716
Channel spectra
For the derivatives with respect to the channel spectrum parameters, a relationship similar to Eq. (4) may be obtained. y
= B(V)b(“) y
.
(8)
Again, b(v) is by far the most time-consuming part of the equation for a computer, but has already been computed for the nominal case using Eq. (1). B(v) has also been calculated for the nominal case, but can be quickly recomputed. The channel spectra are computed using an equation equivalent to that of Niple.3’ C(V) = P(V)=/{1
~(v){Re[b(v)lcos[e(v)l- Im[Wv)lWVv)l),
-Tca,[
1 -COS$(V
-&)]~+{~~~aiSin$(V
(9) -&))‘y
(lo)
(11)
khan
z, aisin
Nchiln 1 ai i= I is the number of channel spectra and ai, Pi and di are the amplitude, period In these equations NChan and phase of channel i. The derivatives with respect to the amplitude of channel i can be readily calculated. 1-
WV) -=~[j?(cosd,aai
l)-sind,]+C(v) 3
{Im[b(v)]sin d, + Re[b(v)][cos di - l]},
(12)
(13) p = tan
8,
(14)
p(v){Re[b(v)]sin 8 + Im[b(v)]cos 0} II=
9
P(1 + P’) p =
(15)
1 - C a,(1 - cos di).
(16)
i=l
The derivatives with respect to the period of channel i can also be determined.
at(v) diai
-=F api
, 1
rl(cosdi+/Isindi)+
sin di - 5 cos di)
5=%Z4” C %(P
ai sin di.
1,
(17)
(18)
i=l
The derivative with respect to the phase of channel i is particularly simple if the derivative with respect to the period of the same channel was calculated.
acw pi aqv) -= V+aP,' ah
(19)
Instrumental line shape
A derivative with respect to a parameter such as Fourier transform spectrometer phase error involved in the convolution with the instrument profile or for a spectral line parameter is again analogous to Eq. (4).
as(v) = a4
-
B(v)C(v)
WV)
-.
a4
(20)
117
Nonlinear least squares
B(v) and C(v) are simple enough that they can be easily recomputed in case they are not saved from the evaluation of the nominal case. The function b(v) can be expressed as a convolution of the line profile, w(v), and the monochromatic spectrum, J(v). Jc w(v)J(v) dv. s --3cI
b(v)=
(21)
For parameters contributing to w(v) there are two approaches possible for determining the derivative. First, note that the partial derivative of b(v) may be expressed as a modified convolution with the partial derivative of w(v) replacing w(v).
s
WV) a9
= +
-=
J(v)-
wv) a4
dv
(22)
.
If the partial derivative of w(v) is known, this convolution is accurate and, because the values of J(v) can be stored from the nominal case, faster than Eq. (2). If the partial derivative is not known, Eq. (2) may be used. Considerable computation time may be saved even here since J(v) need not be recomputed. Monochromatic Spectrum For parameters derived for b(v).
contributing
to the monochromatic
spectrum,
a similar integral can be
dJ(v>dv, w(v)---co 84 s
db(v) -= a4
=
(23)
This leaves only the determination of aJ(v)/aq. In most cases this can be handled without resorting to the method of finite differences [Eq. (4)]. Except in complex cases (for example, line mixing cases not simplified by too many approximations), the monochromatic spectrum can be expressed in terms of a line shape for each individual spectral line using Eq. (24). (24) A,, = mass path of gas i in path j, ZGk= intenstiy of line k of gas i in path j, rcVk = line shape factor for line k of gas i in path j. Here the Cpntity in the solution.
A,,z#
Kuk(v)
is the cross section. For laboratory spectra only ZVkand icOkare varied
Intensity The intensity Z,,,,can be expressed as a function of the line intensity Zoa at temperature T,,, the lower state energy Ei, the temperature q, the second radiation constant C, and several other factors (partition function, Honl-London factors, etc.), R,, which do not depend upon either Zoik or E;;.
The derivative with respect to the intensity of a line at temperature
aJ(v> -J(v)~~~A,Q~(v)$, -= azoik
aJ(v)
-=
azOik
,rA,hjkK,k(“).
J(v) N ‘11, Olk
(26) Olk
j=l
-r
TO can now be calculated.
(27)
, - I
Since J(v) and Key are calculated for the nominal case, the right side of this equation may be readily evaluated. The lower state energy Ei in the evaluation of S(v) only contributes to Zijk. The partial derivatives of J(v) with respect to E;; may then be calculated similarly.
718
D. Chris Benner et al
aJ(v) -= aE;;.
- J(v)C, ‘f” AijIUk
(28)
j=l
In the special case that there is only one path or that all paths have the same temperature T, these two derivatives are related closely. (29) Placing this relation into Eqs. (20) and (23), the following computationally found.
rapid relationship is
(30) Line shape
Some spectral line parameters (such as position and width parameters) contribute to only the line shape r$(v). aJw _=
-J(v)JkNfh/l,y. j=
a%k
I
(31)
I
A combination of the real and imaginary parts of the complex error function is used for most line shapes. For example, the Voigt function is simply the real part of the complex error function. The complex error function depends on only two parameters, x and y. The functions K(x, y) and L(x, y) are then defined in terms of the complex error function E(x,y). E(x, y) = KG, y) + iL(x, y).
(32)
Schreier46 reviews many of the techniques that have been proposed for calculation of E(x, y) and shows that the partial derivatives of K(x, y) and L(x, y) with respect to x and y can be expressed in terms of K(x, y), L(x, y), x and y only. The following relations can be derived from his equations.
W-T Y) ax
aw, Y 1 =2 ay
(33)
=~[YW,Y)-~W,Y)I,
yK(x, y) + xL(x, y) - -!-
1
Jn '
(34)
(35) aw,~)
ay
aK(x,~)
=ax’
(36)
Line position and pressure shift
The line position Vline is expressed in one of the models of Smith et alZSin terms of the zero pressure position, vo, pressure shift at temperature To, 6, temperature dependence of the pressure shift, 6’, gas sample pressure P and temperature T. Vline= ~0 + [6 + 6 ‘(T
- To)]P.
(37)
For the Voigt profile, r$(V) is just K(x, y) for the appropriate gas, path and spectral line. The parameter y is independent of vline9 vo, 6, and 6’ while x depends upon v, Vline and the Doppler width of the spectral line, an. x =
Jln 2(VaD
Vline)
(38)
719
Nonlinear least squares
The derivatives of +.(v) with respect to v,,, 6 and 6’ are then expressible in terms of K(x,, , yijk) and L(x,, , yijk). Consider, for example, the form for vO.
aKuk (V> aK(Xvkhk 1
---=
aVO,L
)
axok av+
axuk
m%k
3 Yijk) aY,ik
04
ayr/k avO,k
(39) ’
The second term on the right-hand side of Eq. (39) is zero. This equation can be expressed in a more convenient form by using Eqs. (33) (37) and (38). -2JiZ aKvk(v) -----= (4) b,k L (xijk 7 Yyk > $k K(-%,k , yi/k )I. a%, %,k Similar derivatives may be found for 6 and 6 ‘. General relations among these three derivatives hold. aKvk(v)
=
aslk
p
aKvk(V> ’
dvo,i
(41) ’
Under certain conditions (such as a spectrum obtained with only one path containing the absorbing gas), even more general relationships hold. (43) (44) If pressure shifts induced by different gases are separated in the solution, there are similar relationships that can be formed with these three derivatives that will involve the mixing ratios of the corresponding gases. Lorentz
halfwidth
A similar procedure may be followed with the Lorentz halfwidth, CI~,and the variables contained in it, namely, the Lorentz halfwidth at temperature To, the pressure P and the temperature dependence of the halfwidth n.
For the Voigt profile an equation analogous to Eq. (39) may be obtained. Since x does not depend upon CX~or n, considerable simplification is possible. The relationship between y and these two parameters is given by: (46)
The derivative of the line shape factor again derived.
Key
with respect to tx[ for the Voigt profile case can be
Jn
The derivative with respect to the temperature derivative.
1 .
I’,,kK(x,k,,y,,)+x,L(x,,y,)-l
(47)
dependence of the halfwidth is related to this
(48)
If only one path contains the absorbing gas or all paths have the same temperature fundamental relation can be derived.
T, a more
(49)
D. Chris Benner et al
720
Relations similar to Eqs. (47)-(49) may be derived for cases in which broadening due to different gases are handled separately. These equations involve the mixing ratios of the gases. Doppler
width
The derivative with respect to Doppler width a: at temperature TO, is somewhat more complicated since both x and y depend upon a:, but this can be simplified by relating it to previously derived derivatives. (50) If all paths are at the same temperature T and pressure P (or the pressure shift is zero at T) or if there is only one path, a more general relation holds. as(v) I = -c aa “D,, 1[
In our implementation,
(V-
Vline,k > y
+
a& 91.
O,k
(51) f-d
we solve for the molecular mass, M, rather than for the Doppler width. a0
D
=v”
(52)
c
The derivatives with respect to a: and A4 can be related.
a~,,,(4 ah(v) aa& aM,=x abf,'
aa; _=-_ aM
(53)
4
(54)
2~’
An analogous expression is derivable for the partial derivative of S(v). Nonvoigt
spectral lines
Nonvoigt line shapes can be treated in the same way if Key can still be expressed in terms of the complex error function. For example, Pine and Looney36 define a line profile that includes pressure narrowing. K(V) = Re
E(x, Y + e)
1
(55) [ 1 - JneE(x, y + e) ’ E was defined in terms of K and L in Eq. (32). The parameter e defines the pressure narrowing. This function can be differentiated with respect to x and y to obtain the same derivatives previously defined for the Voigt function. In addition, the derivative with respect to e may also be determined. All these derivatives are expressible in terms of K, L, x, y, e and v only. Similarly, the line shape for the Rosenkrantz approximation to line coupling is given in terms of E(x,y) by Gentry and Strow.37 rc(v) = Re[E(x, v)] + PYIm[E(x, r)].
(56)
Again this can be differentiated with respect to x and y to obtain the derivatives derived above for the Voigt case or be differentiated with respect to Yik,,the line coupling coefficient. To combine collisional narrowing and line coupling, a first approximation is to use the complex form of Key from Eq. (55) for E(x, y) in Eq. (56). All derivatives can still be expressed in terms of K, L, x, y, e, Y and v. NOMENCLATURE
A-Mass path of gas a-Amplitude of the channel spectrum B-Continuum &-Continuum polynomial coefficient b--Gas spectrum convolved with instrument line shape C-Channel spectrum
Nonlinear least squares
C,-Second radiation constant c--Speed of light d-Auxillary channel spectrum variable E-Complex error function E “-Lower state energy e-Collisional narrowing coefficient of spectral line Z-Spectral line intensity Z&Spectral line intensity at standard temperature i-Square root of - 1 [Eq. (32) only] i-Index variable J-Monochromatic gas spectrum j-Index variable (over physical paths) K-Real part of the complex error function k-Index variable (over spectral lines) k-Boltzman’s constant [Eq. (52) only] L--Imaginary part of the complex error function M-Mass of the absorbing molecule m-Continuum polynomial order N-Polynomial order of continuum Nchan-Number of channel spectra Npath-Number of physical paths in the solution n-Temperature dependence of the halfwidth exponent P-Period of channel spectrum [Eqs. (10-19) only] P-Gas pressure q-Parameter to be determined by the least squares fit R-Collection of miscellaneous factors S-Spectral amplitude T-Gas temperature TO-Standard temperature w-Instrumental line profile x-Voigt x parameter Y-Spectral line coupling coefficient (Rosenkrantz approximation) y-Voigt y parameter z-Zero level offset cm-Doppler halfwidth of spectral line cc,-Lorentz halfwidth of spectral line c$-Lorentz halfwidth coefficient of spectral line B-Auxillary channel spectrum variable b-Pressure shift coefficient 6’-Temperature dependence of the pressure shift coefficient c-Small finite number q-Auxillary channel spectrum variable 8-Auxillary channel spectrum variable k--Spectral line shape function p-Auxillary channel spectrum variable v-Wavenumber v,i,,-Wavenumber position of spectral line V ,,,,,-Wavenumber of start of fitted interval v,-Zero pressure wavenumber position of spectral line <-Channel spectrum auxillary variable X-Ratio of circumference to diameter of a circle p-Channel spectrum auxillary variable b--Wavenumber of a channel spectrum peak
721