CN117872449A - Nuclide identification method in complex environment - Google Patents

Abstract

The invention relates to a nuclide identification method in a complex environment, which comprises the following steps: s1, background subtraction; s2, solving the content of natural nuclides: calculating the content value of the natural nuclide to be removed according to the gamma energy spectrum data to be analyzed after background subtraction and the inverse matrix of the response coefficient of the natural nuclide to be removed; s3, removing natural nuclides: according to the count on each track address and the content value of the natural nuclide to be removed, gamma energy spectrum data to be analyzed after the natural nuclide is removed is obtained; s4, smooth spectral lines; s5, peak searching; s6, determining peak positions and nuclide types. According to the invention, under the condition of unknown nuclides, the interference of the multi-energy peak nuclides on nuclide identification is reduced by removing the multi-energy peak nuclides, and the accuracy of nuclide identification is improved; the smooth times of the spectral lines are fully automatically adjusted, so that the smooth effect of the spectral lines is improved; the combined peak searching method combining the derivative peak searching method and the symmetrical zero area peak searching method has the advantages of high peak searching speed, high peak searching precision, short peak searching time and light load on the CPU.

Description

Nuclide identification method in complex environment

Technical Field

The invention relates to a nuclide identification method, in particular to a nuclide identification method in a complex environment.

Background

With the continuous development of society, nuclear technology application and nuclear radiation monitoring technology play an increasingly important role in the fields of daily life, social and economic construction, environmental protection, national defense construction, nuclear anti-terrorism and the like, and the generated radioactive hazard is increasingly valued by more and more people. On the one hand, the threat of nuclear material diffusion and nuclear terrorism is increasingly serious, and new requirements are put on nuclear security and nuclear detection technology; on the other hand, nuclear energy is increasingly used in human society, and although production can be promoted, the frequency of nuclear accidents is also increasing. The catastrophic consequences of nuclear radiation are alarming, so that the concern about nuclear safety is also increasing, and the monitoring and identification of radionuclides is a hotspot problem in current research. Therefore, the national natural science foundation committee is also very important for the development of nuclear technology, and in the "research report of the development strategy of nuclear technology science", the research of nuclear analysis methods and applications thereof is listed as a direction of preferential support.

As one of the most important means in the nuclear detection and analysis method, the nuclide identification technology can realize detection, positioning and qualitative identification of radioactive substances, can also alarm radioactive pollution, and can realize monitoring of radioactive substances in all directions. The method can rapidly give out the property of nuclide, is favorable for commanders to make targeted accurate decisions, and has important significance for improving the accident handling efficiency. In addition, the nuclide identification device can rapidly, flexibly and conveniently perform radioactive measurement, and is widely applied to the aspects of nuclear power stations, food hygiene inspection, environmental radiation pollution, nuclear waste monitoring, geoscience, nuclear medicine instrument radioactive detection, radioactive source theft prevention, nuclear terrorism and the like.

In the conventional analysis of gamma energy spectrum, as shown in fig. 1, the nuclide identification method is generally implemented by searching a nuclide library to perform gamma energy spectrum characteristic peak matching. As a key ring of the nuclide identification algorithm, after a peak searching algorithm finds a plurality of peak positions from the gamma energy spectrum, the peak positions are converted into energy values according to a set energy scale formula, and the identification algorithm matches the obtained peak position energy values with the characteristic peak energies of all nuclides in the library one by one, and the nuclide with the best matching degree is regarded as the identified nuclide.

At present, most of common nuclide identification algorithms lack nuclide identification capability under uranium series Ra-226 and thorium series Th-232 backgrounds, and nuclide identification accuracy is low.

Disclosure of Invention

The invention aims to provide a nuclide identification method in a complex environment, so as to solve the problem of inaccurate nuclide identification in uranium series Ra-226 and thorium series Th-232 backgrounds.

The purpose of the invention is realized in the following way:

the nuclide identification method in the complex environment comprises the following steps:

s1, background subtraction: background subtraction is carried out on gamma energy spectrum data to be analyzed by utilizing an SNIP algorithm;

s2, solving the content of natural nuclides: calculating the content value of the natural nuclide to be removed according to the inverse matrix of the response coefficient of the natural nuclide to be removed and the count of gamma energy spectrum data to be analyzed after background subtraction on each channel address, wherein the characteristic peak of the natural nuclide to be removed is a multipotent peak;

s3, removing natural nuclides: calculating the contribution value of the natural nuclide to be removed on each channel site according to the content value of the natural nuclide to be removed, subtracting the contribution value of the natural nuclide to be removed on each channel site from the count of the gamma energy spectrum data to be analyzed after background subtraction on each channel site, and obtaining the gamma energy spectrum data to be analyzed after the natural nuclide is removed;

s4, smooth spectral line: calculating the gamma total count of gamma energy spectrum data to be analyzed after natural nuclides are removed, calculating the smoothing times according to the gamma total count after the natural nuclides are removed and the relation between the smoothing times and the gamma total count, and carrying out smoothing treatment on the gamma energy spectrum data to be analyzed after the natural nuclides are removed based on the smoothing times;

s5, peak searching: carrying out peak searching on gamma energy spectrum data to be analyzed after smooth processing by utilizing a derivative peak searching method and a symmetrical zero area peak searching method to obtain a characteristic peak array;

s6, determining peak positions and nuclide types: converting the addresses of characteristic peaks in the characteristic peak array into characteristic peak energy, and determining the nuclide types according to the converted characteristic peak energy.

Further, in step S2, the specific way of calculating the inverse matrix of the response coefficient of the natural nuclide to be removed is as follows:

s2a-1, calculating a response coefficient matrix of natural nuclides to be rejected based on the total count of known mixed nuclides on each track site and the content of each nuclide in the known mixed nuclides, wherein the known mixed nuclides comprise the natural nuclides to be rejected;

s2a-2, calculating a response coefficient inverse matrix of the natural nuclide to be removed according to the response coefficient matrix of the natural nuclide to be removed.

Further, the specific way of calculating the content value of the natural nuclide to be removed in the step S2 is as follows:

and multiplying the inverse matrix of the response coefficient of the natural nuclide to be removed by the count of the gamma energy spectrum data to be analyzed after background subtraction on each channel address, and calculating the content value of the natural nuclide to be removed.

