CN102841063B - Method for tracing and identifying charcoal based on spectrum technology - Google Patents

Abstract

The invention discloses a method for tracing and identifying charcoal based on a spectrum technology, which comprises the following steps of: (1) building a charcoal modeling spectrum library; (2) carrying out principal component analysis for the charcoal modeling spectrum library, and taking the front three principal components according to the variance contribution rate of principal components; (3) calculating a three-dimensional central coordinate of the charcoal variety according to the score values of the front three principal components; (4) scanning the spectrum of a sample to be identified; (5) carrying out principal component analysis for the spectrum of the sample to be identified, extracting the front three principal components and calculating the mahalanobis distances between the front three principal components and the central coordinates of the varieties according to the variance contribution rate of the principal components; and (6) classifying the sample to be identified to the corresponding charcoal varieties according to the minimal mahalanobis distance between the sample to be identified and the central coordinate of each variety, thereby finishing the identifying and tracing of the charcoal. The invention belongs to a method for monitoring machine learning, achieves the aim of quickly tracing and identifying the charcoal through the optical indirect testing and numerical value processing ways, so that the method can be applied to the laboratories or the job sites.

Description

A kind of charcoal based on spectral technique discrimination method of tracing to the source

Technical field

The present invention relates to a kind of charcoal based on spectral technique discrimination method of tracing to the source.

Background technology

Charcoal is the product that living beings form after rough burning thermal cracking under anoxia condition, and solubility is extremely low, and factor of porosity and specific surface area are large, and absorption affinity, resistance to oxidation and antibiont capacity of decomposition are strong.Charcoal increases the aspects such as remittance reduction of discharging and contaminated environment reparation at soil fertility improvement, the solid carbon of soil and has remarkable effect.Charcoal physicochemical property is not only subject to the impact of the technological parameters such as temperature in carbonization process, time, intensity, also relevant with biological material.According to (Xu R such as Xu, Xiao S, Yuan J, et al.Adsorption of methyl violet from aqueous solutions by the biochars derived from crop residues.Bioresource Technology, 2011,102 (22): 10293-10298.), there is difference to the adsorptive power of methyl violet in aqueous solution (being commonly called as gentian violet, a kind of coloring agent) in the charcoal taking castor-oil plant stalk, peanut stalk, soybean stalk and rice husk as raw material.Charcoal uses department often also to need, to the charcoal discriminating of tracing to the source, therefrom to filter out the charcoal with specific material properties.The living beings in various sources are prepared into after charcoal through charring process processing, are all aterrimus on apparent, the naked eyes examination that is not easily distinguishable.And traditional physical and chemical testing complex steps, the test duration is long, is not suitable for the quick discriminating of charcoal.Therefore be necessary to carry out trace to the source fast discrimination method research of charcoal.

In recent years, infrared spectrum technology because its test is quick, reproducible, measuring accuracy compared with high, be easy to on-the-spotly use, the advantages such as agent of being excused from an examination obtain extensive concern.Chinese scholars is to orange silk tree (Chia C H, Gong B, Joseph S D, et al.Imaging of Mineral-enriched biochar by FTIR, Raman and SEM-EDX.Vibrational Spectroscopy.2012.In Press.), pecan shell (Novak J M, Busscher W J, Watts D W, et al.Short-term CO ₂mineralization after additions ofbiochar and switchgrass to a Typic Kandiudult.Geoderma.2010, 154 (3-4): 281-288.), switchgrass (Kumar S, Loganathan V A, Gupta R B, et al.An Assessment of U (VI) removal from groundwater using biochar produced from hydrothermal carbonization.Journal of Environmental Management.2011, 92 (10): 2504-2512.), beet root (Dong X, Ma L Q, Li Y.Characteristics and mechanisms of hexavalent chromium removal by biochar from sugar beet tailing.Journal of Hazardous Materials.2011, 190 (1-3): 909-915.) etc. charcoal carried out diffuse reflectance infrared spectroscopy analysis.But utilizing spectral technique charcoal to be traced to the source fast in discriminating research, yet there are no bibliographical information.

