A method for assessing nitrate treatment in water bodies based on multispectral imagery

Abstract

This invention discloses a method for assessing nitrate treatment in water bodies based on multispectral imagery, specifically relating to the fields of image processing and multispectral image analysis. The method includes acquiring multispectral images before and after treatment, extracting water surface texture from the multispectral images, performing edge detection on the extracted water surface texture to obtain disturbance boundary points, and recording the positions of these disturbance boundary points at different acquisition times to form a disturbance boundary point sequence. This invention replaces the traditional image matching and image correction system with a fluid coordinate grid based on water disturbance boundary points, disturbance trajectories, and flow direction fields. This allows for the establishment of reliable pixel correspondences between multi-temporal multispectral images before and after treatment, even when the reference is destroyed. This solves the problem in existing technologies where the effectiveness of nitrate treatment cannot be reliably assessed using multispectral images.

Description

Technical Field

[0001] This invention relates to the field of image processing and multispectral image analysis technology, and more specifically, to a method for assessing nitrate treatment in water bodies based on multispectral images. Background Technology

[0002] In current water nitrate treatment assessment technologies based on multispectral images, the mainstream methods usually rely on relatively stable spatial reference targets such as shorelines, riverbed textures, and fixed structures. First, basic image correction is completed through geometric registration and radiometric calibration. Then, image matching algorithms such as feature point extraction, regional correspondence, and deformation model fitting are combined to unify multispectral images from different time phases into the same spatial reference frame, which supports the comparative analysis of nitrate concentration distribution before and after treatment. However, in actual governance projects, nitrate reduction is often accompanied by strong intervention measures such as bottom sediment dredging, ecological base laying, river channel shaping and slope reinforcement. These operations systematically reshape the micro-topography and sedimentary structure of the riverbed, causing large-scale destruction or replacement of the original bottom texture. The composition of bottom particles, reflective roughness and turbidity field distribution may be fundamentally changed. In such a scenario of overall reconstruction of bottom reference, the stable control points and optical response models on which traditional image correction is based are no longer effective. Image matching algorithms also have difficulty in determining whether pixels in two temporal images still correspond to the same water body unit due to the lack of traceable spatial features. Consequently, they cannot distinguish whether the spectral differences are due to changes in nitrate concentration or to the reconstruction of the reflective field caused by the evolution of bottom structure. Therefore, the assessment of the effectiveness of nitrate treatment using multispectral images must rely on accurate image matching and correction to ensure the spatial comparability of images from different times. However, the nitrate treatment process significantly alters the riverbed structure, causing the traditional registration and correction system based on bottom texture and fixed features to lose its physical basis. Ultimately, this results in the inability to establish a reliable correspondence between multi-temporal images, which severely weakens the credibility of nitrate treatment assessment. This is precisely the technical problem that urgently needs to be solved. Summary of the Invention

[0003] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a method for assessing nitrate treatment in water bodies based on multispectral images. By constructing a fluid coordinate grid based on water disturbance boundary points, disturbance trajectories, and flow direction fields to replace the traditional image matching and image correction system with fixed reference objects, a reliable pixel correspondence can still be established between multi-temporal multispectral images before and after treatment, even when the reference objects are destroyed. This solves the problem mentioned in the background art that the effect of nitrate treatment cannot be reliably assessed using multispectral images.

[0004] To achieve the above objectives, the present invention provides the following technical solution: a method for assessing nitrate treatment in water bodies based on multispectral imagery, comprising: S1. Perform acquisition operations on multispectral images before and after treatment, extract water surface texture from the multispectral images, perform edge detection on the extracted water surface texture to obtain disturbed boundary points; record the positions of the disturbed boundary points at different acquisition times to form a disturbed boundary point sequence; S2. Calculate the position difference between adjacent time points for each point in the disturbance boundary point sequence to form a disturbance displacement sequence; perform cumulative processing on the disturbance displacement sequence in chronological order to solve for the disturbance trajectory; S3. Calculate the tangential direction segment by segment for the disturbance trajectory, and perform directional statistics on the tangential direction according to the spatial region to form the flow direction field; arrange directional constraint lines in the image space according to the flow direction field, and solve the fluid coordinate grid. S4. Perform coordinate mapping on the multispectral images before and after treatment according to the fluid coordinate grid to form two sets of coordinate mapping point sets; perform point-by-point comparison on the disturbance trajectory morphology of the two sets of coordinate mapping point sets, calculate the morphological difference of each comparison point, and rearrange the coordinate mapping point sets according to the morphological difference to form corresponding position point sets; S5. Calculate the spectral difference curve pixel by pixel for the multispectral image after adjustment of the corresponding location point set, perform regional stability screening on the spectral difference curve, screen out the spectral difference results of unstable regions, and perform spatial aggregation on the screened spectral difference results to form a nitrate change map.

