Abstract:
An apparatus and method for non-contact temperature measurement of an object, using least-squares-based multiwavelength pyrometry techniques. Radiances from an object are detected by a spectrograph/detector apparatus and are converted into electronic signals readable by a computer. The computer then operates on these signals as data to be curve-fit, using least squares analysis, to a predetermined theoretical function for the dependence of the radiance on the wavelength. When the computer has minimized the least-squares difference function, the computer identified a parameter representing the temperature and reports this value to the user, along with a collaterally calculated maximum error in the temperature estimate.

Description:
The Government has rights in this invention pursuant to Contract Number N00014-80-C-0384 awarded by the Office of Naval Research. 
    
    
     This is a continuation of copending application Ser. No. 07/320,888, filed on Mar. 9, 1989, now abandoned, which is a continuation-in-part of copending application Ser. No. 296,538, filed on Jan. 12, 1989, now abandoned. 
    
    
     BACKGROUND OF THE INVENTION 
     This invention relates to non-contact temperature pyrometry, and more particularly to least-squares-based multiwavelength pyrometry techniques. 
     Pyrometric non-contact temperature measurement of high-temperature sources has been known for at least sixty years. Such measurement makes use of the Planck Radiance Formula: 
     
         N(λ)=ε(λ)C.sub.1 λ.sup.-5 [exp/(C.sub.2 /λT)-1].sup.-1 
    
     where: 
     N is the spectral radiance emission, ε(λ) is the emissivity of the material at wavelength λ, T is the temperature, and C 1  and C 2  are Planck radiation constants. Strictly speaking, many pyrometry techniques use the Wien Approximation to the Planck Radiance Formula, which is: 
     
         N(λ)=ε(λ)C.sub.1 λ.sup.-5 exp(-C.sub.2 /λT) 
    
     By measuring the spectral radiance emission, N, at a wavelength λ, and by supplying the appropriate values for ε, C 1 , and C 2 , an estimate of the temperature of the source can thus be calculated. 
     Historically there have been two distinct techniques for such calculations. The older technique is sometimes known as ratio pyrometry. This technique involves measurement of radiance emission at a number of different wavelengths in the attempt to eliminate the emissivity term by making ratios of the measured radiances. The other, newer, method of pyrometry is known as multiwavelength pyrometry. In this method, again, measurements of the spectral radiance emissions of the object are taken at several wavelengths. The processing of this data is then performed by a variety of techniques, the most accurate of which have proven to be least-squares-based multiwavelength techniques. These involve fitting the radiance emission data to an assumed emissivity functional form. Ratio techniques have not, in general, provided adequately accurate temperature estimates for broad industrial usage. The often large inaccuracies of the ratio techniques have been attributed to the fact that they require unrealistic assumptions to be made about the nature of the emissivity, ε, in the Planck formula. Ratio techniques assume that both (in the case of two-color) or certain (in the cases of three-color or four-color) of the emissivities at the measured wavelengths be equal. Multiwavelength, particularly a least-squares-based, techniques have been somewhat more successful, largely because they more reasonably assume a wavelength-dependent emissivity function, rather than that all or some of the emissivities are equal. In fact, with certain materials, these techniques have proven to be accurate to within one percent. With other materials, however, results have been unsatisfactory. It has been assumed by previous investigators that the unsatisfactory results have been due to two sources: first, incorrect form or lack of sufficient complexity of the assumed emissivity model function; and second, so-called &#34;correlation effects&#34; due to the inability of curve-fitting routines to distinguish in certain circumstances between changes in emissivity and changes in temperature. 
     A need has therefore been felt for a method and device which can yield an accurate estimate of the temperature of a source object, along with providing a determination of the largest possible error in the temperature estimate. 
     SUMMARY OF THE INVENTION 
     In response to this longfelt need, the applicant has developed an apparatus and method for non-contact temperature measurement of an object, using least-squares-based pyrometry techniques. The invention includes a detector and a computer. 
     The detector apparatus, in the preferred embodiment, includes a spectrograph which responds to incoming radiances from the object, and separates these incoming radiances into small wavelength bands. These bands are then converted into electronic radiance signals by photosensitive devices, one for each wavelength band. The photosensitive devices provide electronic radiance signals of strength corresponding to the intensity of the incoming radiances at their respective wavelength bands. In this embodiment, the spectrograph is internally calibrated to scatter or eliminate undesirable or extraneous incoming radiances. 
     In another embodiment, a baffle and filter apparatus, operatively connected to the spectrograph, is required to eliminate detection of undesirable stray light and undesirable wavelengths of light from the object. 
     The computer responds to the electronic radiance signals generated by the photosensitive devices of the spectrograph. The computer operates on the raw data (i.e. the electronic radiance signals), along with a predetermined emissivity model function, to generate an estimated value of the temperature of the object, along with an estimate of the largest possible error in the value for the temperature. In one embodiment, the emissivity model function is input directly by the user. In another embodiment, the emissivity model function is selected by the user from a bank of potential emissivity model functions included within the computer. In yet another embodiment, the computer itself selects the appropriate emissivity model function from a bank of potential emissivity model functions. The emissivity model function is a function which models the wavelength dependence of the emissivity for the material from which the object is made. The computer is adapted to first multiply the chosen emissivity model function by the Plank Radiance Formula (or the Wien Approximation thereto, in the case of linear techniques) in order to generate a radiance/wavelength function. This radiance/wavelength function constitutes a model for the theoretical dependence of the incoming radiance on the wavelength of detection. The computer then curve-fits the data (i.e. the electronic radiance signals) with the theoretical model, radiance/wavelength function. The object is to minimize the difference between the respective radiance values predicted by the radiance/wavelength function and those detected electronic radiance signals incoming from the object. The computer is adapted to generate a least-squares difference function, which is related to the difference between the detected electronic radiance signals and the respective radiance values predicted by the radiance/wavelength function. The computer is then programmed to minimize this difference function by finding values for the unknown parameters in the radiance/wavelength function which provide a minimum for the difference function. The computer is then programmed to identify an estimated value of the temperature from a preselected one of the unknown parameters of the radiance wavelength function. Lastly, the computer is programmed to calculate an estimate of the largest possible error in the value for the temperature from the now determined values of the unknown parameters in the radiance/wavelength function. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a schematic representation of the apparatus configuration of the preferred embodiment, wherein the self-filtering spectrograph/detector apparatus is shown connected to the computer. 
     FIG. 2 is a schematic representation of the baffle/filter apparatus, operatively connected to the detector, which comprises a spectrograph and a photodetector housing. The detector is further shown connected to the computer. 
     FIG. 3 is a flow-chart, demonstrating the steps by which this invention determines the estimate of the temperature along with the largest possible error in the temperature determination. 
    
    
     DESCRIPTION OF THE PREFERRED EMBODIMENT 
     The preferred embodiment contains two devices: a detector and a computer. The detector, in turn, includes a spectrograph and a series of photodetectors. 
     The apparatus configuration of the preferred embodiment is represented in FIG. 1. A Photo-Research Spectra Scan PR-700 PC® or similar equipment is used as the spectrograph/detector (01). This device automatically filters out undesirable light and converts the incoming radiances into electronic radiance signals usable by the computer (02). FIG. 2 demonstrates the configuration of an alternative embodiment wherein the baffle/filter apparatus is operatively connected to a non-self-filtering spectrograph. Spectral radiance emissions emanate from the source object (10). These are then passed through limiting optical filters (11) and baffles (12) to scatter and eliminate undesirable wavelengths of light. A number of filters and baffles may be included and arranged in any order. The radiances then pass through focusing optics (13) of either reflective or refractive form. The focused beam then passes through dispersive optics (14). These can take the forms of either reflective optics (i.e., a spectrograph), or transmissive (i.e., narrow band filters). The beam may then be passed (though it is optional) through channeling optics (15) into the photodetector housing (16), inside of which are the photosensitive detectors. Lastly, the photodetectors convert the radiances into electronic radiance signals, which are conveyed to the computer (17) for processing. 
     FIG. 3 is a flow-chart describing the steps by which the preferred embodiment estimates the temperature of the source object from the spectral radiance it emits. Steps 25 and 26 are performed by the spectrograph in the preferred embodiment, and the baffle/filter apparatus in conjunction with the spectrograph in the alternate embodiment. Steps 27 through 32 are performed by the computer. 
     After the computer receives the electronic radiance signals (ERS), an appropriate emissivity model function (EMF) must be chosen (27). By the various embodiments, this function may be supplied to the computer by the user, or the user may select the function from a bank of functions predetermined within the computer, or the computer itself may be programmed to select the function from the bank. After the function has been chosen, the computer multiplies the EMF by the radiance/temperature relation (RTR), which for non-linear least-squares analysis is the Planck Radiance Formula, while for linear analysis the RTR is the Wien Approximation to the Plank Radiance Formula. The multiplication results in the Radiance/Wavelength Function (RWF) (28). 
     The EMF has a functional form of dependence on the wavelength, with unknown parameters as coefficients. The EMF can assume many different forms. Examples of the form might be linear, exponential, quadratic, or a Fourier series in the wavelength. Since the wavelength values will be known by the computer, effectively it is these unknown parameters which are the unknown variables in the function. Thus, for the EMF 
     
         EMF=ε(α.sub.0, α.sub.1, . . . , α.sub.m-1) 
    
     where ε is the emissivity, and the α&#39;s are the unknown parameters. This leads to 
     
         N.sub.i =ε(λ.sub.i, α.sub.0, α.sub.1, . . . , α.sub.m-1)C.sub.1 λ.sub.i.sup.-5 exp(-C.sub.2 /λ.sub.i T), 
    
     where N i  is defined as the theoretical value of the radiance at the ith wavelength, as predicted by the radiance/wavelength function; λ i  is the value of the ith wavelength band; T is the temperature; the α&#39;s are the unknown parameters; and C 1  and C 2  are the Planck Radiation constants. 
     The computer then establishes a difference function (DF), formed as the difference between the RWF and the input ERS values from the detector (29): ##EQU1## where n is defined as the total number of wavelength bands at which measurements of the radiance of the object were taken; m is the total number of unknown parameters; f i  is the measured value of the radiance at the ith wavelength band (i.e., the ERS corresponding to the ith wavelength); and dλ i  is the uncertainty in the value of the ith wavelength (i.e. the width of the band). 
     After establishing the DF, the computer is programmed to curve-fit the RWF to the ERS by finding those values for the unknown parameters, α 0 , . . . ,α m-1 , which minimize the DF (30). At this point, the user may have an option. If the EMF chosen was a linear exponential function of the wavelength, the user may opt for the faster, though less flexible, linear least-squares method of curve-fitting. Otherwise, the non-linear method is used. Both of these methods are now described. 
     I 
     Linear Least-Squares Method 
     The linear least-squares method uses linear algebra techniques to solve for the values of the unknown parameters. Consider the problem of fitting some data (b 1 ,b 2 . . . ,b n ) to a polynomial of degree m. If n&gt;m, we say that the problem is &#34;overdetermined&#34;. If, for example, n=3 and m=1 then we will have to solve the following system of equations 
     
         a.sub.0 +a.sub.1 Y.sub.1 =b.sub.1 
    
     
         a.sub.0 +a.sub.1 Y.sub.2 =b.sub.2 
    
     
         a.sub.0 +a.sub.1 Y.sub.3 =b.sub.3 
    
     Here, the unknown parameters are the a i  coefficients; because the equations are linear with respect to these parameters, we can derive the following equations. The overspecified system of equations is inconsistent and in general has no solution. One can get around this problem by defining an error function (E) and then minimizing it; more precisely, we will minimize the square root of the sum of the squares of the error. 
     
         E=SQRT[(a.sub.0 +a.sub.1 Y.sub.1 -b.sub.1).sup.2 +(a.sub.0 +a.sub.1 Y.sub.2 -b.sub.2).sup.2 +(a.sub.0 +a.sub.1 Y.sub.3 -b.sub.3).sup.2 ] 
    
     or in matrix form: ##EQU2## The minimum of E 2  will be found by solving the system of equations ##EQU3## and in matrix form: 
     
         A.sup.T AX-A.sup.T B=0 
    
     In general the method of Least Squares involves finding the minimum of ||AX-B||. Thus, 
     
         ρ=||AX-B|| 
    
     and at the minimum 
     
         ρ.sub.LS =||AX.sub.LS -B||. 
    
     From trigonometry, we can define ##EQU4## Then it can be shown that in the worst case ##EQU5## where u s  (t) is called the &#34;condition number&#34; of t with respect to s, and ##EQU6## A +  =Pseudo-Inverse of A, σ 1  =Largest Singular Value of A, 
     σ n  =Smallest Singular Value of A. 
     Equations (1) and (2) describe the effects of perturbations in the data B and the coefficient matrix A on the solution X. In the simplest terms the condition numbers can be thought of as amplification factors for the translation of a perturbation in B or A on to a perturbation in X LS . Whereas θ is a measure of how closely the equations AX can match the data B. 
     For the case being discussed here, perturbations in B would arise from noise in the measured data. Perturbations in A would be caused primarily by roundoff error during the LS calculations and would not be very significant unless A was terribly conditioned (i.e. had a very large condition number). 
     So far we have made no mention of the number of data points n that we are dealing with. In general the actual condition number for perturbations in the data is lower by an approximate multiple of n -1/2 , i.e., ##EQU7## It should be kept in mind that this is only an order of magnitude approximation. 
     The measured data (f i ) for Multiwavelength Pyrometry (assuming perfect measurements) is the product of two functions: 
     
         f.sub.i =ε(λ.sub.i)C.sub.1 λ.sub.i.sup.-5 [exp(C.sub.2 /λ.sub.i T)-1].sup.-1 
    
     where 
     λ i  =ith wavelength, 
     ε(λ)=Emissivity. 
     In general we do not know much about ε(λ), but we can make some simplifying assumptions: 
     i) ε(λ),ε&#39;(λ),ε&#34;&#34;(λ) are 
     continuous over some sufficiently small (λ 1 , λ 2 ) 
     ii) ε(λ) is single valued over some 
     sufficiently small (λ 1 , λ 2 ) 
     Given these assumptions we can approximate the emissivity with any suitable analytic function such as a polynomial in λ or an exponential raised to a polynomial in λ. This function f is to be modeled by a fitting function g. At this point, however, f is nonlinear in at least one of the unknowns, the temperature T. So we rearrange it, after making the Wien approximation, to get 
     
         ln(C.sub.1.sup.-1 f.sub.i λ.sub.i.sup.5)=lnε-C.sub.2 λ.sub.i.sup.-1 T.sup.-1 
    
     All the quantities on the left side of this equation are known and on the right side the unknowns are ε and temperature. If we define ##EQU8## the emissivity is then modeled by ##EQU9## and the unknown parameters are the m-1 coefficients (α 0 , α 1 , . . . , α m-2 ) of the m-2 degree polynomial and the temperature T. We now have g as a linear function of the unknown parameters. 
     Translating this into the terms used to describe the linear LS technique we get ##EQU10## 
     Using this method, the values of the unknown parameters X can be evaluated. As shown, the value for the temperature is inversely proportional to the value of α m-1  : 
     
         T=C.sub.2 /α.sub.m-1 
    
     The error is determined from the condition number and the closeness of fit as determined by θ. The larger the condition number, the more inaccurate the estimate of the temperature will be for a given noise level. 
     II 
     Non-Linear Least Squares Method 
     If the EMF is not a linear exponential function of the wavelength, non-linear least-squares methods are used. The technique is as follows. 
     Consider the case where one is trying to fit the data f with the function N. From the above ##EQU11## 
     One now wishes to evaluate the parameters α j  at the minimum of χ 2 . If N is nonlinear in the α j , then the system of equations resulting from differentiating χ 2  with respect to the α j  will be nonlinear; and linear techniques such as Gaussian elimination, cannot be used to solve this system of equations. One would then have to look for some sort of iterative technique to locate the minimum in χ 2 . 
     The solution technique used here is described by Bevington. The algorithm used is termed the Gradient Expansion algorithm and was first presented by Marquardt. This algorithm searches the surface defined by χ 2  over the m dimensions for a minimum. When χ 2  is large this method follows the gradient of χ 2 . As χ 2  gets smaller, the algorithm linearizes χ 2  by a Taylor expansion. Finally when χ 2  is very small the algorithm does a linear fit to the data. How easy this task will be depends upon the shape of the surface defined by χ 2 . 
     To study the behavior of χ 2 , we expand it in a Taylor Series around the minimum. We start by defining the following ##EQU12## J&#39;=Gradient at β=β&#39;, H&#39;=Hessian at β=β&#39;, 
     
         δβ=β-β&#39;. 
    
     Assuming χ 2  is quadratic in a sufficiently small neighborhood around β=β&#39;, then, expanding χ 2  in a Taylor series, we can write 
     
         χ.sup.2 (β)˜χ.sup.2 (β&#39;)+J&#39;.sup.T δβ+1/2δβ.sup.T H&#39;δβ. 
    
     If β&#39; is an unconstrained minimum and the function N is a reasonable model for the data f, then J=0 and χ 2  (β&#39;)≈0 by definition. The above equation then reduces to 
     
         χ.sup.2 (β)≧1/2δβ.sup.T H&#39;δβ. 
    
     This inequality describes an m dimensional ellipsoid for each value ν of χ 2  (β). The ellipsoids are concentric and centered around δβ=0. The ellipsoid of interest covers a volume of uncertainty in the parameter space, i.e., we cannot determine any of the parameters with an error less than L multiplied by the smallest axis of the ellipsoid. The quantity L is a constant on the order of n 1/2 . We can predict L only in order of magnitude terms. 
     Our interest, however, is more in the magnitude of the largest possible error, i.e. we wish to locate the points on the ellipsoid farthest from the center. Thus we wish to find the vector δβ which maximizes 
     
         δβ.sup.T H&#39;δβ=2ν. 
    
     It can be shown, using the method of Lagrange multipliers, that the desired vector is an unnormalized eigenvector of H&#39; with an eigenvalue η. So, 
     
         δβ.sup.T δβ=2ν/η. 
    
     Each one of these eigenvectors forms an axis of the ellipsoid, and so the worst error will lie in the direction corresponding to the eigenvector with the smallest eigenvalue. 
     
         ||δβ||=[2ν/η.sub.smallest ].sup.1/2 
    
     In general, 
     
         ||δβ||≦[2ν/η.sub.smallest ].sup.1/2 
    
     and in particular 
     
         ||δβ|| is of the order of [2ν/nη.sub.smallest ].sup.1/2. 
    
     If the data f i  are measured with an uncertainty Δ 1 , then we can tolerate a change of ##EQU13## In χ 2  without having any reason to prefer one value of δβ over another as long as they all lie in the region define by ##EQU14## This means that the m-dimensional ellipsoid described above is now defined by ##EQU15## rather than just by ν. So the error in the parameter estimation is now given by ##EQU16## In general, ##EQU17## and in particular, 
     
         ||δβ|| is of the order of ##EQU18## 
    
     Eigenvalue calculations are difficult and sometimes unstable, so we wish to avoid them if possible. Fortunately, for the cases we deal with here, H&#39; is a Hermitian matrix. Hermitian matrices have the property that the singular values are equal to the absolute value of the eigenvalues. Additionally the singular vectors are equal to the eigenvectors if the matrix is positive semidefinite, otherwise, the singular and eigenvectors agree to within a factor of (-1) 1/2 . So, we can write the equation as: 
     
         ||δβ|| is of the order of ##EQU19## where σ.sub.smallest is the smallest singular value of H. 
    
     In the non-linear technique, the user is allowed to determine which of the unknown parameters contained in β is to serve as the parameter corresponding to the temperature estimate. In turn, the worst possible error, presented in percentage form is ##EQU20## 
     After the computer has determined the values for the temperature and error, these quantities are displayed to the user. 
     Other embodiments are within the following claims.