Summary of the invention

In order to reach the object of easy quick discriminating charcoal kind, the invention provides a kind of charcoal based on spectral technique discrimination method of tracing to the source.

(1) set up charcoal kind library of spectra: the biological carbon materials of collecting known not homology matter kind, every kind of material 5-10 sample, utilize spectrometer to gather the reflected spectrum data of each sample, the spectroscopic data of recording wavelength scope 350-1100nm, sampling interval 1nm, totally 751 wavelength points, record the light intensity value of each wavelength points, and make baseline offset and proofread and correct, set up charcoal modeling library of spectra.

Baseline offset updating formula is as follows:

F(x)＝f(x)-minf(X)

In formula, f (x) is the light intensity of wavelength x, and X represents all wavelengths variable, and F (x) represents the light intensity value after wavelength x proofreaies and correct, and min (X) represents the minimum value of the light intensity value of all wavelengths point; .

Charcoal modeling library of spectra is made up of the spectrum record of all samples.The spectrum record of each sample comprises the light intensity value (751 field) after charcoal variety name (1 field), different wave length point calibration, conventionally also can comprise and sort out code (1 field).

(2) charcoal modeling library of spectra is carried out to principal component analysis (PCA), by major component variance contribution ratio size, get front 3 major components of described contribution rate maximum;

The calculation procedure of described principal component analysis (PCA) is as follows:

For n sample, each sample, with carrying out descriptive study object with p variable, is used respectively X ₁, X ₂x _prepresent X _prepresent the vector that the light intensity value of p all samples of wavelength points forms in order, the p dimension random vector that this p variable forms is X=(X ₁, X ₂x _p) ^t.If the average of random vector X is μ, correlation matrix is R, and X is carried out to linear transformation, considers the linear combination of original variable:

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=u_{21}{X}_1+u_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ \·\\ \·\\ \·\\ Z_p=u_{p1}{X}_1+u_{p2}{X}_2+\·\·\·+u_{\mathrm{pp}}{X}_p\end{array}\right.$

Major component is incoherent linear combination Z ₁, Z ₂z _p, and Z ₁x ₁, X ₂x _plinear combination in variance the maximum, Z ₂be and Z ₁variance the maximum in incoherent linear combination, Z _pbe and Z ₁, Z ₂z _p-1all variance the maximum in incoherent linear combination.

Calculate correlation matrix R, R can reflect degree of correlation between data,

$R=\left[\begin{array}{cccc}{R}_{11}& R_{12}& \·\·\·& R_{1p}\\ R_{21}& R_{22}& \·\·\·& R_{2p}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ R_{p1}& R_{p2}& \·\·\·& R_{\mathrm{pp}}\end{array}\right]$

Wherein R _ij(i, j=1,2 ... p) be original variable X _iwith X _jrelated coefficient.Because R is that real symmetric matrix (is R _ij=R _ji), only need to calculate triangle element or lower triangle element on it, R _ijcomputing formula be:

$R_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}{\sqrt{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)}^2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}^2}}$

N represents the sample number of all samples;

Separate secular equation according to correlation matrix R

|λE-R|=0

Obtain eigenvalue λ _i(i=1,2 ... p).Because R is positive definite matrix, so its eigenvalue λ _ibe all positive number, by λ _iorder sequence by size,

λ ₁≥λ ₂≥…≥λ _i≥0

Obtain respectively corresponding to eigenvalue λ _iproper vector e _i(i=1,2 ..., p), require || e _i||=1,

$\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{e}_{\mathrm{ij}}^2=1$

Wherein e _ijrepresent vectorial e _ij component;

Eigenwert is the variance of each major component, and its size has reflected the influence power of each major component.Major component Z _icontribution rate W _icomputing formula be