[0005] In a preferred embodiment, S1 further includes acquiring multispectral image data by performing synchronous acquisition operations on multispectral images before and after treatment, performing channel-by-channel reading on the multispectral image data, locking each spectral channel sequentially by spectral channel index, reading the reflection intensity value of the water surface area pixel by pixel in each spectral channel, summarizing the read water surface reflection intensity values ​​to generate a channel reflection set of water surface reflection, and rearranging the channel reflection set according to the channel order to form a channel reflection sequence. The pixel difference distribution is obtained by performing inter-channel difference calculation on the channel reflection sequence, and the neighborhood difference change calculation is performed on the pixel difference distribution. The difference change is obtained by calculating the difference change of adjacent pixels point by point, and then the gradient is reconstructed according to the pixel index to form the water surface texture. By performing difference calculation on the neighboring pixels of each pixel in the water surface texture, each difference is compared with the adjacent difference point by point to solve the difference mutation position, and the difference mutation position is constructed into a mutation point set according to the pixel index. The adjacent mutation points in the mutation point set are expanded point by point to form a continuous mutation chain, and the pixels in the continuous mutation chain are marked as perturbation boundary points. The spatial coordinates of the disturbance boundary points are obtained by reading their positions point by point at different acquisition times. The spatial coordinates are then bound to a time index. Time series records are formed by writing the corresponding indexes according to the acquisition times. The time series records are then arranged in order to form a disturbance boundary point sequence and output.

[0006] In a preferred embodiment, S2 further includes performing point-by-point reading of the perturbation boundary point sequence to obtain the coordinates of the perturbation boundary points at adjacent times, performing subtraction operations on the horizontal and vertical coordinates of adjacent perturbation boundary points to obtain the horizontal difference and the vertical difference, and recording the obtained horizontal difference and vertical difference in the order of position index to form a position difference record; By writing the position difference records in the order of the acquisition time of the disturbance boundary point sequence, the horizontal and vertical differences corresponding to each time point are kept consistent with the time order of the original disturbance boundary point sequence, forming a position difference set organized by time. The horizontal and vertical differences of the first acquisition moment in the location difference set, arranged in chronological order of acquisition time, are written as the horizontal and vertical cumulative values ​​for that acquisition moment, respectively. In subsequent acquisition times, the horizontal difference at the current acquisition time is added to the horizontal cumulative value written at the previous acquisition time in the cumulative value set formed in the acquisition time order to obtain the horizontal cumulative value at the current acquisition time. The vertical difference at the current acquisition time is added to the vertical cumulative value written at the previous acquisition time in the cumulative value set formed in the acquisition time order to obtain the vertical cumulative value at the current acquisition time. Then, the horizontal and vertical cumulative values ​​at the current acquisition time are written in the acquisition time order to form a cumulative value set arranged by time.

[0007] In a preferred embodiment, S2 further includes, during the process of performing time-by-time accumulation on the position difference set, before writing the horizontal and vertical accumulated values ​​at the current acquisition time into the accumulated value set formed in the order of acquisition time, determining whether the difference between the horizontal accumulated value at the current acquisition time and the horizontal accumulated value written at the previous acquisition time in the set, and whether the difference between the vertical accumulated value at the current acquisition time and the vertical accumulated value written at the previous acquisition time in the set exceeds the allowable range of continuous change: If none of them exceed the allowable range of continuous change, then the horizontal and vertical cumulative values ​​at the current acquisition time are written into the cumulative value set in the order of acquisition time. If any item exceeds the allowable range of continuous change, the horizontal and vertical cumulative values ​​to be written at the current acquisition time will be set to the horizontal and vertical cumulative values ​​written at the previous acquisition time in this set, respectively, and written into the cumulative value set in the order of acquisition time. After the accumulated values ​​at all acquisition times are written, the horizontal and vertical accumulated values ​​arranged in the acquisition time order in the accumulated value set are organized to form a disturbance displacement sequence. By performing point-by-point connection on the disturbance displacement sequence, the horizontal and vertical cumulative values ​​at adjacent acquisition times in the disturbance displacement sequence are connected in chronological order to form a continuous point series, and it is determined whether there is a positional break between adjacent connection segments: If no positional break occurs, the continuous point series is connected in the order of acquisition time to form a connection result, and the connection result is organized into a disturbance trajectory in the order of acquisition time. If a positional break occurs, the horizontal and vertical cumulative values ​​corresponding to the acquisition time at which the positional break occurs will be reset to the horizontal and vertical cumulative values ​​of the previous acquisition time arranged in the order of acquisition time. Based on the reset horizontal and vertical cumulative values, the disturbance displacement sequence will be calibrated once. Then, the disturbance trajectory will be generated by connecting the points one by one according to the calibrated disturbance displacement sequence and organizing them in the order of acquisition time.

[0008] In a preferred embodiment, S3 further includes performing segment-by-segment reading of the disturbance trajectory, calculating the difference between the spatial coordinates corresponding to adjacent acquisition times in the disturbance trajectory to obtain the difference between the horizontal coordinates and the difference between the vertical coordinates, calculating the ratio between the difference between the horizontal coordinates and the difference between the vertical coordinates to obtain the direction ratio, and then writing the direction ratios in the order of acquisition time to form a direction ratio sequence. The direction ratios of each segment in the direction ratio sequence are decomposed into lateral direction components and longitudinal direction components. Normalized direction components are obtained by performing proportional normalization on the lateral and longitudinal direction components, and written in the order of acquisition time to form tangential directions. Perform a direction consistency comparison on the tangential directions within the same region to determine whether the direction difference between any tangential direction within the region and other tangential directions within the same region exceeds a preset direction deviation threshold. If the direction difference between all tangential directions in the region does not exceed the preset direction deviation threshold, then the first tangential direction in the region arranged in the order of collection time is taken as the region's direction statistics result. If there is a tangential direction whose direction difference exceeds the preset direction deviation threshold, then the tangential direction whose direction difference exceeds the threshold is proportionally corrected, and the corrected tangential direction is used as the direction statistics result of the region. The directional statistical results of each region are then written into the flow direction field in regional order. The directional statistical results of each region in the flow direction field are used as the directional constraint lines to form directional constraint lines within the directional constraint execution area, and the continuity of line segment connections between adjacent directional constraint lines is judged. If the direction difference between adjacent direction constraint lines does not exceed the preset direction continuity threshold, the direction constraint lines are written in spatial coordinate order; if the direction difference between adjacent direction constraint lines exceeds the preset direction continuity threshold, the direction statistics results for this area are corrected once, and the corrected direction statistics results are re-connected within the area before being written into the direction constraint lines. By sequentially connecting the directional constraint lines according to the image space, and indexing and organizing the connections between the directional constraint lines according to spatial coordinates, a fluid coordinate grid is formed.

