Food nondestructive testing method and system based on hyperspectral imaging

Abstract

The invention provides a food nondestructive testing method and system based on hyperspectral imaging, and relates to the technical field of food detection.Position information of a leaf wrinkle area is determined by analyzing and judging hyperspectral image data, and a leaf vegetable surface pesticide residue distribution diagram is calculated according to a pesticide characteristic absorption waveband range; according to the method, the technical problem of low detection precision caused by diffuse reflection interference caused by leaf wrinkles in the prior art is solved, so that whether pesticide residues exist or not is judged in an auxiliary manner, the situation that vegetables with the pesticide residues flow into the market is reduced, and food poisoning events caused by the pesticide residues are further reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of food detection, and more specifically, to a non-destructive food detection method and system based on hyperspectral imaging. Background Art

[0002] As an advanced technology integrating spectral analysis and image processing, hyperspectral imaging technology has received extensive attention and applications in recent years in fields such as agriculture, food detection, and environmental monitoring. Compared with traditional spectral analysis methods, hyperspectral imaging can not only provide rich spectral information but also combine spatial information to form a three-dimensional hyperspectral data cube with spatial and spectral dimensions. Hyperspectral imaging technology has significant advantages in food quality and safety detection, especially in the non-destructive detection of pesticide residues in fruits and vegetables. Through hyperspectral imaging technology, not only can the types and distributions of food surface contaminants be identified, but also the changes in its internal components can be analyzed, avoiding the destructiveness and time-consuming nature of traditional chemical detection. However, in the detection of foods with complex surface structures such as leafy vegetables, due to the natural folds of the leaves, their surface reflection characteristics are easily affected by factors such as the illumination angle and surface curvature, resulting in deviations in the spectral reflectance in the hyperspectral data. This deviation will directly affect the extraction accuracy of spectral features and further reduce the reliability of the detection results. Therefore, how to effectively eliminate the interference of leaf folds on hyperspectral data has always been one of the research difficulties in the related technical fields.

[0003] Currently, some studies have attempted to reduce the impact of folds on hyperspectral data through methods such as global mean compensation and spectral correction based on empirical models. However, these methods usually rely on specific experimental conditions or assumptions and are difficult to be widely promoted in practical applications. For example, global mean compensation ignores the local differences in spectral characteristics between the folded area and the non-folded area, easily causing excessive or insufficient spectral correction; while the spectral correction method based on empirical models requires a large amount of sample data to be obtained in advance, and the applicability of the model is often limited by the types of leaves and environmental conditions. In addition, when dealing with hyperspectral data, the existing technology usually does not accurately mark and analyze the spatial domain of the folded area, but directly processes the overall spectral data, which will lead to inaccuracies in the spectral feature extraction results, especially in application scenarios such as pesticide residue detection that highly rely on spectral accuracy. It can be seen that there is still much room for improvement in the spatial domain analysis and spectral domain compensation of hyperspectral data in the current technology. Summary of the Invention

[0004] To solve the above technical problems, the present invention provides a non-destructive food detection method and system based on hyperspectral imaging, which can solve the problem of rapid non-destructive detection of organophosphorus pesticide residues on the surface of leafy vegetables (such as lettuce, leeks, etc.) to a certain extent, especially solve the technical problem of low detection accuracy caused by diffuse reflection interference due to leaf folds in the prior art.

[0005] According to one aspect of the present invention, there is provided a non-destructive food detection method based on hyperspectral imaging, which includes:

[0006] Obtain hyperspectral image data of the leafy vegetable to be measured, and obtain a three-dimensional hyperspectral data cube including a spatial dimension and a spectral dimension;

[0007] Identify the folded area of the three-dimensional hyperspectral data cube, use spatial domain gray gradient analysis to determine the position information of the leaf fold area, and obtain the position marking data of the leaf fold area;

[0008] According to the position marking data, perform adaptive spectral compensation on the folded area in the three-dimensional hyperspectral data cube, establish a compensation factor by calculating the spectral reflectance difference between adjacent non-folded areas in the folded area, and obtain the compensated spectral data after eliminating the influence of folds;

[0009] Based on the compensated spectral data, within the characteristic absorption band range of organophosphorus pesticides, extract the second derivative characteristics of the spectral curve, and calculate the distribution map of organophosphorus pesticide residues on the surface of the leafy vegetable to be measured.

[0010] Further, identifying the folded area of the three-dimensional hyperspectral data cube includes:

[0011] Extract a gray-scale image with a fixed wavelength from the three-dimensional hyperspectral data cube as a reference image, and perform noise reduction processing;

[0012] Calculate the gray-scale gradient amplitude of the noise-reduced gray-scale image, and mark the pixel points with a gray-scale gradient greater than the gray-scale gradient threshold as candidate fold edge points;

[0013] Perform morphological processing on the candidate fold edge points to obtain an initial folded area;

[0014] Optimize the initial folded area to obtain the position marking data of each folded area.

[0015] Further, performing morphological processing on the candidate fold edge points includes multiple morphological processings, determining whether to perform dilation operations according to adjacent connected regions, and finally performing boundary smoothing processing to obtain the boundary of the folded area after morphological processing.

[0016] Further, based on the morphological characteristics of the blade folds, an adaptive elliptical structural element is constructed. The dilation operation is performed by translating the rotated elliptical structural element point by point in the image, only processing the internal area of the image, and recording the pixel position information after dilation. Finally, a position index table of the fold area is established to accurately mark the fold area.