Further, the specific mode of smoothing the data after removing the natural nuclides is as follows:

and carrying out five-point three-time polynomial smooth calculation on each point in the gamma energy spectrum data to be analyzed according to the smooth times to obtain the gamma energy spectrum data to be analyzed after being smoothed.

Further, the relation between the number of smoothness and the gamma total count value is:

Q=-0.00798x+24.44

wherein,Qrepresenting the number of repetitions of the smoothing process,xrepresenting the total gamma count per 1 second.

Further, the specific mode of peak searching for the smooth gamma energy spectrum data to be analyzed in the step S5 is as follows:

s5-1, carrying out peak searching on the smooth gamma energy spectrum data to be analyzed by adopting a derivative peak searching method to obtain a first characteristic peak address, and storing the first characteristic peak address into a first characteristic peak array until the last address of the gamma energy spectrum data to be analyzed is traversed;

s5-2, carrying out secondary peak searching on the gamma energy spectrum data to be analyzed after being smoothed according to a symmetrical zero-area peak searching method and the first characteristic peak array, and determining a second characteristic peak array.

Further, the specific mode of carrying out peak searching on the smooth gamma energy spectrum data to be analyzed in the step S5-1 is as follows:

s5-1-1, calculating a first derivative value of each channel address in the smooth gamma energy spectrum data to be analyzed to obtain a first derivative spectrum;

s5-1-2. If the address isiThe first derivative value of (2) is positive, adjacent track addressi+1The first derivative value of (2) is negative, then at the track addressiSum of addressesi+1Track address with characteristic peak in betweenxWherein, the method comprises the steps of, wherein,x=i+0.5the method comprises the steps of carrying out a first treatment on the surface of the Traversing the whole first derivative spectrum to obtain an initial characteristic peak arrayX[x ₁ ,x ₂ …x _n ]；

S5-1-3. If the address isjThe first derivative value of (2) is negative, adjacent addressesj+1The first derivative value of (2) is positive, then at the track addressjSum of addressesj+1With a characteristic peak boundary betweenzWherein, the method comprises the steps of, wherein,z=j+0.5the method comprises the steps of carrying out a first treatment on the surface of the Traversing the whole first derivative spectrum to obtain a characteristic peak boundary arrayZ[z ₁ ,z ₂ …z _n ]；

S5-1-4. Taking initial characteristic peak arrayX[x ₁ ,x ₂ …x _n ]Is a track address of a characteristic peak of (a)x _i Array at characteristic peak boundaryZFind its nearest left and right boundariesz _m Andz _n the method comprises the steps of carrying out a first treatment on the surface of the If the characteristic peak width isWFailure to meet the condition FWHM<W<(4 XFWHM), thenx _i A trace address that is not a characteristic peak; if the condition FWHM is satisfied<W<(4 xfwhm), then go to the next step, where,W=z _n -z _m FWHM is full width at half maximum;

s5-1-5. If the characteristic peak is located at the trackx _i Less than or equal to 2 times the background per second countx _i A trace address that is not a characteristic peak; if the track address of the characteristic peakx _i Is greater than 2 times the background, then the trace address of the characteristic peakx _i Is the first characteristic peak address and willx _i And storing the first characteristic peak array.

Further, the specific way of determining the second characteristic peak array in step S5-2 is as follows:

s5-2-1. For each first characteristic peak addressx _i Setting a calculation interval to obtain a calculation interval of each first characteristic peak address, wherein the calculation interval is as follows:x _i -2×FWHM～x _i +2×FWHM；

s5-2-2 for each first characteristic peak address calculation intervalx _i -2×FWHM～x _i +2 xfwhm) by performing a convolution transform for each address:wherein (1)>Is the address of the trackiValues after convolution transformation +.>Is the address of the tracki+jIs used for counting the number of the counts,C _j is peak-like function, ++>，/>；

S5-2-3, judging gamma energy spectrum data to be analyzed in the calculation interval of each first characteristic peak address, if the gamma energy spectrum data is positive, determining that the gamma energy spectrum data is positiveR _i >fWill thenR _i Corresponding track addressiDetermining the first characteristic peak address as a second characteristic peak address, storing the second characteristic peak address into a second characteristic peak array until the calculation interval of each first characteristic peak address in the first characteristic peak array is traversed, wherein,，/>is->Root mean square of>，fIs a peak-finding threshold.

Further, in step S6, the specific way of converting the address of the characteristic peak in the characteristic peak array into the characteristic peak energy is as follows:

converting the track address of the characteristic peak into the energy of the characteristic peak according to an energy scale formula, wherein the energy scale formula is thatT= 3.05x-2.5Wherein, the method comprises the steps of, wherein,xrepresenting the location of the track(s),Trepresenting energy.

Further, in step S6, the specific manner of determining the nuclide species according to the converted characteristic peak energy is as follows:

selecting any nuclide, calculating the absolute value of the difference between the energy of the converted characteristic peak and the standard energy of the characteristic peak of the selected nuclide, and if the calculated absolute value is smaller than or equal to the preset energy window width, the gamma energy spectrum data to be analyzed comprise the selected nuclide.

The invention completes the verification test of the nuclide identification algorithm, and the performance index is superior to the related calibration standard. The invention completes the related test at the nuclear industry radioactive investigation metering station, and can identify radionuclides such as Am-241, ba-133, co-60, cs-137, K-40, ra-226, th-232 and the like; under the natural background condition, the intensity of a single nuclide is in the range of 0.1-0.7 mu Sv/h, the nuclide identification rate is 100%, and the confidence coefficient is more than or equal to 97.8%; the mixed nuclide identification accuracy is 100%, the confidence coefficient is more than or equal to 96.7%, and the performance index is superior to the related calibration standard. The invention adopts SNIP algorithm to carry out background subtraction, reduces the interference of cosmic rays and instruments, and improves the accuracy of nuclide identification. By adopting a unique overlapping peak decomposition method and removing the multi-energy peak nuclide, the interference of the multi-energy peak nuclide on nuclide identification is reduced, and the nuclide identification rate reaches 100%. The smooth times of the spectral lines are fully automatically adjusted, so that the smooth effect of the spectral lines is improved; the combined peak searching method combining the derivative peak searching method and the symmetrical zero area peak searching method has the advantages of high peak searching speed, high peak searching precision, short peak searching time, light load on a CPU and the like, and is particularly suitable for portable nuclide identification equipment.

