Method and system for identifying strip intercropping of orchard and grain

Abstract

The present application relates to the technical field of remote sensing image recognition, in particular to a fruit and grain strip compound planting recognition method and system, a multi-year multi-temporal vegetation index time sequence is constructed pixel by pixel and a land unit is divided; a multi-year woody component and an annual herbaceous component as a residual are separated by annual cycle harmonic fitting with the same set of harmonic coefficients in each year, without crop calendar normalization or end member unmixing; a peak phase and a rising and decaying scale of each crop of the herbaceous component are analyzed by a non-symmetrical Gaussian pulse mixing model and a crop mouth overlap degree is constructed; a cross-year phase consistency is obtained from a unit complex vector average module length of the woody component year by year; a directional modulation significance and a strip period are obtained by one-dimensional power spectrum analysis of a multi-directional integral field of the herbaceous component in a growth period, and a grain strip width ratio is obtained by inverse cosine inversion of a ratio of a second harmonic and a fundamental frequency amplitude; three criteria are normalized and a product is taken to determine the compound planting and distinguish intercropping, relay cropping and rotation.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image recognition technology, specifically to a method and system for identifying strip-shaped intercropping of medicinal herbs, fruit trees, grain crops, and cash crops in medicinal herb orchards. Background Technology

[0002] In mountainous and high-altitude areas, intercropping grain crops in strips between rows of medicinal herb gardens and orchards is an important way to simultaneously ensure food security and cash crop yields. Taking Liangshan Prefecture as an example, potatoes, corn, and early rice are intercropped in strips between rows in Artemisia argyi gardens; early rice, chili peppers, and potatoes are intercropped in pomegranate gardens; and potatoes, buckwheat, and radishes are intercropped in alpine blueberry gardens. Year-round production is achieved through precise crop rotation. This model can alleviate competition for land between grain and cash crops and the conversion of arable land to non-grain crops. However, its promotion and evaluation depend on a crucial prerequisite: the ability to accurately determine the distribution of these intercropping crops at a regional scale, whether they are intercropping, relay cropping, or crop rotation, the proportion of grain crops in the strips, and the crop rotation methods. On-site garden-by-garden surveys are difficult to scale up in mountainous areas, making remote sensing identification the only feasible method.

[0003] However, existing remote sensing crop identification methods almost entirely fail in medicinal herb and fruit orchards, and the reason can be attributed to a common technical approach. Patent application CN115661554A discloses a crop classification method based on time-series remote sensing imagery. This method's crop classification model consists of a learnable location encoding module, a spatial encoding module, a temporal encoding module, and a classification decoder. It uses gated recurrent units to encode calendar time and fuse it with spatial features, employs cross-domain temporal encoding to eliminate inconsistencies in the temporal distribution of the dataset, and first performs time-shifted enhancement training on the source dataset, followed by semi-supervised training with time-shifted enhancement attenuation in both the source and target domains. All the processing of the time dimension in this method aims to align the time series from different sources and years to a unified time reference, i.e., common phenological normalization. Patent application CN108399400A discloses a method and system for early crop identification based on high-resolution remote sensing data. This method collects Sentinel-1 and Sentinel-2 data as well as ground data, obtains synthetic time series of a specified time length using a data synthesis method, extracts classification features, evaluates feature contribution using random forest and JM distance methods to obtain optimal features, and then performs crop identification based on an immune system network algorithm. An adaptive processing mechanism is designed for mixed pixels caused by farmland fragmentation, where phenological features are only one of the parallel classification features. Patent application CN110751019A discloses a method and device for automatic crop extraction from high-resolution images based on deep learning. This method trains convolutional neural networks according to phenological stages to obtain multiple phenological classification models, matches the corresponding model according to the input data type, and extracts crop distribution, explicitly stating that it does not rely on multi-temporal series.

[0004] Although the three methods differ in their classifiers, they share the same preprocessing paradigm: treating sowing date differences as interference that needs to be aligned and eliminated, treating mixed pixels as pathological problems that require endmember spectral unmixing, treating radiation differences caused by mountain slope aspect as systematic errors that must be corrected first, and treating sub-pixel scale stripes as high-frequency texture noise that should be filtered out. This paradigm is suitable for contiguous monoculture farmland, but when used in strip intercropping of grain and cash crops in medicinal orchards, it destroys the information it seeks: Firstly, crop rotation is a controlled parameter deliberately designed and precisely executed by agronomists, and it is the defining temporal characteristic that distinguishes the intercropping pattern from monoculture. Aligning the time series to a unified benchmark is tantamount to erasing this characteristic. As a result, the phenology of fruit trees and grain crops within mixed pixels is forcibly aligned to the same phase, and the mixed time series is compressed into a blunted single peak, almost indistinguishable from the time series of monoculture orchards. Firstly, the planting of crops in medicinal orchards is divisible; secondly, the row spacing in medicinal orchards is typically 3.0m to 4.0m, and the width of the grain strips between rows is only 1.5m to 2.5m, far smaller than the pixel size of commonly used medium-to-high resolution remote sensing images. After the stripes are filtered out in the preprocessing stage, the structural attribute of the strip and the grain-to-cash crop bandwidth ratio can no longer be recovered; thirdly, the topographic radiation correction in mountainous and steep slope areas itself has a large error. The distortion introduced by the correction residual in the temporal amplitude is of the same order of magnitude as the difference in the actual phenological amplitude, making any criterion based on amplitude unreliable. The combined result of these three factors is that the grain-to-cash crop intercropping in the park is degraded into a single category of orchard in the identification results of existing methods, and its area, pattern category, and grain-to-cash crop bandwidth ratio are all unavailable. These three items are precisely the quantities required for food security assessment and monitoring of non-grain conversion of arable land.

[0005] Therefore, the core technical problem facing this field is: how to restore the phenological misalignment and sub-pixel stripes that are eliminated in the preprocessing stage as interference in medicinal and fruit orchards where the strip width is much smaller than the pixel size, mixed pixel endmembers are unavailable, mountain radiation distortion is significant, and field samples are extremely scarce, to a criterion that can directly determine the category of grain-cash crop strip intercropping pattern and invert the grain-cash crop bandwidth ratio. Summary of the Invention

[0006] Addressing the core bottleneck in existing remote sensing identification of strip-shaped intercropping of medicinal herbs, orchards, and grain-crop crops—namely, the degradation of intercropping into a single orchard category due to phenological alignment erasing of crop rotation phase differences and the filtering out of sub-pixel fringes as high-frequency noise, making it impossible to determine the pattern category and grain-crop bandwidth ratio—this invention provides a method and system for identifying strip-shaped intercropping of medicinal herbs, orchards, and grain-crop crops. This method transforms phenological misalignment from an interference to be eliminated into a defining criterion for the intercropping pattern; it transforms mixed pixels from an abundance problem to be unmixed into a two-component problem that can be naturally separated in the temporal spectrum; and it transforms sub-pixel fringes from texture noise to be filtered out into a directional modulation signal that can invert the grain-crop bandwidth ratio. Under the constraints of not increasing image resolution, not acquiring endmember spectra, and not performing topographic radiometric correction, this method achieves the determination of the intercropping pattern category and the quantitative inversion of the grain-crop bandwidth ratio from the principle level of agronomic crop rotation design and sub-pixel modulation.