[0017] Further, the compensation factor determines whether to use the distance weighting method or the average spectrum of the global reference to calculate the compensation factor according to the Euclidean distance between each pixel point in the fold area and the sampling point of the nearest reference area.

[0018] Further, the distance weighting method includes:

[0019] Measure the Euclidean distance between the pixel point to be processed in the fold area and the reference sampling point, select multiple reference points with the smallest distance, and calculate the weight coefficient according to their distances;

[0020] Combine the spectral reflectance difference and the weight coefficient, calculate the weighted spectral difference and sum it up to obtain the spectral compensation factor of the pixel point to be processed;

[0021] If the pixel point is located in the fold edge area, calculate the attenuation coefficient based on its shortest distance to the fold boundary, and multiply it by the spectral compensation factor to obtain the final compensation value;

[0022] Repeat the above process to calculate the spectral compensation factors at each wavelength position.

[0023] Further, calculating the compensation factor through the average spectrum of the global reference includes:

[0024] Extract the annular reference area outside the fold, divide it into upper and lower sub-areas, and calculate the spectral mean value of each wavelength;

[0025] When the spectral difference between the two sub-areas is less than the preset value, take the average value as the reference spectrum, otherwise select the spectral of the sub-area closer to the pixel point to be processed as the reference;

[0026] Calculate the difference between the pixel to be processed and the reference spectrum, perform clipping processing on the over-limit part to generate the initial compensation factor, and calculate the position weight according to the distance between the pixel and the fold center line, and multiply the two to obtain the final compensation factor. Repeat the process to generate the compensation sequence of each wavelength.

[0027] Further, extracting the second derivative feature of the spectral curve includes:

[0028] Preprocess the compensated spectral data;

[0029] Based on the smoothed spectral data, use the five-point derivative formula to calculate the first derivative, and then use the same method to calculate the second derivative for the first derivative curve.

[0030] Further, the five-point derivative formula calculates the first derivative by selecting two adjacent wavelength points on both sides of the current wavelength point to form a five-point derivative window, and calculating the first derivative of the center point using fixed coefficients; while the three-point derivative formula is used for correction at both ends of the spectral curve; finally, the complete first derivative curve is generated by translating the derivative window point by point.

[0031] According to another aspect of the present invention, there is provided a non-destructive food detection system based on hyperspectral imaging, which includes:

[0032] An acquisition module, configured to acquire hyperspectral image data of the to-be-detected leafy vegetables, and obtain a three-dimensional hyperspectral data cube including a spatial dimension and a spectral dimension;

[0033] A position marking module, configured to identify the wrinkled area of the three-dimensional hyperspectral data cube, determine the position information of the leaf wrinkled area by using spatial domain gray scale gradient analysis, and obtain the position marking data of the leaf wrinkled area;

[0034] A compensation module, configured to adaptively perform spectral compensation on the wrinkled area in the three-dimensional hyperspectral data cube according to the position marking data, establish a compensation factor by calculating the spectral reflectance difference between the adjacent non-wrinkled areas of the wrinkled area, and obtain the compensated spectral data after eliminating the influence of wrinkles;

[0035] An extraction module, configured to extract the second derivative features of the spectral curve within the characteristic absorption band range of the organophosphorus pesticide based on the compensated spectral data, and calculate and obtain the distribution map of the organophosphorus pesticide residue concentration on the surface of the to-be-detected leafy vegetables.

[0036] Compared with the prior art, the non-destructive food detection method and system based on hyperspectral imaging provided by the present invention determine the position information of the leaf wrinkled area by analyzing and judging the hyperspectral image data, and then calculate the pesticide residue distribution map on the surface of the leafy vegetables through the characteristic absorption band range of the pesticide. In this way, the technical problem of low detection accuracy caused by the diffuse reflection interference due to leaf wrinkles in the prior art can be solved, so as to help judge whether there is pesticide residue, reduce the inflow of vegetables with pesticide residue into the market, and further reduce the food poisoning incidents caused by pesticide residue. Description of the Drawings

[0037] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the following drawings are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts. In the drawings:

[0038] Figure 1Flow chart of a non-destructive food detection system based on hyperspectral imaging according to an embodiment of the present invention. Detailed implementation manners

[0039] Hereinafter, exemplary embodiments of the present invention will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments of the present invention. It should be understood that the present invention is not limited by the exemplary embodiments described herein.

[0040] Figure 1 Flow chart of a non-destructive food detection method based on hyperspectral imaging according to an embodiment of the present invention. As Figure 1 shown, in the non-destructive food detection method based on hyperspectral imaging, it includes:

[0041] S1: Obtain hyperspectral image data of the to-be-detected leafy vegetables. Scan the to-be-detected leafy vegetables in a direction perpendicular to the leaf surface through the hyperspectral imaging device to obtain a three-dimensional hyperspectral data cube including a spatial dimension and a spectral dimension.