Drawings

FIG. 1 is a flow chart of a generic species identification algorithm.

Fig. 2 is a flow chart of the present invention.

FIG. 3 is a characteristic peak diagram of Ra-226, th-232 and Ba-133.

FIG. 4 is a decomposition spectrum of a mixed spectrum of Ba-133 containing polynuclein after removal of natural nuclides.

FIG. 5 is a decomposition spectrum of a mixed spectrum of polynuclein containing Cs-137 after removal of natural nuclides.

FIG. 6 is a graph showing the effect of smoothing Cs-137 a different number of times.

FIG. 7 is a graph showing the effect of smoothing a general office place K-40 for various times.

FIG. 8 is a K-40 original spectrum and its first derivative spectrum.

FIG. 9 is a K-40 smooth 20 th order spectrum and its first derivative spectrum.

Fig. 10 is a flow chart of the combined peak finding method.

FIG. 11 is a flow chart of a derivative peak finding method.

FIG. 12 is a Cs-137 smooth spectrum and its first derivative spectrum.

FIG. 13 is a graph showing the effects of a smooth spectrum of Cs-137 and a symmetrical zero-area peak finding method on Cs-137.

FIG. 14 is an Am-241 peak finding test chart.

FIG. 15 is a Co-60 peak finding test chart.

Fig. 16 is an energy scale graph.

Fig. 17 is a graph of mixed radionuclide identification test results.

Detailed Description

The present invention will be described in further detail below.

As shown in fig. 2, the nuclide identification method in the complex environment of the present invention includes the following steps:

s1, background subtraction: and performing background subtraction on gamma energy spectrum data to be analyzed by utilizing an SNIP algorithm.

Under the influence of radioactive factors of cosmic rays and instrument materials, certain radioactivity exists in the environment where us is located, namely, the environment background. Therefore, background subtraction is required for gamma-spectroscopy data to be analyzed, and interference of the instrument itself and cosmic rays is eliminated. SNIP (Simple Non-equivalent Peak) algorithm is used for background subtraction, and interference of factors such as natural background count, cosmic rays and the like is removed.

S2, solving the content of natural nuclides: and calculating the content value of the natural nuclide to be removed according to the inverse matrix of the response coefficient of the natural nuclide to be removed and the count of gamma energy spectrum data to be analyzed after background subtraction on each channel address.

Wherein, the characteristic peak of the natural nuclide to be removed is a multipotent peak.

In nuclide identification, a detector is usually NaI (Tl) crystal, which has the characteristics of high quality, low cost, rapidness and high efficiency, but the energy resolution is not high, and is generally about 8%. Thus, a plurality of adjacent characteristic peaks are overlapped, and the detection and identification of the radionuclide are more difficult.

As shown in Table 1, uranium and thorium are multipotent peaks, so that a plurality of adjacent characteristic peaks are overlapped, and aiming at the problem that a plurality of adjacent characteristic peaks are overlapped in the nuclide identification process, the invention provides a gamma energy spectrum overlapping peak analysis method.

TABLE 1 characteristic peaks of common nuclides

For example, a 609.3keV peak for Ra-226, a 583.2keV peak for Th-232, and a 661.7keV for artificial nuclide Cs-137, will have characteristic peak overlaps; characteristic peak overlap occurs between characteristic peak 1120.3keV of Ra-226 and characteristic peak 1173.2keV of artificial nuclide Co-60; as shown in FIG. 3, characteristic peak overlap occurs for characteristic peak 351.9keV for Ra-226, characteristic peak 338.3keV for Th-232, and characteristic peak 356.0keV for artificial nuclide Ba-133.

The energy resolution of the detector affects the count of one address in the gamma energy spectrum data to be analyzed, which not only has the contribution of the nuclide corresponding to the energy, but also may have the contribution of the peak similar to the energy or the superposition of Compton platforms generated by high-energy gamma rays, wherein the count on each address is the sum of the counts of each single nuclide on the address, and can be expressed as:

（11）

wherein,c _i to mix nuclides at the firstiThe total count on the track is calculated,mto measure the species of nuclides in the gamma energy spectrum,jis the firstjThe species of the species,x _j is the first in the mixed sourcejThe value of the content of the species nuclide,a _ij is the firstjThe seed components areiResponse coefficient of the address.

Using a model body source containing natural nuclides to be removed, and using the content value of each nuclide in known mixed nuclidesx _j And the measured count at each track addressc _i And calculating a response coefficient matrix of each nuclide according to the formula (11), wherein the response coefficient matrix comprises a response coefficient matrix of the natural nuclide to be removed. The content of each species in the different mixed species may not be the same, but the response coefficient matrix of the species is unchanged, thus solving the unknown mixtureWhen nuclides are used, the response coefficient matrix of the nuclides can be used for solving the content of the nuclides.

For example, the response coefficient matrix is calculated using the measured counts at each trace site and the nuclide content of the "response coefficient calculation model sources (including radium model source, thorium model source, potassium model source, and background model source)".

Transforming equation (11) according to the response coefficient matrixa _ij Obtaining a response coefficient inverse matrix, and obtaining the content value of the natural nuclide to be removed according to the response coefficient inverse matrix and the count on each channel address of the gamma energy spectrum data to be analyzed after background subtractionx _j ：

（12）

S3, removing natural nuclides: calculating the contribution value of the natural nuclide to be removed on each channel site according to the content value of the natural nuclide to be removed, subtracting the contribution value of the natural nuclide to be removed on each channel site from the count of the gamma energy spectrum data to be analyzed after background subtraction on each channel site, and obtaining the gamma energy spectrum data to be analyzed after the natural nuclide is removed.

And multiplying the content value of the natural nuclide to be removed by the response matrix of the natural nuclide to be removed by using a formula (11), wherein the multiplied result is the count (contribution value) of the natural nuclide to be removed on each track address.