[0007] The technical solution of this invention is as follows: a method for identifying strip-shaped intercropping of medicinal herbs, fruit trees, grain crops, and cash crops in a medicinal herb orchard, comprising the following steps: acquiring multi-year, multi-temporal remote sensing images covering the medicinal herb orchard area; calculating the vegetation index time series pixel by pixel and dividing the area into plot units according to plot boundaries; performing least-squares fitting on the vegetation index time series of each pixel using an annual periodic harmonic model that uses the same set of harmonic coefficients for each year, taking the fitted value as the perennial woody component, and taking the difference between the vegetation index time series and the fitted value as the annual herb component; fitting the annual herb component with an asymmetric Gaussian pulse mixture model to obtain the peak phase, rising scale, and decay scale of each crop, and obtaining the crop overlap by dividing the difference between the sum of the scales of two adjacent crops and the difference between their peak phases; fitting the first-order harmonics of the perennial woody component year by year to obtain the annual phase, and taking the annual phase... The average modulus of the unit complex vector corresponding to the position is used as the interannual phase consistency; one-dimensional power spectrum analysis is performed on the growth period integral field of the annual herb component in multiple directions within the plot unit, and the maximum value of the ratio of the peak value of each direction spectrum to the mean value of the spectrum in that direction is taken as the directional modulation significance, and the reciprocal of the corresponding peak frequency is taken as the strip period. The grain-economic crop bandwidth ratio is obtained by dividing the inverse cosine of the ratio of the amplitude of the second harmonic to the fundamental frequency in that direction by pi; the product of the interannual phase consistency, the directional modulation significance and the normalized crop overlap is taken as the joint discriminant value to determine the grain-economic crop strip intercropping and distinguish between intercropping, relay cropping and rotation, and output the mode category, grain-economic crop bandwidth ratio and crop rotation period; anisotropic Gaussian kernels are constructed with the direction corresponding to the spectral peak as the strip direction to perform neighborhood weighted smoothing on the vegetation index time series and iterate alternately.

[0008] This invention also provides a system for identifying strip-shaped intercropping of medicinal herbs, fruit trees, and cash crops, comprising: a time series construction and plot unit division module for constructing a vegetation index time series and dividing plot units; a woody herbaceous component decomposition module for separating perennial woody components and annual herbaceous components by fitting a shared harmonic coefficient across years; a crop rotation analysis module for analyzing the parameters of each crop rotation and constructing the crop rotation overlap using an asymmetric Gaussian pulse mixture model; a cross-year phase consistency verification module for determining the affiliation of woody crop parks based on the unit complex vector average modulus of the phases each year; a directional modulation analysis module for outputting directional modulation significance, strip period, strip direction, and inverting the ratio of cash crop to cash crop bandwidth; and a joint discrimination and feedback output module for jointly determining the pattern category based on three criteria and iteratively feeding back the anisotropy kernel with the strip direction. Each of the above modules corresponds one-to-one with each step of the aforementioned method.

[0009] Firstly, this invention makes the strip intercropping of grain and cash crops identifiable. This effect stems from two technical features in the independent claims: separating the two components using a least-squares fitting of an annual cycle harmonic model that uses the same set of harmonic coefficients across years; and constructing the crop overlap degree by the difference between the sum of the scales of two adjacent crops and the peak phase. The mechanism lies in the fact that perennial woody plants recur in the same phase over many consecutive years, while annual herbs vary from year to year and exhibit seasonal pulses. The two are naturally separated in the temporal spectrum, thus requiring no end-member spectral analysis for decomposition. After decomposition, the agronomic design parameter of crop overlap directly manifests as the relationship between the phase difference of adjacent pulses and the sum of scales, which is precisely the defining characteristic that distinguishes the intercropping pattern from monoculture.

[0010] Secondly, this invention makes the grain-to-wave bandwidth ratio, a sub-pixel scale quantity, measurable. This effect stems from the technical feature of obtaining the grain-to-wave bandwidth ratio by dividing the inverse cosine of the ratio of the second harmonic to the fundamental frequency by pi. The mechanism is that for a periodic strip whose width is the product of the strip period and the duty cycle, the magnitude of each order of the Fourier coefficients of its spatial waveform is proportional to the sine of the product of the duty cycle and the order. The ratio of the second order to the first order is exactly equal to the cosine of the product of the duty cycle and pi. Therefore, the duty cycle can be directly analytically obtained from this ratio using the inverse cosine. This relationship does not require the strip to be spatially resolvable, only that the strip period falls within the frequency band allowed by the sampling theorem.

[0011] Third, this invention exhibits a significant synergistic effect. A single cross-year harmonic decomposition simultaneously generates three criterion branches: the woody component provides cross-year phase consistency to determine orchard affiliation; the herbaceous component's pulse sequence provides crop overlap to determine intercropping, relay cropping, and rotation; and the herbaceous component's growth period integral field provides directional modulation significance and the ratio of grain to cash crop bandwidth to determine zonal structure. These three are combined through product rather than weighted sum, corresponding to the definition of compound planting where the three conditions of woody coexistence, zonal existence, and crop overlap must simultaneously be met. Furthermore, the zonal direction obtained from directional modulation analysis is fed back for anisotropic smoothing of the time series, synchronously improving the signal-to-noise ratio of all three branches. The existence of any single component alone is insufficient.

[0012] Fourth, this invention conveniently eliminates the need for topographic radiation correction in mountainous areas. This effect stems from the structural property that all the aforementioned criteria are phase quantities or ratios of the same dimension. The mechanism lies in the fact that the radiation distortion caused by the slope aspect of mountains is approximately represented as a multiplicative factor in the vegetation index time series, while the harmonic coefficients of the least squares fitting scale proportionally with this factor, the phase given by the four-quadrant arctangent remains unchanged, the peak phase and rise decay scale of the impulse fitting remain unchanged, and the peak-to-average power ratio and harmonic amplitude ratio of the power spectrum remain unchanged; therefore, all the criteria of this invention are strictly insensitive to this multiplicative factor. This is particularly important in steep mountainous areas because the residual of the topographic radiation correction itself is of the same order of magnitude as the actual phenological differences, and existing amplitude-based methods must be corrected first, yet remain unreliable even after correction. This invention is suitable for applications such as verifying the area of ​​grain-cash crop intercropping in mountainous medicinal herb gardens and orchards, surveying model categories, monitoring the conversion of arable land to non-grain crops, and evaluating the effectiveness of demonstration base promotion. Attached Figure Description

[0013] Figure 1 This is a flowchart of the method for identifying strip-shaped intercropping of medicinal herbs, orchards, grain crops, and cash crops according to an embodiment of the present invention.

[0014] Figure 2 This is an architectural diagram of the strip-shaped intercropping identification system for medicinal herbs, orchards, grains, and cash crops according to an embodiment of the present invention. Detailed Implementation

[0015] The invention will be further described below with reference to the accompanying drawings and embodiments. This embodiment focuses on three specialty economic crop parks within Liangshan Yi Autonomous Prefecture, Sichuan Province: an Artemisia argyi park at an altitude of 1650m with a row spacing of 3.2m, where potatoes and corn are intercropped in strips between rows; a pomegranate park at an altitude of 1180m with a row spacing of 4.0m, where early rice, chili peppers, and potatoes are intercropped between rows; and an alpine blueberry park at an altitude of 2450m with a row spacing of 3.0m, where potatoes, buckwheat, and radishes are intercropped between rows. The slopes of the three parks range from 8 to 26 degrees, encompassing both sunny and semi-shaded slopes. Remote sensing data utilizes multispectral imagery from four complete years (2021-2024), with a spatial resolution of 2.0m, a revisit period of 5 days, and 32 to 41 effective observations per year after cloud masking.