[0009] In a preferred embodiment, S4 further includes performing pixel-by-pixel reading on the multispectral images before and after treatment, writing the grid index of each pixel in the fluid coordinate grid, and binding the grid index with the spatial coordinates of the corresponding pixel to form a coordinate mapping point set before and after treatment; The coordinate mapping points with the same grid index in the coordinate mapping point set before and after the treatment are read point by point. The difference between the horizontal and vertical coordinates of the corresponding points is calculated to form the morphological difference. The morphological difference is then compared point by point to determine whether the morphological difference exceeds the preset difference threshold. If the difference threshold is not exceeded, the corresponding mapping point is written into the original index to form the initial corresponding point record; if the difference threshold is exceeded, the spatial coordinates of the mapping point are adjusted once, and the adjusted spatial coordinates are written into the original index to form the initial corresponding point record. For the initial corresponding point record, perform coordinate continuity detection within the region, calculate the difference between the spatial coordinates of adjacent corresponding points within the same region, and determine whether the difference exceeds the continuity threshold: If the continuity threshold is not exceeded, the corresponding points in the region are written into the corresponding point set of the region in grid order; if the continuity threshold is exceeded, the corresponding points in the region whose difference exceeds the continuity threshold are adjusted once, and the adjusted corresponding points are written into the corresponding point set of the region in grid order. All corresponding point sets in all regions are indexed and organized according to the spatial order of the fluid coordinate grid to form corresponding location point sets.

[0010] In a preferred embodiment, S5 further includes performing pixel-by-pixel reading on the pre-treatment multispectral image and the post-treatment multispectral image after adjustment of the corresponding location point set, writing the reflection intensity values ​​of the pre-treatment pixels at the same corresponding location point in each spectral channel in spectral channel order to form a pre-treatment channel reflection sequence, and writing the reflection intensity values ​​of the post-treatment pixels at the same corresponding location point in each spectral channel in spectral channel order to form a post-treatment channel reflection sequence. The channel-by-channel difference was calculated between the pre-treatment channel reflection sequence and the post-treatment channel reflection sequence, and the channel difference was recorded according to the spectral channel index to form an initial spectral difference curve. The initial spectral difference curve is read point by point within the region. The spectral differences of adjacent pixels within the same region are calculated in pixel order to obtain the spectral difference value. It is then determined whether the spectral difference value exceeds a preset stability threshold. If the preset stability threshold is not exceeded, the spectral differences of adjacent pixels in the region are written into the region stability difference record in the region order. If the preset stability threshold is exceeded, the spectral difference of the pixel will be adjusted in one step, and the adjusted spectral difference will be written into the region stable difference record in the order of the regions. All regional stable difference records are sorted in order of execution region to form a stable spectral difference set.

[0011] In a preferred embodiment, S5 further includes performing spatial reading on the stable spectral difference set, connecting the spectral differences within the same spatial neighborhood point by point according to spatial coordinate order, calculating the difference between the spectral differences between adjacent connected points, and determining whether the difference between adjacent connected points exceeds a preset aggregation threshold. If the preset aggregation threshold is not exceeded, the adjacent connection points will be written in spatial coordinate order to form the initial aggregation record; If the preset aggregation threshold is exceeded, the spectral difference of the connection point will be aggregated and corrected once, and the corrected connection point will be written into the initial aggregation record in spatial coordinate order. All initial aggregated records are spatially ordered to form aggregated spectral difference records, and then the aggregated spectral difference records are indexed according to the spatial order of the images to form a nitrate variation map.

[0012] The technical effects and advantages of this invention are as follows: This invention replaces the solid ground object reference system that has been destroyed after treatment by constructing a fluid coordinate grid based on perturbation boundary points, perturbation trajectories and flow direction fields. This enables multi-temporal multispectral images to still establish stable spatial correspondences even when reference objects fail, thus solving the problem that existing image matching and image correction failures lead to unreliable assessment of nitrate changes. This invention reads and reconstructs the channel reflection sequence channel by channel, obtains the water surface texture and locates the disturbance boundary points by the difference change, so that the dynamic water texture becomes a spatial reference that can be traced across time phases, ensuring that image registration no longer depends on riverbed texture or fixed structures, and improving the usability of feature extraction in governance scenarios. This invention constructs a disturbance trajectory by performing accumulation, continuity judgment, and one-time correction on disturbance boundary points, so that water disturbance presents a continuous structure in the time dimension, providing stable directional information for subsequent fluid coordinate construction; This invention calculates the tangential direction segment by segment of the disturbance trajectory, performs regional direction statistics, and generates direction constraint lines, enabling the fluid coordinate grid to express the flow direction structure of water in different regions and improving the stability of image spatial mapping under local complex disturbance conditions. This invention performs coordinate mapping and constructs corresponding point sets on multispectral images before and after treatment under a fluid coordinate grid, and performs stability screening and spatial aggregation on the spectral difference results, so that the formation of the nitrate change map is based on the spectral differences of real homologous pixels, ensuring that the treatment effect assessment has sufficient data reliability. Attached Figure Description