For example, when the multi-energy peak nuclide in the gamma energy spectrum data to be analyzed is radium element and thorium element, multiplying the response coefficient matrix of radium element and the content value of radium element to obtain the contribution value of radium element on each channel address of the gamma energy spectrum data to be analyzed, obtaining the contribution value of thorium element on each channel address of the gamma energy spectrum data to be analyzed by the method, subtracting the contribution value of radium element and thorium element on each channel address from the total count of the gamma energy spectrum data to be analyzed on each channel address, and obtaining the gamma energy spectrum data to be analyzed after eliminating radium and thorium interference, thereby completing the decomposition of overlapping peaks under the background interference of uranium system Ra-226 and thorium system Th-232.

（13）

Wherein,to at the firstiTotal count on track, +.>Is nuclidemIn the first placeiContribution value on track->Is nuclidenIn the first placeiContribution value on track->Is the firstiThe remaining counts on the tracks.

As shown in fig. 4, when the gamma energy spectrum data to be analyzed comprises radium element, thorium element and barium element, after eliminating the interference of Ra-226 and Th-232, the characteristic peak of Ba-133 is very obvious, and the uniformity is high, similar to the original spectrum form of single nuclide Ba-133; as shown in fig. 5, when the gamma energy spectrum data to be analyzed comprises radium element, thorium element and cesium element, after eliminating the interference of Ra-226 and Th-232, the characteristic peak of Cs-137 is very obvious, and the uniformity is high, similar to the original spectrum form of single nuclide Cs-137; the method achieves the aim of decomposing the overlapped peaks, can easily identify nuclides through the following steps, and improves the accuracy of nuclide identification.

S4, smooth spectral line: and calculating the gamma total count of gamma energy spectrum data to be analyzed after removing natural nuclides, calculating the smoothing times according to the gamma total count after removing the natural nuclides and the relation between the smoothing times and the gamma total count, and carrying out smoothing treatment on the gamma energy spectrum data to be analyzed after removing the natural nuclides based on the smoothing times.

In the nuclear radiation detection process, factors such as noise, detector materials, radioactive decay and the like of a hardware circuit can influence the measurement of gamma energy spectrum, so that measured gamma energy spectrum data have certain statistical fluctuation, namely, the gamma energy spectrum data can deviate from a theoretical result.

The statistical fluctuation is represented in the peak searching process, and the false peak and the weak peak can be covered to cause the loss of the peak, so that the accurate identification of nuclides is affected, and the problems of missing identification, misidentification and the like are caused. Although the gamma-ray energy is discrete, there is some logical correlation between the measured gamma-ray energy as the gamma-ray energy is continuous. The mathematical statistics method is used for processing each data, so that most of the influence of statistical fluctuation can be eliminated, and the important characteristics of the original spectrum can be still maintained.

Therefore, the actually measured gamma spectrum data must be smoothed. The smoothing of the gamma-energy spectrum data requires the preservation of important features of the measured gamma-energy spectrum data, the smoothing of the gamma-energy spectrum data essentially removes noise from the gamma-energy spectrum signal, and the smoothing process typically filters the gamma-energy spectrum data using a time domain least squares method.

The basic principle of the time domain least square method is: left and right of the original spectrum data at a certain pointmTrack data, i.e. co-fetch (2m+1) original spectrum data, the corresponding track address coordinates are%-m,-m+1,…-2,-1,0,1,2,…m-1,m) The corresponding count is%y _-m ,y _-m+1 ,…y _-2 ,y _-1 ,y ₀ ,y ₁ ,y ₂ ,…y _m-1 ,y _m ). By using one of the followingmThe polynomial with the point as the center carries out least square fitting on the acquired spectrum, and the calculated value is the firstmData after point smoothing. And then moving to the right point by point, and calculating each spectrum data in the same way to obtain a complete smooth spectrum. The polynomial used to fit the spectral data can be expressed as: n≤2m+1（14）

wherein,iis the firstiA point of the light-emitting diode is located,is the firstiIndividual pointsIs a fit of (a) to the fitting value of (b).

According to the least squares principle, i.e. letting the fitting valuesAnd actual measured valuesy _i The sum of squares of the differences is minimal, namely:（15）

obtaining coefficients according to equation (15)b _nk And willb _nk The value of (2) is brought into formula (14) to obtain a fitting value of each pointAnd obtaining corresponding fitting smooth data. In practical application, a smooth formula with different points is adopted according to the actual measurement spectrum data condition, and a general formula of the smooth spectrum data with the points is as follows: />（16）

Wherein,is normalized constant, ++>Is smooth and is->And->Can be calculated from the corresponding filters. For a Savitzky-Golay filter, the coefficient calculation formula of the quadratic or cubic polynomial spectral smoothing formula is:

（17）

wherein,w=5,7,…,2m+1。

in actual measurement, the width of the track value is fittedwIt is not preferable to be greater than FWHM. Such asFruit fit trace value widthwToo wide, the smooth spectral line can lead the characteristic peak voltage to be low and flat, and lead the peak boundary to be lifted, so that the spectral data is distorted and a certain error is brought to the following peak searching; if the width of the fitting channel valuewToo narrow will reduce the link between the spectral data, the effect of removing statistical fluctuations is not obvious, and good smoothing effect is not achieved.

According to the formula (17), a five-point cubic polynomial smooth calculation formula can be obtained:

（18）

in order to effectively reduce the influence of noise and statistical fluctuation in the gamma energy spectrum data to be analyzed, repeated smoothing processing is performed on the gamma energy spectrum data to be analyzed for a plurality of times. However, as the number of smoothness increases, the CPU load increases; and causes distortion in the shape of the spectrum, annihilation of weak peaks, enlargement of peak widths, reduction of peak heights, and the like.

As shown in fig. 6, the resulting curve is smoother as the number of smoothness increases, but the peak position remains substantially unchanged. Therefore, the gamma energy spectrum data to be analyzed is processed smoothly, so that the noise and the interference of statistical fluctuation can be reduced, the difficulty of peak searching is reduced, the analysis of the energy spectrum is not greatly influenced, and the nuclide identification failure is caused.

However, the number of smoothness is greatly dependent on the intensity of the radiation source. As shown in FIG. 7, since K-40 has a low radioactivity, the smoothness effect is poor, and a significant fluctuation can be seen at 10 times of smoothness. When the number of smoothness reaches 20, a satisfactory smoothness is achieved, in sharp contrast to the satisfactory effect achieved by only 5 times of smoothness of Cs-137 in fig. 6.