[0016] Before proceeding with each step, let's clarify the correspondence between several terms used in this paper and standard terms in this field. The perennial woody component and annual herbaceous component referred to in this paper correspond to the distinction between evergreen background signals and seasonal vegetation signals in remote sensing phenology. The difference lies in the fact that this paper presents an explicit decomposition method using the fitting of shared harmonic coefficients across years, rather than relying on endmember spectra. The interannual phase consistency referred to in this paper, i.e., the average composite vector length in circular statistics, is a standard statistical quantity. The crop overlap referred to in this paper is a remote sensing measurable expression of the concept of crop rotation in agronomy, explicitly defined by pulse parameters. The directional modulation saliency referred to in this paper, i.e., the peak-to-average power ratio of the spatial power spectrum, is a standard signal processing quantity. The grain-to-earth bandwidth ratio referred to in this paper corresponds to the duty cycle of a periodic square wave in signal processing. The one-dimensional power spectrum, fast Fourier transform, Bayesian information criterion, and Levenberg-Marquardt method referred to in this paper all use their standard meanings in their respective fields.

[0017] like Figure 1 As shown, the method in this embodiment includes steps S1 to S6.

[0018] Step S1: Vegetation index time series construction and plot subdivision. After radiometric calibration, atmospheric correction, and geometric fine correction of each image, the normalized vegetation index is calculated pixel by pixel to obtain the vegetation index time series. ,in and Let be the row and column numbers of the pixels, both of which are positive integers. The observation time is expressed in days per year (d). Masking of clouds, cloud shadows, and snow is determined using a combination of mass band and spectral thresholding. Masked pixels are recorded as missing measurements in the time series without interpolation filling, because subsequent harmonic and impulse fitting are solved using least squares on the effective observations; missing measurements only affect the sample size and do not introduce spurious values. It should be noted that topographic radiation correction is not performed in this step; the reason for this is explained in the invariance argument at the end of step S2.

[0019] The plot units were divided using multi-scale segmentation based on the composite image of the highest annual vegetation index. The segmentation scale parameter was set to 60, the shape factor to 0.3, and the compactness to 0.5. This set of parameters ensured that the average area of ​​the segmented objects was comparable to the measured average area of ​​0.62 hm² in the local park. The segmentation results were then overlaid with slope and aspect layers derived from the digital elevation model. When the aspect difference between adjacent segmented objects exceeded the aspect tolerance of 45 degrees, they were split into different plot units. This is because different aspects result in different radiation conditions, and the park's row orientation typically follows contour lines and changes with slope aspect. Consistent row orientation must be maintained within the same plot unit; otherwise, two spurious spectral peaks will appear in subsequent directional spectral analysis. Segmented objects with an area smaller than the minimum plot area of ​​0.05 hm² were merged into the nearest adjacent plot unit in terms of spectral distance to avoid insufficient spectral resolution due to small fragments.

[0020] Regarding the boundary conditions for plot unit size: For the one-dimensional power spectrum to form a discernible peak at the strip frequency, the plot unit must have a span perpendicular to the strip direction of at least three strip periods. In this embodiment, the strip period is 3.0m to 4.0m, and the pixel size is 2.0m; therefore, the span must be at least 12.0m, or 6 pixels. Plot units that do not meet this condition will not undergo directional spectral analysis and will be directly labeled as plots with insufficient spatial scale. Only a temporal criterion will be provided as a reference conclusion, which will be listed separately in the results and not included in the area statistics to avoid unreliable determinations of total pollution.

[0021] Step S2: Decomposition of woody and herbaceous components sharing harmonic coefficients across years. This step splits the vegetation index time series of mixed pixels into two branches, without phenological normalization or endmember spectral unmixing. The rationale is that the phenology of perennial woody crops is determined by species and latitude, recurring in phase on the same day of the year over consecutive years; the sowing date of annual food crops is determined by agricultural arrangements and crop rotation design, varying year by year, and its signal exhibits a distinct seasonal pulse pattern. Therefore, the two are naturally separable in the temporal spectrum: the former is a strictly annual cycle component, and the latter is a non-repeating pulse superimposed on it across years. Based on this, an annual cycle harmonic model is constructed. , And make continuous All valid observations for each year are combined into a system of equations, and the same set of harmonic coefficients is used for all years to solve the equations using least squares: , in: Let be the perennial woody component of this pixel, and be a scalar with a value range of . , dimensionless, is calculated from this formula, and is the sole source for determining the affiliation of the park; For constant terms, and For the first The coefficients of the cosine and sine are both scalars and dimensionless, and are obtained by the least squares solution described above. The harmonic order is set to 2, which retains both the annual and semi-annual cycles. The value is based on the fact that the phenological curve of woody crops can only show two transitions, leaf expansion and leaf fall, within an annual cycle. A harmonic order of 1 cannot depict the summer plateau of evergreen woody plants, while a harmonic order of 3 or higher will begin to absorb the energy of herbaceous pulses, resulting in incomplete separation of the two branches. In this embodiment, the measured pulse amplitude of the herbaceous component decreased by 21% when the order was 3 compared to when the order was 2. The cycle is annual, taking 365.25 days, with the unit being days. For the first The number of years of a single effective observation, expressed in days (d). For the first The effective observation set for the year; The number of years participating in the joint fitting is set to 4. The key point of this formula is that all years share the same set of coefficients: mathematically, it forces only components that repeat across years to be classified as woody components, thus achieving separation without any endmember spectra.

[0022] The weight of annual herbs is taken as residuals: ,in: This represents the annual herbaceous component of the pixel, which is a scalar quantity. In this embodiment, the value range is... to Dimensionless, calculated by this formula, it is the common input for both stubble analysis and directional modulation analysis. The following demonstrates that all criteria of this invention are insensitive to multiplicative radiation distortion caused by mountain slope aspect, which is also the basis for not performing topographic radiation correction in this step. Let the time series of the true vegetation index of a certain pixel be... The observation timeline after being affected by slope aspect is as follows ,in Let be a multiplicative distortion factor, be a positive scalar, dimensionless, and approximately constant within the same plot unit. Due to the linearity of least squares, for... The fitted coefficients satisfy That is, all coefficients are scaled proportionally; thus, the phase is given by the first-order coefficients through the arctangent of the four quadrants. Remains unchanged. Accordingly, the herbaceous component satisfies The peak amplitude obtained from its pulse fitting is scaled proportionally, while the peak phase and the rising decay scale remain unchanged; the overlap of the stubble formed by the peak phase and scale therefore remains unchanged; the integral field of the growing season is scaled proportionally, and its power spectrum is adjusted accordingly. Due to scaling, the ratio of peak-to-average power ratio (PAPR) to the amplitude of each harmonic remains unchanged. In summary, the four criteria used in this invention—cross-year phase consistency, crop overlap, directional modulation significance, and grain-to-economic bandwidth ratio—are either phase quantities or ratios of the same dimension, and are applicable to... Strictly insensitive. It should be added that when... When the power spectrum is not constant within the plot unit but varies gradually along the slope, its spatial scale is much larger than the strip period. It only affects the low-frequency end of the power spectrum and does not affect the ratio of the fundamental frequency to the second harmonic at the strip frequency. Therefore, the above conclusion still holds.