[0013] Figure 1 This is a flowchart of the method steps of the present invention. Detailed Implementation

[0014] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0015] Refer to the instruction manual appendix Figure 1 An embodiment of the present invention provides a method for assessing nitrate treatment in water bodies based on multispectral imagery, comprising: S1. Perform acquisition operations on multispectral images before and after treatment, extract water surface texture from the multispectral images, perform edge detection on the extracted water surface texture to obtain disturbed boundary points; record the positions of the disturbed boundary points at different acquisition times to form a disturbed boundary point sequence; S2. Calculate the position difference between adjacent time points for each point in the disturbance boundary point sequence to form a disturbance displacement sequence; perform cumulative processing on the disturbance displacement sequence in chronological order to solve for the disturbance trajectory; S3. Calculate the tangential direction segment by segment for the disturbance trajectory, and perform directional statistics on the tangential direction according to the spatial region to form the flow direction field; arrange directional constraint lines in the image space according to the flow direction field, and solve the fluid coordinate grid. S4. Perform coordinate mapping on the multispectral images before and after treatment according to the fluid coordinate grid to form two sets of coordinate mapping point sets; perform point-by-point comparison on the disturbance trajectory morphology of the two sets of coordinate mapping point sets, calculate the morphological difference of each comparison point, and rearrange the coordinate mapping point sets according to the morphological difference to form corresponding position point sets; S5. Calculate the spectral difference curve pixel by pixel for the multispectral image after adjustment of the corresponding location point set, perform regional stability screening on the spectral difference curve, screen out the spectral difference results of unstable regions, and perform spatial aggregation on the screened spectral difference results to form a nitrate change map.

[0016] In S1, the process also includes acquiring multispectral image data by performing synchronous acquisition operations on multispectral images before and after treatment, reading the multispectral image data channel by channel, locking each spectral channel sequentially by spectral channel index, reading the reflection intensity value of the water surface area pixel by pixel in each spectral channel, summarizing the read water surface reflection intensity values ​​to generate a channel reflection set of water surface reflection, and rearranging the channel reflection set according to the channel order to form a channel reflection sequence. The pixel difference distribution is obtained by performing inter-channel difference calculation on the channel reflection sequence, and the neighborhood difference change calculation is performed on the pixel difference distribution. The difference change is obtained by calculating the difference change of adjacent pixels point by point, and then the gradient is reconstructed according to the pixel index to form the water surface texture. By performing difference calculation on the neighboring pixels of each pixel in the water surface texture, each difference is compared with the adjacent difference point by point to solve the difference mutation position, and the difference mutation position is constructed into a mutation point set according to the pixel index. The adjacent mutation points in the mutation point set are expanded point by point to form a continuous mutation chain, and the pixels in the continuous mutation chain are marked as perturbation boundary points. The spatial coordinates of the disturbance boundary points are obtained by reading their positions point by point at different acquisition times. The spatial coordinates are then bound to a time index. Time series records are formed by writing the corresponding indexes according to the acquisition times. The time series records are then arranged in order to form a disturbance boundary point sequence and output.

[0017] In S2, the process also includes reading the perturbation boundary point sequence point by point to obtain the coordinates of the perturbation boundary points at adjacent time points, performing subtraction operations on the horizontal and vertical coordinates of adjacent perturbation boundary points to obtain the horizontal and vertical differences, and recording the obtained horizontal and vertical differences in the order of the position index to form a position difference record. By writing the position difference records in the order of the acquisition time of the disturbance boundary point sequence, the horizontal and vertical differences corresponding to each time point are kept consistent with the time order of the original disturbance boundary point sequence, forming a position difference set organized by time. The horizontal and vertical differences of the first acquisition moment in the location difference set, arranged in chronological order of acquisition time, are written as the horizontal and vertical cumulative values ​​for that acquisition moment, respectively. In subsequent acquisition times, the horizontal difference at the current acquisition time is added to the horizontal cumulative value written at the previous acquisition time in the cumulative value set formed in the acquisition time order to obtain the horizontal cumulative value at the current acquisition time. The vertical difference at the current acquisition time is added to the vertical cumulative value written at the previous acquisition time in the cumulative value set formed in the acquisition time order to obtain the vertical cumulative value at the current acquisition time. Then, the horizontal and vertical cumulative values ​​at the current acquisition time are written in the acquisition time order to form a cumulative value set arranged by time.

[0018] In S2, before writing the horizontal and vertical accumulated values ​​at the current acquisition time into the accumulated value set formed in the order of acquisition time, it is also necessary to determine whether the difference between the horizontal accumulated value at the current acquisition time and the horizontal accumulated value written at the previous acquisition time in the set, and whether the difference between the vertical accumulated value at the current acquisition time and the vertical accumulated value written at the previous acquisition time in the set exceeds the allowable range of continuous change, during the process of performing time-by-time accumulation on the set of position difference values. If none of them exceed the allowable range of continuous change, then the horizontal and vertical cumulative values ​​at the current acquisition time are written into the cumulative value set in the order of acquisition time. If any item exceeds the allowable range of continuous change, the horizontal and vertical cumulative values ​​to be written at the current acquisition time will be set to the horizontal and vertical cumulative values ​​written at the previous acquisition time in this set, respectively, and written into the cumulative value set in the order of acquisition time. After the accumulated values ​​at all acquisition times are written, the horizontal and vertical accumulated values ​​arranged in the acquisition time order in the accumulated value set are organized to form a disturbance displacement sequence. By performing point-by-point connection on the disturbance displacement sequence, the horizontal and vertical cumulative values ​​at adjacent acquisition times in the disturbance displacement sequence are connected in chronological order to form a continuous point series, and it is determined whether there is a positional break between adjacent connection segments: If no positional break occurs, the continuous point series is connected in the order of acquisition time to form a connection result, and the connection result is organized into a disturbance trajectory in the order of acquisition time. If a positional break occurs, the horizontal and vertical cumulative values ​​corresponding to the acquisition time at which the positional break occurs will be reset to the horizontal and vertical cumulative values ​​of the previous acquisition time arranged in the order of acquisition time. Based on the reset horizontal and vertical cumulative values, the disturbance displacement sequence will be calibrated once. Then, the disturbance trajectory will be generated by connecting the points one by one according to the calibrated disturbance displacement sequence and organizing them in the order of acquisition time.