$W_i=\frac{{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{j=1}^p{\mathrm{\λ}}_i}$

According to major component variance contribution ratio size, get front 3 major component Z of described contribution rate maximum ₁, Z ₂, Z ₃.

Arrange the λ of first three according to size order ₁, λ ₂, λ ₃, calculate according to the following formula the major component load of three major components

$u_{\mathrm{ij}}=p(Z_i,X_j)=\sqrt{{\mathrm{\λ}}_i}{e}_{\mathrm{ij}}(i,j=\mathrm{1,2},\·\·\·,p)$

After obtaining the load of each major component, can calculate according to the following formula the score of each major component

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=l_{21}{X}_1+l_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ Z_3=u_{31}{X}_1+u_{32}{X}_2+\·\·\·+u_{3p}{X}_p\end{array}\right.$

Z ₁, Z ₂, Z ₃represent 3 new variables indexs of n sample; [Z ₁z ₂z ₃] be the new variables matrix of n × 3;

(3) according to front 3 major component Z ₁, Z ₂, Z ₃, the data of the different samples of same charcoal kind are got to average, calculate the three-dimensional center coordinate of each kind of charcoal.

${\stackrel{\‾}{x}}_{i,k}=\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{x}_{i,j,k}}{m},(k=\mathrm{1,2,3}),$

K represents k kind major component, namely k dimension;

charcoal modeling library of spectra i kind sample spectrum k dimension major component mean value;

X _{i, j, k}: charcoal modeling library of spectra i kind j sample spectrum k dimension principal component scores value;

M: the sample number of charcoal modeling library of spectra i kind;

The three-dimensional center coordinate of biological carbon i kind is

(4) sample spectral scan to be identified: utilize spectrometer to gather the reflected spectrum data of each sample, the spectroscopic data of recording wavelength scope 350-1100nm, record the light intensity value of 751 wavelength points identical with step (1), and make baseline offset and proofread and correct, the same step of method (1), sets up sample library of spectra to be identified.

(5) principal component analysis (PCA) is done in sample light spectrum to be identified storehouse, by principal component contributor rate size, extract front 3 major components of contribution rate maximum, the same step of method (2), calculates the score value that obtains of each major component, thereby obtains the three-dimensional center coordinate of each sample to be identified; .And calculate the mahalanobis distance between the three-dimensional center coordinate of sample to be identified and each kind.To appointing a sample x=(x to be identified ₁, x ₂, x ₃) ^t, the centre coordinate of i kind mahalanobis distance formula is:

$D(x,\stackrel{\‾}{x})=\sqrt{{(x-\stackrel{\‾}{x})}^t{S}^{-1}(x-\stackrel{\‾}{x})}$

In formula, S ^-1be the inverse matrix of S, S is the covariance matrix of the principal component scores matrix of all samples of i kind in charcoal modeling library of spectra,

$S=\left[\begin{array}{ccc}{S}_{11}& S_{12}& S_{13}\\ S_{21}& S_{22}& S_{23}\\ S_{31}& S_{32}& S_{33}\end{array}\right]$

Wherein S _kl(k, l=1,2,3) are the principal component scores value vector Z of all sample k dimensions in i kind _kprincipal component scores value vector Z with all sample l dimensions _lbetween covariance, and S _kl=S _lk, its computing formula is:

$S_{\mathrm{kl}}=\frac{\underset{t=1}{\overset{m}{\mathrm{\Σ}}}(Z_{\mathrm{tk}}-{\stackrel{\‾}{Z}}_k)(Z_{\mathrm{tl}}-{\stackrel{\‾}{Z}}_l)}{m-1}$

M: the sample number of charcoal modeling library of spectra i kind;

(6) according to sample to be identified to the minimum mahalanobis distance between each variety central coordinate, sample to be identified is referred to corresponding charcoal kind, thereby completes the charcoal discriminating of tracing to the source.

The biological carbon materials of not homology matter kind of the present invention can be from agriculture and forestry organic waste material, house refuse, charcoal base manure, farm animal excrement etc.

Particularly, below 14 kinds of charcoals spectroscopic data by the inventive method detect modeling, its three-dimensional major component centre coordinate is as shown in the table:

