Real Time Magnetic Resonance Analysis With Non-Linear Regression Means

This application is closely related to the U.S. Pat. No. 5,015,954 issued on 14 May 1991 to Dechene et al., and to U.S. Pat. No. 5,049,819 issued Sep. 17, 1991 to Dechene et al., both entitled "Magnetic Resonance Analysis in Real Time, Industrial Usage Mode"; and to U.S. patent application Ser. No. 07/794,931 filed Nov. 20, 1991 entitled "Improved Magnetic Resonance Analysis in Real Time, Industrial Usage Mode"; and to U.S. patent application Ser. No. 07/885,633 filed May 19, 1992 entitled "NMR Analysis of Polypropylene in Real Time". All of these patents and patent applications are of common assignment with this application, and the disclosures of all are hereby incorporated herein by reference, as though set out at length herein.

The present invention relates to the analytic means relating a Nuclear Magnetic Resonance (NMR) free induction decay (FID) curve to the physical quantities of the target nuclei in the samples under test. More particularly, the present invention relates to regression means that provides correlation between the function equations derived from the FID and the types, properties and quantities of target nuclei in the sample under test.

Pulsed NMR techniques are used in instruments for the measurement of the type, property and quantity of lattice bound and free, magnetically active nuclei within a sample. Some of the substances and properties that have been measured by NMR techniques are: moisture, polymers and copolymers, oils, fats and crystalline materials.

Pulsed NMR uses a burst or pulse of energy that is designed to excite the nuclei of a particular nuclear species of a sample being measured (the protons, or the like, of such sample having first been precessed in an essentially static magnetic field); in other words the precession is modified by the pulse. After the application of the pulse there occurs a free induction decay (FID) of the magnetization associated with the excited nuclei. That is, the transverse magnetization associated with the excited nuclei relaxes back to its equilibrium value of zero. This relaxation produces a changing magnetic field which is measured in adjacent pickup coils. A representation of this relaxation is the FID curve.

The analysis method described herein and in the above related patents and applications is to decompose the FID waveform into a group of separate time function equations. The coefficients of these equations are derived from the FID by use of a Marquardt-Levenberg (M-L) iterative approximation that minimizes the Chi-squared function--a technique well known in the art. Some of the time function equations found useful are: Gaussians, exponentials, Abragams, and trigonometric. From these time functions a set of parameters is calculated. Some of these parameters are ratios of the y-axis intercepts, squares and cross products of these ratios, and decay times for each of the time curves. In addition the sample temperature may form the basis for another parameter.

But, relating these parameters, quantitatively and qualitatively, back to the species of target nuclei is required. In the above referenced patent applications, the system is calibrated with known samples, and a `regression line` is generated which relates the parameters to the types, properties and quantities of the target nuclei. An unknown sample is introduced and the time functions are derived via the M-L iteration, and the parameters are calculated. The parameters are "regressed" via the "regression line" to yield the types, properties and quantities of target nuclei in the unknown sample. That is, the measured parameters from the unknown FID are used with the "regression line", and the types, properties and quantities in the unknown sample are determined. It is to be understood that the multidimensional "regression line" may not be graphically represented. As a simple regression technique example, consider that the grade point average of each of the students at a college were related to that student's SAT score and high school standing (forming a three dimensional space). The line formed is a "regression line" (which may be graphed). A new student's grade point average may be predicted by placing the student's SAT and high school standing on the "regression line,, and "reading" the grade point average.

It is a principal object of the present invention to relate the type, property and quantity of target nuclei of interest accurately and precisely.

The above object is met in an NMR system that effects a reliable extraction of free induction decay data in a way that is practical in a wide variety of applications, including industrial and medical. The NMR system is calibrated by measuring known samples of target nuclei and, from the FIDs generated, forming a multi-dimensional, non-linear regression relationship to the types, properties and quantities of target nuclei. The FIDs are decomposed or transformed into a set of equations for the calibration samples from which a set of parameters is generated. From these parameters a non-linear regression function is calculated relating the type, property and quantity of target nuclei to the parameters. An unknown sample FID is decomposed or transformed as were the known samples, the parameters are calculated and these parameters are used with the non-linear regression function to determine the type, property and quantity of target nuclei in the unknown sample. In a preferred embodiment, the FID is decomposed into multiple time equations via M-L processes and parameters are calculated for each of these time equations. In another preferred embodiment the parameters are non-dimensional in order to eliminate concentrations and the like from the measurements.