[0019] In S3, the process also includes reading the disturbance trajectory segment by segment, calculating the difference between the spatial coordinates corresponding to adjacent acquisition times in the disturbance trajectory to obtain the difference between the horizontal and vertical coordinates, calculating the ratio between the difference between the horizontal and vertical coordinates to obtain the direction ratio, and then writing the direction ratios in the order of acquisition time to form a direction ratio sequence. The direction ratios of each segment in the direction ratio sequence are decomposed into lateral direction components and longitudinal direction components. Normalized direction components are obtained by performing proportional normalization on the lateral and longitudinal direction components, and written in the order of acquisition time to form tangential directions. Perform a direction consistency comparison on the tangential directions within the same region to determine whether the direction difference between any tangential direction within the region and other tangential directions within the same region exceeds a preset direction deviation threshold. If the direction difference between all tangential directions in the region does not exceed the preset direction deviation threshold, then the first tangential direction in the region arranged in the order of collection time is taken as the region's direction statistics result. If there is a tangential direction whose direction difference exceeds the preset direction deviation threshold, then the tangential direction whose direction difference exceeds the threshold is proportionally corrected, and the corrected tangential direction is used as the direction statistics result of the region. The directional statistical results of each region are then written into the flow direction field in regional order. The directional statistical results of each region in the flow direction field are used as the directional constraint lines to form directional constraint lines within the directional constraint execution area, and the continuity of line segment connections between adjacent directional constraint lines is judged. If the direction difference between adjacent direction constraint lines does not exceed the preset direction continuity threshold, the direction constraint lines are written in spatial coordinate order; if the direction difference between adjacent direction constraint lines exceeds the preset direction continuity threshold, the direction statistics results for this area are corrected once, and the corrected direction statistics results are re-connected within the area before being written into the direction constraint lines. By sequentially connecting the directional constraint lines according to the image space, and indexing and organizing the connections between the directional constraint lines according to spatial coordinates, a fluid coordinate grid is formed.

[0020] In S4, the process also includes performing pixel-by-pixel readings on the multispectral images before and after the treatment, writing the grid index of each pixel in the fluid coordinate grid, and binding the grid index with the spatial coordinates of the corresponding pixel to form a coordinate mapping point set before and after the treatment. The coordinate mapping points with the same grid index in the coordinate mapping point set before and after the treatment are read point by point. The difference between the horizontal and vertical coordinates of the corresponding points is calculated to form the morphological difference. The morphological difference is then compared point by point to determine whether the morphological difference exceeds the preset difference threshold. If the difference threshold is not exceeded, the corresponding mapping point is written into the original index to form the initial corresponding point record; if the difference threshold is exceeded, the spatial coordinates of the mapping point are adjusted once, and the adjusted spatial coordinates are written into the original index to form the initial corresponding point record. For the initial corresponding point record, perform coordinate continuity detection within the region, calculate the difference between the spatial coordinates of adjacent corresponding points within the same region, and determine whether the difference exceeds the continuity threshold: If the continuity threshold is not exceeded, the corresponding points in the region are written into the corresponding point set of the region in grid order; if the continuity threshold is exceeded, the corresponding points in the region whose difference exceeds the continuity threshold are adjusted once, and the adjusted corresponding points are written into the corresponding point set of the region in grid order. All corresponding point sets in all regions are indexed and organized according to the spatial order of the fluid coordinate grid to form corresponding location point sets.

[0021] In S5, the method also includes performing pixel-by-pixel reading on the pre-treatment multispectral image and the post-treatment multispectral image after adjustment of the corresponding location point set, writing the reflection intensity values ​​of the pre-treatment pixels at the same corresponding location point in each spectral channel in the order of the spectral channel index to form a pre-treatment channel reflection sequence, and writing the reflection intensity values ​​of the post-treatment pixels at the same corresponding location point in each spectral channel in the order of the spectral channel index to form a post-treatment channel reflection sequence. The channel-by-channel difference was calculated between the pre-treatment channel reflection sequence and the post-treatment channel reflection sequence, and the channel difference was recorded according to the spectral channel index to form an initial spectral difference curve. The initial spectral difference curve is read point by point within the region. The spectral differences of adjacent pixels within the same region are calculated in pixel order to obtain the spectral difference value. It is then determined whether the spectral difference value exceeds a preset stability threshold. If the preset stability threshold is not exceeded, the spectral differences of adjacent pixels in the region are written into the region stability difference record in the region order. If the preset stability threshold is exceeded, the spectral difference of the pixel will be adjusted in one step, and the adjusted spectral difference will be written into the region stable difference record in the order of the regions. All regional stable difference records are sorted in order of execution region to form a stable spectral difference set.