Thus, the number of smoothness should be related to the intensity of the radiation source, i.e. the number of smoothness should be related to the gamma total count.

Fitting a straight line by a least square method, and calculating an insertion point to obtain the relation between the smooth times and the gamma total count:

Q=-0.00798x+24.44（19）

wherein,Qrepresenting repetitive lightThe number of times of the sliding is counted,xrepresenting the total gamma count per 1 second.

And calculating the gamma total count of the gamma energy spectrum data to be analyzed after removing the natural nuclide within 1 second, and calculating the smooth times of the gamma energy spectrum data to be analyzed according to the gamma total count of the gamma energy spectrum data to be analyzed after removing the natural nuclide.

S5, peak searching: and carrying out peak searching on the gamma energy spectrum data to be analyzed after the smoothing treatment by using a derivative peak searching method and a symmetrical zero area peak searching method to obtain a characteristic peak array.

The peak searching process of the nuclide identification algorithm can be regarded as searching a Gaussian peak meeting the condition (usually, the characteristic peak of the gamma energy spectrum meets Gaussian distribution) from the function point of view, and the position of the highest value of the Gaussian characteristic peak is the peak position of the characteristic peak.

In function analysis, as a data feature which can well reflect the change trend of a function, a derivative is a very important mathematical tool. The characteristic of the variation of the derivative of the gamma energy spectrum around the characteristic peak position can be used to determine the characteristic peak position. Aiming at the condition of single characteristic peak, the first derivative has better recognition efficiency; however, for more complex cases, the use of the first derivative alone is not effective in locating the characteristic peaks, and is required in combination with a symmetrical zero area peak finding method.

The invention combines the derivative peak searching method and the symmetrical zero area peak searching method, thereby improving the accuracy of nuclide identification.

First, derivative peak finding is used. According to the characteristics of the Gaussian function, the derivative of the gamma energy spectrum has a specific change near the peak area of the characteristic peak, particularly when passing through the peak top of the characteristic peak, and the characteristic peak can be found by utilizing the characteristic. Since the first derivative is greater than zero at the leading edge of the characteristic peak; at the trailing edge of the characteristic peak is less than zero; at the top of the characteristic peak, the first derivative is equal to zero. Therefore, when the first derivative is obviously changed from positive to negative gradually, the existence of the characteristic peak and taking the peak position of the characteristic peak as the track address with the first derivative equal to zero can be considered, which is the principle of the derivative peak finding method.

The first derivative formula adopted by the invention is as follows:

（20）

wherein,is the address of the trackiFirst derivative of>Is the address of the tracki-2Count of->Is the address of the tracki-1Count of->Is the address of the tracki+1Count of->Is the address of the tracki+2Counting.

The accuracy of the peak searching result has a great relation with the effect of spectral line smoothness, and when the gamma energy spectrum data to be analyzed is not subjected to spectral line smoothness, a characteristic peak can not be found almost; when the smooth effect of the spectral line is poor, deviation of peak searching results, such as missing peak identification, misplugging peak identification and the like, can be caused. As shown in fig. 8, when the smoothing effect is not ideal, the first derivative of the gamma-energy spectrum data to be analyzed is cluttered, and a characteristic peak cannot be found; as can be seen from fig. 9, when the smoothing effect of the gamma energy spectrum data to be analyzed is remarkable, the first derivative of the gamma energy spectrum data also has a process from positive number to zero and then to negative number, and the trace address with the zero derivative is the peak position of the characteristic peak.

When searching for the peak position of the characteristic peak by using the first derivative peak searching method, some limiting conditions are added to strengthen the accuracy of peak searching and eliminate the interference of radioactive fluctuation or other factors: first derivative two zero value distances from negative to positiveWShould be between 1 and 4 times FWHM. FWHM can be determined from instrument energy resolution, for a 1024-channel gamma spectrometer of NaI crystals, the FWHM typically takes a value of 18. By the limiting condition, the misjudgment of the characteristic peak when the smooth effect of the spectral line is poor can be eliminated. Features (e.g. a character)The peak-to-peak address count should be greater than 2 times the background value per second. By this limitation, the influence of the radioactive fluctuation and the weak peak can be eliminated.

As shown in fig. 10, the steps of the combined peak finding method are as follows:

s5-1, rapidly screening a first characteristic peak address by adopting a derivative peak searching method, and storing the first characteristic peak address into a first characteristic peak array until the last address of gamma energy spectrum data to be analyzed is traversed.

As shown in fig. 11, the derivative peak finding method specifically includes the following steps:

s5-1-1, calculating a first derivative spectrum in the gamma energy spectrum data to be analyzed after the smoothing treatment by using a formula (20) according to the track address from small to large.

S5-1-2, determining the channel address (namely peak position) of the characteristic peak by using a zero value method. On-road addressiIs positive and adjacent addressesi+1When the first derivative value is negative, the address is thatiSum of addressesi+1Track address with a characteristic peak in betweenxWherein, the method comprises the steps of, wherein,x=i+0.5. Traversing the whole first derivative spectrum to obtain an initial characteristic peak arrayX[x ₁ ,x ₂ …x _n ]；

S5-1-3, determining peak areas by using a zero value method. On-road addressjThe first derivative value of (2) is negative, adjacent addressesj+1The first derivative value of (2) is positive, then at the track addressjSum of addressesj+1With a characteristic peak boundary betweenzWherein, the method comprises the steps of, wherein,z=j+0.5the method comprises the steps of carrying out a first treatment on the surface of the Traversing the whole first derivative spectrum to obtain a characteristic peak boundary arrayZ[z ₁ ,z ₂ …z _n ]；

S5-1-4, judging the peak width. Taking initial characteristic peak arrayX[x ₁ ,x ₂ …x _n ]Is a track address of a characteristic peak of (a)x _i Array at characteristic peak boundaryZFind its nearest left and right boundariesz _m Andz _n the method comprises the steps of carrying out a first treatment on the surface of the If the characteristic peak width isWFailure to meet the condition FWHM<W<(4 XFWHM), descriptionx _i A trace address that is not a characteristic peak; if the condition is satisfied, go to step S5-1-5, in whichW=z _n -z _m . By the limiting condition, the misjudgment of the characteristic peak when the smooth effect of the spectral line is poor can be eliminated.