The present invention may be used to advantage with any number or type of time or frequency functions derived from an FID waveform, including Fourier transform functions.

Other objects, features, and advantages will be apparent from the following detailed description of preferred embodiments taken in conjunction with the accompanying drawing(s) in which:

FIGS. 1-2 is a block/schematic drawing of a pulsed NMR system suitable for measuring a range of industrial materials,

FIG. 3 is a graphical representation of an FID and its component curves,

FIG. 4 is a flow chart of a preferred embodiment of the present invention, and

FIG. 5 is a flow chart of the steps to establish an effective industrial measurement.

FIG. 1 shows transverse and cross sections with block diagram inserts of an NMR apparatus and method where the present invention may be used to advantage. An industrial process line IPL has material flowing as indicated by arrow A. Some of the material is captured by a probe P and fed through an inlet line LI to a sample region S1. The region is defined by a tube 98 typically about 30 cm long made of an essentially non-magnetic, nonconducting material which does not itself generate substantially interfering FID signals (glass, certain ceramics, certain plastics or hybrids may be used). The sample region is defined between inlet and outlet valves V1 and V2. Gas jets J are also provided. These are pulsed on/off repeatedly to agitate fluent sample materials during sample admission and expulsion. The region S2 is the critical portion of the sample. It is surrounded by a sample coil 100 tuned to resonance and driven by a tuning circuit 102 and related transmitter/receiver controller 104. Grounded loops 101 are Lenz Law shields which are provided above and below coil 100 to help shape the field of coil 100--i.e., contain the field established by an excitation pulse. The controller 104 includes an on-board microprocessor and required power supply elements, memory, program and I/O decoding suitable to interconnect to the hardware shown and to an external microcomputer 106 with keyboard 108, monitor (or other display) 110, recorder 112 and/or process controller 114 (to control the process at IPL). The operator initiates and controls operation from the display keyboard 108 and the resulting data and signals are subsequently shown on the display 110 and utilized in 112 and/or 114. The computer 106 also controls instrument operation conditions.

The region S2 of tube 98 and coil 100 are in a static, but adjustable, crossing magnetic field defined by a magnetic assembly 116 which comprises a yoke 118, pole pieces 120, surrounding Helmholtz coils 124, and a coil current generator 117. The critical sample region S2 of the tube 98 and magnet are contained in a metallic (but non-ferromagnetic) box 126 with highly thermally conductive face-plates 128 and internal partitions 130 and over-all mass related to each other to minimize harmonics and other interferences with a signal emitted from coil 100 to a sample and/or returned from the sample for pick-up by coil 100 and its tuned circuit 102 and transmit/receive controller 104.

The magnetic assembly 116 including yoke 118, and other parts therein as shown on FIGS. 1-2, is in turn contained in an environmental control chamber 132 with optional inert gas fill and purge controls (not shown), an internal gas heater 134, a motor M driving fan 136, and a temperature sensor 138 which can be applied to the yoke or other detection region whose temperature is reflective of the temperature at pole pieces 120 and in the sample region therebetween. A thermal controller 140 processes temperature signals from 138 to adjust heating/circulation at 134/136 as a coarse control and to adjust current through the Helmholtz coils 124 at magnet pole pieces 120 as a sensitive and fast fine control, as well as implementing general control instructions of computer 106. Further thermal stabilization may be provided by a closed loop heat exchanger 142 having pump 144 and coils 146 attached to yoke 118 and coils 148 attached to the plates 128 of box 126.

The strength, consistency and constancy of the magnetic field between poles 120 in the region S2 of the sample is thus controlled by a uniform base magnetic field in the entire region S2. The Helmholtz coils 124 are energized by the coil current controller 117 to accurately trim the final magnitude of the field in which the sample is placed. This field is the vector addition of the fields due to the magnet poles 120 and the Helmholtz coils 124. The controller 117 sets the current through the Helmholtz coils 124 using current generators. The coils 124 are wound around the magnet pole pieces such that the magnetic field created by the current in the coils 124 can add to or subtract from the field created by the magnet pole pieces. The magnitude of the current through the coils 124 determines the strength of the field added to or subtracted from the field due to the magnet pole pieces (and related yoke structure) alone.