[0022] S5 also includes performing spatial readings on a stable spectral difference set, connecting the spectral differences within the same spatial neighborhood point by point according to spatial coordinates, calculating the difference between the spectral differences between adjacent connected points, and determining whether the difference between adjacent connected points exceeds a preset aggregation threshold. If the preset aggregation threshold is not exceeded, the adjacent connection points will be written in spatial coordinate order to form the initial aggregation record; If the preset aggregation threshold is exceeded, the spectral difference of the connection point will be aggregated and corrected once, and the corrected connection point will be written into the initial aggregation record in spatial coordinate order. All initial aggregated records are spatially ordered to form aggregated spectral difference records, and then the aggregated spectral difference records are indexed according to the spatial order of the images to form a nitrate variation map.

[0023] Existing multispectral image assessment methods generally require acquiring multi-temporal images before and after remediation. These images are then aligned to the same coordinate system through geometric image correction, radiometric image correction, and image matching based on fixed features. Spectral differences are then analyzed to infer nitrate variation trends. However, in actual remediation projects, nitrate reduction is often accompanied by significant construction activities such as sediment dredging, riverbed reshaping, bank reinforcement, and ecological substrate laying. These activities directly alter the reflective structure beneath the water body, causing irreversible changes to fixed features, bottom texture, and sedimentary morphology. This results in the original image correction model lacking a reference basis. Furthermore, image matching algorithms are unable to find homologous points with continuous spectral significance and spatial location in different temporal images, completely breaking the traditional pixel-consistency-based matching chain. In this context, the spectral differences between pre-remediation and post-remediation images no longer clearly correspond to changes in nitrate concentration, but are likely merely reflection field reconstructions caused by changes in riverbed structure. Therefore, it is necessary to construct a completely new image matching and image correction system that does not rely on fixed references, enabling multispectral images to maintain spatial comparability even under drastic environmental changes after remediation. In this invention, firstly, in step S1, channel reflection sequences are read and generated channel by channel to deduce water surface texture from channel difference changes. Then, perturbation boundary points are identified through the difference mutation chain, making the perturbation boundary points the only dynamic reference objects in this invention that can maintain consistency across time phases without relying on the riverbed or structures. Subsequently, in step S2, the perturbation boundary points are reconstructed into perturbation trajectories through position difference, cumulative value sequence, continuity judgment, and one-time correction. This trajectory is essentially a deformation record of the water body in the temporal dimension, with continuous, detectable, and cumulative dynamic attributes, thus constructing a flowing reference framework for cross-temporal image matching. Based on this, step S3 decomposes the disturbance trajectory into tangential directions segment by segment, and then forms a flow direction field through direction statistics, so that the dynamic pattern of the water body itself becomes the core basis to replace the traditional geometric image correction. The directional constraint line is directly derived from the flow direction field, so that the image space is re-divided according to the real flow structure of the water body. The fluid coordinate grid formed in the end constitutes a new spatial correction system, replacing the traditional spatial transformation model that relies on fixed ground objects. Through this structure, even if the riverbed is completely reshaped, the images of different time phases can still maintain a strict one-to-one correspondence in the fluid coordinate system, realizing image correction that does not rely on solid reference objects. Step S4 performs coordinate mapping on the images before and after the treatment under the fluid coordinate grid. By judging the morphological difference, adjusting the coordinates, and detecting the regional continuity, a stable set of corresponding location points is constructed. This allows image matching to be based on a dynamic reference frame rather than a static land cover reference, fundamentally avoiding matching errors caused by riverbed changes. Finally, step S5, under the completely reconstructed image matching and image correction system, performs stability screening and spatial aggregation on the spectral difference, ensuring that the nitrate change map reflects the true change in nitrate content, rather than the image spurious changes caused by the treatment. Therefore, the reason why this invention can solve the problems that traditional methods cannot solve is that it does not make minor repairs to the existing image matching and image correction system, but starts from the actual governance scenario and completely reconstructs a dynamic spatial registration system driven by water disturbance. This allows multispectral images to maintain cross-temporal comparability even under the condition that governance causes changes in land features, and ultimately ensures the effectiveness and credibility of nitrate governance assessment.

[0024] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for evaluating nitrate treatment of a water body based on multi-spectral images, characterized in that, include: S1. Perform acquisition operations on multispectral images before and after treatment, extract water surface texture from the multispectral images, perform edge detection on the extracted water surface texture to obtain disturbed boundary points; record the positions of the disturbed boundary points at different acquisition times to form a disturbed boundary point sequence; S2. Calculate the position difference between adjacent time points for each point in the disturbance boundary point sequence to form a disturbance displacement sequence; perform cumulative processing on the disturbance displacement sequence in chronological order to solve for the disturbance trajectory; S3. Calculate the tangential direction segment by segment for the disturbance trajectory, and perform directional statistics on the tangential direction according to the spatial region to form the flow direction field; arrange directional constraint lines in the image space according to the flow direction field, and solve the fluid coordinate grid. S4. Perform coordinate mapping on the multispectral images before and after treatment according to the fluid coordinate grid to form two sets of coordinate mapping point sets; perform point-by-point comparison on the disturbance trajectory morphology of the two sets of coordinate mapping point sets, calculate the morphological difference of each comparison point, and rearrange the coordinate mapping point sets according to the morphological difference to form corresponding position point sets; S5. Calculate the spectral difference curve pixel by pixel for the multispectral image after adjustment of the corresponding location point set, perform regional stability screening on the spectral difference curve, screen out the spectral difference results of unstable regions, and perform spatial aggregation on the screened spectral difference results to form a nitrate change map.