On the whole, the present invention carries out differential test according to charcoal spectral signature to unknown sample, is a kind of Supervised machine learning method.This method utilizes optics indirectly testing and numerical value to process means, reaches the charcoal discriminating object of tracing to the source fast.

The present invention has multiple advantages compared with conventional physical and chemical analysis technology:

(1) fast detecting.The single charcoal test sample time is only several seconds.

(2) on-the-spot use.Spectrometer is easy to carry, not limited by place.

(3) batch testing.Can carry out analyzing and testing to great amount of samples, the agent of being excused from an examination, pollution-free.

Brief description of the drawings

Fig. 1 is the mirror method for distinguishing schematic diagram of tracing to the source fast of the charcoal based on spectral technique.

The typical emission spectra figure of 14 kinds of charcoals in Fig. 2 embodiment.

Embodiment

Collect following 14 provenance matter charcoals, wherein 7 kinds of agricultural wastes (corn, rice chaff, peanut shell, pecan shell, rice straw, straw), 4 kinds of forestry waste (mao bamboon, edible fungi residue, wood chip, caragana microphylla), 1 kind of house refuse, a kind of charcoal composite fertilizer, a kind of excreta of cultivation (pig manure).For simplicity, table 1 has been listed these charcoal Chinese and English title and initialisms.

Table 1

Charcoal is levigate with mortar, obtains particle diameter and is less than 0.1mm sample.Every kind of charcoal is prepared 10 samples, and totally 140 samples are used for analytical test.

The measurement mechanism of charcoal sample comprises the compositions such as the Maya2000Pro spectrometer, 10W halogen light source, Transmission Fibers, sample rack, computing machine, power supply that marine optics company produces.Wavelength coverage is selected 350-1100 wave band.Fibre-optical probe is vertically positioned over sample top 1cm left and right place, and source light is connected to probe by optical fiber, and light, after sample reflection, is received through Optical Fiber Transmission to spectrometer by probe again, and spectroscopic data sends on computing machine through USB interface.Use computer software to process spectroscopic data, and provide result of calculation.

(1) set up 14 kinds of charcoal modeling library of spectra: to above 14 kinds of biological carbon materials of homology matter kind not, 8 samples of every kind of material, utilize spectrometer to gather the reflected spectrum data of each sample, the spectroscopic data of recording wavelength scope 350-1100nm, sampling interval 1nm, totally 751 wavelength points, record light intensity value, and make baseline offset and proofread and correct, set up charcoal modeling library of spectra.Canonical biometric charcoal emission spectrum as shown in Figure 2.

Baseline offset updating formula is as follows:

F(x)＝f(x)-minf(X)

In formula, f (x) is the light intensity value of wavelength x, and X represents all wavelengths variable, and F (x) represents the light intensity value after wavelength x proofreaies and correct, and minf (X) represents the minimum value of the light intensity value of all wavelengths point;

Charcoal modeling library of spectra is made up of the spectrum record of all samples.The spectrum record of each sample comprises the light intensity value (751 field) after charcoal variety name (1 field), different wave length point calibration, sorts out code (1 field).

(2) to charcoal modeling library of spectra, utilize Unscrambler X10.1 software principal component analysis (PCA) functional module to calculate, by principal component contributor rate size, get front 3 major components of contribution rate maximum.

The Computing Principle of described principal component analysis (PCA) is as follows:

For n sample, each sample carrys out descriptive study object with p variable, uses respectively X ₁, X ₂x _prepresent X _prepresent the vector that the light intensity value of p all samples of wavelength points forms in order, the p dimension random vector that this p variable forms is X=(X ₁, X ₂x _p) ^t;

Calculate correlation matrix R,

$R=\left[\begin{array}{cccc}{R}_{11}& R_{12}& \·\·\·& R_{1p}\\ R_{21}& R_{22}& \·\·\·& R_{2p}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ R_{p1}& R_{p2}& \·\·\·& R_{\mathrm{pp}}\end{array}\right]$