[0042] Obtaining the hyperspectral image data of the to-be-detected leafy vegetables specifically includes: placing the to-be-detected leafy vegetables on a sample stage with a black non-reflective background, and the sample stage is arranged parallel to the horizontal plane; controlling the light source system of the hyperspectral imaging device to emit a continuous spectrum with a wavelength range of 250-1000 nm, and the light source system is located at a position 60°±5° directly above the to-be-detected leafy vegetables; adjusting the spectral camera of the hyperspectral imaging device so that the optical axis of the spectral camera is perpendicular to the surface of the to-be-detected leafy vegetables, and the height of the spectral camera from the surface of the to-be-detected leafy vegetables is 50 cm±2 cm; setting the scanning rate of the hyperspectral imaging device to 30 frames per second, and under the illumination of the light source system, performing line-by-line push-broom scanning on the to-be-detected leafy vegetables through the spectral camera to obtain a three-dimensional hyperspectral data cube with a spatial resolution of 1024×1024 pixels and a spectral resolution of 2.5 nm, wherein the three-dimensional hyperspectral data cube includes two spatial dimension information and one spectral dimension information, the spatial dimension information characterizes the two-dimensional spatial distribution characteristics of the surface of the to-be-detected leafy vegetables, and the spectral dimension information characterizes the spectral reflection characteristics of each point on the surface of the to-be-detected leafy vegetables.

[0043] S2: Identify the wrinkled area of the three-dimensional hyperspectral data cube, and use spatial domain gray gradient analysis to determine the position information of the leaf wrinkled area to obtain the position marking data of the leaf wrinkled area.

[0044] Identify the wrinkled areas in the three-dimensional hyperspectral data cube, which specifically includes: extracting the grayscale image at a wavelength of 550 nm from the three-dimensional hyperspectral data cube as the reference image. The reference image corresponds to the reflection peak of chlorophyll and is used to highlight the leaf wrinkling characteristics; performing Gaussian filtering on the reference image with a 5×5 pixel window to obtain the denoised filtered image; calculating the horizontal grayscale gradient value and the vertical grayscale gradient value of each pixel point in the filtered image, and obtaining the grayscale gradient amplitude of the pixel point through the square sum and square root operations on the horizontal grayscale gradient value and the vertical grayscale gradient value; setting the grayscale gradient amplitude threshold to 1.5 times the average value of the grayscale gradient amplitude of the current image, and marking the pixel points greater than the grayscale gradient amplitude threshold as candidate wrinkled edge points; performing morphological processing on the candidate wrinkled edge points to obtain the initial wrinkled area; calculating the grayscale mean and standard deviation of the pixel points in the initial wrinkled area, and removing the pixel points whose grayscale values deviate from the grayscale mean by more than 2 times the standard deviation to obtain the optimized wrinkled area; recording the pixel position information and the corresponding area information of the optimized wrinkled area as the position marking data of the leaf wrinkled area, where the position marking data includes the starting pixel coordinates, the ending pixel coordinates, and the area size of each wrinkled area.

[0045] Perform morphological processing on the candidate wrinkled edge points, which specifically includes: performing a closing operation on the candidate wrinkled edge points with a circular structuring element with a radius of 3 pixels to obtain the result of the first morphological processing, where the closing operation includes performing a dilation operation first and then an erosion operation; when there are connected regions with an area less than 30 pixels in the result of the first morphological processing, deleting these small-area connected regions to obtain the result of the second morphological processing; extracting the contour of each connected region in the result of the second morphological processing and calculating the ratio of the perimeter to the area of the contour. When the ratio is greater than 0.8, it is determined as a noise point and deleted to obtain the result of the third morphological processing; if the distance between adjacent connected regions in the result of the third morphological processing is less than 5 pixels, perform a dilation operation to connect these adjacent connected regions into a complete wrinkled area; when the ratio of the major axis to the minor axis of the connected wrinkled area is greater than 5, divide the area along the major axis direction, and set the division point at the midpoint position of the major axis; perform boundary smoothing processing on each divided wrinkled area, and perform an opening operation with a circular structuring element with a radius of 2 pixels, where the opening operation includes performing an erosion operation first and then a dilation operation, and finally obtain the boundary of the wrinkled area after morphological processing.

[0046] In view of the characteristics of the wrinkled areas on the surface of leafy vegetables, the present invention specifically improves the dilation operation, which specifically includes:

[0047] Based on the morphological characteristics of leaf wrinkles, an adaptive elliptical structural element is constructed. The long axis direction of the elliptical structural element is parallel to the local wrinkle trend. The length of the long axis of the elliptical structural element is 5 pixels, and the length of the short axis is 3 pixels. Select the current pixel to be processed in the original image, and calculate the gray gradient direction within the 9×9 pixel neighborhood window centered on this pixel. The gray gradient direction is used as the local wrinkle trend. After the gray gradient direction is determined, rotate the elliptical structural element to be parallel to the gray gradient direction to obtain the rotated structural element. When the pixel to be processed is located at the edge of the wrinkle area, based on the gray mean and variance of the neighborhood of this pixel, dynamically adjust the ratio of the long axis to the short axis of the elliptical structural element. When the gray variance is greater than the preset threshold, increase the ratio of the long axis to the short axis to 2:1. When the gray variance is less than the preset threshold, reduce the ratio of the long axis to the short axis to 3:2. Perform a dilation operation on the rotated structural element. When the center point of the structural element coincides with any pixel point in the neighborhood window, set the values of all pixel points within the coverage area of this structural element to 1. Translate the structural element point by point along the row direction and column direction of the original image, and repeat the above operations for each pixel point. When the neighborhood window intersects with the image boundary, only perform dilation processing on the internal area of the image. Record the pixel position information after dilation and establish a position index table for the wrinkle area.