2. The method for assessing nitrate treatment in water bodies based on multispectral imagery according to claim 1, characterized in that: In S1, the process also includes acquiring multispectral image data by performing synchronous acquisition operations on multispectral images before and after treatment, reading the multispectral image data channel by channel, locking each spectral channel sequentially by spectral channel index, reading the reflection intensity value of the water surface area pixel by pixel in each spectral channel, summarizing the read water surface reflection intensity values ​​to generate a channel reflection set of water surface reflection, and rearranging the channel reflection set according to the channel order to form a channel reflection sequence. The pixel difference distribution is obtained by performing inter-channel difference calculation on the channel reflection sequence, and the neighborhood difference change calculation is performed on the pixel difference distribution. The difference change is obtained by calculating the difference change of adjacent pixels point by point, and then the gradient is reconstructed according to the pixel index to form the water surface texture. By performing difference calculation on the neighboring pixels of each pixel in the water surface texture, each difference is compared with the adjacent difference point by point to solve the difference mutation position, and the difference mutation position is constructed into a mutation point set according to the pixel index. The adjacent mutation points in the mutation point set are expanded point by point to form a continuous mutation chain, and the pixels in the continuous mutation chain are marked as perturbation boundary points. The spatial coordinates of the disturbance boundary points are obtained by reading their positions point by point at different acquisition times. The spatial coordinates are then bound to a time index. Time series records are formed by writing the corresponding indexes according to the acquisition times. The time series records are then arranged in order to form a disturbance boundary point sequence and output.

3. The method for assessing nitrate treatment in water bodies based on multispectral imagery according to claim 2, characterized in that: In S2, the process also includes reading the perturbation boundary point sequence point by point to obtain the coordinates of the perturbation boundary points at adjacent time points, performing subtraction operations on the horizontal and vertical coordinates of adjacent perturbation boundary points to obtain the horizontal and vertical differences, and recording the obtained horizontal and vertical differences in the order of the position index to form a position difference record. By writing the position difference records in the order of the acquisition time of the disturbance boundary point sequence, the horizontal and vertical differences corresponding to each time point are kept consistent with the time order of the original disturbance boundary point sequence, forming a position difference set organized by time. The horizontal and vertical differences of the first acquisition moment in the location difference set, arranged in chronological order of acquisition time, are written as the horizontal and vertical cumulative values ​​for that acquisition moment, respectively. In subsequent acquisition times, the horizontal difference at the current acquisition time is added to the horizontal cumulative value written at the previous acquisition time in the cumulative value set formed in the acquisition time order to obtain the horizontal cumulative value at the current acquisition time. The vertical difference at the current acquisition time is added to the vertical cumulative value written at the previous acquisition time in the cumulative value set formed in the acquisition time order to obtain the vertical cumulative value at the current acquisition time. Then, the horizontal and vertical cumulative values ​​at the current acquisition time are written in the acquisition time order to form a cumulative value set arranged by time.

4. The water nitrate treatment assessment method based on multispectral imagery according to claim 3, characterized in that: In S2, before writing the horizontal and vertical accumulated values ​​at the current acquisition time into the accumulated value set formed in the order of acquisition time, it is also necessary to determine whether the difference between the horizontal accumulated value at the current acquisition time and the horizontal accumulated value written at the previous acquisition time in the set, and whether the difference between the vertical accumulated value at the current acquisition time and the vertical accumulated value written at the previous acquisition time in the set exceeds the allowable range of continuous change, during the process of performing time-by-time accumulation on the set of position difference values. If none of them exceed the allowable range of continuous change, then the horizontal and vertical cumulative values ​​at the current acquisition time are written into the cumulative value set in the order of acquisition time. If any item exceeds the allowable range of continuous change, the horizontal and vertical cumulative values ​​to be written at the current acquisition time will be set to the horizontal and vertical cumulative values ​​written at the previous acquisition time in this set, respectively, and written into the cumulative value set in the order of acquisition time. After the accumulated values ​​at all acquisition times are written, the horizontal and vertical accumulated values ​​arranged in the acquisition time order in the accumulated value set are organized to form a disturbance displacement sequence. By performing point-by-point connection on the disturbance displacement sequence, the horizontal and vertical cumulative values ​​at adjacent acquisition times in the disturbance displacement sequence are connected in chronological order to form a continuous point series, and it is determined whether there is a positional break between adjacent connection segments: If no positional break occurs, the continuous point series is connected in the order of acquisition time to form a connection result, and the connection result is organized into a disturbance trajectory in the order of acquisition time. If a positional break occurs, the horizontal and vertical cumulative values ​​corresponding to the acquisition time at which the positional break occurs will be reset to the horizontal and vertical cumulative values ​​of the previous acquisition time arranged in the order of acquisition time. Based on the reset horizontal and vertical cumulative values, the disturbance displacement sequence will be calibrated once. Then, the disturbance trajectory will be generated by connecting the points one by one according to the calibrated disturbance displacement sequence and organizing them in the order of acquisition time.

5. The method for assessing nitrate treatment in water bodies based on multispectral imagery according to claim 4, characterized in that: In S3, the process also includes reading the disturbance trajectory segment by segment, calculating the difference between the spatial coordinates corresponding to adjacent acquisition times in the disturbance trajectory to obtain the difference between the horizontal and vertical coordinates, calculating the ratio between the difference between the horizontal and vertical coordinates to obtain the direction ratio, and then writing the direction ratios in the order of acquisition time to form a direction ratio sequence. The direction ratios of each segment in the direction ratio sequence are decomposed into lateral direction components and longitudinal direction components. Normalized direction components are obtained by performing proportional normalization on the lateral and longitudinal direction components, and written in the order of acquisition time to form tangential directions. Perform a direction consistency comparison on the tangential directions within the same region to determine whether the direction difference between any tangential direction within the region and other tangential directions within the same region exceeds a preset direction deviation threshold. If the direction difference between all tangential directions in the region does not exceed the preset direction deviation threshold, then the first tangential direction in the region arranged in the order of collection time is taken as the region's direction statistics result. If there is a tangential direction whose direction difference exceeds the preset direction deviation threshold, then the tangential direction whose direction difference exceeds the threshold is proportionally corrected, and the corrected tangential direction is used as the direction statistics result of the region. The directional statistical results of each region are then written into the flow direction field in regional order. The directional statistical results of each region in the flow direction field are used as the directional constraint lines to form directional constraint lines within the directional constraint execution area, and the continuity of line segment connections between adjacent directional constraint lines is judged. If the direction difference between adjacent direction constraint lines does not exceed the preset direction continuity threshold, the direction constraint lines are written in spatial coordinate order; if the direction difference between adjacent direction constraint lines exceeds the preset direction continuity threshold, the direction statistics results for this area are corrected once, and the corrected direction statistics results are re-connected within the area before being written into the direction constraint lines. By sequentially connecting the directional constraint lines according to the image space, and indexing and organizing the connections between the directional constraint lines according to spatial coordinates, a fluid coordinate grid is formed.

6. The method for assessing nitrate treatment in water bodies based on multispectral imagery according to claim 5, characterized in that: In S4, the process also includes performing pixel-by-pixel readings on the multispectral images before and after the treatment, writing the grid index of each pixel in the fluid coordinate grid, and binding the grid index with the spatial coordinates of the corresponding pixel to form a coordinate mapping point set before and after the treatment. The coordinate mapping points with the same grid index in the coordinate mapping point set before and after the treatment are read point by point. The difference between the horizontal and vertical coordinates of the corresponding points is calculated to form the morphological difference. The morphological difference is then compared point by point to determine whether the morphological difference exceeds the preset difference threshold. If the difference threshold is not exceeded, the corresponding mapping point will be written into the original index to form the initial corresponding point record. If the difference exceeds the preset threshold, the spatial coordinates of the mapping point will be adjusted once, and the adjusted spatial coordinates will be written into the initial corresponding point record according to the original index. For the initial corresponding point record, perform coordinate continuity detection within the region, calculate the difference between the spatial coordinates of adjacent corresponding points within the same region, and determine whether the difference exceeds the continuity threshold: If the continuity threshold is not exceeded, the corresponding points in the region are written into the corresponding point set of the region in grid order; If the continuity threshold is exceeded, a one-time position correction is performed on the corresponding points in the region whose difference exceeds the continuity threshold, and the corrected corresponding points are written into the corresponding point set of the region in grid order; All corresponding point sets in all regions are indexed and organized according to the spatial order of the fluid coordinate grid to form corresponding location point sets.

7. The method for assessing nitrate treatment in water bodies based on multispectral imagery according to claim 6, characterized in that: In S5, the method also includes performing pixel-by-pixel reading on the pre-treatment multispectral image and the post-treatment multispectral image after adjustment of the corresponding location point set, writing the reflection intensity values ​​of the pre-treatment pixels at the same corresponding location point in each spectral channel in the order of the spectral channel index to form a pre-treatment channel reflection sequence, and writing the reflection intensity values ​​of the post-treatment pixels at the same corresponding location point in each spectral channel in the order of the spectral channel index to form a post-treatment channel reflection sequence. The channel-by-channel difference was calculated between the pre-treatment channel reflection sequence and the post-treatment channel reflection sequence, and the channel difference was recorded according to the spectral channel index to form an initial spectral difference curve. The initial spectral difference curve is read point by point within the region. The spectral differences of adjacent pixels within the same region are calculated in pixel order to obtain the spectral difference value. It is then determined whether the spectral difference value exceeds a preset stability threshold. If the preset stability threshold is not exceeded, the spectral differences of adjacent pixels in the region are written into the region stability difference record in the region order. If the preset stability threshold is exceeded, the spectral difference of the pixel will be adjusted in one step, and the adjusted spectral difference will be written into the region stable difference record in the order of the regions. All regional stable difference records are sorted in order of execution region to form a stable spectral difference set.

8. The method for assessing nitrate treatment in water bodies based on multispectral imagery according to claim 7, characterized in that: S5 also includes performing spatial readings on a stable spectral difference set, connecting the spectral differences within the same spatial neighborhood point by point according to spatial coordinates, calculating the difference between the spectral differences between adjacent connected points, and determining whether the difference between adjacent connected points exceeds a preset aggregation threshold. If the preset aggregation threshold is not exceeded, the adjacent connection points will be written in spatial coordinate order to form the initial aggregation record; If the preset aggregation threshold is exceeded, the spectral difference of the connection point will be aggregated and corrected once, and the corrected connection point will be written into the initial aggregation record in spatial coordinate order. All initial aggregated records are spatially ordered to form aggregated spectral difference records, and then the aggregated spectral difference records are indexed according to the spatial order of the images to form a nitrate variation map.

Images