S5-1-5, judging the peak height. Track address of characteristic peakx _i Should be greater than 2 times the background, if this condition is not met,x _i a trace address that is not a characteristic peak; if the condition is satisfied, the track address of the characteristic peakx _i Is the first characteristic peak address and willx _i Storing the first characteristic peak arrayCHN[]Until the last track address of the gamma energy spectrum data to be analyzed is traversed. By this limitation, the influence of the radioactive fluctuation and the weak peak can be eliminated.

As shown in FIG. 12, the derivative peak finding method has better recognition efficiency, can easily recognize the characteristic peak of Cs-137, and has the ordinate of Cs-137 counted per second and the unit of s ^-1 。

S5-2, according to a symmetrical zero area peak finding method and a first characteristic peak arrayCHN[]And carrying out peak searching again on the smooth gamma energy spectrum data to be analyzed, and determining a second characteristic peak array.

The symmetrical zero area peak finding method refers to: the "window" function of zero area and the gamma energy spectrum data to be analyzed complete the convolution transformation, and the "window" function is a symmetric function, such as a gaussian function, a square wave function, etc. The convolution transform for the linear basis will be zero, with only where the peak is present being non-zero.

S5-2-1. For the first characteristic peak arrayCHN[]Each first characteristic peak addressx _i Setting a calculation interval to obtain a calculation interval of each first characteristic peak address, wherein the calculation interval is as follows:x _i -2×FWHM～x _i +2×FWHM。

s5-2-2 for each first characteristic peak address calculation intervalx _i -2×FWHM～x _i Convolving (i.e., an expression of a "window" function) each track address within +2 xFWHM:（21）

wherein,is the address of the trackiValues after convolution transformation +.>Is the address of the tracki+jIs used for counting the number of the counts,C _j as a function of the peak-like shape,，/>，2m+1in the form of a width, the width,σ=fwhm/2.355 is a peak width parameter;

s5-2-3, judging gamma energy spectrum data to be analyzed in the calculation interval of each first characteristic peak address, whenR _i >fWhen it willR _i Corresponding track addressiDetermining the second characteristic peak address and storing the second characteristic peak address into a second characteristic peak arrayR[]Traversing the calculation interval of each first characteristic peak address in the first characteristic peak array; when (when)R _i ≤fDescription of the embodimentsR _i The corresponding track address is not a characteristic peak, wherein,，/>is->Is set at the root mean square of (c),，fis a peak-finding threshold.

When (when)R _i Exceeding a predetermined peak finding thresholdfWhen a characteristic peak is found. According toR _i The peak position can be determined by the channel address corresponding to the positive electrode value of (a), and the peak boundary can be defined byR _i The distance between two adjacent minima on both sides of the positive peak of (c). By the limiting condition, weak peaks can be found on high background, or peaks close to the weak peaks can be found, and the peak-finding device has better high background inhibition capability and weak peak identification capability.

The symmetrical zero-area peak searching method has better high background inhibition capability and weak peak identification capability, but has larger calculated amount, and as shown in figure 13, the symmetrical zero-area peak searching method can identify the characteristic peak of Cs-137, the ordinate is the count per second of Cs-137, and the unit is s ^-1 。

Due to the influence of statistical fluctuation and the complexity of gamma energy spectrum, finding all characteristic peaks and precisely determining peak positions are important tasks in gamma energy spectrum analysis. However, it is difficult to find weak peaks on a high background, or to find peaks that are very close together. The first derivative peak searching method is superior in determining the boundary of the peak area, but the derivative peak searching method has insufficient resolution capability as a whole, fully considers the advantages of each peak searching method, and adopts the derivative peak searching method and the symmetrical zero-area peak searching method to perform combined peak searching so as to meet the requirements of high peak searching speed and high peak searching precision.

The derivative peak searching method and the symmetrical zero area peak searching method are combined to perform peak searching, and the method has the advantages of high peak searching speed, high peak searching precision, short peak searching time consumption, light load on a CPU (central processing unit) and the like, and is particularly suitable for a portable nuclide identifier.

As shown in FIG. 14 and FIG. 15, the peak searching capability test is carried out by using the artificial nuclides Am-241 and Co-60, and the characteristic peaks of each artificial nuclide can be easily positioned by combining the peak searching method, so that the performance is excellent.

S6, determining peak positions and nuclide types.

The combined peak searching method only accurately obtains the peak position array of the characteristic peakR[]This is not sufficient. The peak address of the characteristic peak is converted into the energy of the characteristic peak, so that the species of the gamma energy spectrum data to be analyzed can be determined according to the energy of the characteristic peak. The energy scale refers to converting peak addresses of characteristic peaks into characteristic peak energies.

As shown in table 2, the energy corresponding to the characteristic peak address can be determined by means of the energy scale. The energy scale formula is different for different detector crystals and different multi-channel analyzers. Table 2 shows the relationship between the peak energy and the channel address of each nuclide characteristic measured for the existing NaI crystal and the self-developed multi-channel analyzer.

TABLE 2 relation table of energy and address of characteristic peaks of various nuclides

According to the invention, a straight line is fitted through a least square method, and the interpolation point calculation is carried out, so that the relation between energy and the track address can be obtained, namely an energy scale formula:

T=3.05x-2.5（22）

wherein,xrepresenting the location of the track(s),Trepresenting energy (keV).

As shown in fig. 16, the characteristic peak energy is obtained from the characteristic peak address.

Because the resolution of the detector itself and the error generated by the energy scale formula, the energy corresponding to the peak address found by the combined peak finding method cannot be completely consistent with the energy of the characteristic peak in the nuclide library, so that a proper energy window width is required to be setW _E . If equation (23) is satisfied, it can be identified as the characteristic energyE _r And calculating the nuclides corresponding to the energies of all the characteristic peaks.

（23）

Wherein,E _r the standard energy of the characteristic peak of the nuclide (table 1),Ethe energy of the characteristic peak after the channel address conversion of the second characteristic peak,W _E is the window width.

That is, if peak positions obtained by the peak finding method are combinedR[]After conversion by the energy scale formula, the result isEAnd core(s)Characteristic peak energy of certain nuclide in element libraryE _r Phase difference is smaller than window widthW _E The nuclide is considered to be present.