Wherein R _ij(i, j=1,2 ... p) be original variable X _iwith X _jrelated coefficient, and R _ij=R _ji, its computing formula is:

$R_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}{\sqrt{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)}^2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}^2}}$

Separate secular equation according to correlation matrix R,

|λE-R|=0

Obtain eigenvalue λ _i(i=1,2 ... p); By λ _iorder sequence by size,

λ ₁≥λ ₂≥…≥λ _i≥0

Obtain respectively corresponding to eigenvalue λ _iproper vector e _i(i=1,2 ..., p), require || e _i||=1,

$\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{e}_{\mathrm{ij}}^2=1$

Wherein e _ijrepresent vectorial e _ij component;

Calculate major component Z _icontribution rate W _i, computing formula is

$W_i=\frac{{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{j=1}^p{\mathrm{\λ}}_i}$

According to principal component contributor rate size, get front 3 major component Z of described contribution rate maximum ₁, Z ₂, Z ₃;

Arrange the λ of first three according to size order ₁, λ ₂, λ ₃, calculate according to the following formula the major component load of three major components

$u_{\mathrm{ij}}=p(Z_i,X_j)=\sqrt{{\mathrm{\λ}}_i}{e}_{\mathrm{ij}}(i=\mathrm{1,2,3},j=\mathrm{1,2},\·\·\·,p)$

After obtaining the load of each major component, calculate according to the following formula the score of each major component

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=l_{21}{X}_1+l_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ Z_3=u_{31}{X}_1+u_{32}{X}_2+\·\·\·+u_{3p}{X}_p\end{array}\right.$

Z ₁, Z ₂, Z ₃represent 3 new variables indexs of n sample; [Z ₁z ₂z ₃] be the new variables matrix of n × 3;

(3) according to front 3 major components, the data of the different samples of same charcoal kind are got to average, calculate the three-dimensional center coordinate of each kind of charcoal.

${\stackrel{\‾}{x}}_{i,k}=\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{x}_{i,j,k}}{m},(k=\mathrm{1,2,3})$

K represents k kind major component, namely k dimension;

charcoal modeling library of spectra i kind sample spectrum k dimension major component mean value;

X _{i, j, k}: charcoal modeling library of spectra i kind j sample spectrum k dimension major component value;

M: the sample number of charcoal modeling library of spectra i kind; In the present embodiment, m=8.

The three-dimensional center coordinate of biological carbon i kind is

Acquired results is as table 2.

Table 2

(4) sample spectral scan to be identified: utilize spectrometer to gather the reflected spectrum data of each sample, the spectroscopic data of recording wavelength scope 350-1100nm, sampling interval 1nm, totally 751 wavelength points, record light intensity value, and make baseline offset and proofread and correct, the same step of method (1), 2 samples of each kind, totally 28 samples, set up sample library of spectra to be identified.

(5) by step (2) same method, principal component analysis (PCA) is done in sample light spectrum to be identified storehouse, by major component variance contribution ratio size, extract front 3 major components of contribution rate maximum, calculate the score value that obtains of each major component, thereby obtain the three-dimensional center coordinate of each sample to be identified; As shown in table 3.

Table 3

And calculate the mahalanobis distance between the three-dimensional center coordinate of each sample to be identified and each kind, the results are shown in Table 4.

To arbitrary sample x=(x to be identified ₁, x ₂, x ₃) ^t,, the centre coordinate of i kind mahalanobis distance formula is:

$D(x,\stackrel{\‾}{x})=\sqrt{{(x-\stackrel{\‾}{x})}^t{S}^{-1}(x-\stackrel{\‾}{x})}$

In formula, S ^-1be the inverse matrix of S, S is the covariance matrix of the principal component scores matrix of all samples of i kind in charcoal modeling library of spectra,

$S=\left[\begin{array}{ccc}{S}_{11}& S_{12}& S_{13}\\ S_{21}& S_{22}& S_{23}\\ S_{31}& S_{32}& S_{33}\end{array}\right]$