[0023] Step S3: Asymmetric Gaussian Pulse Mixture Decomposition and Crop Overlap. Herbaceous components are often not unimodal within a growing season: potato and corn intercropping in Artemisia argyi fields, early rice and radish intercropping in pomegranate fields, and potato and buckwheat intercropping in blueberry fields can all generate two or even three growth pulses in the same year. Under intercropping conditions, the decay phase of the preceding crop and the rising phase of the following crop coexist, causing partial overlap between the two pulses. If a simple extreme value detection is used to extract the phenological point, the overlapping two crops will be merged into one. Therefore, the herbaceous component is represented as the sum of several asymmetric Gaussian pulses: , in: For the first The peak phase of a crop crop, i.e., the accumulated days of the crop's growth reaching its maximum, is a scalar quantity. In this embodiment, the value ranges from 95 days to 285 days, with the unit being days. It is obtained by fitting this formula and is the direct source of the crop rotation period. and These are the rising and falling scales for the crop, respectively. They are scalars, ranging from 8d to 40d, with units of d. Different values ​​of these two scales constitute the asymmetry of the pulse. This asymmetry is necessary because the vegetative growth period of crops is usually longer than the maturity and harvest period. If forced symmetry is used, the fitting will sacrifice the peak phase to accommodate the width on both sides. In this embodiment, the root mean square error of the peak phase increased from 3.6d to 11.2d when forced symmetry was used. This represents the peak amplitude of this crop, which is dimensionless. The number of crop rotations is a positive integer, automatically selected between 1 and 4 by the Bayesian information criterion. The penalty term coefficient is calculated as the product of the number of parameters and the logarithm of the number of observations; this criterion is well-known in the field. All parameters are solved using the Levenberg-Marquardt method, with initial values ​​given by the local maxima of the herbaceous component. The overlap between adjacent crop rotations is constructed using the following formula: ,in: For the first Stubble and First The overlap between stubble joints is a scalar quantity with a value range of [value range missing]. The dimensionless value is calculated from this formula and serves as the discriminant for the three crop rotation patterns: intercropping, relay cropping, and crop rotation. The numerator is the sum of the decay scale of the previous crop and the rise scale of the subsequent crop minus the difference in the peak phase of the two crops. Its sign directly reflects whether the growth periods of the two crops overlap. The denominator is the sum of the same scale, which serves as a normalization function to make this value independent of the length of the crop growth period.

[0024] The discrimination rule is: A value greater than zero indicates intercropping with a symbiotic period, meaning the subsequent crop has already emerged before the previous crop has been harvested; Tolerance range near zero The inner judgment is that it is an intercropping with a stubble immediately following the first stubble; Crops smaller than the lower bound of this interval are considered to be in successive crop rotation. The width of the tolerance interval is determined by the field survey plots: in this embodiment, among the 46 known pattern plots, the intercropping plots... Distribution concentrated in to The sample plots are concentrated in to The rotation plots are concentrated in to The tolerance range is calculated by taking the 95% coverage area of ​​the intercropping distribution and leaving a 0.02 margin on both sides. If the range is too wide, a large number of relay cropping and crop rotation crops will be classified as intercropping; if it is too narrow, the same pattern will jump between the three categories in different years due to slight adjustments in the sowing date. When adding new field plots, they are classified according to their... The quantiles of the distribution are recalculated for this interval.

[0025] The crop rotation period is taken as the difference between the peak phase of two adjacent crops and the sum of the scales of the two crops, i.e. The unit is d. When the value is negative, it indicates the number of days that two crops coexist; when it is positive, it indicates the number of days that there is no crop rotation between the two crops. It can be directly used to assess the tightness of crop rotation.

[0026] Step S4: Cross-year phase consistency test of perennial woody components. A pixel containing multiple cropping of herbaceous plants is insufficient to prove that the pixel is located in an intercropping area; the same applies to double-cropping rotation on pure cultivated land. The distinction should ideally be based on whether the pixel simultaneously contains perennial woody components. However, directly judging by the amplitude of woody components can lead to an unforeseen reverse failure: newly planted alpine blueberries have a canopy closure of less than 0.25 in the first three years after planting, and their woody component amplitude is only one-quarter to one-third of that in mature areas, making them difficult to distinguish from the stubble background of rotated cultivated land. Therefore, newly established demonstration areas are more likely to be missed as cultivated land, and these are precisely the areas most in need of monitoring coverage.

[0027] The solution is to focus on the phase's repetitiveness across years, rather than the amplitude: For perennial woody plants, regardless of canopy closure, leaf expansion and leaf fall occur around the same day of the year; the sowing date for annual crops varies from year to year. Based on this, a first-order harmonic is fitted to each year, and its phase is determined: , in: For this pixel in the th The phase of the woody component in a given year is a scalar with a range of values ​​of [value range missing]. The unit is rad, which is calculated by this formula. Its physical meaning is the annual cycle phase corresponding to the peak value of woody phenology in that year. and They are respectively only using the first The first-order cosine coefficients and first-order sine coefficients obtained from annual observation fitting are both dimensionless. This is a four-quadrant arctangent function, with the first independent variable being the ordinate and the second independent variable being the abscissa, and its return range is... Here, a four-quadrant arctangent must be used instead of a regular arctangent; otherwise, the phase will fold in the opposite direction in the second and third quadrants, resulting in falsely low inter-year consistency values. Inter-year phase consistency is calculated using the average composite vector length from circular statistics. ,in: The inter-year phase consistency of this pixel is denoted by , which is a scalar and has a value range of . , dimensionless, is calculated by this formula. A value close to 1 indicates that the phase of each year is highly concentrated, that is, the pixel contains perennial woody plants, while a value close to 0 indicates that the phase of each year is diffuse, that is, it does not contain perennial woody plants. The imaginary unit; The module length is 1 and the argument angle is 1. unit complex vector; The modulus of a complex number; The definition is the same as the aforementioned formula. The key property of this formula is that the contribution each year is normalized to a complex vector with a magnitude of one, thus this quantity is completely independent of the amplitude of the woody component. This is precisely why newly planted areas with low canopy closure can still be correctly identified. The measured distribution in this embodiment is: mature pomegranate orchards. The average annual variability (AGR) ranged from 0.88 to 0.96 in mature Artemisia argyi orchards, 0.79 to 0.91 in newly planted blueberry orchards (two years old), and 0.74 to 0.86 in pure cultivated land with a two-crop-a-year rotation system, with a significant gap between the two distributions. Based on this, a normalized scaling parameter for interannual phase consistency was established. Taking 0.55 means making The normalized value of this criterion is 0.5; this parameter is calibrated by the field sample plots according to the equal error rate points of the subject operating characteristic curve. If the value is too high, the newly planted area is excluded, and if it is too low, the rotation farmland is included.

[0028] Boundary cases regarding year numbers: Phase consistency requirements across years The average modulus of the unit complex vector must be at least 3; otherwise, the average modulus of the two vectors is statistically insufficient to distinguish between concentration and dispersion. When the angle between the two vectors is less than 120 degrees, the average modulus is always greater than 0.5, and the criterion loses its discriminative power. For demonstration bases established less than three years ago, this embodiment does not use this criterion, but instead uses a combined determination of directional modulation significance and crop overlap, as described later. The resulting plot is marked as a short-sequence determination plot. Simultaneously, if there are fewer than 12 valid observations in a given year, that year is not included in the annual phase fitting. Subtract one and record it accordingly.

[0029] Step S5: Directional sub-pixel modulation spectrum analysis and grain-crop bandwidth ratio inversion. At this point, one key attribute remains unutilized: the strip structure of the grain-crop intercropping. In the medicinal herb orchard, the row spacing is 3.0m to 4.0m, and the width of the grain strips between rows is only 1.5m to 2.5m. However, the pixel size of the image in this embodiment is 2.0m, meaning a single strip does not span two pixels and is indistinguishable in the spatial domain. But the indistinguishability of the strips in the spatial domain does not mean their information has been lost: the periodicity of the strips causes the area ratio of grain crops in adjacent pixels to fluctuate periodically. This fluctuation, as directional sub-pixel modulation, is completely preserved in the local spatial spectrum.

[0030] First, the modulated carrier field is constructed. The integral of the annual herbaceous component over its growth period is taken and normalized using an annual cycle: ,in: Let be the integral field of the growth period of this pixel. Let be a scalar, and in this embodiment, its value ranges from 0.004 to 0.096. It is dimensionless and is calculated by this formula. Its value is proportional to the product of the area ratio of grain crops in this pixel and their growth status. Therefore, it is the spatial sampling of the strip duty cycle. and The starting and ending dates of the reproductive period integral are calculated by taking the 5% quantile of the first rise starting point and the 95% quantile of the last decay ending point of all pixels in the plot unit, respectively, in days. The definition is the same as the formula mentioned above. The purpose of normalizing with an annual period is to make the field dimensionless, so that the subsequent peak-to-average power ratio and harmonic amplitude ratio do not introduce dimensions. Dimensional verification: the integrand is dimensionless, the integral variable is in d, the unit of the integral result is d, and it is dimensionless after dividing by the annual period in d.

[0031] One-dimensional power spectrum analysis along the direction: Within the plot unit, traverse from 0 degrees to 178 degrees with a directional step size of 2 degrees, and perform one-dimensional power spectrum analysis along the direction for each direction. ,Will A one-dimensional sequence is obtained by projecting and averaging along the orthogonal direction of this direction. , For along direction The distance coordinates are in meters; a fast Fourier transform is performed on the sequence, and the square of the modulus is taken to obtain the one-dimensional power spectrum in that direction. Spatial frequency, unit: Projective averaging along orthogonal directions is necessary: ​​the strip is theoretically uniform along its own direction, and projective averaging suppresses random noise in that direction by the square root of the number of pixels, while the periodic signal of the strip remains unaffected. The directional modulation saliency is given by the following formula: ,in: Here, represents the directional modulation saliency, and is a scalar with a value range of . , dimensionless, is calculated by this formula, and its value is the lifting factor of the band spectral peak in the optimal direction relative to the spectral background in that direction, which is the criterion for the existence of band structure; The search interval for the band frequency is determined by the reciprocal of the minimum and maximum row spacing in the medicinal herb orchard area. In this embodiment, we take... The corresponding row spacing range is 2.5m to 4.0m. The purpose of limiting this range is to exclude low-frequency fluctuations caused by plot boundaries, field ridges and terrain shadows, as well as high-frequency fluctuations caused by sensor noise. This represents the arithmetic mean over that interval. The direction in which the maximum value is taken is denoted as the strip direction. The corresponding spectral peak frequencies are denoted as strip period The unit is meters (m). The inversion of the grain bandwidth ratio is based on the harmonic structure of a periodic square wave. Assume the strip has a periodicity in the direction perpendicular to the strip. Duty cycle is The square wave distribution (high values ​​are taken at the grain stripe and low values ​​at the fruit tree row), then its th The modulus of the order Fourier coefficients satisfies Therefore, the ratio of the second order to the first order is... This ratio can be directly read from the power spectrum, since the power spectrum is the square of the amplitude. ,in: The grain-to-cash crop band ratio, i.e., the proportion of the grain crop band width to the band period, is a scalar quantity with a value range of [value missing]. The dimensionless value is obtained by calculation using this formula. Its product with the strip period is the strip width of the grain crop. and These are the power spectral values ​​at the fundamental frequency and second harmonic, respectively, along the strip direction; the square root sign converts the power ratio back to the amplitude ratio. It is an inverse cosine function, whose independent variable is the dimensionless amplitude ratio and whose range is... Dividing by pi gives Duty cycle within the interval. There is an ambiguity here that must be addressed: duty cycle. and Given identical second-order and first-order amplitude ratios, the above formula can only determine the branch less than or equal to 0.5. The disambiguation is based on the level of the integral field during the growth period; herbaceous growth is high in the grain strip and low in the fruit tree row. Therefore, in the one-dimensional sequence along the strip direction, the proportion of segments with values ​​higher than the sequence mean is the estimate of the true duty cycle; if this proportion is greater than 0.5, then the true duty cycle is taken as... In this embodiment, the actual duty cycles of the three zones are 0.44, 0.38, and 0.47, respectively, all below 0.5. The disambiguation step did not trigger inversion, but this step is necessary for the novel wide-row, densely planted pattern. Furthermore, when the signal-to-noise ratio of the fundamental spectrum peak is insufficient, i.e. Normalized scaling parameter below directional modulation significance When the value is 4.0, the spectral value at the second harmonic is dominated by noise. At this time, the bandwidth ratio is not output, but only the strip period is output, and it is marked in the results as bandwidth ratio cannot be inverted.

[0032] Step S6: Joint discrimination based on three criteria, pattern classification, and strip direction feedback. At this point, three independent criterion branches have been obtained: inter-year phase consistency. To determine whether a pixel is located within a woody crop area, the directional modulation saliency is... To answer whether the pixel is within a banded structure, the degree of overlap of the cut edges. The question asks whether multiple grain crops are grown between rows and how these crops are connected. The definition of strip intercropping of grain and cash crop requires all three conditions to be met simultaneously; failure to meet any one condition disqualifies the crop from this type: a pure orchard with only woody plants and no striping is a strip-sown monoculture; and crop rotation with only multiple crops and no coexistence of woody plants is a crop rotation. Based on this, a product-type joint discriminant value is constructed: ,in: Let be the joint discriminant value of this pixel, and let be a scalar with a value range of . , dimensionless, is calculated by this formula, and it is the final basis for determining the strip intercropping of grain and cash crops; It is the maximum value among all adjacent stubble pairs of this pixel, and is dimensionless; The normalized scaling parameter for cross-year phase consistency is set to 0.55; The normalized scaling parameter for directional modulation saliency is set to 4.0, where... The reason for subtracting one is that the theoretical lower limit of the peak-to-average power ratio (PAPR) without modulation is 1, not 0. As the discriminant center for the degree of overlap of the stubble, take This refers to the lower limit of the intercropping tolerance range, which means that intercropping and relay cropping are included in the compound planting, while pure crop rotation with successive cropping is not included. The discriminant softening width for crop overlap is set to 0.06, dimensionless, and taken as half the gap width between the intercropping and rotation quadrats. If the value is too small, the discriminant value of samples within the gap will jump between 0 and 1, becoming sensitive to interannual fluctuations in sowing date; if it is too large, rotation will be largely included in the intercropping calculation. The three terms are multiplied rather than weighted summed because weighted summing allows a high value in one term to compensate for a zero value in another. For example, if a pure orchard accidentally shows a high value due to plot boundary effects... Even with weighted summation, it may still exceed the threshold, while with product, it may exceed the threshold. Approaching zero causes the entire discriminant value to approach zero. Greater than the discrimination threshold The pixel with a value of 0.42 was identified as a grain-cash crop strip intercropping. Along with the aforementioned three normalized scaling parameters, the data was calibrated using field quadrats: A total of 46 quadrats with confirmed planting patterns were established in the three parks (28 intercropping quadrats, 9 pure orchards, 5 strip-seeded monoculture fields, and 4 crop rotation quadrats). The parameters were determined based on the 3D feature distribution of pixels within the quadrats according to the equal error rate points of the receiver operating characteristic (ROC) curve. In this embodiment, the equal error rate was 6.8%. When adding new field quadrats, the discrimination threshold was determined according to... Perform an index-weighted update, where The threshold obtained solely from the calibration of newly added quadrats. To update the step size, we set it to 0.25, which is dimensionless. This value makes the equivalent memory length approximately four batches of quadrats, which is comparable to the frequency of the annual survey.

[0033] The pattern category is determined by a combination of crop overlap and strip period: A value greater than 0.08 is considered a plagiarism, while a value falling within... Internal judgment is intermittent, less than However, this pixel still meets the first two criteria for being classified as a crop rotation composite crop; simultaneously, it is determined whether the strip period falls within the nominal row spacing of the cash crops in the park. Within the interval, cross-validation is performed; those falling outside the interval are marked as having abnormal row spacing and prompted for review. The output includes the mode category, the grain-to-cash crop bandwidth ratio, the crop rotation period obtained from step S3, and the grain crop band width given by the product of the grain-to-cash crop bandwidth ratio and the band period.

[0034] Finally, feedback is performed along the strip direction. The signal-to-noise ratio (SNR) of the aforementioned directional modulation analysis is limited by timing noise. The conventional approach to improving SNR, isotropic neighborhood smoothing, inadvertently mixes adjacent strips into the same smoothing window, thus erasing the modulation being measured. The larger the smoothing window, the more thorough the modulation attenuation. The correct approach is to perform smoothing only along the strip direction: , in: Let be the weights of the anisotropic Gaussian kernel, and let be a scalar with a range of values. , dimensionless, is obtained by calculation using this formula; and These represent the neighboring pixels relative to the center pixel along the strip direction. The distance component perpendicular to this direction, in meters; The core bandwidth along the strip direction is taken as 1 to 3 times the strip period, i.e., 3.0m to 12.0m, in meters. Its upper limit is limited by the length of the plot unit along the strip direction. The core bandwidth is perpendicular to the strip direction and is no greater than one-quarter of the strip period, i.e., no more than 1.0m. The unit is m. This upper limit is the key to this formula. When it exceeds one-quarter of the period, the adjacent strips begin to enter the smoothing window. The modulation amplitude decays according to the sine function of the ratio of the core bandwidth to the period. When it is half the period, the modulation is completely canceled. The vegetation index time series is smoothed and is dimensionless. For the set of neighboring pixels, take the central pixel as the center and the semi-major axis as the semi-major axis. semi-short shaft All pixels within the ellipse. Steps S2 to S5 are re-executed with the smoothed time sequence to obtain new strip directions; the process terminates when the absolute value of the difference between the strip directions obtained from two adjacent iterations is less than the direction convergence threshold of 2 degrees, or when the number of iterations reaches the upper limit of 5 iterations, and the result at termination is used as the final output. In this embodiment, the typical number of convergence iterations for the three parks is 2 to 3. After iteration, the directional modulation significance increases from 5.7 in the first iteration to 9.4, and the estimated standard deviation of the grain-to-bandwidth ratio decreases from 0.061 to 0.023. It should be noted that if the first iteration yields... That is, lower than If the direction of the strip is unreliable, then anisotropic smoothing is not performed and the strip is treated directly as a stripless structure to avoid artificially creating false strip modulation by using the wrong direction for directional smoothing.

[0035] To verify the effectiveness of this invention, three control groups were established within a 1200-mu (approximately 80 hectares) area encompassing three orchards: Artemisia argyi, pomegranate, and alpine blueberry. Control group 1 was implemented according to the technical route of CN115661554A, employing a time-aligned spatiotemporal coding classification model; Control group 2 was implemented according to the technical route of CN108399400A, employing multi-source synthetic temporal sequence, feature optimization, and a classifier; the experimental group was implemented according to this invention. All three groups used identical imagery and the same field quadrats, with 60% of the quadrats used for training or calibration and 40% used for accuracy verification. Control groups 1 and 2 underwent topographic radiometric correction according to their respective technical routes, while the experimental group did not. The results are shown in the table below.

[0036]

[0037] As shown in the table, the overall accuracy of binary classification in the experimental group improved from 68.9% to 71.3% in the control group to 92.6%, and the Kappa coefficient improved from 0.39 to 0.43 to 0.85. Almost all errors in the two control groups stemmed from misclassifying intercropping as pure orchards. This was rooted in the fact that aligning the time dimension eliminated the phase difference between crop rotations, resulting in a dulled single-peak mixing sequence. The three classifications of intercropping, relay cropping, and crop rotation, as well as the three outputs of grain-cash crop bandwidth ratio and crop rotation period, were completely unavailable in the two control groups; their output granularity remained at the land type level. The missed detection rate in newly planted blueberry orchards two years after planting decreased from 58.2% to 63.6% to 9.1%, confirming the necessity of the design that uses phase rather than amplitude to determine the woody plant classification. The accuracy difference between sunny and semi-shaded slopes narrowed from 12.1 to 14.8 percentage points to 2.3 percentage points. Furthermore, the experimental group did not perform topographic radiation correction, which experimentally verified the theoretical argument that the criteria of this invention are insensitive to multiplicative radiation distortion. In terms of computational overhead, the experimental group consumed approximately three times more time per square kilometer than the control group, mainly due to two nonlinear fitting operations per pixel and multi-directional spectral scanning per plot. However, the entire process for a 1200-acre area took less than 12 minutes on a conventional workstation. Moreover, this processing was an annual survey task rather than a real-time task, and therefore did not constitute an application bottleneck.

[0038] Furthermore, this invention produced an unexpected effect in its application: the crop-to-cash crop bandwidth ratio obtained from the inversion showed a systematic variation along the slope within the same planting area, with the bandwidth ratio at the upper part of the slope being on average 0.07 lower than that at the lower part. Field verification revealed that this was because farmers actively narrowed the crop strips at the upper part of the slope to reduce soil erosion; this subtle agronomic difference had never been recorded by any monitoring methods. This indicates that the crop-to-cash crop bandwidth ratio output by this invention is not only a verification indicator but can also serve as an observation of agronomic behavior, used to evaluate the actual implementation of intercropping technical regulations under different site conditions.

[0039] like Figure 2 As shown, this embodiment also provides a strip-shaped intercropping identification system for medicinal orchards, grain and cash crops, including a time-series construction and plot unit division module, a woody herbaceous component decomposition module, a crop rotation analysis module, a cross-year phase consistency test module, a directional modulation analysis module, and a joint discrimination and feedback output module. The six modules correspond one-to-one with steps S1 to S6 of the aforementioned method embodiment. Their internal algorithms and parameter values ​​are the same as those of the aforementioned method embodiment and will not be repeated.

[0040] The temporal series construction and plot unit division module takes as input multi-year, multi-temporal, multispectral imagery covering the medicinal herb and fruit orchard, a digital elevation model, and the orchard boundary vector. After radiometric calibration, atmospheric correction, geometrical fine correction, and cloud / snow masking, this module calculates the normalized vegetation index (NDI) pixel-by-pixel, outputting a vegetation index time series with annual day as the time coordinate. Simultaneously, it performs multi-scale segmentation using composite images of the highest annual vegetation index values, overlays them with slope and aspect layers, splits them according to aspect tolerance, and merges them according to the smallest plot area, outputting plot unit vectors and their spans along the vertical strip direction. This information is then used by the directional modulation analysis module to determine whether the spatial scale requirements for spectral analysis are met. This module explicitly excludes topographic radiometric correction.

[0041] The woody herbaceous component decomposition module is connected to the time series construction and plot unit division module. For each pixel, it merges all valid observations from consecutive years into a set of equations, solves it using a second-order annual cycle harmonic model that uses the same set of harmonic coefficients for all years, and outputs the perennial woody component and the annual herbaceous component as the residual. Simultaneously, this module fits the first-order harmonics for each year and outputs the first-order cosine and sine coefficients for each year, which are used by the inter-year phase consistency verification module. The two outputs of this module constitute the common upstream of the subsequent three criterion branches.

[0042] The crop rotation analysis module is connected to the woody herbaceous component decomposition module. After determining the number of asymmetric Gaussian pulses using the Bayesian information criterion, it fits the peak amplitude, peak phase, rising scale, and decay scale of each crop using the Levenberg-Marquardt method, outputs the crop rotation overlap of each adjacent crop pair and its maximum value, and gives a preliminary judgment of intercropping, crop rotation, or crop rotation according to the positive or negative value and tolerance range of the overlap. The module also outputs the crop rotation connection period obtained by subtracting the difference in peak phase from the sum of the scales.

[0043] The inter-year phase consistency test module is connected to the woody herbaceous component decomposition module. It calculates the phase of the first harmonic coefficient of each year using the arctangent of the four quadrants, maps the phase of each year to a unit complex vector, takes the magnitude of the sum of the vectors and divides it by the number of years, and outputs the inter-year phase consistency. When the number of valid years is less than three years, the module sets the short sequence flag and does not output the criterion. Instead, the joint discrimination and feedback output module uses two criteria to determine the consistency.

[0044] The directional modulation analysis module is connected to the woody herbaceous component decomposition module and the time series construction and plot unit division module. Within the plot unit that meets the spatial scale requirements, the growth period integral field of the annual herbaceous component is traversed by projection averaging and fast Fourier transform with a step size of 2 degrees. The spectral peak is taken within the frequency band defined by the row spacing range of the park, and the directional modulation significance, strip period and strip direction are output. The bandwidth ratio is obtained by inverting the square root of the ratio of the second harmonic and the fundamental frequency power spectrum values ​​through inverse cosine. At the same time, the duty cycle is disambiguated according to the proportion of the growth period integral field that is higher than its mean value.

[0045] The joint discrimination and feedback output module is connected to the aforementioned three criterion modules. It calculates the joint discrimination value based on the product of the three criteria after normalization and compares it with the discrimination threshold. It outputs the pixel-level determination of grain-crop strip intercropping, pattern category, grain-crop bandwidth ratio, grain crop strip width and crop rotation period. It also calculates the area of ​​various intercropping types according to pattern category and generates a distribution thematic map. This module also constructs an anisotropic Gaussian kernel based on the strip direction given by the directional modulation analysis module. It feeds back to the input of the woody herbaceous component decomposition module to perform neighborhood weighted smoothing of the vegetation index time series, forming an alternating iterative closed loop of direction, smoothing and spectral analysis. It terminates when the direction converges or the number of iterations reaches the upper limit.

[0046] The system also includes a parameter calibration module connected to the joint discrimination and feedback output module. This module determines three normalized scaling parameters and a discrimination threshold based on the three-dimensional criterion feature distribution of pixels within the field survey quadrats, according to the equal error rate points of the receiver operating characteristic (ROC) curve. The discrimination threshold is updated exponentially when a new quadrat is added. The system runs on a workstation equipped with a multi-core processor, and both the spectral analysis between plot units and the two nonlinear fitting operations between pixels are distributed in parallel.

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

Claims

1. A method for identifying strip-planted intercropping of medicinal herbs, orchards, grains, and cash crops, characterized in that... include: Acquire multi-year, multi-temporal remote sensing images covering the medicinal herb and fruit orchard, calculate the vegetation index time series pixel by pixel, and divide the plots into plot units according to the plot boundaries; For the vegetation index time series of each pixel, the least squares fit is performed using the annual cycle harmonic model that uses the same set of harmonic coefficients for each year. The fitted value is taken as the perennial woody component, and the difference between the vegetation index time series and the fitted value is taken as the annual herb component. The annual herb component is fitted with an asymmetric Gaussian pulse mixture model to obtain the peak phase, rising scale, and decay scale of each crop. The crop overlap is obtained by dividing the difference between the sum of the scales of two adjacent crops and the difference between their peak phases by the sum of the scales. The first harmonics of the perennial woody components are fitted to each year to obtain the annual phase, and the average magnitude of the unit complex vector corresponding to each year's phase is taken as the inter-year phase consistency. Within the plot unit, a one-dimensional power spectrum analysis is performed on the growth period integral field of the annual herb component along multiple directions. The maximum value of the ratio of the peak value of each direction spectrum to the mean value of the direction spectrum is taken as the directional modulation significance, and the reciprocal of the corresponding peak frequency is taken as the strip period. The grain-earth bandwidth ratio is obtained by dividing the inverse cosine of the amplitude ratio of the second harmonic to the fundamental frequency in that direction by pi. The product of the cross-year phase consistency, the directional modulation salience, and the normalized crop overlap is taken as the joint discrimination value. Pixels with a value greater than the discrimination threshold are judged as grain-cash crop strip intercropping. Intercropping, relay cropping, and crop rotation are distinguished according to the sign of the crop overlap and the strip period. The output mode category, grain-cash crop bandwidth ratio, and crop transition period determined by the difference of the peak phase are output. Using the direction corresponding to the maximum value as the strip direction, an anisotropic Gaussian kernel with a kernel bandwidth greater than that in the vertical direction is constructed. After performing neighborhood weighted smoothing on the vegetation index time series, the above steps are repeated until the change in the strip direction between two adjacent changes is less than the direction convergence threshold.

2. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grain crops, and cash crops according to claim 1, is characterized in that, The least-squares fitting of the annual cycle harmonic model that uses the same set of harmonic coefficients for all years includes: An annual harmonic model containing constant terms and the first two sine and cosine terms is constructed with the accumulated days of the year as the independent variable and one year as the basic cycle. All observations over several consecutive years are combined into a set of equations, and the same set of harmonic coefficients is shared for each year and solved by least squares. In the least squares solution, the vegetation index time series crop phenology for each year is not normalized, nor is the vegetation index time series endmember spectral unmixing performed. Among them, the perennial woody component, the annual herbaceous component, and the peak phase, the interannual phase consistency, the grain-economic bandwidth ratio and the cropping overlap derived from the two are all phase quantities or ratios of the same dimension. When the vegetation index time series is scaled as a whole by the multiplicative factor caused by the mountain slope aspect, its value remains unchanged.

3. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grain crops, and cash crops according to claim 2, is characterized in that... The fitting of annual herb components using an asymmetric Gaussian pulse mixture model includes: The annual herb component is represented as the sum of several asymmetric Gaussian pulses, each of which is determined by four parameters: peak amplitude, peak phase, rising scale before the peak phase, and decay scale after the peak phase. The number of asymmetric Gaussian pulses is determined using the Bayesian information criterion, and all parameters are solved using the Levenberg-Marquardt method. Based on the tolerance range of the overlap degree of the cropping, which is greater than zero, equal to zero, or less than zero, the two adjacent crops are respectively determined to be intercropping with a symbiotic period, intercropping with closely following cropping, and crop rotation with successive cropping.

4. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grain crops, and cash crops according to claim 3, is characterized in that... The method of taking the average magnitude of the unit complex vector corresponding to the phase of each year as the cross-year phase consistency includes: For each year, the first-order sine and cosine terms of the annual periodic harmonic model are individually fitted with all observations within that year, and the inverse of the coefficients of the first-order sine term and the arctangent of the coefficients of the first-order cosine term in the four quadrants are taken as the phase of that year. The phase of each year is mapped to a complex vector with a magnitude of one. The magnitude of the sum of the complex vectors of all years is divided by the number of years to obtain the phase consistency across years. The interannual phase consistency is independent of the amplitude of the perennial woody component. Newly planted areas with a canopy closure lower than the lower limit of canopy closure are still determined to be woody crop areas according to the interannual phase consistency when the amplitude of the perennial woody component is lower than the lower limit of amplitude.

5. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grain crops, and cash crops according to claim 4, is characterized in that, The one-dimensional power spectrum analysis along multiple directions includes: Within the plot unit, the directions are traversed at a set directional step size. For each direction, the fertility period integral field is projected and averaged along the orthogonal direction of that direction to obtain a one-dimensional sequence. The one-dimensional sequence is then subjected to a fast Fourier transform and the square of the modulus is taken to obtain the one-dimensional power spectrum of that direction. The range of the one-dimensional power spectrum is limited to the spatial frequency interval corresponding to the minimum and maximum row spacing of the medicinal herb orchard, and the peak value of the spectrum is taken within this interval. In the direction corresponding to the peak frequency, the ratio of the spectral amplitude at twice the peak frequency to the spectral amplitude at the peak frequency is taken, and after inverse cosine operation and divided by pi, the bandwidth ratio is obtained.

6. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grain crops, and cash crops according to claim 5, is characterized in that... The construction of anisotropic Gaussian kernels for neighborhood-weighted smoothing of vegetation index time series includes: A two-dimensional anisotropic Gaussian kernel is constructed with the strip direction as the main axis, and the kernel bandwidth along the strip direction is set to be one to three times the strip period, and the kernel bandwidth perpendicular to the strip direction is set to be no more than one-quarter of the strip period. The anisotropic Gaussian kernel is used to perform neighborhood weighted smoothing on the vegetation index images at the same observation time. Then, the least squares fitting, the asymmetric Gaussian pulse mixture model fitting, and the one-dimensional power spectrum analysis are re-executed on the smoothed vegetation index time series. The iteration terminates when the absolute value of the difference between the strip directions obtained from two adjacent iterations is less than the direction convergence threshold, or when the number of iterations reaches the upper limit of the number of iterations, and the result at the time of termination is used as the output.

7. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grains, and cash crops according to claim 1, is characterized in that, The division of the medicinal herb and fruit orchard into plot units according to plot boundaries includes: Multi-scale segmentation is used to segment the composite image of the vegetation index maximum value within the year in the remote sensing image, and it is overlaid with the slope and aspect layers derived from the digital elevation model. Adjacent segmented objects whose aspect changes exceed the aspect tolerance are split into different plot units. For subdivided objects with an area smaller than the minimum plot area, they are merged into the adjacent plot unit with the closest spectral distance. For plots spanning more than two slope aspect zones, the one-dimensional power spectrum analysis is performed independently within each slope aspect zone.

8. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grain crops, and cash crops according to claim 1, is characterized in that, Also includes: Within the medicinal herb and fruit park, sample plots with planting patterns confirmed through field investigations were selected. The three-dimensional feature distribution of the pixels within the sample plot, consisting of the inter-year phase consistency, the directional modulation salience, and the crop overlap, was used to determine the normalized scaling parameters and the discrimination thresholds for each plot based on the equal error rate points of the receiver operating characteristic curve. When a new field survey plot is added, the discrimination threshold is updated using an exponentially weighted moving average. When the proportion of pixels within a plot that is identified as being planted in a strip of grain and cash crop intercropping is lower than the lower limit, the entire plot is identified as not being planted in intercropping and is marked as a plot with inconsistent pixel-level identification in the results.

9. The method for identifying strip-planted intercropping of medicinal herbs, orchards, grain crops, and cash crops according to claim 1, characterized in that, Also includes: The number of pixels in the strip-shaped intercropping of grain and cash crops in the medicinal herb and fruit orchard is counted according to the aforementioned pattern category, and then multiplied by the pixel area to obtain the intercropping area of ​​each pattern category. The product of the grain-to-cash crop band width ratio and the band period is used as the grain crop band width, and its mean and standard deviation are statistically analyzed by plot unit. The pattern category, the area of ​​the intercropping, the width of the grain crop strip, and the crop rotation period are written into the thematic output file, and a thematic map of the distribution of grain and cash crop strip intercropping is generated.

10. A strip-planting identification system for medicinal herb orchards and grain-economic crop intercropping, used to implement the strip-planting identification method for medicinal herb orchards and grain-economic crop intercropping as described in any one of claims 1-9, characterized in that, include: The temporal sequence construction and plot unit division module is used to acquire multi-year, multi-temporal remote sensing images covering the medicinal herb and fruit orchard, calculate vegetation index per pixel to obtain the vegetation index temporal sequence, and divide the medicinal herb and fruit orchard into plot units according to the plot boundaries. The woody herbaceous component decomposition module is used to perform least squares fitting on the vegetation index time series of each pixel using an annual cycle harmonic model that uses the same set of harmonic coefficients for each year. The fitted value is taken as the perennial woody component, and the difference between the vegetation index time series and the perennial woody component is taken as the annual herbaceous component. The crop rotation analysis module is used to fit the annual herb components with an asymmetric Gaussian pulse mixture model to obtain the peak phase, rising scale and decay scale of each crop rotation, and to construct the crop rotation overlap between adjacent crops accordingly. The inter-year phase consistency test module is used to fit the first harmonic of the perennial woody component to obtain the phase of each year, and take the average magnitude of the unit complex vector corresponding to the phase of each year as the inter-year phase consistency. The directional modulation analysis module is used to perform one-dimensional power spectrum analysis on the growth period integral field of the annual herbaceous component in multiple directions within the plot unit, outputting directional modulation significance, strip period, strip direction, and inverting the grain bandwidth ratio from the ratio of second harmonic amplitude to fundamental frequency amplitude. The joint discrimination and feedback output module is used to take the product of the cross-year phase consistency, the directional modulation significance and the normalized crop overlap as the joint discrimination value to determine the strip intercropping of grain and cash crop and to distinguish between intercropping, relay cropping and rotation patterns. It outputs the pattern category, the grain-cash crop bandwidth ratio and the crop transition period, and constructs an anisotropic Gaussian kernel with the strip direction to feed back to the woody herbaceous component decomposition module for neighborhood weighted smoothing and alternating iteration.

Images