The actual determination of the current through the Helmholtz coils is accomplished by carrying out the magnetic energy and resonance techniques hereinafter described in preliminary runs and adjusting Helmholtz current until the maximum sensitive resonance is achieved, and then setting the Helmholtz current off resonance by a given offset, of about 0.1-3 KHz.

The major elements of electrical controls are tuner 102, including coils 100 and 101 and variable capacitors 102-1 and 102-2, resistor 102-3 and diodes 102-4 and constructed for tuning to Q of twenty to sixty to achieve coil 100 resonance, and control 104 including a transmit/receive switch 104-1, a transmitter 104-2 and receiver 104-3, a crystal oscillator 104-4, gated pulse generator (PPG) 104-5, and phase shifter 104-6. The crystal provides a nominal twenty Megahertz carrier which is phase modulated or demodulated by the MOD, DEMOD elements of transmitter 104-2 and receiver 104-3. The receiver includes variable gain amplifier elements 104-31 and 104-32 for operation. The analog signals received are fed to a high speed at least 12 bit flash A/D converter 105-1 and internal (to the instrument) CPU element 105-2, which provides data to an external computer 106 which has a keyboard 108, monitor 110, modem 109, recording elements 112 and process controller elements 114, e.g., for control of valves V1, V2 via valve controls 115 and/or to coil current controls 117, all via digital-analog converters (not shown).

The analog signal FID curve is conditioned by a Bessel filter which acts as a prefilter and an anti-aliasing filter as the subsequent sampling is done at 10 MHz. After digitization the signal may be time smoothed by a fast Fourier transform filter program. The combination of these filters produces a relative improvement in signal to noise ratios which enhances the accuracy of the system.

The excitation of coil 100 and excitation-precession of the sample's proton content and subsequent relaxation/decay produces a received FM signal that, after demodulation, controlled gain amplification, and A/D conversion produces the free induction decay (FID) curve.

Referring to FIG. 3, the digitized FID curve data are stored in the external computer 106 where a program finds the best component curves to fit each stored FID curve. In this preferred embodiment there are three component curves, a fast Gaussian, a slow modified Gaussian and an exponential. Other preferred embodiments have more or less than three component curves and other curve types. The determination of the types of curves which make up the FID curve is important because, once the curves are known, they can be extended back to a time origin (shown as A.sub.O, B.sub.O and E.sub.O at tO, i.e., excitation of a Cycle 1), which is close to the center of the transmitted burst signal. This is important since there are saturation effects of the instrument's electronic gear which occur during and immediately after the excitation burst signal. During this time, measurements cannot be accurately taken, yet the area of interest under the curve, which is a measure of the number of nuclei in the sample, extends from the immediate end of the excitation burst to where the curve is too small to be digitized or it is in the noise.

The entire curve is decomposed into component curves and these curves are fitted to the data by an iterative process based upon the Marquardt-Levenberg (M-L) approximation technique applied automatically through a structured realization in software. This technique is used to determine the magnitude of all the parameters, constants, frequencies, etc. which best fit the FID curve. This is an iterative technique where the entire curve is determined at once. The M-L iteration process performs the curve fitting by attempting to minimize the Chi-Squared error function (the sum of the squared differences between the measured data points and the data points from the derived equation). The results of the M-L approximation are accepted if the Chi Squared error is small enough, if not, the M-L fitting procedure may be reapplied with a different set of starting guesses. If this process also fails, the sample is discarded and a new sample obtained. The M-L technique is documented in the following references: Ind. Appl. Math., vol. 11, pp. 431-441 by D. W. Marquardt, 1963; Data Reduction and Error Analysis for the Physical Sciences (New York, McGraw Hill), Chapter 11 by Philip R. Bevington, 1969; and The State of the Art in Numerical Analysis (London: Academic Press, David A. H. Jacobs, ed 1977), chapter III.2 by J. E. Dennis. As applied to the measurement regime of interest herein, in a preferred embodiment of the present invention, the selected parameters taken from the derived curves are the y-axis intercept ratios, time constants, frequency terms and other parameters described below.