Wherein S _kl(k, l=1,2,3) are the principal component scores value vector Z of all sample k dimensions in i kind _kprincipal component scores value vector Z with all sample l dimensions _lbetween covariance, and S _kl=S _lk, its computing formula is:

$S_{\mathrm{kl}}=\frac{\underset{t=1}{\overset{m}{\mathrm{\Σ}}}(Z_{\mathrm{tk}}-{\stackrel{\‾}{Z}}_k)(Z_{\mathrm{tl}}-{\stackrel{\‾}{Z}}_l)}{m-1}$

M: the sample number of charcoal modeling library of spectra i kind;

Table 4

(6) according to sample to be identified to the minimum mahalanobis distance between each variety central coordinate, sample to be identified is referred to corresponding charcoal kind, as shown in table 4, thereby complete the charcoal discriminating of tracing to the source.According to the method for the invention, above-mentioned 28 samples to be identified are sorted out to test, wherein 27 samples are sorted out correctly, and 1 F class sample (house refuse charcoal) is mistaken for L class (Caragana Microphylla biological charcoal), and driscrimination error is 3.6%.

Claims (1)

1. the discrimination method of tracing to the source of the charcoal based on spectral technique, is characterized in that the step of the method is as follows:

(1) set up charcoal kind library of spectra: the biological carbon materials of collecting known not homology matter kind, every kind of material 5-10 sample, utilize spectrometer to gather the reflected spectrum data of each sample, the spectroscopic data of recording wavelength scope 350-1100nm, sampling interval 1nm, record altogether the light intensity value of 751 wavelength points, and make baseline offset and proofread and correct, set up charcoal modeling library of spectra;

Baseline offset updating formula is as follows:

F(x)=f(x)-minf(X)

In formula, f (x) is the light intensity value of wavelength x, and X represents all wavelengths variable, and F (x) represents the light intensity value after wavelength x proofreaies and correct, and minf (X) represents the minimum value of the light intensity value of all wavelengths point;

Charcoal kind library of spectra is made up of the spectrum record of all samples, and the spectrum record of each sample comprises the light intensity value after charcoal variety name, different wave length point calibration;

(2) charcoal modeling library of spectra is carried out to principal component analysis (PCA), by major component variance contribution ratio size, get front 3 major components of described contribution rate maximum;

The calculation procedure of described principal component analysis (PCA) is as follows:

For n sample, each sample carrys out descriptive study object with p variable, uses respectively X ₁, X ₂x _prepresent X _prepresent the vector that the light intensity value of p all samples of wavelength points forms in order, the p dimension random vector that this p variable forms is X=(X ₁, X ₂x _p) ^t;

Calculate correlation matrix R,

$R=\left[\begin{array}{cccc}{R}_{11}& R_{12}& \·\·\·& R_{1P}\\ R_{21}& R_{22}& \·\·\·& R_{2p}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ R_{p1}& R_{p2}& \·\·\·& R_{\mathrm{pp}}\end{array}\right]$

Wherein R _ijfor original variable X _iwith X _jrelated coefficient, i, j=1,2 ... p, and R _ij=R _ji, its computing formula is:

$R_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}\left(X_{\mathrm{ki}}{\stackrel{\‾}{X}}_i\right)(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}{\sqrt{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{{(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)}^2}^{}{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}}^{}{{(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}_{}}^2}}$

N represents the sample number of all samples;

Separate secular equation according to correlation matrix R,

|λE-R|＝0

Obtain eigenvalue λ _i, i=1,2 ... p; By λ _iorder sequence by size,

λ ₁≥λ ₂≥…≥λ _i≥0

Obtain respectively corresponding to eigenvalue λ _iproper vector e _i, i=1,2 ..., p, requires || e _i||=1,

$\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{e}_{\mathrm{ij}}^2=1$

Wherein e _ijrepresent vectorial e _ij component;

Calculate major component Z _icontribution rate W _i, computing formula is