Width of energy windowW _E Directly related to the success rate of nuclide identification. If the window width is wideW _E The value is too narrow, and the error generated by calculation of the energy scale and peak position finally leads to the energy used for nuclide identificationEAnd characteristic peak energyE _r Is greater than the window widthW _E Resulting in a reduction of the species identified and omission of some species.

If the window width is wideW _E If the value is too wide, the probability of occurrence of coherent nuclides (nuclides with similar or same energy characteristic peaks are called coherent nuclides) during nuclide identification is greatly increased, and the false identification rate of nuclide identification is improved.

The invention introduces the concept of confidence to evaluate the confidence level of nuclide identification, and the greater the confidence level is, the energy of nuclide identification is illustratedEThe closer the characteristic energy isE _r The method comprises the steps of carrying out a first treatment on the surface of the Conversely, the smaller the confidence, the energy of nuclide identification is accounted forEAnd characteristic energyE _r The phase difference is large.

The confidence coefficient is calculated according to the formula:

（24）

wherein alpha is the confidence coefficient,ΔEis the energy identified by the nuclideEAnd characteristic peak energy in nuclide libraryE _r A difference between; ETOL is the energy window range of a nuclear species.

S7, performance test

The invention utilizes radioactive reference sources Am-241, ba-133, co-60, cs-137, K-40, ra-226 and Th-232 to carry out nuclide identification performance test, and verifies the feasibility and accuracy of the invention.

S7-1. Identification test of single radionuclides.

As shown in table 3, first, positions of each radioactive reference source at which the surrounding dose equivalent rate of the background is subtracted and which is generated at the detection point is respectively 0.1 [ mu ] Sv/h, 0.3 [ mu ] Sv/h, 0.5 [ mu ] Sv/h and 0.7 [ mu ] Sv/h are determined. And secondly, adjusting the positions of the reference sources to change the intensity of each radioactive reference source, and testing the nuclide identification capability of the instrument under different radioactive intensities. The nuclide measurement and identification live time was 60 seconds, and for each radionuclide, 4 sets of tests were performed with different source intensities (i.e., ambient dose equivalent rates), 10 times per set of measurements, i.e., 40 times per nuclide test. The identification probability is calculated according to the ratio of the number of times that the radionuclide can be correctly identified by the nuclide identification algorithm to the total number of tests.

According to the invention, the radionuclide identification performance test is carried out by utilizing radioactive reference sources Am-241, ba-133, co-60, cs-137, K-40, ra-226 and Th-232, statistical results are obtained, the identification rate of a single radioactive reference source is 100%, the average confidence coefficient of each radioactive reference source is more than 97%, the excellent radionuclide identification effect is obtained, and the requirements of the calibration standard are met.

TABLE 3 Single radionuclide identification test

S7-2, identification test of mixed radionuclides.

The surrounding dose equivalent rate generated at the detector position for each radionuclide was measured to be 0.5 mu Sv/h above background by combining Am-241, ba-133, co-60, cs-137, ra-226 and Th-232 radionuclides.

According to the invention, spectrum analysis is carried out on the mixed nuclide, the identification capability of a nuclide identification algorithm is tested, the nuclide measurement and identification living time is controlled to be 60 seconds, 40 times of tests are carried out, and the identification probability is calculated according to the ratio of the number of times that the nuclide identification algorithm can correctly identify the group of radionuclides to the total number of times of the tests.

As shown in fig. 17, the right side shows the identified individual species, and the confidence of the respective species correspondence. The left side shows the serial numbers respectively: representing the current 8 th measurement data, and the total measurement times is 40; electric quantity: the current electric quantity of the measuring instrument is 72%; high pressure: the current working high-voltage value of the measuring instrument is 552.7V; time: the active time under the current measurement state is 60.0s, and the total duration is 61.6s; dose rate: indicating that the measured dose rate value is 2.56 mu Sv/h in the current environment state; track address: the number of addresses at the cursor positioning position is 1005, and the energy size is 3063keV; counting: the count of the track address where the cursor is located is 4.

As shown in Table 4, the present invention has excellent nuclide recognition effect on mixed radionuclides Am-241, ba-133, co-60, cs-137, ra-226 and Th-232, with a nuclide recognition rate of 100% and an average confidence of each nuclide of more than 96%.

TABLE 4 statistical table of mixed radionuclide identification capability test results

/>

Claims (10)

1. The nuclide identification method in the complex environment is characterized by comprising the following steps:

s1, background subtraction: background subtraction is carried out on gamma energy spectrum data to be analyzed by utilizing an SNIP algorithm;

s2, solving the content of natural nuclides: calculating the content value of the natural nuclide to be removed according to the inverse matrix of the response coefficient of the natural nuclide to be removed and the count of gamma energy spectrum data to be analyzed after background subtraction on each channel address, wherein the characteristic peak of the natural nuclide to be removed is a multipotent peak;

s3, removing natural nuclides: calculating the contribution value of the natural nuclide to be removed on each channel site according to the content value of the natural nuclide to be removed, subtracting the contribution value of the natural nuclide to be removed on each channel site from the count of the gamma energy spectrum data to be analyzed after background subtraction on each channel site, and obtaining the gamma energy spectrum data to be analyzed after the natural nuclide is removed;

s4, smooth spectral line: calculating the gamma total count of gamma energy spectrum data to be analyzed after natural nuclides are removed, calculating the smoothing times according to the gamma total count after the natural nuclides are removed and the relation between the smoothing times and the gamma total count, and carrying out smoothing treatment on the gamma energy spectrum data to be analyzed after the natural nuclides are removed based on the smoothing times;

s5, peak searching: carrying out peak searching on gamma energy spectrum data to be analyzed after smooth processing by utilizing a derivative peak searching method and a symmetrical zero area peak searching method to obtain a characteristic peak array;

s6, determining peak positions and nuclide types: converting the addresses of characteristic peaks in the characteristic peak array into characteristic peak energy, and determining the nuclide types according to the converted characteristic peak energy.

2. The method for identifying nuclides in a complex environment according to claim 1, wherein the specific way of calculating the inverse matrix of the response coefficient of the natural nuclides to be removed in step S2 is as follows:

s2a-1, calculating a response coefficient matrix of natural nuclides to be rejected based on the total count of known mixed nuclides on each track site and the content of each nuclide in the known mixed nuclides, wherein the known mixed nuclides comprise the natural nuclides to be rejected;

s2a-2, calculating a response coefficient inverse matrix of the natural nuclide to be removed according to the response coefficient matrix of the natural nuclide to be removed.

3. The method for identifying nuclides in a complex environment according to claim 2, wherein the specific way of calculating the content value of the natural nuclides to be removed in step S2 is as follows:

and multiplying the inverse matrix of the response coefficient of the natural nuclide to be removed by the count of the gamma energy spectrum data to be analyzed after background subtraction on each channel address, and calculating the content value of the natural nuclide to be removed.

4. The method for identifying nuclides in a complex environment according to claim 2, wherein the specific way of smoothing the data after removing natural nuclides is as follows:

and carrying out five-point three-time polynomial smooth calculation on each point in the gamma energy spectrum data to be analyzed according to the smooth times to obtain the gamma energy spectrum data to be analyzed after being smoothed.

5. The method for identifying nuclides in a complex environment according to claim 4, wherein the relationship between the number of smoothness and the total gamma count value is:

Q=-0.00798x+24.44

wherein,Qrepresenting the number of repetitions of the smoothing process,xrepresenting the total gamma count per 1 second.

6. The method for identifying nuclides in a complex environment according to claim 1, wherein the specific mode of peak searching for the smooth gamma-energy spectrum data to be analyzed in step S5 is as follows:

s5-1, carrying out peak searching on the smooth gamma energy spectrum data to be analyzed by adopting a derivative peak searching method to obtain a first characteristic peak address, and storing the first characteristic peak address into a first characteristic peak array until the last address of the gamma energy spectrum data to be analyzed is traversed;

s5-2, carrying out secondary peak searching on the gamma energy spectrum data to be analyzed after being smoothed according to a symmetrical zero-area peak searching method and the first characteristic peak array, and determining a second characteristic peak array.

7. The method for identifying nuclides in a complex environment according to claim 6, wherein the specific method for carrying out peak searching on the smooth gamma energy spectrum data to be analyzed in step S5-1 is as follows:

s5-1-1, calculating a first derivative value of each channel address in the gamma energy spectrum data to be analyzed after the smoothing treatment to obtain a first derivative spectrum;

s5-1-2. If the address isiThe first derivative value of (2) is positive, adjacent track addressi+1The first derivative value of (2) is negative, then at the track addressiSum of addressesi+1Track address with characteristic peak in betweenxWherein, the method comprises the steps of, wherein,x=i+0.5the method comprises the steps of carrying out a first treatment on the surface of the Traversing the whole first derivative spectrum to obtain an initial characteristic peak arrayX[x ₁ ,x ₂ …x _n ]；

S5-1-3. If the address isjThe first derivative value of (2) is negative, adjacent addressesj+1The first derivative value of (2) is positive, then at the track addressjSum of addressesj+1With a characteristic peak boundary betweenzWherein, the method comprises the steps of, wherein,z=j+0.5the method comprises the steps of carrying out a first treatment on the surface of the Traversing the whole first derivative spectrum to obtain a characteristic peak boundary arrayZ[z ₁ ,z ₂ …z _n ]；

S5-1-4. Taking initial characteristic peak arrayX[x ₁ ,x ₂ …x _n ]Is a track address of a characteristic peak of (a)x _i Array at characteristic peak boundaryZFind its nearest left and right boundariesz _m Andz _n the method comprises the steps of carrying out a first treatment on the surface of the If the characteristic peak width isWFailure to meet the condition FWHM<W<(4 XFWHM), thenx _i A trace address that is not a characteristic peak; if the condition FWHM is satisfied<W<(4 xfwhm), then go to the next step, where,W=z _n -z _m FWHM is full width at half maximum;

s5-1-5. If the characteristic peak is located at the trackx _i Less than or equal to 2 times the background per second countx _i A trace address that is not a characteristic peak; if the track address of the characteristic peakx _i Is greater than 2 times the background, then the trace address of the characteristic peakx _i Is the first characteristic peak address and willx _i And storing the first characteristic peak array.

8. The method for identifying nuclides in a complex environment according to claim 7, wherein the determining the second characteristic peak array in step S5-2 is performed in the following specific manner:

s5-2-1. For each first characteristic peak addressx _i Setting a calculation interval to obtain a calculation interval of each first characteristic peak address, wherein the calculation interval is as follows:x _i -2×FWHM～x _i +2×FWHM；

s5-2-2 for each first characteristic peak address calculation intervalx _i -2×FWHM～x _i +2 xfwhm) by performing a convolution transform for each address:wherein (1)>Is the address of the trackiValues after convolution transformation +.>Is the address of the tracki+jIs used for counting the number of the counts,C _j is peak-like function, ++>，/>；

S5-2-3, judging gamma energy spectrum data to be analyzed in the calculation interval of each first characteristic peak address, if the gamma energy spectrum data is positive, determining that the gamma energy spectrum data is positiveR _i >fWill thenR _i Corresponding track addressiDetermining the first characteristic peak address as a second characteristic peak address, storing the second characteristic peak address into a second characteristic peak array until the calculation interval of each first characteristic peak address in the first characteristic peak array is traversed, wherein,，/>is->Root mean square of>，fIs a peak-finding threshold.

9. The method for identifying nuclides in a complex environment according to claim 1, wherein the specific way of converting the addresses of the characteristic peaks in the characteristic peak array into the energy of the characteristic peaks in step S6 is as follows:

converting the track address of the characteristic peak into a special according to an energy scale formulaThe energy of the characteristic peak and the energy scale formula are as followsT=3.05x- 2.5Wherein, the method comprises the steps of, wherein,xrepresenting the location of the track(s),Trepresenting energy.

10. The method for identifying nuclides in a complex environment according to claim 1, wherein the specific manner of determining the nuclide type according to the converted characteristic peak energy in step S6 is as follows:

selecting any nuclide, calculating the absolute value of the difference between the energy of the converted characteristic peak and the standard energy of the characteristic peak of the selected nuclide, and if the calculated absolute value is smaller than or equal to the preset energy window width, the gamma energy spectrum data to be analyzed comprise the selected nuclide.