Other known in the art iterative techniques which may be applied instead of or with the Marquardt-Levenberg, include: Gauss-Newton and "steepest descent" (found in the above J. E. Dennis reference), Newton-Raphson (known in the art), or like techniques, including combinations of these techniques.

One of the major difficulties in making use of iterative curve fitting techniques (such as Marquardt-Levenberg) is their tendency to reach incorrect solutions. Such solutions frequently (but not always) contain parameters which would imply a negative quantity of protons or an exponential "decay" which grows with time. These incorrect solutions lead to serious errors in the result found for a physical sample, for example, the density in polyethylene or the quantity of xylene solubles in polypropylene.

The usual methods of handling these difficulties have been:

(1) have a human evaluate the result and eliminate those solutions that are ridiculous, and/or

(2) put a series of upper and lower bounds on each parameter beyond which the fitting procedure is forbidden to go.

In an on-line situation where readings are generated every few minutes, the first approach obviously cannot be used, and in the case of polyolefins the second approach fails because the bounds for each parameter depend on the actual values of the other parameters (recall that for polypropylene and polyethylene the model equations involve ten or more parameters).

We have evolved a Marquardt Reference Ratio (MRR) to handle this difficulty. MRR is a ratio although other techniques (differences, for example) could be used.

As discussed herein, the techniques to find a property of an unknown sample include calibration by applying the M-L technique to reach solutions for a group of FIDs of samples with known properties. The various amplitudes and time constants in the solutions are combined to produce a number of ratios, cross products and higher order parameters. These parameters may undergo various non-linear transformations and are finally regressed multi-dimensionally to obtain the coefficients of the regression equation to use in predicting a property of an unknown sample, say, for example, density. Each of the parameters contributes to the overall prediction of density. Moreover, in the nature of things, these parameters tend to be relatively highly correlated among themselves, e.g., a large crystalline content must necessarily correspond to a small amorphous content (comparing the modified Gaussian to the exponential in the polyethylene FID solution). This means that overlapping density information is contained in many of the parameters used in the regression equation. Similar arguments apply to other properties, such as xylene solubles in polypropylene.

To make use of this high correlation (continuing the density example), the parameters are divided into subgroups (two roughly equal groups in a preferred embodiment) and each of these groups is regressed on density to obtain two further predictions of density based on each subgroup, as follows:

D1 (density)=F(subgroup 1)

D2 (density)=G(subgroup 2)

Because of the correlation, discussed above, among the parameters, the functions F and G (above) result in predictions D1 and D2 which are only slightly less accurate than the density prediction based on the entire set of variables. The ratio (MRR) or the difference (MDR) are formed as follows:

MRR=D1/D2=F/G

MDR=D1-D2=F-G

MRR has a nominal value of one and MDR zero. MRR and MDR are sensitive measures of whether or not a particular proposed M-L solution for an unknown sample belongs to the set of (calibrated) data from which the functions F and G were derived. If the calculated ratio or difference of D1 and D2 for a proposed M-L solution for fitting the FID of an unknown sample lies outside reasonably well-defined limits (usually +/-3 sigma), the proposed M-L solution may be assumed to be bad and is discarded.

Once the equation of the FID curve is known, each component curve can be extrapolated back to the mid point of the excitation signal to establish the intercept of each said component curve.

The resulting data utilized in the computer 106 (FIGS. 1-2) is the equation for the FID curve as composed of a number of component curves. Each of these curves (and their intercepts) has been experimentally and theoretically related to particular nuclei of interest. In particular, when the FID curve equation is determined, the ratios of the y-axis intercepts, the cross product and squares of these ratios and the decay times for each of the curve components, the product temperature and a cosine term form a multidimensional model.

Calibration of the system is accomplished by measuring a number of known samples and using the M-L technique to derive the model equation constants associated with each known sample. Various non-linear transforms may then be applied to these constants, usually with the goal of linearizing their relationship to the dependent (i.e., predicted) parameter. Useful non-linear functions include exponential, logarithmic, powers and cross products of the independent (i.e., measured) parameters.

FIG. 4 is a flow chart of the steps used in a preferred embodiment of the present invention. The first step is to measure samples with known types and quantities of target nuclei. The FID curve is digitized via a flash converter of at least 12 ands preferably 14 bits accuracy and stored in computer memory. The next step is to apply the M-L iterative process to derive curve coefficients from the stored FIDs to a given Chi Squared error. In step three the ratios of Y-axis intercepts, squares and cross products of these ratios, decay times and temperatures are calculated. In the next step, the various non-linear transformations to be used are determined, and the types, properties and quantities of target nuclei in the known samples are related to the constants by a regression against these transformed parameters--the "regression function." Step five is to record, digitize and store the FID for an unknown sample and derive the curve coefficients. The parameters are calculated for the unknown, and these parameters with desired non-linear transforms are used in the regression equation to determine the actual type, property and quantity of target nuclei in the unknown sample. Ratios are used since constants with dimensions of weight would require the samples to be carefully weighed before measuring.

FIG. 5 is a flow chart showing the steps of measurement to establish effective industrial measurement. A single FID curve is established to see if the sample area is clear (Quick FID) in an abbreviated cycle of attempting to establish an FID curve. If the sample region is not clear (N), measurement is interrupted to allow valve V2 to open and jets J and gravity to clear the region. A new Quick FID step establishes clearance. Then another sample is admitted by closing valve V2, opening valve V1 and making such adjustments of probe P and line L1 as may be necessary to assure sample acquisition. Jets J adjust and stabilize the new sample.

Temperature controls 134-138 and 142-146, described above, may be used to establish very coarse and less coarse thermal controls countering sample and ambient temperature variations.

An electronic signal processing apparatus baseline is established in 3-4 cycles (each having (+) and (-) subcycles with addition of (C+) and (C-) to detect a baseline offset and compensate for it). It would be feasible to avoid this baseline step determination and simply deal with it as an additional parameter (i.e. eleventh dimension in the M-L analysis, but this would increase iteration time).

Further adjustment is established by coils 124 to adjust HO (i.e., resonance) and this is enabled by ten to twenty field check cycles of FID curve generation. The (C-) FID is subtracted from the (C+) FID, (this process eliminates small baseline offsets) to obtain a workable digitized FID signal--which has a maximum value at resonance. HO is adjusted via coil current generator 117 and coils 124 until such maximum is achieved, and then HO is changed to offset the system a known amount from resonance. Adequate field adjustment is usually made in less than seven cycles.

Then five to one hundred cycles are conducted to obtain a useable measurement. Each of those five to one hundred cycles involves a modulated transmission/reception/flash A-D conversion, and storage of data. The curves are then averaged for M-L curve fitting, and the above listed intercepts and ratios are established. Similar cycles, but somewhat abbreviated can be applied for Quick FID, field check and baseline correction purposes. Each of the sub-cycles [(+) and (-)] of each such cycle involves a capture and utilization of thousands of FID points in data reduction. Each sub-cycle occurs on the order of a second and the number of such sub-cycles employed depends on the desired smoothing and signal to noise ratio (S/N); generally S/N improves in a square root relationship to the number of cycles.

As noted in above cited Dechene et al. references, in requiring greater accuracy and reliability, sample tube composition can distort readings. If glass is not used (and it is preferred to avoid glass in industrial usage), then the replacement should not be a hydrocarbon plastic. But fluorocarbons can be effective in several applications since signals from fluorine appear far from resonance. These signals can be distinguished from hydrogen at the levels of sensitivity required and if desired can be filtered (or distinguished). In other cases of higher sensitivity measurements, e.g., for gauging relative proportions of amorphous and crystalline species in mixtures thereof, the sample container should be glass or non-protonic ceramic. In some instances, however, fluorocarbon or reinforced fluorocarbon can be used acceptably for polymer measurements. In all such cases the point is to avoid sample containers with species that can couple with transmitted energy and generate a FID decay curve mimicking the samples.

It will now be apparent to those skilled in the art that other embodiments, improvements, details, and uses can be made consistent with the letter and spirit of the foregoing disclosure and within the scope of this patent, which is limited only by the following claims, construed in accordance with the patent law, including the doctrine of equivalents.