$W_i=\frac{{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{j=1}^p{\mathrm{\λ}}_i}$

According to major component variance contribution ratio size, get front 3 major component Z of described contribution rate maximum ₁, Z ₂, Z ₃;

Arrange the λ of first three according to size order ₁, λ ₂, λ ₃, calculate according to the following formula the major component load of three major components

$u_{\mathrm{ij}}=p(Z_i,X_j)=\sqrt{{\mathrm{\λ}}_i}{e}_{\mathrm{ij}}$

Wherein i=1,2,3, j=1,2 ..., p

After obtaining the load of each major component, calculate according to the following formula the score of each major component

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=u_{21}{X}_1+u_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ Z_3=u_{31}{X}_1+u_{32}{X}_2+\·\·\·+u_{3p}{X}_p\end{array}\right.$

Z ₁, Z ₂, Z ₃represent 3 new variables indexs of n sample; [Z ₁z ₂z ₃] be the new variables matrix of n × 3;

(3) according to front 3 major component Z ₁, Z ₂, Z ₃, the data of the different samples of same charcoal kind are got to average, calculate the three-dimensional center coordinate of each kind of charcoal;

${\stackrel{\‾}{x}}_{i,k}=\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{x}_{i,j,k}}{m}$

K=1,2,3, k represents k kind major component, namely k dimension;

charcoal modeling library of spectra i kind sample spectrum k dimension major component mean value;

X _{i, j, k}: charcoal modeling library of spectra i kind j sample spectrum k dimension principal component scores value;

M: the sample number of charcoal modeling library of spectra i kind;

The three-dimensional center coordinate of biological carbon i kind is

(4) sample spectral scan to be identified: utilize spectrometer to gather the reflected spectrum data of each sample, the spectroscopic data of recording wavelength scope 350-1100nm, record the light intensity value of 751 wavelength points identical with step (1), and make baseline offset and proofread and correct, the same step of method (1), sets up sample library of spectra to be identified;

(5) principal component analysis (PCA) is done in sample light spectrum to be identified storehouse, by major component variance contribution ratio size, extract front 3 major components of contribution rate maximum, the same step of method (2), calculate the score value that obtains of each major component, thereby obtain the three-dimensional center coordinate of each sample to be identified; Calculate the mahalanobis distance between the three-dimensional center coordinate of sample to be identified and each kind, to arbitrary sample x=(x to be identified ₁, x ₂, x ₃) ^t, the centre coordinate of i kind mahalanobis distance formula is:

$D(x,\stackrel{\‾}{x})=\sqrt{{(x-\stackrel{\‾}{x})}^t{S}^{-1}(x-\stackrel{\‾}{x})}{}^{}$

In formula, S ^-1be the inverse matrix of S, S is the covariance matrix of the principal component scores matrix of all samples of i kind in charcoal modeling library of spectra,

$S=\left[\begin{array}{ccc}{S}_{11}& S_{12}& S_{13}\\ S_{21}& S_{22}& S_{23}\\ S_{31}& S_{32}& S_{33}\end{array}\right]$

Wherein S _klit is the principal component scores value vector Z of all sample k dimensions in i kind _kprincipal component scores value vector Z with all sample l dimensions _lbetween covariance, k, l=1,2,3, and S _kl=S _lk, its computing formula is:

$S_{\mathrm{kl}}=\frac{\underset{t=1}{\overset{m}{\mathrm{\Σ}}}(Z_{\mathrm{tk}}-{\stackrel{\‾}{Z}}_k)(Z_{\mathrm{tl}}-{\stackrel{\‾}{Z}}_l)}{m-1}$

M: the sample number of charcoal modeling library of spectra i kind;

(6) according to sample to be identified to the minimum mahalanobis distance between each variety central coordinate, sample to be identified is referred to corresponding charcoal kind, thereby completes the charcoal discriminating of tracing to the source.2, the method for claim 1, is characterized in that: the three-dimensional major component centre coordinate of spectrum of following 14 kinds of charcoals is as shown in the table: