System for Identifying Known Materials Within a Mixture of Unknowns

This invention relates to methods and systems for identifying known materials or artifacts from within an unknown mixture. More particularly it can relate to systems and methods for identifying chemicals or biological constituents within a mixture of known and unknown substances; it can also relate to systems and methods for identifying shapes or patterns from a noise filled background. Still more particularly, it relates to an expert system which analyses the interaction of a first set of wavelengths of electromagnetic radiation on an unknown mixture in relation to a database of responses to that radiation of known constituents to provide concentrations of the known constituents in the mixture. Still more particularly, the expert system has two components to analyse the information received from the mixture describing the interaction of the radiation, a first subsystem utilizing a multivariate multidimensional patch algorithm trained by an evolutionary algorithm to provide the concentrations and a second subsystem utilizing a neural network also trained by an evolutionary algorithm to provide the concentration information.

Searching for the presence or absence of a known material within an unknown mixture is a classic problem with a variety of methods of solution. Most college students who took physical chemistry classes will remember the linear process-of-elimination techniques used in the laboratory to determine the unknown chemical from the sample presented to them for analysis. With the advent of the availability of increasingly powerful and inexpensive computing power, certain computer implemented processes for analyzing the problem have been brought to bear. These include various numerical analysis techniques (i.e. multivariate analysis, monte carlo techniques, etc.), genetic algorithms, and neural networks. Unfortunately, the nature of the problem presented herein is so computationally intense that solutions using these computer implemented processes are still beyond of the capability of readily available computers. What is needed is a new approach based upon but significantly more powerful than the presently known computer-implemented techniques to provide a fast and accurate solution to the problem.

This invention is an expert system wvhich employs either or both of two computer-implemented processes to analyze data taken which describes the responses of an unknown mixture to multiple wavelengths of electromagnetic radiation to determine the amount of a particular known constituent in the unknown mixture. Both processes employ a new type of genetic algorithm, hereinafter called an evolutionary algorithm, to train and optimize either a numerical analysis process or a neural network to identify the known constituent.

The evolutionary algorithm is more powerful and flexible than its antecedents. It reworks the standard mating and mutation techniques to achieve, very efficiently and rapidly, an optimized solution set (genetic sequence). For the numerical analysis, a new multidimensional multivariate analysis method has been created which, when combined with the evolutionary algorithm, can achieve heretofore unattainable solutions to complex identification problems. For the neural network, the flexibility the network has been increased in order that it can be efficiently optimized by the evolutionary algorithm.

In the numerical analysis techniques the evolutionary algorithm trains and optimizes the technique using a first set of wavelengths of electromagnetic radiation to interact with a known sample. The responses of the known sample are recorded at a second set of wavelengths. An optimized solution set is then created for each known constituent for which a search is desired. The unknown mixture is then irradiated with the same first set of wavelengths and data are recorded at the second set of wavelengths. The optimized solution sets for the various known constituents are then run on the data from the unknown mixture to determine the amount, if any, of each of the knowns.

In the neural network analysis, the neural network is trained for a given set of one or more known constituents where the concentrations of the members of the set change from generation to generation, but the members of the set remain the same. After some number of generations, the now-trained neural net has had its neurons and interconnection weights reconfigured by the evolutionary algorithm to provide a very fast and efficient measure of the concentrations of any of the members of the training set that are present in the unknown.

The expert system combines these two subsystems to provide two different pathways to the desired concentration information from the unknown mixture with very useful cross checking information as to the accuracy of the information.

FIGS. 1A and 1B are three dimensional graphs of the uv fluorescence spectra of acetone and benzene respectively with horizontal scales of 41 excitation wavelengths and 1024 fluorescence (return) wavelengths and a vertical scale for magnitude of the returns.

FIG. 2 is a chart showing the decomposition of a two dimensional patch into a one dimensional string, with the attendant loss of important information associated with the grouping with neighboring pixels.

FIG. 3 is a diagram showing how the summation of a number of fluorescence return patches for known chemicals multiplied by their concentrations plus noise add up to the fluorescence return from the unknown mixture.

FIG. 4 is a diagram showing the various constituents of the uncertainty of a concentration estimate.

FIG. 5 is a diagram showing several different patch strategies with a variety of allowed pixel groupings.

FIGS. 6A-6E are a series of graphs comparing the use of one dimensional and two dimensional multivariate algorithms.

FIG. 6A is group of three arbitrary chemicals whose gaussian-shaped spectral signatures are placed on a diagonal without noise.

FIG. 6B adds a contaminant whose concentration is 100 times higher to the graph of FIG. 6A.

FIG. 6C shows the results of running a one-dimensional multivariate analysis of the data from FIG. 6B at varying degrees of noise.

FIG. 6D shows the results of running a two-dimensional multivariate analysis for three patches of the data from FIG. 6C at varying degrees of noise.

FIG. 6E is a pixel plot of the three patches used in the analysis for FIG. 6D.

FIG. 7A is a graph showing the spectral response of a mixture of methanol, xylene and toluene without noise.

FIG. 7B is a graph which takes the information from FIG. 7A and adds up to 100% noise.

FIG. 7C is a plot of the output of the patch multivariate algorithm showing the successful identification of the three known constituents through the entire noise range.

FIG. 8 is a diagram of the steps of the evolutionary algorithm as applied to the multivariate patch technique with the upper left comer showing a genetic sequence containing at least information regarding pixels within a patch, the next step showing reproduction (mating) probability as a function of the fitness of a genetic sequence(determined by its success in identifying a chemical), the next step showing segments of the sequences being exchanged between mating sequences with offspring with altered sequences appearing in the next generation, in the next step some of the offspring also receive mutations to random genes within the sequence with the new generation repeating the process.

FIG. 9 is a diagram showing that portion of a gene sequence which controls the patching strategy by describing the patch boundaries.

FIGS. 10a-10F are a series of plots showing the evolution of a series of 24 patches on toluene, xylene and methanol.

FIG. 10A shows fitness level as a function of the generation for the average (darker shaded) and best (lighter shaded) sequences in each generation.

FIG. 10B shows the average patch size history over 500 generations.

FIG. 10C shows the patch locations on the pixel array for the average genetic sequence.

FIG. 10D shows the patch locations on the pixel array for a "best fitness" genetic sequence after 500 generations.

FIG. 10E is a plot of the training set with the magnitude returns associated with the pixels being represented by the shading.

FIG. 10F is a plot of the fitnesses of all 100 genetic sequences in the 500th generation with the raw fitness appearing as dots at the top edge of the plot and the weighted fitnesses appearing as asterisks.

FIG. 11 is a diagram of a typical neural network used herein with each pixel of the complex mixture return being uniquely connected to one neuron in the input layer, the various interconnections between the neurons in the various layers and the single outputs from the neurons in bottom layer, with a detailed inset of one of the bottom layer neurons with the weights associated with the input connections to it and its single output as the summation of the inputs multiplied by their respective weights.

FIGS. 12A-12F are a series of plots which show the evolution of a four layer neural network at the time of the 50th generation.

FIG. 12A is a plot of the training set used for this generation.

FIG. 12B is a plot of the average and best fitness scores for the genetic sequences in each of the 50 generations.

FIG. 12C is a plot of the raw and weighted fitnesses of the 100 genetic sequences at the 50th generation.

FIG. 12D is a plot of the weights associated with the pixel inputs in the top layer of the net at generation 50.

FIG. 12E is a plot of the weights associated with each of the inputs to each of the. neurons in the second layer of the net.

FIG. 12F is a plot of the weights associated with each of the inputs to each of the neurons in the bottom layer of the net.

FIGS. 13A-13F correspond to the plots of FIGS. 12A-12F but the data come from the neural network genetic sequences at the 3000th generation.

FIG. 14 is a flow chart showing the composition and operation of the elements of an expert system incorporating the multivariate patch algorithms and neural networks trained by the evolutionary algorithm.

FIG. 15 is a diagram illustrating how the concentrations of the various constituents add to form the mixture returns in a one-dimensional pixel array.

FIG. 16 is a schematic diagram of the hardware layout of the system.

This detailed description of the intention will begin with a discussion of one preferred embodiment involving recognition of chemicals within a distant plume on the basis of UV fluorescence when irradiated with a multispectral uv laser pulse. This was the first application developed for this invention. Following this will be a discussion of other uses of the invention and a concluding section with a detailed discussion of the mathematics involved in this invention.

One purpose of this invention is to estimate the individual chemical concentrations in a multicomponent mixture using a premeasured chemical catalog of in calibrated spectral signatures. It should be able to reject noise and optical contaminants such as scattered light, and be able to reject unknown or uninteresting chemicals.

A multispectral laser remote sensing system generates a range-resolved 2-D spectral signature of a remote object or plume. The 2-D spectra can be visualized as a surface: the surface height representing return signal intensity, with the excitation wavelength from the laser source along one axis and emission (or fluorescence) wavelength from the interrogated object along the other axis. Two typical uv spectra taken by a uv fluorometer are shown in FIGS. 1A and 1B, for acetone and benzene respectively. Note that, at least in the ultraviolet, spectral features can be broad and featureless: this is not necessarily so in the visible and infrared. The system of this invention can use both broad and featureless as well as sharply defined spectra.

When the object or plume is made up of many components, the returning 2-D mixture signal is assumed to be a linear superposition of the contributions from each component. The signal contribution from each component depends on its unique spectral shape amplified by its relative concentration.

Historically, algorithms designed to estimate concentrations of components in mixtures used a chemometric approach involving 1-D least-square minimizations (see "Numerical Recipes: The Art of Scientific Computing", W. H. Press et al., Cambridge University Press, 1986, p. 509) and in some cases normal-component analysis (see "Factor Analysis in Chemistry", E. R. Malinowsky and D. G. Howery, Krieger Publishing Co., 1989, p. 141). These techniques do tend to work well for many applications. Nonetheless these conventional methods do not utilize all the available information in spectra, such as the correlation of one pixel's intensity with another. A pixel here is defined as a unique point representing the intersection one of the wavelength bands from the first set of wavelengths from the laser source and one of the wavelength bands in the second set of wavelengths for the returns from the irradiated mixture. As an aside, it should be noted that while for this preferred embodiment the first wavelengths will differ from the second wavelengths, for measurements involving elastic backscattering effects this may not necessarily be the case.

For example when using a conventional 1-D method, the pixels of the mixture can be arbitrarily scrambled without changing the resulting concentration estimates, as long as the identical pixels in the database catalog are scrambled in the same way. In 2-D (and in higher dimensions) the spectra can be converted into a 1-D problem by stringing the pixels out, as shown in FIG. 2. In higher dimensions, such as in 2-D, the pixel-to-pixel correlations are even more structured to the human eye, and more complex in the sense that pixels along diagonals, and well-separated pixels can be correlated with each other.

In order to take advantage of the additional information available in spectral signatures I have developed new techniques of analyzing spectral signatures, including new ways of adaptively extracting minute chemical concentrations in complex multicomponent mixtures. The chemical recognition and analysis research has concentrated on multivariate methods, neural nets and genetic optimization. Other techniques can be productive as well, for example massively parallel implementations, accurate uncertainty estimates on concentrations, and intelligent techniques to search for chemicals when the database catalog becomes very large (.about.hundreds). These techniques have important uses that go well beyond laser remote sensing that will be discussed further below.

The Multivariate Approach

A core algorithm in the chemical recognition method and system is the multivariate `patch` algorithm. This algorithm computes the chemical concentrations and estimates the uncertainty of those concentrations. Conceptually the multivariate method is like a least squares fit, where the fitted parameters (the concentrations) are chosen to minimize the residuals between the measured unknown mixture and the fitted solution. The `patch` algorithm extends this approach in the sense that it minimizes the residuals on pixel sets which collectively contain the most important features in a particular chemical spectrum. Since neighboring pixels often contain valuable correlated information, this approach enables better concentration estimation and better noise rejection. The algorithm incorporates conventional 1-D chemometrics as a limiting case. Visually these correlated sets of pixels appear as patches when overlaid onto a 2-D spectrum, hence the name. This is illustrated in FIG. 3, where the mixture is on the left, and is equal to a sum of chemical samples multiplied by their respective concentrations, and may be contaminated with noise. One patch occupies the same pixels in the mixture and in the chemical components.

Mathematically, the `patch` algorithm computes the residuals between the mixture and a hypothetical solution, for all pixels in each patch, then minimizes this residual, or a function of this residual, for every patch independently. The user can chose to minimize the squared sum, the sum of the squared (as in least squares), or the summed absolute value (as in robust estimation) of the residuals in every patch. There are trade-offs associated with these choices. The system of equations is inverted by computing the pseudo-inverse of the patch-sum matrix using singular value decomposition. The concentrations are then computed when the pseudo inverse is multiplied by the mixture and component-sum matrix. The mathematical details are discussed in greater detail below.

After the concentrations are estimated, the algorithm then computes an uncertainty estimate. This uncertainty estimate is derived from the first three terms of the Taylor series expansion of the rate of change of concentration from both mixture 100, chemical database 130 and mixture-chemical (cross product) uncertainty 120. This is shown in FIG. 4. Conventional algorithms usually assume only the mixture has measurement uncertainty (term 100 in the figure), with the curious result that the uncertainty in the concentrations is computed from a covariance matrix that is independent of the mixture itself. In other words all mixtures would have the same uncertainty in the concentrations. The covariance matrix is the first term in the Taylor expansion, with the other new terms exhibiting the expected loss of confidence as the concentrations become relatively small or when there is a chemical missing from the database. I have also implemented a Monte Carlo uncertainty estimate, so that, if the mixture measurements do not obey a normal statistical distribution, I can still estimate concentration errors with a real system response.

Patches are lists of correlated pixels, but how are these lists chosen? This question of patching strategy is a complicated one. Patches can overlap, be sparse, or widely separated. In fact a patch does not have to made of contiguous pixels at all. Pixels outside of a patch are ignored. Although any reasonable patching strategy will yield very good results, for many problems the optimization of the patch algorithm is too labor intensive and too subjective. Some of the possibilities are shown in FIG. 5. The optimization of the patches will be discussed in the next section on genetic optimizations by the evolutionary algorithm of this invention.

As an example of the benefits of the flexibility of the patch algorithm, consider the following problem. A set of three gaussian-shaped `computer synthesized` chemicals are created with their peaks on a diagonal, and the spectra of these three chemicals are stored in the multivariate algorithm's chemical database. The image is made up of 40.times.40 pixels. Mixtures are computer created by combining these three chemicals with concentrations 0.5, 1.0, and 1.5 respectively. Then a fourth chemical is added to the mixtures, off diagonal, with a very large concentration of 100.0. Both the conventional 1-D algorithm and the `patch` algorithm are given this problem, with increasing (0% to 100%) random noise added. The results are shown in FIGS. 6A-6E. The 1-D algorithm is unable to reject the contaminate as noise is increased.

To demonstrate the algorithm on `real` data, FIGS. 7A, 7B, and 7C show the results of the `patch` algorithm identifing methanol, xylene, and toluene with relative concentrations of 0.5, 1.0, and 1.5 respectively, as up to 100% noise is added. The spectrum for each chemical was taken with a uv fluorometer, and the mixtures were generated by the computer by adding the spectra according to the concentrations, and then adding additional random noise. The random noise was uniformly distributed from--1.0 to 1.0 multiplied by the percentage. The size of the images was 41.times.1024 pixels. No attempt was made to remove the scattered light from the raw data. The code achieves good concentration estimates (accurate to a few percent) even at 100% noise.

Genetic Optimization of the Patch Algorithm

To overcome the need for trial and error optimization of the `patch` algorithm I have developed a novel Genetic Algorithm (GA) to optimize the lists of correlated pixels for any set of candidate chemicals. In a genetic algorithm a set of genetic-like sequences is created in which each sequence can completely describe a possible solution to a problem. Every genetic sequence, or `bug`, must compete with its peers on the problem, and the most fit are allowed to generate offspring for the next generation, with occasional mutation to introduce new genetic material into the population. Generation after generation the population relentlessly improves its fitness. Genetic algorithms are exceptionally good at search and optimization when applied to problems with very large multidimensional solution spaces. GA's, used carefully, are not easily trapped in local minima or maxima, a problem that can plague hill climbing methods or variations on Newton's method. A schematic diagram of a genetic algorithm is shown in FIG. 8. A good background article on genetic algorithms can be found in the July 1992 issue of Scientific American entitled "Genetic Algorithms" by John H. Holland, pp. 66-72.

Applied to the patch algorithm problem, the GA's genetic sequences are the patchlists, which describe which pixels are in a given patch. The GA ferrets out which pixels contribute to the solution and organizes them optimally. FIG. 9 shows a GA controlling rectangular patches. I have also written GA's not restricted to rectangular patches, and are currently evaluating the trade-offs.

In the training sessions I have conducted so far, the genetically optimized code can frequently achieve an order of magnitude improvement in speed and accuracy, and frequently comes up with novel and better solutions to complex problems than those anticipated by the inventor. The speed increase comes from the fact that the genetically optimized `patch` algorithm typically uses only a small fraction of the total pixels in the solution, greatly reducing the number of operations required.

FIGS. 10A, 10B, 10C, 10D, 10E and 10F illustrate a GA optimizing 24 patches on toluene, xylene, and methanol. The upper left panel, FIG. 10A, is a historical plot over 500 generations of the average fitness (lower curve) and superbug fitness (upper curve). A superbug is the most fit bug in a generation. The fitness is defined so that a perfect bug gets a 10. The middle upper panel, FIG. 10B, is a history of the patch size. The fact that the curve flattens out is a good indication the GA had found the optimal solution after about 200 generations. The right upper panel, FIG. 10C, shows the location of the 24 patches corresponding to an `average` bug. The patches are not really as distorted as the figure indicates. There is a severe aspect ratio distortion since there are 41 rows and 1024 columns of pixels. There are 100 bugs competing in each generation. The lower left panel, FIG. 10D, shows the superbug patches after 500 generations. Experience has taught that it is usually better to use the average solution rather than the super solution since the average is adapted to a wide range of conditions whereas the super solution is adept only on a special case. The lower middle panel, FIG. 10E, shows a snapshot of the training set: a linear combination of toluene, xylene, and methanol added together using random concentrations and 10% noise added as well. The right bottom panel, FIG. 10F, is a snapshot of the fitness of all 100 bugs taken at the end of the run. The actual raw fitness is indicated by the dots at the top of the panel. The asterisks are the weighted fitnesses. Weighted fitnesses are a tool used in GA's to force competition very early and very late in the run when all bugs are almost equally fit.

It is important to note that the GA is run only once for a given application. Once the patches are found that optimally solve for a given set of chemicals, that patch set can be used over and over. Also note that a variety of strategies can be used to train the GA. For example patches can be optimized to recognize certain chemicals while at the same ignoring others and so on. It is possible to train the GA using the real system response to optimize concentration estimates while rejecting systemic distortion, noise and contaminates. In other words, the actual operating parameters of the instruments used to take the measurements can themselves be optimized as well as the configuration of the data patches.

Restating some of the major points on the evolutionary algorithm, the algorithm begins as a statement of everything necessary to solve the identification/measurement problem. Each thing or "gene" within this statement or "genetic sequence" is assigned an allowable range of numbers. A generation of given number, say 100 or 1000, of genetic sequences or "bugs" is then compared against a data set for a known material or set of known materials, and the bugs with lowest residuals are assigned higher fitness numbers. Those bugs below a certain fitness level become extinct. Of the remaining bugs, those with the higher fitnesses are allowed to reproduce and do so by exchanging gene sequences. Also, mutations are introduced into each generation. Generation by generation the bugs evolve to a point where there are no longer improving to a significant degree. At this point the optimized solution sets are determined.

This evolutionary algorithm has certain characteristics which form a unique combination of capabilities otherwise unknown in the art. As opposed to the bit stream comparisons made in the Scientific American article, this algorithm assigns a humber with a variable allowed range for each measured characteristic. For the example above, three pieces of data are associated with each pixel--radiating wavelength, return wavelength and intensity of return wavelength. For the purpose of this discussion, this creates a one-dimensional mixture array, with one value associated with each pixel in the array. When the pixels are considered in patches, the array becomes two-dimensional. If one were to record more information, i.e. strength of radiating wavelength, polarization of radiating and return wavelengths, rate of attenuation of return wavelengths, etc., one will build up an n-dimensional matrix with regard to the radiation interaction with the database of known constituents. The patch algorithm is designed to handle these higher dimensions, and the evolutionary algorithm is designed to optimize it. As discussed immediately above, for each of these pieces of information or genes there is number with a assigned allowable range. This allowable range determines how much genetic mutation is allowed for that gene. This range can vary depending upon how far along the evolution of the bug has progressed. A broader range allowing greater mutations is often desirable at the beginning of the evolution while a narrower range towards the end of the evolution is often better.

Specifically, the evolutionary algorithm modifies the patches by evolving their boundaries and included pixels with their associated return signal magnitudes on the basis of a particular fitness parameter. The best patches have the highest fitnesses and are included as members of the patch database of known chemicals. The actual patch database or library of patches for known materials is made up of normalized values for the magnitudes of the measurements associated with each chemical or constituent. Determination of the actual concentrations of the known materials from the returns from the unknown mixture requires carefully calibrated information about the transmitted signal intensities, the system response characteristics and other factors to be able to translate the raw return information into concentration information. The procedure is complicated but routine.

This algorithm also provides for the optimization of more than one fitness criteria. For example, it may be desirable to have one patch set optimized for speed of determination with another optimized for accuracy. The algorithm will split the population into two populations which are then separately optimized for speed and accuracy respectively.

Chemical Recognition Using Neural Nets

As the chemical database grows one will need to intelligently choose which chemicals in the database will be actually used in the quantitative analysis. Some chemicals of interest are mission determined, but an intelligent neural net scanning the input could automate this selection, note unusual occurrences, and provide a useful double check on concentration estimates. Neural nets are attractive in the sense that they can be very fast to evaluate. A feed forward neural net has an easy to follow structure, has favorable scaling with number of chemicals (linear) and can be evaluated in a straightforward series of multiplies and adds. Neural nets do not require an iterative solve or a matrix inversion as do multivariate methods. The disadvantage of a neural net is that it is difficult to construct (train), and does not provide a rigorous uncertainty estimate. These deficiencies may someday be relaxed as research in neural nets continues at a rapid pace. A neural net designed for chemical recognition is potentially quite different from neural nets used in pattern recognition, which are basically classifiers. In the chemical recognition problem similar spectra are superimposed, so instead of asking `what letter is this?`, one is asking `how much of every letter?`, is superimposed on top of every other in the image.

I have developed a new genetically trained neural net for chemical recognition, as shown in FIGS. 11 and 12. The most common training technique of neural nets is back propagation. The novel use of a genetic algorithm to train the net has advantages in that a GA can also design the net. There may well be speed and accuracy advantages compared to the conventional back propagation training method.

FIG. 11 illustrates a typical neural network for this application. It has inputs of the return data of the response of a complex unknown mixture to the input radiation. Each of the pixels from the apparatus receiving the return information is connected to one of the neurons in the top layer with a single input connection with each such input being multiplied by a weight. Each input layer neuron has a separate connection to each of the neurons in the next layer down, again with each connection having a weight multiplier associated with it at the receiving neuron. This intermediate or hidden layer could be only one layer or could be a stack of n layers. This intermediate stage of layers (or just a single intermediate layer) is similarly connected to an output layer of neurons. Each neuron in this layer has only a single output which, as in the outputs of the neurons in the layers above is a function of the sum of all the individual inputs to a particular neuron multiplied by the weights associated with each of the individual inputs. A final layer, not shown, computes the concentrations or presence of various known chemicals.

In FIGS. 12A, 12B, 12C, 12D, 12E and 12F, a 4 layer neural net is training on three gaussian-shaped chemicals along the diagonal, similar to the test case displayed in the earlier FIGS. 6A, 6B, 6C, 6D and 6E. The upper left panel, FIG. 12A, is a snapshot of the training set used during generation 50. The training set is recomputed every generation, using a different set of concentrations for each of the three chemicals, chosen randomly between 0.0, and 1.0. The second panel, FIG. 12B, on the top row is a history of the raw fitness throughout the first fifty generations. The lower curae is the average fitness of all the 100 nets that are competing in each generation, while the top points are the fitness of the superbug at each generation. The gap between the superbugs and the average population is an indication that the general population has a lot of learning to do. The good performance of a superbug compared to the average on these early training sets indicates that it is too overspecialized for a particular special case. The next panel, FIG. 12C, is a snapshot of the raw and scaled fitness at generation 50. The right most panel, FIG. 12D, is a plot of the weights on all 400 neurons in the topmost layer. The weights are initially random, and are not showing much structure this early in the training session. The middle panel, FIG. 12E, plots the 400 weights of each of the 12 neurons in the second layer for a total of 4800 weights. The bottom panel, FIG. 12F, plots the 12 weights of each of the 12 neurons in the bottom layer for a total of 144 weights. The bottom 4th layer is not plotted; its output is the net's concentration estimates and is used to compute the fitness.

FIGS. 13A, 13B, 13C, 13D, 13E and 13F show the same net after 3000 generations. The fitness of both the superbugs and the overall population are achieving fitnesses over 8, which is very good. The population is still learning at this point, and its ultimate capability can improve still further. Notice the plots of the weights now exhibit considerable structure; the information in certain pixels being amplified and differentiated by the net, while other pixels are multiplied by zero and are effectively discarded.

The evolutionary algorithm as applied to the neural net is also somewhat of a departure from the standard practice. This evolutionary algorithm has the capability to alter the number of neurons in a layer and also the interconnects between layers by "pruning" out those interconnection with zero weight multipliers. It also has the capability to mutate the neural net by adding in neurons and interconnections with initially random weights. While the neural network is initially assigned a structure with a given number of layers and neurons, the evolutionary algorithm takes control over its structure once it starts running. The net is trained by supplying a different training set (the same combination of known samples but at different concentrations) to the net at each generation and evolving the weights and interconnections of the net on the basis of the fitness of the various bugs in each generation in much the same fashion as for the patch technique discussed above.

The genetically trained multivariate analysis patch technique can be combined with the genetically trained neural net to form a powerful expert system organized as in FIG. 14. The data received from the sensors would be fed initially into the neural net which will quickly identify the most important constituents in the mixture and estimate their concentrations. The net would hand off the candidate list from the chemical database to the multivariate processor which then utilizes this list in addition to any user-specified chemicals for analysis and uncertainty estimates. The pre-screen by the neural net cuts down the workload for the multivariate analysis data processing significantly. Also a comparison of the concentration estimates from the neural net with the output of the multivariate processing is a very useful consistency check.

The technique of the first preferred embodiment has also been applied to the problem of recognition of microorganisms in a biologic sample. As with the example above, uv laser induced fluorescence spectra were analyzed for E. coli bacteria in a bovine serum albumen mixture. The system was able to recognize the E. coli bacteria spectral signature and estimate the concentration. This has immediate applications for detection of food contamination and also for diagnostic use by medical personnel for the presence of pathogens in or on a patient. A suitably packaged system could be used to diagnose a streptococcus infection in a patient's throat. It could also be used to search for biologic contamination in hospital operating rooms and on sterile equipment. Indeed a suitably trained database could recognize subtle genetic variations between samples of DNA or could provide an early detection method for cancerous or otherwise abnormal cells on the basis of any characteristic of the subject cells that can be distinguished on the basis of the interaction of the radiation with those cells.

Other uses of the techniques disclosed herein will find use in a variety of areas. Petrochemical industries have a particular need to continuously monitor the composition of materials in pipelines and feedstocks. As agriculture refines its usage of fertilizers and pesticides further, the invention could be used as a real time monitor of the fertilizer/pesticide content of a field in the pathway of an application vehicle and to direct a system on the vehicle to apply the fertilizer/pesticide onto only those portions of the field that need it.

Returning to the hardware and the physics of the invention, the first preferred embodiment utilized uv fluorescence spectroscopy to analyze the mixtures. This inelastic process is useful only with relatively complicated molecules. As was discussed above, this technique will work with any radiative interaction with a material that provides a return that somehow serves to distinguish the material of interest. Both elastic and inelastic interactions can be measured and used herein. Elastic interactions occur with all materials when irradiated with electromagnetic waves. However, these effects are often more subtle than the inelastic measurements and are sometimes nothing more distinct than a small shift in the lifetime of the return pulse at various return wavelengths as the difference between two materials. More robust effects are absorption, changes in polarization, and resonance Raman spectroscopy. For the uv spectroscopy, two different types of sources can be used. A flash lamp with various bandpass filters can be used as a low resolution source. A tunable laser operating between 200 and 400 nm is a good high resolution source. It is not necessary to use the very same source for the creation of the chemical database and for the actual measurement of the unknown mixture so long as the electromagnetic waves are provided at the same frequencies. As was mentioned above, very careful attention must be paid to calibrating all of the data generation equipment. "The hardware used in the system is shown schematically in FIG. 16. An illumination/excitation source 100 is used to irradiate the sample 102. The energy coming out of the sample is collected by a detector 104. The source and detector are controlled by and the data are collected by the data collection and controller 106. The collected data are sent to the processor 110. If the known samples are being run for the first time, the catalog of their optimized responses are stored in the catalog of known samples 108 as a product of the processors operations on the data. However, if this step has already been done, the catalog of known samples is referred to by the processor when determinations about an unknown sample are being made.

The remainder of this discussion is devoted to the mathematics underlying the invention.

1.1 `Traditional` 1-D Chemometrics

In this section we develope a formalism for a algorithm to estimate chemical concentrations contained in a complex mixture. We assume that each chemical sample has a unique 1-D spectral signature, and we have a set of these spectral signatures saved in a reference database set. We assume the mixture, or unknown is made up of a linear combination of the chemicals in the database. The purpose of this section is to solve the problem in a conventional way for comparison to the more sophisticated methods in 1-D, 2-D, and 3-D developed in subsequent sections. For now we assume the mixture, or unknown, is a 1-D linear array of values which represent a measurement.

FIG. 15 represents how the magnitudes associated with the pixels for each sample in the database multiplied by their concentrations add up to the measurements from what is shown here as a mixture of the three samples with no noise.

In our application the measurement will be a UV fluorescence intensity taken at a variety of emission wavelengths, but the general approach described here is valid for any measurement curve that is a mixture of other curves.

For convenience we will frequently refer to the elements of the mixture and sample arrays as pixels. Each pixel has a value or intensity. The 1-D mixture array M with n pixels m.sub.1 m.sub.2 m.sub.3 . . . m.sub.n will be represented by a matrix with one row,

Each chemical sample S.sub.l with n pixels s.sub.1l s.sub.2l s.sub.3l . . . s.sub.nl will also be represented by a matrix with one row, as in

The second subscript l (l=0, 1, 2, . . . , N.sub.s) is the chemical identifier. There are a total of N.sub.s chemicals in the database.

Since ›M! is assumed to be made up of a linear combinations of the samples, ›M! can be written schematically as:

The residuals can be considered a combination of measurement uncertainties such as calibration uncertainties, signal and instrument effects, noise, filtering, and the presence of chemicals in the mixture which are not in the database. Writing out equation 1.1-3 fully results in the equation,

Now to solve for the concentrations C.sub.1, C.sub.2, C.sub.3 . . . C.sub.l . . . C.sub.N, we will use a least squares minimization approach. We define ##EQU1## where .delta..sub.t.sup.2 is the weighted squared residual for pixel i. The weight .sigma..sub.i can be interpreted in many ways. For our applications the weights .sup..sigma. i will frequently be the uncertain estimates for mixture measures of each pixel, i.e. ##EQU2## n.sub.l =# measurements of m.sub.l. Another interpretation of .sup..sigma. i has to do with the `value` of that pixel. We will explore the use of this interpretation in follow on reports dealing with the genetic algorithm optimizations. In many cases for simpler problems, and for interesting limiting cases, the weights can be simply set to one and ignored, with excellent results.

Now to solve for the concentrations, we minimize the sum of the squared weighted residuals, ##EQU3## It is worth noting here the choice of this functions to minimize is arbitrary, and many other choices exist. We will explore this more later in this report. In the literature minimizing the sum of the squares is the classical least squares method, while minimizing the sum of the absolute values results in a `robust estimation` method that is less sensitive to out lying points, but the robust method is frequently solvable only by iterative means. There are many other interesting modifications to the minimization function "Numerical Recipes: The Art of Scientific Computing", W. H. Press et al., Cambridge University Press, 1986.

To minimize .sub.x.sup.2, equation 1.1-6 is differentiated with respect to each of the concentrations; ##EQU4## Equations 1.1-7a to 1.1-7e can be rewritten as (with the factor 2 divided out), ##EQU5## Now equation 1.1-8a to 1.1-8e can be rewritten to facilitate a matrix solution, ##EQU6## Now equations 1.1-9a to 1.1-9e become the matrix equation, ##EQU7## Now if we define ##EQU8##

Finally equation 1.1-10 can be written in the very simple form,

Or explicitly, ##EQU9## Or even using the `repeated subscript implies summation` convention,

We are not quite done yet. The concentrations must be solved for. In this exsection the matrix .alpha..sub.mn is square with size N.sub.s by N.sub.s. Given the typical restrictions, such as a nonzero determinant an inverse to ›.alpha.!, namely ›.alpha.!.sup.-1 exits with the property ›.alpha.! ›.alpha.!.sup.-1 =1 ›.alpha.!.sup.-1 can be found in several ways. To maintain consistency with later sections of this report, SVD (Singular Value Decomposition) methods are used as follows. ›.alpha.! can be written as the product of three special matrices, ##EQU10## Where ›U! and ›V! have inverses that are their transposes. SVD algorithms are readily available see Press, et al., supra, and they provide the right hand side of 1.1-17 given ›.alpha.!. The advantages of using SVD methods are very important and have broad relevance to this problem which go beyond just providing the inverse. This will be discussed more later. Hence we get, ##EQU11## So that ##EQU12## which is of the form,

This is the solution. The concentrations are found. Note that in many real applications, only ›.beta.! changes with new mixtures. Therefore the inverse of .alpha. needs to be computed only once for a given set of sample in the edatabase. This makes repetitive estimates of concentrations extremely efficient, involving only a matrix multiply.

The question of uncertainty in concentration estimation will be deferred until later in this report. A good discussion can be found in Numerical Recipes see Press, et al., supra. Basically, assuming a normal distribution of errors on the mixture the uncertainty in concentration estimates is given by the SVD inverse matrix which can be shown to be the covariance matrix. In cases were there are uncertainties in the database the situation is much more complicated, and will be discussed later.

1.2 Basic `Patch` Algorithm in 1-D

In the traditional approach outlined above, all pixels were used in an order-insensitive way. The contribution from each pixel in the sums can be added in any order without changing the vale of the sum. this means that the pixels could even be scrambled randomly as long as the mixture and database sample are scrambled identically, without changing the answer. To address this we have developed a new algorithm, described below

1.2.1 Minimization of residuals for a single patch

The sum of the squared differences .delta..sub.L.sup.2 between pixels located in common patches in both the mixture and its constituents can be expressed as: ##EQU13##

To minimize .delta..sub.L.sup.2, differentiate with respect to each of the concentrations and get: ##EQU14##

Simplifying, and using the simpler notation ##EQU15##

Which can be written as a matrix equation: ##EQU16##

1.2.2 Multiple patches in 1-D

Multiple patches are handled by stacking the equations for the patches on top of each other vertically: Simplifying, and using the simpler notation ##EQU17##

In the previous section the minimization of .delta..sub.L.sup.2 for a single patch was demonstrated. In our case we our more interested in the general case of many patches. We want to minimize the sum of all .delta..sub.L.sup.2 as given by: ##EQU18##

Notice there is no requirement for patches to be sized similarly or located in any special place. The sums are for all elements in a patch, which can be any size or shape. ##EQU19## 1.3 Basic `Patch` Algorithm in 2-D

In this section the mathematics for a 2-D image-recognition algorithm is described. The problem to be solved is: given a 2-D `unknown`, what are the concentrations of it's 2-D constituents? The unknown mixture is assumed to be in the form a 2-D array and is assumed to be made up of a linear superposition of many 2-D elemental constituent images. Such 2-D arrays can be from any source, from a set of pixels describing objects in a picture, to the application of primarily interested: 2-D multispectral UV LIDAR images. In the LIDAR application multispectral 2-D images of the UV florescence of many chemicals are measured and stored individually in a database. When the algorithm is given an unknown, it estimates the concentrations of the constituents in the database which are present in the unknown.

The arrays need not be square, but for simplicity we assume the image is rectangular. Even though this requirement is not really necessary, it is expedient to limit 2-D images to be rectangular when using some computer languages such as FORTRAN. The unknown mixture will be denoted by ›M! with dimensions m.times.n: ##EQU20##

Subscripts are used to designate array elements. It is convenient to think of each m.sub.1l as an intensity at pixel location (i,j).

Now suppose the mixture ›M! is made up of the linear superposition of N.sub.s constituents. ›M! may also contain constituents not in the data base, plus some uncertainty or noise. Each constituent, or sample ›S.sub.l !, (where l.ltoreq.N.sub.s) is also a 2-D m.times.n array. Superscripts are used to identify different samples, ##EQU21##

Since ›M! is assumed to be made up of a linear combinations of the samnples, ›M! can be written as:

The problem is to estimate the concentrations C.sub.l,k=1,N.sub.s ##EQU22##

One method to solve this problem is to minimize the sum of the squared differences between pixels in the unknown and a hypothetical set of sample concentrations. This least-squares method is usually well deternuned and straight forward but does not take advantage of the 2-D feature of the problem. In a pixel-by-pixel minimization, ordering is not important and information can be under utilized. The pixel-by-pixel minimization is a special case of the patch algorithm described next. To take advantage of the 2-D nature of the problem it is possible to minimize the squared error between a group of near-neighbor pixels in bothi the unknown mixture and its hypothetical constituent makeup. A set of near-neighbor pixels in a group is defined here to be a `patch` of the array.

Different patching strategies are possible, as illustrated in FIG. 5. Bad data can be ignored while emphasizing regions of strong pixel-to-pixel correlation.

1.3.1 Minimization of residuals for a single patch

The sum of the squared difference .delta..sub.p.sup.2 between pixels located in common patches in both the mixture and its constituents can be expressed as: ##EQU23##

To minimize .delta..sub.P.sup.2, differentiate with respect to each of the concentrations and get: ##EQU24##

Simplifying, and using the simpler notation ##EQU25##

Which can be written as a matrix equation: ##EQU26##

1.3.2 Multiple patches in 2-D

Multiple patches are handled by stacking the equations for the patches on top of each other vertically: Simplifying, and using the simpler notation ##EQU27## In the previous section the minimization of .delta..sub.P.sup.2 for a single patch was demonstrated. In our case we our more interested in the general case of many patches. We want to minimize the sum of all .delta..sub.P.sup.2 as given by: ##EQU28##

Notice there is no requirement for patches to be sized similarly or located in any special place. The sums are for all elements in a patch, which can be any size or shape. ##EQU29## 1.4 Basic `Patch` Algorithm in 3-D

1.4.1 Minimization of residuals for a single patch

The sum of the squared differences .delta..sub.P.sup.2 between pixels located in common patches in both the mixture and its constituents can be expressed as: ##EQU30##

To minimize .delta..sub.P.sup.2, differentiate with respect to each of the concentrations and get: ##EQU31##

Simplifying, and using the simpler notation ##EQU32##

Which can be written as a matrix equation: ##EQU33##

1.4.2 Multiple patches in 3-D

Multiple patches are handled by stacking the equations for the patches on top of each other vertically: Simplifying, and using the simpler notation ##EQU34## In the previous section the minimization of .delta..sub.V.sup.2 for a single patch was demonstrated. In our case we our more interested in the general case of many patches. We want to minimize the sum of all .delta..sub.V.sup.2 as given by: ##EQU35##

Notice there is no requirement for patches to be sized similarly or located in any special place. The sums are for all elements in a patch, which can be any size or shape. ##EQU36## 1.5 Error Analysis

1.5.1 `Traditional` 1-D Chemometrics

Since ›M! is assumed to be made up of a linear combinations of the samples, ›M! can be written as: ##EQU37##

To minimize .sub.x.sup.2, differentiate with respect to each of the concentrations and get: ##EQU38## which is of the form,

If we define ##EQU39##

Then using SVD, . . . ##EQU40##

So that ##EQU41## which is of the form,

Now what is the uncertainty of C? three terms ##EQU42##

Error analysis derivation based on the assumption that ##EQU43## the first term is commonly found and the others are commonly ignored. This is because it is usually assumed that the fitting functions are exact, or in the case of curve fitting the locations x at which measurements y are taken are known exactly. This leads to the curious result if the uncertainty of all measurements are equal then the variances of the fitted parameters (the intercept and slope in line fitting) depend only the x's and not the y's. In our case this is not true. both the mixture and the database samples have measurement uncertainties associated with them.

Definitions . . . ##EQU44##

To simplify the derivation, consider the basic equation (repeated subscripts indicate summation) ##EQU45## in the non-square patch case it will look like ##EQU46## for patch a and b . . . Note the matrix on the left does not need to be square, . . . it will have a pseudo-inverse as long as the problem is well posed. Then the concentrations can be solved for, by computing ›.alpha.!, the pseudo inverse of ›a!. Here we define a pseudo inverse to be the non-square matrix ›.alpha.! which has the property ##EQU47## (a square matrix with zero's everywhere except 1's on the diagonals). ##EQU48## For the purpose of this derivation we assume the pseudo inverse can computed at will. ln practice we will use singular value decomposition, which not only provides the inverse, but does so in a form that is particularly useful and remarkably free from numerical problems.

Now to compute the derivative terms . . . ##EQU49##

Note the rule (a and .alpha. are pseudo inverses of each other) ##EQU50## with the result ##EQU51## this enables us to take the derivative functionally ##EQU52## recognizing that the product of the two central matrixes is just the concentration matrix, and bringing the right hand term into the product, ##EQU53##

Now generalize a little bit, so we don't have to repaint three equations over and over, explicitly writing out the sum on j for now. The kronecker delta is used. ##EQU54##

Consider the product, for l=1 ##EQU55## and for l=2 ##EQU56## so that in general ##EQU57## and therefore ##EQU58## now note the presence of the residual r ##EQU59## or for a particular concentration ##EQU60## now to compute the squared value, ##EQU61## in summation notation . . . ##EQU62## summing over all pixels . . . ##EQU63##

The first term is rather remarkable, it is just 1/c.sup.2 times the classical result, for the first time we see the variance getting bigger as the concentration gets smaller. Also note the scaling for as the residual grows the variance grows rapidly. This will be a great indication that chemical are missing from the database.

this equation can now be plugged into the first equation in this chapter.

For completeness, . . .

summation implied . . . ##EQU64##

Now to compute the `classical result` . . . the Rate of change of concentration estimate with respect to mixture pixel variance.

Using a double subscript i to indicate a sum over pixels different from p . . . ##EQU65## using the summation convention on i and j and dropping the subscript p for now... ##EQU66## using the summation convention . . . ##EQU67## this equation can now be plugged into the first equation in this chapter, due to ##EQU68## this is the classical result. The variance associated with this term depends on the s's, the samples in the database. You could give the algorithm any mixture and this term doesn't change.

Covariant or terms (for the same pixel p therefore dropping subscript p for now) . . . ##EQU69## this equation can now be plugged into the first equation in this chapter, notice that in many cases these terms will be negative, this is good.