[0048] Preferably, the existing dilation operation cannot adapt to the diverse characteristics of leaf wrinkles. The present invention adopts an adaptive elliptical structural element based on the local gradient direction. The direction of the structural element is aligned with the wrinkle trend in real time, and the ratio of the long axis to the short axis is dynamically adjusted according to the local gray characteristics, significantly improving the detection accuracy of the wrinkle edge. In particular, the recognition accuracy of irregular wrinkles and slight wrinkles is increased by more than 30%. And by calculating the local gradient direction and gray statistical characteristics at one time, directly guiding the shape adjustment of the structural element, avoiding the repeated feature extraction and parameter optimization process, while ensuring the detection accuracy, the calculation efficiency is increased by about 50%.

[0049] More specifically, the determination of the preset threshold is described in detail: Based on the statistical analysis of 500 samples of leafy vegetables, an adaptive threshold determination method is adopted, which specifically includes: obtaining the gray variance distribution within a 9×9 pixel neighborhood window for each sample, and sorting the gray variance distribution according to the numerical size; statistically analyzing the gray variance distribution characteristics of the wrinkled areas and non-wrinkled areas in the 500 samples, calculating the average gray variance of the wrinkled area and recording it as M1, and the average gray variance of the non-wrinkled area and recording it as M2; when the overall brightness level of the image is within the normal exposure range, setting the preset threshold T to (M1 + M2) / 2; if the image shows underexposure, that is, the average value of the overall gray value of the image is less than 128, adjusting the preset threshold T to 0.8×(M1 + M2) / 2; when the image shows overexposure, that is, the average value of the overall gray value of the image is greater than 200, adjusting the preset threshold T to 1.2×(M1 + M2) / 2; for each sample to be detected, automatically select the corresponding preset threshold based on its overall brightness level.

[0050] On the other hand, for the characteristics of the wrinkled areas on the surface of leafy vegetables, an erosion operation is performed on the wrinkled edge candidate points, which specifically includes:

[0051] Based on the morphological characteristics of the wrinkled edges, an adaptive-scale erosion structural element is constructed. When processing the main body area of the wrinkles, a circular structural element with a radius of 3 pixels is used. When processing the details of the wrinkled edges, the structural element is reduced to a radius of 2 pixels; in the binary image to be eroded, select the current pixel to be processed, and establish a 7×7 pixel detection window centered on this pixel; when the structural element is completely located within the area with a pixel value of 1, keep the value of the current central pixel as 1. If the structural element is partially or completely located within the area with a pixel value of 0, then change the value of the current central pixel to 0; for the pixel points at the edge of the wrinkled area, calculate the pixel value distribution density within their neighborhood. When the density is greater than 80%, a large-scale structural element is used. When the density is less than 80%, a small-scale structural element is used to achieve adaptive erosion; move the structural element pixel by pixel along the row direction and column direction of the binary image, and repeat the above erosion operation; if a partial area of the detection window exceeds the image boundary, treat the pixel values of the exceeded part as 0 values; when the erosion processing of all pixel points is completed, count the positions of the pixel points with a remaining pixel value of 1, and establish a new contour index table for the wrinkled area.

[0052] It should be noted that although the above description outlines the specific steps for identifying wrinkled areas in a three-dimensional hyperspectral data cube, in practical applications, appropriate adjustments may be required according to different varieties of leafy vegetables and their wrinkled characteristics. For example:

[0053] In terms of wavelength selection:

[0054] For dark green leafy vegetables (such as spinach, lettuce etc.), the reflection peak at 550 nm wavelength is the most obvious, and the wrinkle features are easy to identify; while for light green leafy vegetables (such as lettuce, cabbage etc.), it may be necessary to select images within the wavelength range of 500 - 520 nm as the reference to obtain better wrinkle detection effect.

[0055] In terms of filtering processing:

[0056] When processing fresh leafy vegetables, Gaussian filtering with a 5×5 pixel window is sufficient to remove noise; but for leafy vegetables with a longer storage time, due to possible slight shrinkage on the surface, it may be necessary to expand the filtering window to 7×7 pixels to better distinguish natural wrinkles and shrinkage caused by storage.

[0057] In terms of morphological processing:

[0058] For large - leaf vegetables (such as celery leaves, Chinese cabbage leaves etc.), the wrinkled areas are usually large and continuous. At this time, a larger structuring element (such as 5×5 pixels) can be used for morphological processing; while for small - leaf vegetables (such as coriander, shepherd's purse etc.), smaller structuring elements (such as 3×3 pixels) are needed to retain fine wrinkle features.

[0059] In terms of threshold setting:

[0060] Leafy vegetables picked in spring have more regular leaf wrinkles due to better growing environments. At this time, the threshold of the gray - level gradient amplitude can be set to 1.5 times the mean value; while leafy vegetables picked in winter may have less regular wrinkles due to changes in the growing environment, and it may be necessary to reduce the threshold to 1.3 times the mean value to ensure the integrity of detection.

[0061] In terms of region connection:

[0062] For tightly packed leafy vegetables such as cabbage, the threshold for judging the distance between adjacent wrinkled areas can be set smaller (such as 3 pixels); while for loosely packed leafy vegetables such as loose - leaf lettuce, the distance threshold can be appropriately increased (such as 7 pixels) to avoid over - connection.

[0063] These specific examples illustrate that algorithm parameters need to be flexibly adjusted according to the actual application scenarios to adapt to the characteristic differences of different types of leafy vegetables. No matter what parameter configuration is adopted, the ultimate goal is to accurately identify the wrinkled areas on the leaf surface and provide reliable spatial position information for subsequent spectral analysis and quality evaluation. In practical applications, a parameter template library for different leafy vegetable varieties can be established to achieve fast and accurate identification of wrinkled areas.

[0064] S3: According to the position marking data, perform adaptive spectral compensation on the wrinkled areas in the three-dimensional hyperspectral data cube. By calculating the spectral reflectance difference between the adjacent non-wrinkled areas of the wrinkled areas, establish a compensation factor to obtain the compensated spectral data after eliminating the influence of wrinkles.

[0065] Performing adaptive spectral compensation on the wrinkled areas in the three-dimensional hyperspectral data cube according to the position marking data specifically includes:

[0066] For each marked wrinkled area, determine its annular non-wrinkled area with a 5-pixel width on the periphery as the reference area; within the reference area, evenly select 12 sampling points along the wrinkle edge, and extract complete hyperspectral data for each sampling point; calculate the mean spectral reflectance of the 12 sampling points at each wavelength as the standard spectrum of the normal leaf tissue around the wrinkled area; for each pixel point in the wrinkled area, determine the Euclidean distance between it and the nearest reference area sampling point; when the Euclidean distance is less than the preset threshold, calculate the spectral compensation factor of the pixel point in a distance-weighted manner, and the weight is inversely proportional to the distance; if the Euclidean distance is greater than the preset threshold, then calculate the compensation factor using the average spectrum of the global reference area; when the spectral reflectance difference is greater than 30% in the visible light band (400 - 700 nm), correct the compensation factor using piecewise linear interpolation to avoid overcompensation; for the near-infrared band (700 - 1000 nm), limit the compensation factor within the range of ±20% based on the optical properties of the leaf tissue; according to the corrected compensation factor, perform spectral correction on each pixel point in the wrinkled area at all wavelengths to obtain the compensated spectral data; to ensure the continuity of the spectrum, use Gaussian smoothing transition at the edge of the wrinkled area, and the size of the smoothing window is 3 pixels.

[0067] Among them, calculating the spectral compensation factor of the pixel point to be processed in a distance-weighted manner specifically includes:

[0068] Measure the Euclidean distance from the pixel point to be processed in the wrinkled area to all the sampling points in the reference area. When the distance is less than the threshold of 12 pixels, select the 4 reference sampling points with the smallest distance; after selecting the 4 nearest reference sampling points, calculate the weight coefficient of each reference sampling point respectively. The weight coefficient is equal to the reciprocal of the distance of this reference point divided by the sum of the reciprocals of the distances of the 4 reference points; for each reference sampling point, calculate the spectral reflectance difference between it and the pixel point to be processed at the same wavelength position. The spectral reflectance difference is equal to the spectral reflectance of the reference sampling point minus the spectral reflectance of the pixel point to be processed; multiply the spectral reflectance difference obtained for each reference sampling point by its corresponding weight coefficient to obtain the weighted spectral difference corresponding to this reference point; sum up the weighted spectral differences of the 4 reference sampling points to obtain the spectral compensation factor of the pixel point to be processed at this wavelength; when the pixel point to be processed is located in the wrinkled edge area, measure the shortest distance from this pixel point to the wrinkled boundary, and divide this distance by the preset transition zone width to obtain the attenuation coefficient; multiply the attenuation coefficient by the aforementioned spectral compensation factor to obtain the final spectral compensation factor; repeat the above steps to calculate the spectral compensation factor of the pixel point to be processed at each wavelength position.

[0069] On the other hand, the compensation factor is calculated using the average spectrum of the global reference area, which specifically includes: extracting the spectral data of all pixel points in the annular reference area with a width of 5 pixels outside the wrinkled area; for the annular reference area, divide it into upper and lower sub-areas along the wrinkled direction; calculate the average spectral reflectance of all pixel points in the upper and lower sub-areas at each wavelength position respectively to obtain two reference spectra; when the reflectance difference between the reference spectra of the upper and lower sub-areas at the same wavelength position is less than 5%, take the average of the spectral reflectances of the two sub-areas as the standard reference spectrum at this wavelength; if the reflectance difference between the reference spectra of the upper and lower sub-areas at a certain wavelength position is greater than 5%, select the spectral reflectance of the sub-area closer to the pixel point to be processed as the standard reference spectrum at this wavelength; for the pixel point to be processed in the wrinkled area, calculate the spectral reflectance difference between it and the standard reference spectrum at each wavelength position; when the spectral reflectance difference is greater than the preset threshold, use a piecewise linear function to limit the difference to avoid over-compensation; take the limited spectral reflectance difference as the initial compensation factor at this wavelength position; based on the relative position of the pixel point to be processed and the wrinkled center line, calculate the position weight coefficient, which is maximum at the wrinkled center line and decreases from 1 to 0.2 towards the edge; multiply the position weight coefficient by the initial compensation factor to obtain the final compensation factor at this wavelength position; repeat the above steps to obtain the compensation factor sequence of the pixel point to be processed at all wavelength positions.

[0070] It should be noted that the method for determining the 5% reflectance difference threshold includes: First, analyze the hyperspectral data of a large number of healthy vegetation leaves, and statistically analyze the spectral reflectance variation range of different parts of the same leaf under natural conditions; when measuring different regions of the same leaf, in the visible light band (400 - 700 nm), the normal spectral reflectance natural variation is usually between 3% - 4%; in the near-infrared band (700 - 1000 nm), due to the differences in the internal structure of the leaves, the normal variation can reach 4% - 6%; comprehensively considering the natural variation characteristics of different bands, select 5% slightly higher than the upper limit of natural variation as the reference threshold; when dealing with vegetation of different species or different growth periods, the threshold can be appropriately adjusted according to the pre-experiment results, and the specific adjustment range is between 4% - 7%; if large natural variations are observed in certain characteristic bands, higher thresholds can be set separately for these bands, with a maximum of no more than 10%; when significant changes occur in environmental conditions (such as light intensity, measurement angle, etc.), the threshold needs to be recalibrated.

[0071] S4: Based on the compensated spectral data, within the characteristic absorption band range of organophosphorus pesticides (290 - 400 nm), extract the second derivative characteristics of the spectral curve, and calculate the distribution map of the organophosphorus pesticide residue concentration on the surface of the to-be-detected leafy vegetables.

[0072] Based on the compensated spectral data, extract features and calculate the pesticide residue concentration distribution within the characteristic absorption band range of organophosphorus pesticides, specifically including:

[0073] Preprocess the compensated spectral data. Use the Savitzky-Golay filtering algorithm to smooth the spectrum in the wavelength range of 290 - 400 nm, set the window width to 9 wavelength points, and the polynomial order to 3; based on the smoothed spectral data, calculate the first derivative using the five-point derivative formula, and then calculate the second derivative of the first derivative curve using the same method to obtain the second derivative characteristics of the spectral curve; when obvious troughs appear at the three characteristic wavelength points of 315 nm, 340 nm, and 375 nm on the second derivative curve, it indicates the presence of organophosphorus pesticide residues; calculate the trough depths of the second derivative at these three characteristic wavelength points, and convert the trough depths into pesticide residue concentration values according to the pre-established concentration calibration curve; if the second derivative value at a certain wavelength point is lower than the detection limit threshold, then record the concentration value at this point as zero; for each detected pixel point, calculate the weighted average based on the concentration values at the three characteristic wavelength points, and the weight coefficients are determined according to the detection sensitivities of each wavelength point, which are 0.4, 0.35, and 0.25 respectively; establish the correspondence between the image resolution and the actual spatial scale, and normalize the concentration value of each pixel point to the unit of mg / kg; use the bilinear interpolation method to smoothly transition the concentration values between adjacent pixel points to generate a continuous concentration distribution map; when the concentration difference between adjacent pixel points is greater than 1.5 mg / kg, add sampling points in the transition area to improve the interpolation accuracy; perform false color enhancement display on the generated concentration distribution map, establish the correspondence between the concentration value and the color, display the high-concentration area as red, and the low-concentration area as blue.

[0074] The method for calculating the first derivative of the spectral curve using the five-point derivative formula is as follows:

[0075] First, select two adjacent wavelength points on both sides of the wavelength point to be differentiated, forming a differentiation window containing five equally spaced wavelength points. Denote the current wavelength point as the center point, the two points on its left as the first left point and the second left point, and the two points on its right as the first right point and the second right point. Multiply the spectral reflectance of the second left point by the coefficient -1, multiply the spectral reflectance of the first left point by the coefficient -8, multiply the spectral reflectance of the first right point by the coefficient 8, and multiply the spectral reflectance of the second right point by the coefficient -1. Divide the algebraic sum of the above four products by twelve times the wavelength interval to obtain the first derivative value at the center point. When dealing with the wavelength points at both ends of the spectral curve, use the three-point differentiation formula for calculation. For the starting endpoint, multiply the spectral reflectance of the current point by the coefficient -3, the first right point by the coefficient 4, and the second right point by the coefficient -1, and divide the algebraic sum of the products by twice the wavelength interval. For the ending point, multiply the spectral reflectance of the current point by the coefficient 3, the first left point by the coefficient -4, and the second left point by the coefficient 1, and divide the algebraic sum of the products by twice the wavelength interval. Translate the differentiation window point by point and perform the same differentiation operation on each wavelength point to obtain a complete first derivative curve. If the intervals between wavelength points are unequal, the calculation coefficients need to be corrected according to the actual interval sizes to ensure the accuracy of the derivative calculation. Check the obtained first derivative curve to ensure that the continuity and smoothness of the curve meet the requirements of subsequent analysis. More specifically, the process of calculating the first derivative of the spectral curve using the five-point differentiation formula can be expressed by the following formula:

[0076]

[0077] where \(R'(\lambda)\) is the first derivative value of the spectral reflectance at the wavelength point \(\lambda\), \(R(\lambda)\) is the spectral reflectance value at the wavelength point \(\lambda\), \(\Delta\) is the interval between adjacent wavelength points, \(\lambda\) is the current wavelength point to be differentiated, and \(R(\lambda\pm\Delta)\) and \(R(\lambda\pm2\Delta)\) represent the spectral reflectance values at one interval and two intervals away from the current wavelength point, respectively.

[0078] The trough depths of the second-order derivative at three characteristic wavelength points (315nm, 340nm and 375nm) were calculated, specifically including: taking 5 wavelength points on the left and right of the target wavelength point to form a calculation window of 11 wavelength points; performing parabolic fitting on the second-order derivative data in the calculation window to obtain a quadratic function describing the shape of the local trough; when the goodness of fit R2 is greater than 0.95, the minimum value of the quadratic function is used as the trough depth; if the goodness of fit R2 is less than 0.95, the calculation window is expanded to 15 wavelength points and refitted; when determining the trough position, first determine whether the target wavelength point is located at the local minimum value; when the target wavelength point is at the local minimum value, When the second-order derivative value is greater than that of its left and right adjacent points, the calculation window is offset toward the trough until the true local minimum point is found; when calculating the trough depth, the average value of the second-order derivatives of the two end points of the calculation window is used as the baseline, and the vertical distance from the lowest point of the trough to the baseline is measured as the trough depth value; if multiple local minima appear in the calculation window, the trough closest to the target wavelength point and with the greatest depth is selected as the characteristic trough; when the trough is asymmetric, different weight coefficients are used on both sides of the calculation window, and data points close to the center of the trough are given a larger weight; an effective detection threshold for the trough depth is set, and when the calculated trough depth is less than the threshold, the trough depth of the feature point is recorded as zero. In more detail, the quadratic function describing the shape of the local trough can be expressed as follows:

[0079]

[0080] in:

[0081]

[0082] Where D(λ) is the calculated value of the trough depth, λ is the wavelength value, and λ min is the wavelength value at the trough position, λ l and λ r are the wavelength values of the left and right endpoints of the calculation window, f(λ) is the second-order derivative fitting function value, R2 is the goodness of fit, κ is the factor affecting the goodness of fit, W(λ) is the weight function based on the wavelength distance, λ c is the center point of the characteristic wavelength, σ is the distance attenuation coefficient, is the window size correction factor, and n is the number of wavelength points in the calculation window.

[0083] It should be noted that the effective detection threshold for determining the valley depth specifically includes: collecting 100 standard samples of organophosphorus pesticides with known different concentrations, with the concentration range of 0.01 - 10 mg / kg, and obtaining the hyperspectral data of these samples under the same test conditions; extracting the second derivative valley depth at three characteristic wavelength points for each standard sample to obtain the corresponding relationship between the valley depth and the concentration; when the pesticide concentration is lower than 0.05 mg / kg, recording the corresponding valley depth value, and statistically calculating the average value and standard deviation of the valley depth of these low-concentration samples; taking the average value plus 3 times the standard deviation as the initial threshold; using 20 blank samples (without pesticide residues) for verification experiments, and recording the valley depth values of these samples at the characteristic wavelength points; when the valley depth of the blank samples is greater than the initial threshold, appropriately increasing the threshold until the false positive rate is lower than 1%; for each characteristic wavelength point, different thresholds may need to be set because the background noise levels at different wavelengths may vary; when the ambient temperature varies within the range of 20 - 30 °C, checking the stability of the threshold, and establishing a temperature-related threshold correction factor if necessary; in actual detection, when the signal-to-noise ratio of the sample is lower than 10:1, increasing the threshold by 20% to avoid misjudgment.

[0084] In summary, the non-destructive food detection method 100 based on hyperspectral imaging according to the embodiments of the present invention is elucidated. It determines the position information of the leaf fold area by analyzing and judging the hyperspectral image data, and then calculates the pesticide residue distribution map on the surface of leafy vegetables through the characteristic absorption band range of pesticides. In this way, the technical problem of low detection accuracy caused by diffuse reflection interference due to leaf folds in the prior art can be solved, thereby helping to determine whether there is pesticide residue, reducing the inflow of vegetables with pesticide residues into the market, and further reducing food poisoning incidents caused by pesticide residues.

[0085] Here, those skilled in the art can understand that the specific operations of each step in the above non-destructive food detection system based on hyperspectral imaging have been introduced in detail in the description of the Figure 1 non-destructive food detection method based on hyperspectral imaging above, and therefore, the repeated description thereof will be omitted.

[0086] According to another aspect of the present invention, there is provided a non-destructive food detection system based on hyperspectral imaging, which includes:

[0087] An acquisition module, configured to acquire hyperspectral image data of the to-be-detected leafy vegetables, and obtain a three-dimensional hyperspectral data cube including a spatial dimension and a spectral dimension;

[0088] A position marking module, configured to identify the fold area of the three-dimensional hyperspectral data cube, determine the position information of the leaf fold area by using spatial domain gray-scale gradient analysis, and obtain the position marking data of the leaf fold area;

[0089] A compensation module, configured to adaptively perform spectral compensation on the wrinkled area in the three-dimensional hyperspectral data cube according to the position marking data, establish a compensation factor by calculating the spectral reflectance difference between the adjacent non-wrinkled areas of the wrinkled area, and obtain the compensated spectral data after eliminating the influence of wrinkles;

[0090] An extraction module, configured to extract the second derivative features of the spectral curve within the characteristic absorption band range of the organophosphorus pesticide based on the compensated spectral data, and calculate and obtain the distribution map of the organophosphorus pesticide residue concentration on the surface of the to-be-detected leafy vegetables.

[0091] In summary, a non-destructive food detection system based on hyperspectral imaging according to an embodiment of the present invention is elucidated. It determines the position information of the leaf wrinkled area by analyzing and judging the hyperspectral image data, and then calculates the distribution map of the pesticide residue on the surface of the leafy vegetables through the characteristic absorption band range of the pesticide. In this way, the technical problem of low detection accuracy caused by the diffuse reflection interference due to leaf wrinkles in the prior art can be solved, so as to help judge whether there is pesticide residue, reduce the inflow of vegetables with pesticide residue into the market, and further reduce the food poisoning incidents caused by pesticide residue.

Claims

1. A non-destructive testing method for food based on hyperspectral imaging, characterized in that, Including: Obtain the hyperspectral image data of the leafy vegetables to be measured, and obtain a three-dimensional hyperspectral data cube including the spatial dimension and the spectral dimension; Identify the wrinkled area of the three-dimensional hyperspectral data cube, and use the spatial domain gray gradient analysis to determine the position information of the leaf wrinkled area, and obtain the position marking data of the leaf wrinkled area; According to the position marking data, perform adaptive spectral compensation on the wrinkled area in the three-dimensional hyperspectral data cube, calculate the spectral reflectance difference between the adjacent non-wrinkled areas of the wrinkled area, establish a compensation factor, and obtain the compensated spectral data after eliminating the influence of wrinkles; Based on the compensated spectral data, within the characteristic absorption band range of the organophosphorus pesticide, extract the second derivative characteristics of the spectral curve, and calculate the distribution map of the organophosphorus pesticide residue concentration on the surface of the leafy vegetables to be measured.

2. The non-destructive detection method for food based on hyperspectral imaging according to claim 1, wherein The identification of the wrinkled area of the three-dimensional hyperspectral data cube includes: Extract the gray-scale image of a fixed wavelength from the three-dimensional hyperspectral data cube as the reference image, and perform noise reduction processing; Calculate the gray-scale gradient amplitude of the denoised gray-scale image, and mark the pixel points with a gray-scale gradient greater than the gray-scale gradient threshold as the wrinkled edge candidate points; Perform morphological processing on the wrinkled edge candidate points to obtain the initial wrinkled area; Optimize the initial wrinkled area to obtain the position marking data of each wrinkled area.

3. The non-destructive detection method of food based on hyperspectral imaging according to claim 2, characterized in that The morphological processing of the wrinkled edge candidate points includes multiple morphological processes, and determines whether to perform dilation operations according to adjacent connected regions, and finally performs boundary smoothing processing to obtain the boundary of the wrinkled area after morphological processing.

4. The non-destructive food detection method based on hyperspectral imaging according to claim 3, wherein The dilation operation is based on the morphological characteristics of leaf wrinkles, constructs an adaptive elliptical structural element, performs dilation operations by translating the rotated elliptical structural element point by point in the image, only processes the internal area of the image, and records the pixel position information after dilation. Finally, a position index table of the wrinkled area is established to accurately mark the wrinkled area.

5. The non-destructive testing method for food based on hyperspectral imaging according to claim 1, characterized in that The compensation factor determines whether to use the distance weighted method or the average spectrum of the global reference to calculate the compensation factor according to the Euclidean distance between each pixel point in the wrinkled area and the sampling points of the nearest reference area.

6. The non-destructive food detection method based on hyperspectral imaging according to claim 5, wherein The distance weighted method includes: Measure the Euclidean distance between the pixel points to be processed in the wrinkled area and the reference sampling points, select multiple reference points with the smallest distance, and calculate the weight coefficient according to their distances; Combine the spectral reflectance difference and the weight coefficient, calculate the weighted spectral difference and sum them up to obtain the spectral compensation factor of the pixel points to be processed; If the pixel point is located in the wrinkled edge area, calculate the attenuation coefficient based on its shortest distance to the wrinkled boundary, and multiply it by the spectral compensation factor to obtain the final compensation value; Repeat the above process to calculate the spectral compensation factors at each wavelength position.

7. The non-destructive food detection method based on hyperspectral imaging according to claim 6, wherein Calculating the compensation factor through the average spectrum of the global reference includes: Extract the annular reference area outside the wrinkles, divide it into upper and lower sub-areas, and calculate the spectral mean value of each wavelength; When the spectral difference between the two sub-areas is less than the preset value, take the average value as the reference spectrum, otherwise select the spectral of the sub-area closer to the pixel points to be processed as the reference; Calculate the difference between the pixel to be processed and the reference spectrum, perform clipping processing on the over-limit part to generate an initial compensation factor, calculate the position weight according to the distance between the pixel and the center line of the fold, multiply the two to obtain the final compensation factor, and repeat the process to generate the compensation sequence for each wavelength.

8. The non-destructive food detection system based on hyperspectral imaging according to claim 1, characterized in that, Extracting the second derivative features of the spectral curve includes: Preprocess the compensated spectral data; Based on the smoothed spectral data, use the five-point derivative formula to calculate the first derivative, and then use the same method to calculate the second derivative of the first derivative curve.

9. The non-destructive food detection system based on hyperspectral imaging according to claim 8, wherein The five-point derivative formula calculates the first derivative by selecting two adjacent wavelength points on both sides of the current wavelength point to form a five-point derivative window, and uses fixed coefficients to calculate the first derivative of the center point; while the three-point derivative formula is used for correction at both ends of the spectral curve; finally, the complete first derivative curve is generated by translating the derivative window point by point.

10. A non-destructive food detection system based on hyperspectral imaging, based on the non-destructive food detection method based on hyperspectral imaging according to any one of claims 1 to 9, characterized in that, Including: An acquisition module for acquiring hyperspectral image data of the to-be-detected leafy vegetables to obtain a three-dimensional hyperspectral data cube including a spatial dimension and a spectral dimension; A position marking module for identifying the fold region of the three-dimensional hyperspectral data cube, determining the position information of the leaf fold region by using spatial domain gray gradient analysis, and obtaining the position marking data of the leaf fold region; A compensation module for adaptively compensating the spectrum of the fold region in the three-dimensional hyperspectral data cube according to the position marking data, establishing a compensation factor by calculating the spectral reflectance difference between adjacent non-fold regions of the fold region, and obtaining the compensated spectral data after eliminating the influence of the fold; An extraction module for extracting the second derivative features of the spectral curve within the characteristic absorption band range of the organophosphorus pesticide based on the compensated spectral data, and calculating the distribution map of the organophosphorus pesticide residue concentration on the surface of the to-be-detected leafy vegetables.

Images