MATLABGUI-based remote sensing image batch index calculation and trend analysis system and method

Abstract

The invention relates to the technical field of remote sensing image processing and ecological environment monitoring, and discloses a remote sensing image batch index calculation and trend analysis system and method based on MATLABGUI. The system comprises a data acquisition and metadata analysis module, a heterogeneous data adaptive mapping module, a multi-dimensional feature matrix construction module, a time sequence trend mining module and a visual interaction module. Remapping the heterogeneous source data into a standard logic spectrum vector by utilizing a wave band mapping lookup table; calculating ecological characteristic components in parallel by adopting a dynamic calculation graph and a matrix operator, and analyzing and synthesizing a remote sensing ecological index through a principal component corrected by characteristic polarity; and carrying out vectorization Sen slope estimation and a Mann-Kendall significance test. According to the method, through a main thread and background task thread separation architecture, multi-source data compatibility, calculation process self-adaption fault tolerance and efficient mining of a long time sequence trend are achieved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of remote sensing image processing and ecological environment monitoring, in particular to a remote sensing image batch index calculation and trend analysis system and method based on MATLAB GUI. BACKGROUND

[0002] Dynamic change monitoring of ecological environment is the core task of global change research and sustainable development management. Remote sensing technology has become the main means of obtaining surface ecological information due to its macro, rapid and periodic repeated observation advantages. Among them, remote sensing ecological index (RSEI) based on vegetation index, water index, building index and land surface temperature can comprehensively reflect the key ecological indicators such as greenness, wetness, dryness and heat of the region, and is widely used in regional ecological quality assessment and evolution analysis.

[0003] However, in the existing long-term ecological monitoring practice, the data processing and analysis process still faces many technical challenges. First, from the preprocessing of the original remote sensing image to the final time series trend analysis, it usually needs to cross multiple professional software platforms for tedious manual interactive operation. This fragmented processing process not only is low in efficiency, but also is prone to result deviation due to human operation errors. Secondly, the physical properties of multi-source remote sensing data are significantly different. Landsat series satellites have evolved through multiple generations of sensors, with changing band settings and radiometric calibration parameters. MODIS and Sentinel-2 data sources also have essential differences in spatial resolution, revisit period and spectral response function. Traditional processing techniques cannot build a unified mapping mechanism to realize seamless splicing and batch standardization processing across sensors and platforms, which seriously restricts the comparability and analysis efficiency of long-term data.

[0004] In addition, in terms of time series trend mining, Sen slope estimation combined with Mann-Kendall significance test (Sen-MK) is considered as the authoritative method for evaluating the long-term trend of environmental variables due to its robustness to outliers and the assumption of data distribution. However, the algorithm implementation involves complex time series matrix operations and statistical inference logic. For researchers in the fields of ecology or geographic information science without deep programming background, implementing Sen-MK analysis at the full-pixel level has a high technical threshold, making it difficult to quickly apply this method to large-scale ecological monitoring business. SUMMARY

[0005] In view of the deficiencies of the prior art, the present application provides a remote sensing image batch index calculation and trend analysis system and method based on MATLAB GUI, which solves the problems of heterogeneous data sources, low calculation efficiency and lagging interactive experience in long-term remote sensing ecological monitoring.

[0006] To achieve the above object, the present application is realized by the following technical solutions: the present application provides a remote sensing image batch index calculation and trend analysis system based on MATLAB GUI, which comprises a data acquisition and metadata analysis module, a heterogeneous data adaptive mapping module, a multi-dimensional feature matrix construction module, a time series trend mining module and a visualization interaction module. The data acquisition and metadata analysis module adopts a combination of file name regular matching and binary file header reading, iterates all source data files under the target storage path, analyzes and acquires the timestamp and sensor type identifier of each file, and constructs a to-be-processed file index table. This mechanism can automatically identify and index remote sensing data from different sources without manual classification.

[0007] In view of the problem that the physical band arrangement order of multi-source remote sensing data is inconsistent, the heterogeneous data adaptive mapping module is internally provided with a band mapping lookup table mechanism. This module remaps source data with different physical storage structures into a standard logical spectrum vector of a unified format according to the sensor type identifier. The standard logical spectrum vector contains predefined logical channels corresponding to blue light, green light, red light, near-infrared, short-wave infrared and thermal infrared bands. For data with inconsistent spatial resolution, the system performs bicubic interpolation resampling during the mapping process to ensure that subsequent calculations are performed under a unified grid dimension.

[0008] In the ecological index calculation link, the multi-dimensional feature matrix construction module uses a full matrix vectorization operator to replace the traditional per-pixel loop calculation, greatly improving the calculation speed of greenness, wetness, aridity and heat components. This module particularly introduces a dynamic computation graph construction mechanism, which adaptively determines to use a full feature input mode containing a heat component or a default feature input mode excluding the heat component by pre-scanning the effectiveness of the thermal infrared channel in the standard logical spectrum vector. When performing principal component analysis transformation to synthesize the remote sensing ecological index matrix, the system automatically determines the first principal component direction by calculating the positivity of the load coefficient, and performs the necessary numerical inversion operation to ensure that the generated index value correctly represents the pros and cons of ecological quality.

[0009] The time series trend mining module introduces a three-dimensional spatio-temporal data cube structure to reorganize discrete multi-temporal remote sensing ecological index matrices into continuous time series data in memory. On this basis, the module performs vectorized non-parametric statistical operations. For Sen slope estimation, the system uses broadcast difference operation to calculate the slope matrix of all time pairs in parallel, and obtains the trend value through median reduction; for Mann-Kendall significance test, the system constructs a statistic matrix based on the slope sign difference. This processing method avoids loop iteration on the time series of each pixel, significantly reducing the time complexity of long-time series and large-area analysis.

[0010] The visual interaction module is built based on the MATLAB AppDesigner framework, and adopts an asynchronous interaction architecture with separation of main thread and background task thread. The complex calculation logic is encapsulated for execution by the background task thread, the main thread polls the background state through a timing listener, and the progress bar and log information are updated in real time, and the finally generated trend classification coding matrix is rendered as a raster image with geographic information.

[0011] To solve the problem of freezing of the graphical user interface (GUI) main thread caused by large-scale remote sensing matrix operations (especially the construction of spatiotemporal data cubes and per-pixel statistics), the present application designs an asynchronous listening mechanism based on parallel computing pool (Parallel Pool) and data queue (DataQueue). The specific technical implementation logic of the mechanism is as follows: Encapsulation and submission of background task object: in response to the user clicking the start processing instruction, the host program first encapsulates the data mapping, feature calculation and trend mining processes described above into an independent function handle using the parfeval function in the MATLAB parallel computing toolbox. Instead of directly calling the function, the system submits it to the background worker (Worker) and returns a future object (Future Object, such as parallel.FevalFuture) immediately. This object serves as the only handle for the background task and resides in the system memory, while the main UI thread is released to remain responsive.

[0012] Construction of common data bus: before submitting the task, the main thread initializes a parallel data queue (parallel.pool.DataQueue) as a common data bus connecting the main UI thread and the background Worker thread. The handle of this queue object is passed as an input parameter to the background calculation function. During the background operation, whenever the heterogeneous data adaptive mapping module completes processing for a year, or the temporal trend mining module completes a matrix operation, the background Worker calls the send instruction to push the structure data containing the current processing progress (percentage value) and running log (string information) into the data queue.

[0013] Asynchronous listening and callback update: the visual interaction module registers a listener callback function (Listener Callback) in the main thread, which is bound to the above-mentioned data queue through the afterEach instruction.

[0014] Listening mechanism: this mechanism is not a traditional dead-loop polling, but is based on event-driven (Event-Driven). When the background Worker sends data to the queue, the callback function in the main thread is automatically triggered.

[0015] Status Update: The callback function receives the progress data in the queue and directly manipulates the UI component handles to update the Value property of the Gauge and the Value property of the TextArea.

[0016] Task Completion Detection: Meanwhile, the system is configured with a Timer Object to detect the State property of the aforementioned Future Object at a low frequency (e.g., 1 Hz). When the state is detected to be finished, the timer triggers the "processing completed" logic, extracts the final trend classification coding matrix from the OutputArguments of the Future Object, and closes the parallel pool resource.

[0017] Through the cooperation of the above-mentioned parfeval (task distribution), DataQueue (cross-thread communication), and Listener (event response), the invention realizes the complete decoupling of computationally intensive tasks and interface rendering tasks.

[0018] The second aspect of the invention provides a remote sensing image batch index calculation and trend analysis method based on MATLAB GUI. This method is applied to the aforementioned system, realizing the full-process automatic processing from raw data input to spatio-temporal trend visualization.

[0019] The method first receives the target path specified by the user through the graphical user interface and parses the spatio-temporal metadata of the source file. Then, according to the sensor type and the preset lookup table, the heterogeneous physical waveband data is standardized into a logical spectrum vector. In the feature extraction stage, the method adaptively detects the waveband integrity of the input data, dynamically adjusts the input dimension of principal component analysis, and uses matrix operation to generate the remote sensing ecological index matrix of each time node in parallel.

[0020] Then, the method stacks and reorganizes the generated two-dimensional index matrix sequence into a three-dimensional spatio-temporal data cube in memory. Based on this data structure, the method performs Sen trend measurement and Mann-Kendall significance classification in parallel for the entire image, generating a trend classification coding matrix containing categories such as significant improvement, slight improvement, stability, slight degradation, and significant degradation. Finally, the matrix is mapped to a pseudo-color image with geographical coordinate projection, and the spatial distribution pattern of ecological change is displayed intuitively in the graphical interface.

[0021] The invention provides a remote sensing image batch index calculation and trend analysis system and method based on MATLAB GUI. It has the following beneficial effects: 1. The application can automatically reorganize the physical wave bands of different source data such as Landsat and Sentinel into a unified dimension standard logical spectral vector according to the sensor type identification through the heterogeneous data adaptive mapping module and the preset wave band mapping lookup table mechanism, which overcomes the limitation of manual classification and processing of different sensor data in traditional methods, realizes seamless batch processing of mixed source remote sensing images, and reduces the labor cost of data preprocessing.

[0022] 2. The application uses the dynamic calculation graph construction logic and feature polarity automatic correction mechanism in the multi-dimensional feature matrix construction module, so that the system can detect the effectiveness of the thermal infrared band in real time during operation, adaptively switch between full feature and default feature input modes, avoid calculation interruption caused by missing of part of the wave band, and automatically correct the index direction by analyzing the positive and negative nature of the principal component load coefficient, thereby solving the problem of ecological significance reversal caused by the randomness of the feature vector direction in principal component analysis, and ensuring the physical consistency of long-time series monitoring results.

[0023] 3. The application uses the three-dimensional spatio-temporal data cube structure and vectorization matrix operator in the time series trend mining module to optimize the traditional pixel-by-pixel loop iteration to full-frame parallel matrix operation, shortens the time consumption of Sen slope estimation and Mann-Kendall test, and separates the main thread and background task thread of the visualization interaction module, effectively avoiding the interface hang-up phenomenon in high-load calculation process, realizing real-time feedback of calculation progress and smooth visualization display of results. BRIEF DESCRIPTION OF DRAWINGS

[0024] Figure 1 is a system architecture diagram of the application; Figure 2 is a method flow diagram of the application; Figure 3 is a RSEI time series evolution trend diagram of a typical area in the Yangtze River Delta of the application embodiment.

[0025] Among them, 110, data acquisition and metadata analysis module; 120, heterogeneous data adaptive mapping module; 130, multi-dimensional feature matrix construction module; 140, time series trend mining module; 150, visualization interaction module. DETAILED DESCRIPTION

[0026] The technical solutions in the embodiments of the application will be described in detail below with reference to the drawings of the application specification. Obviously, the described embodiments are only part of the embodiments of the application, not all embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the application.

[0027] The MATLAB GUI-based remote sensing image batch index calculation and trend analysis method provided by the embodiment of the present application runs in a computer device. The computer device comprises a processor, a memory, a storage medium and an input / output interface connected through a system bus.

[0028] The processor is configured to execute matrix operations and logic control instructions, specifically including but not limited to a central processing unit (CPU) or a graphics processing unit (GPU), for executing parallel computing tasks for high-dimensional grid data. The memory is used to temporarily store data matrices and program instructions during the operation. The storage medium is used to store the heterogeneous remote sensing data source files to be processed and the generated trend analysis result files. The MATLAB runtime environment or the MATLAB complete development environment is deployed on the computer device to support the execution of the graphical user interface and the underlying algorithm.

[0029] Referring to the accompanying drawings Figure 1 The present application provides a MATLAB GUI-based remote sensing image batch index calculation and trend analysis system, which comprises a data acquisition and metadata analysis module 110, a heterogeneous data adaptive mapping module 120, a multi-dimensional feature matrix construction module 130, a time series trend mining module 140 and a visual interactive module 150.

[0030] The data acquisition and metadata analysis module 110 is configured to receive the target folder path selected by the user through the visual interactive module 150 and traverse all the files under the target folder. The data acquisition and metadata analysis module 110 identifies the format attribute of each source data file through a preset file name matching rule and a file header reading instruction. The format attribute includes timestamp information (Year) and sensor type identification.

[0031] The sensor type identification includes but is not limited to Landsat-5 TM, Landsat-7 ETM+, Landsat-8 OLI, Landsat-9 OLI-2 and Sentinel-2 MSL The data acquisition and metadata analysis module 110 indexes and associates the parsed timestamp information with the corresponding source data file path to generate a list of files to be processed.

[0032] The heterogeneous data adaptive mapping module 120 is configured to receive the list of files to be processed and, according to the sensor type identification, reorganize the spectral band data with inconsistent physical storage order in different source data files into a standard logical spectral vector of a unified format.

[0033] Let any one source data file be represented as an input matrix set , which contains physical band channels. The heterogeneous data adaptive mapping module 120 defines a standard logical spectral vector , the vector is defined as follows: ; wherein, represents blue band reflectance data, represents green band reflectance data, represents red band reflectance data, represents near-infrared band reflectance data, represents short-wave infrared first band reflectance data, represents short-wave infrared second band reflectance data, represents thermal infrared band radiance data.

[0034] The heterogeneous data adaptive mapping module 120 internally stores a band mapping lookup table. For each input source data file, the module executes a mapping function mapping the physical band index to the logical index of the standard logical spectral vector . If the source data file corresponds to a Landsat-8 sensor type, the mapping logic executes as assigning the physical 2nd band to ; If the source data file corresponds to a Landsat-5 sensor type, the mapping logic executes as assigning the physical 1st band to . Through this step, the heterogeneous data adaptive mapping module 120 outputs a logically structured data stream.

[0035] The multi-dimensional feature matrix construction module 130, connected with the heterogeneous data adaptive mapping module 120, is configured to calculate ecological index features based on the standard logical spectral vector . The multi-dimensional feature matrix construction module 130 contains a completeness detection unit for scanning whether there are non-empty values in each component in before calculation.

[0036] When component has valid data, the multi-dimensional feature matrix construction module 130 calculates the remote sensing ecological index according to the full-component mode; when component data is missing, the multi-dimensional feature matrix construction module 130 automatically switches to the default mode, and only uses to components to calculate RSEl without interrupting the batch processing process. The multi-dimensional feature matrix construction module 130 outputs a feature result image matrix corresponding to each time node.

[0037] The time series trend mining module 140 is configured to receive the multi-year feature result image matrix output by the multi-dimensional feature matrix construction module 130, and construct a three-dimensional space-time data cube along the time axis direction. The time series trend mining module 140 performs a non-parametric statistical operation on each pixel position in the three-dimensional space-time data cube, calculates a Sen slope value and a Mann-Kendall test statistic, and generates a trend classification coding matrix accordingly.

[0038] The visualization interaction module 150 is configured to provide a human-computer interaction interface. The interface includes an input path selection control, a year interval setting control, a parameter configuration control, and a processing progress display bar. The visualization interaction module 150 converts the trend classification coding matrix into a raster image file (GeoTIFF) with geographic coordinate information, and renders and displays it in the interface preview area.

[0039] Referring to the drawings Figure 2 The present application also provides a MATLAB GUI-based remote sensing image batch index calculation and trend analysis method, comprising the following steps: The specific operation logic of the data acquisition and metadata analysis module 110 is composed of the following technical steps: The data acquisition and metadata analysis module 110 first initializes a file walker. In response to the received root directory path parameter, the file walker performs a recursive scanning operation to retrieve all file objects under the path and its sub-paths. During the scanning process, the file walker applies a preset extension filter to retain only raster image files with a suffix of.tif or.tiff, and filters out auxiliary metadata files (such as.xml,.txt) and system temporary files to generate a candidate file set.

[0040] For each file object in the candidate file set, the data acquisition and metadata analysis module 110 calls the string regular matching unit to perform naming standard analysis. The string regular matching unit has a regular expression library for standard remote sensing data naming rules loaded in memory.

[0041] In one specific embodiment, for Landsat satellite data, the regular expression is configured to extract the character sequence pattern starting with L. The string regular matching unit performs the following analysis logic: read the leading character segment of the file name, extract the sensor identifier field. This field usually occupies the first 4 characters of the file name. For example, for the file name LC08_L1TP_123032_20200501_... Extract the first 1-4 characters LC08 and identify it as the unique feature code of the Landsat-8 OLI sensor data; Extract the 18th-21st characters 2020 and identify them as the time metadata of the acquisition year.

[0042] If the file name does not conform to the preset regular expression rule, the data acquisition and metadata parsing module 110 triggers a file header reading instruction. By calling the LibTIFF library or a similar underlying binary reading interface, the binary header information of the GeoTIFF file is directly accessed. The system reads the tag data field of TIFFTagID 270 (ImageDescription) or 306 (DateTime) to decode the sensor metadata information and imaging time information from the binary stream. This dual parsing mechanism ensures that the system can accurately obtain metadata even if the file name is changed or not properly named.

[0043] The data acquisition and metadata parsing module 110 internally maintains a dynamic index table structure. After completing the parsing of each file, the system creates a record in the index table. The record contains three key fields: Physical storage path field: stores the absolute path address of the image on the disk; Timestamp index field: stores the parsed year value as the primary key for subsequent time sequence sorting; Source type enumeration field: stores the standardized sensor type code (such as mapping LC08, L8, Landsat8 to the internal enumeration value TYPE_L8).

[0044] After completing the traversal, the data acquisition and metadata parsing module 110 sorts all records in the index table in ascending order according to the timestamp index field. If multiple images are detected under the same timestamp index (i.e., there are multiple scene data in the same year), the system records their corresponding spatial row and column numbers and marks them as spatial stitching to-be-processed items. The final ordered index table is passed to the heterogeneous data adaptive mapping module 120 as a structured metadata stream.

[0045] The core function of the heterogeneous data adaptive mapping module 120 is to build a middleware layer independent of the physical sensor hardware characteristics. The module first initializes a standard logical spectral space in the system memory, which is defined as a three-dimensional matrix container with a predetermined channel depth. The channel depth is fixed at 7 layers, corresponding to logical channel indexes 1 to 7, which are logically mapped to: blue channel (Logic_Blue), green channel (Logic_Green), red channel (Logic_Red), near-infrared channel (Logic_NIR), short-wave infrared first channel (Logic_SWIR1), short-wave infrared second channel (Logic_SWIR2), and thermal infrared channel (Logic_Thermal).

[0046] In order to realize the accurate projection from the source data physical address to the above-mentioned standard logical spectral space, the heterogeneous data adaptive mapping module 120 is internally built with an extensible band mapping lookup table. The lookup table is stored in a key-value pair data structure. The key is the sensor type enumeration value (such as TYPE_L8, TYPE_L5, TYPE_S2) passed by the upper module, and the value is a one-dimensional integer array, and the element value in the array represents the physical storage level index of the corresponding band in the original data file of the sensor.

[0047] In a specific operation, when the system receives a data object marked as Landsat-8 OLI / TIRS, the module queries the lookup table and extracts the corresponding index mapping vector . The system performs a matrix slicing operation, directly reads the 2nd layer data in the source file and writes it to the 1st layer (Logic_Blue) of the standard logical space, reads the 3rd layer data and writes it to the 2nd layer (Logic_Green) of the standard logical space, and so on, and writes the 10th layer data to the 7th layer (Logic_Thermal) of the standard logical space. In this process, the aerosol band (1st layer) specific to Landsat-8 is automatically ignored because it is not defined in the mapping vector, thereby realizing the data cleaning function.

[0048] When the system receives a data object marked as Landsat-5 TM, the module queries the lookup table and extracts the index mapping vector . It is worth noting that in the physical storage structure of Landsat-5, the thermal infrared data is usually located in the 6th layer, and the second short-wave infrared band is located in the 7th layer. The mapping logic of the present application explicitly redirects the 7th layer of the source data to the 6th layer (Logic_SWIR2) of the standard logical space, and redirects the 6th layer of the source data to the 7th layer (Logic_Thermal) of the standard logical space. Through this non-linear address remapping mechanism, the system eliminates the band sequence differences between different generations of sensors.

[0049] In addition, the heterogeneous data adaptive mapping module 120 also integrates a spatial resolution normalization unit. During the data extraction process of band mapping, if it is detected that the spatial resolution of the source data is inconsistent with the preset system reference resolution (e.g. 30 meters), the unit will automatically trigger the Resampling subprogram.

[0050] For example, Sentinel-2 data has 10-meter resolution in the red band and 20-meter resolution in the short-wave infrared band. When the mapping logic fills these bands into the standard logical spectral space, the spatial resolution normalization unit adopts bicubic interpolation algorithm to resample all non-reference resolution band matrices to the reference grid size. This step ensures that the final generated standard logical spectral vector is strictly aligned in the spatial dimension (number of rows and number of columns ), so that subsequent matrix algebra operations can be performed legally, avoiding calculation exceptions caused by mismatched matrix dimensions.

[0051] For some special data sources (such as optical satellites without carrying thermal infrared sensors), the corresponding mapping bit in the lookup table is set to a special identifier (such as -1 or NULL). When the mapping module reads this identifier, it will fill a full zero matrix or NaN mask matrix in the corresponding channel layer of the standard logical spectral space (such as Logic_Thermal). This operation provides a clear data-level judgment basis for the integrity detection in the subsequent feature calculation module.

[0052] The multi-dimensional feature matrix construction module 130 internally integrates a dynamic computation graph construction unit, which uses a directed acyclic graph data structure to abstract and manage the feature calculation process. In the directed acyclic graph, the leaf nodes represent the individual component channels in the input standard logical spectral vector, the intermediate nodes represent the calculation operators of single feature indicators (such as NDVI, Wet, NDBSI, LST), and the root node represents the final remote sensing ecological index (RSEI) synthesis operator. The directed edges between nodes represent data dependency relationships.

[0053] Before performing the calculation on each time node's standard logical spectral vector , the dynamic computation graph construction unit first starts the channel integrity scanning program. This program reads the seventh logical channel (i.e. the thermal infrared channel Logic_Thermal) of the , verifies the validity of its data content. If the data in this channel is identified as the invalid mask (such as NaN or full zero matrix) filled by the aforementioned heterogeneous data adaptive mapping module 120, the system generates a missing state signal of Boolean type.

[0054] In response to the missing state signal, the dynamic computation graph construction unit triggers the graph topology reconstruction logic and performs dynamic pruning operations on the dependent path. Specifically, the system first locates to the land surface temperature (LST) calculation node that depends on the thermal infrared channel, and marks the state of this node as inactive (Inactive). Subsequently, the system traverses downstream along the directed edge, blocking the data transmission path from the LST node to the root node (RSEI synthesis node).

[0055] At the same time, the dynamic computation graph construction unit sends a reconstruction instruction to the configuration interface of the root node. In the default full data mode, the root node is configured to receive four input matrices and construct a 4x4 covariance matrix to perform principal component analysis. When receiving the reconstruction instruction, the root node automatically switches its internal matrix constructor, reduces the number of input ports from 4 to 3 (i.e., only receives the feature matrices corresponding to greenness, humidity, and dryness), and redefines the core operation logic as constructing a 3x3 covariance matrix.

[0056] After completing the adaptive adjustment of the above topology and node parameters, the system generates an optimized execution queue. The computing engine loads the memory data in accordance with the queue order and executes the matrix operators corresponding to each activated node in turn. This dynamic scheduling mechanism based on graph theory enables the system to mix processing of data sources containing thermal infrared bands (such as Landsat-8) and data sources lacking thermal infrared bands (such as Sentinel-2 with only optical bands or damaged data) in the same batch processing task without manual intervention or batch operation. By shielding invalid computation branches, the system not only prevents runtime errors caused by missing input dimensions, but also avoids the computational overhead of invalid data, ensuring the continuity and stability of massive heterogeneous data batch processing jobs.

[0057] The adaptive feature calculation unit discards the traditional pixel-by-pixel traversal processing mode and instead adopts a full matrix vectorization operation strategy. This module treats each channel in the standard logical spectral vector as an independent two-dimensional floating-point matrix . During the calculation process, the system calls the underlying linear algebra library (such as BLAS or LAPACK) to perform single instruction multiple data (SIMD) type parallel operations on these two-dimensional matrices.

[0058] For the greenness component, the computing engine extracts the 4th layer (near-infrared matrix ) and the 3rd layer (red light matrix ) from the logical data stack. The system first allocates two temporary buffers in memory to store the difference matrix and the sum matrix of and . Subsequently, the point division operation of the matrix is performed. To prevent division by zero exceptions, the system performs zero value mask filtering on the sum matrix before performing the division, marking the corresponding positions with a non-numeric (NaN) value when the denominator is zero, thereby obtaining the greenness feature matrix For the wetness component, the adaptive feature computation unit performs a linear projection operation based on a coefficient vector. The system first retrieves the corresponding tasseled cap transformation coefficient vector from the pre-configured physical parameter database according to the sensor type identification delivered by the aforementioned metadata parsing module . For example, for Landsat-8 data, the system loads a specific set of empirical coefficients; for Sentinel-2 data, another set of adapted coefficients is loaded. The computation engine performs a weighted superposition operation between the first 6 optical channel matrix of the logical spectral space and the coefficient vector , i.e. , to generate the wetness feature matrix of the whole scene at one time.

[0059] For the dryness component, the computation involves the synthesis of the bare soil index (SI) matrix and the built-up index (IBI) matrix. The system constructs a multi-stage pipeline: first, it computes two intermediate result matrices in parallel and . When computing , a nonlinear combination of three band ratio matrices is involved, which is accomplished by the system through a chain of matrix operations. Finally, the system performs matrix addition and scalar division on the heart rate interval data and the signal strength data (matrix) , and takes their arithmetic mean as the dryness feature matrix .

[0060] For the heat component, it is activated only when the system is in the full-band computation mode. The computation engine reads the 7th logical channel, the thermal infrared matrix . The system first reads the radiation bias and gain parameters and the thermal conversion constant from the metadata, converts the original digital quantization value (DN) matrix into a radiance matrix, and then solves the brightness temperature matrix through the inverse function matrix transformation of the Planck function. Subsequently, the system introduces the specific radiance correction matrix to perform element-wise physical correction on the brightness temperature matrix to generate the final land surface temperature feature matrix .

[0061] After the computation of the above four (or three, depending on the mode) basic feature matrices is completed, the adaptive feature computation unit immediately starts the interquartile normalization subroutine. This subroutine iterates through each feature matrix, identifies its global maximum scalar and minimum scalar , and performs the following normalization mapping: ; wherein, the normalized new value, The original value you want to process.

[0062] This step maps all feature data with different physical meanings to the dimensionless numerical range [0,1], eliminating the interference of dimensional differences on the weights of subsequent principal component analysis and outputting a standardized multidimensional feature dataset.

[0063] After completing the standardization of the basic feature components, the multidimensional feature matrix construction module 130 triggers the principal component analysis subroutine. The core logic of this subroutine is to compress the multidimensional feature space into a one-dimensional comprehensive index space while retaining the maximum information load.

[0064] First, the system executes the adaptive data stack construction step. The computing engine detects the current computing mode status. If the system is in full-band mode, the computing engine will use the standardized greenness matrix ( ), humidity matrix ( ), dryness matrix ( ) and heat matrix ( Stack along the depth direction to construct a dimension of The four-dimensional feature tensor.

[0065] If the system is in default mode, the calculation engine only extracts the greenness, humidity, and dryness matrices, constructing a dimensionless matrix. The three-dimensional feature tensor.

[0066] To ensure the stability of covariance calculation, the system performs validity masking filtering on the feature tensor, removing pixel vectors containing NaN (non-numerical) values, and reshaping the remaining valid observations. Two-dimensional observation matrix (where The total number of valid pixels, (The feature dimension is either 4 or 3).

[0067] Subsequently, the system performs eigenvector solving and projection steps. The computation engine calculates the covariance matrix of the aforementioned two-dimensional observation matrix. The system calls the eigenvalue decomposition algorithm to solve the problem. eigenvalues and its corresponding eigenvectors The system selects the largest eigenvalue. corresponding feature vector As the transformation coefficients of the first principal component.

[0068] Based on this coefficient, the system constructs the initial RSEI matrix. Its computational logic is represented by a linear weighted combination of each feature component: ; in, For feature vectors The weight value of the first component, The normalized feature matrix of the first component.

[0069] Next, the system performs a feature polarity automatic correction step. Since the eigenvector direction in principal component analysis has mathematical uncertainty (may point to the positive or negative direction of the data distribution), the directly calculated may be positively or negatively correlated with ecological quality. In order to ensure the physical meaning consistency of the output results (i.e. the higher the value, the better the ecology), the system has a logical judgment unit built in.

[0070] The logical judgment unit reads the load coefficient of the greenness component (NDVI) in the first principal component . According to ecological principles, vegetation coverage (NDVI) usually has a positive contribution to ecological quality. Therefore, the logical judgment unit performs the following conditional branching: If , it is determined that the current PC1 direction is correct, and is kept unchanged; If , it is determined that the current PC1 direction is reversed, and the system performs a matrix inversion operation: ; This step forces the semantic direction of high-dimensional feature reduction by algorithm, ensuring the logical consistency of the automatic batch processing results. Finally, the system performs the result normalization and output step. In order to facilitate the horizontal comparison between different years, the system again calls the range standardization operator to map the direction-corrected (or ) to the [0, 1] interval, generating the final remote sensing ecological index matrix . The system encapsulates this matrix as a GeoTIFF file format with geographic reference information, and passes it to the next level of time trend mining module 140.

[0071] The time trend mining module 140 performs a data reorganization mechanism based on memory mapping or direct memory access, aiming to convert discrete distributed two-dimensional feature image files into three-dimensional data entities suitable for performing vectorized statistical operations.

[0072] First, the time trend mining module 140 starts the sequence index construction program. This program scans the feature result storage directory output by the previous level module, extracts the file handle of all feature raster files (i.e. images). The system parses the time attribute (year ) associated with each file, and sorts them in chronological order ) The file handle list is sorted in ascending order. In this process, the system builds a one-dimensional time vector , which records the physical time scale of the data cube along the third dimension (z-axis).

[0073] Subsequently, the system performs a three-dimensional space pre-allocation operation. The system reads the header file information of the first valid file in the sorted list to obtain the number of rows (Height), the number of columns (Width), and the data type (such as 32-bit floating point) of the grid matrix. Based on the obtained dimension parameters and the total length of the file sequence , the system allocates a continuous storage space in memory and initializes a three-dimensional zero matrix or a full NaN (Not-a-Number) matrix with dimensions , denoted as data cube . The pre-allocation strategy is used instead of the dynamic expansion strategy to avoid frequent triggering of memory reallocation and data movement during large-scale data writing, thereby ensuring storage efficiency.

[0074] Next, the system enters the layer-by-layer data injection phase. The time trend mining module 140 traverses the sorted file handle list. For the file at the th time point, the system loads it as a two-dimensional matrix . Before injection, the system performs spatial consistency verification to compare whether the row and column dimensions of are strictly consistent with the bottom surface dimensions of the data cube .

[0075] After verification, the system uses matrix slicing assignment operations to map the two-dimensional matrix to the slice position corresponding to the th time point of the three-dimensional data cube: ; Through this loop operation, the originally discrete spatial data stored on the disk is reorganized into a voxel structure in memory, where each voxel represents the feature value of the spatial position at the th row and the th column at the th time point.

[0076] In addition, the time series trend mining module 140 also extracts and binds the geographic reference information (including the projection coordinate system WKT string, the affine transformation six-parameter matrix) of the first file into the metadata field of the data cube object. This ensures that the dimensionless three-dimensional numerical matrix can be reversely restored to the raster map data with geographic spatial properties after subsequent processing. For cases where the entire data cube cannot be accommodated in the memory at one time (e.g., a super large region or a super long time series), the system enables a block mapping logic to divide the entire map into a plurality of sub-cube blocks, which are sequentially loaded into the memory for subsequent trend analysis operations.

[0077] The time series trend mining module 140 adopts a full matrix vectorization calculation strategy to replace the double nested loop structure for a single pixel in the traditional time series analysis, thereby significantly improving the operation efficiency for massive raster data.

[0078] First, the time series trend mining module 140 performs a dimension reshaping operation on the input three-dimensional spatio-temporal data cube . The system flattens the spatial dimensions (rows and columns ) and converts the three-dimensional cube into a two-dimensional pixel time matrix with dimensions , where represents the total number of pixels, and represents the total length of the dimension parameter and the file sequence. In this matrix, each row represents a time series vector of a specific geographic location.

[0079] Next, the system performs a time series pairing index construction step. In order to calculate the Theil-Sen slope, all possible point pair combinations in the time series need to be obtained, which satisfy . The system generates two index vectors in the memory: a source index vector and a target index vector , both with a length of the combination number . These two vectors define the positions of all time point pairs that need to be calculated for the slope.

[0080] Subsequently, the system performs a broadcast difference operation. Using the above index vectors, the system extracts two extended matrices from the matrix : a post-phase matrix with dimensions , whose th column corresponds to the data column of the later observation phase in the th combination; a pre-phase matrix with dimensions​ , the first column corresponds to the first observation in the first combination.

[0081] Meanwhile, the system computes a time interval constant vector , whose elements are the difference values corresponding to the paired years . The system invokes a parallel floating-point operation instruction to perform the following matrix arithmetic operation: ; where is the result of broadcasting the time interval vector to the same dimension as . The output of the operation is a full slope matrix, which stores all possible paired rate of change at each pixel location.

[0082] Finally, the system performs a robust statistical reduction operation. For the matrix, the system performs a row-level median operation along the combination number dimension (i.e., the second dimension). This step automatically filters out the interference of outliers and extracts the Sen slope matrix at each pixel location. The system inversely reshapes the computed trend value vector (length ) into a two-dimensional spatial distribution matrix, i.e., the Sen slope image.

[0083] During the computation, if there are NaN (Not a Number) values in , the system enables the NaN-ignoring flag in the difference operation and median reduction steps, ensuring that the trend is calculated based only on valid observation points, and ensuring the robustness of the algorithm to partial-year data loss.

[0084] The statistical inference unit in the time series trend mining module 140, after completing the Sen slope estimation, immediately performs a Mann-Kendall (MK) test based on the same set of time series data structures. This test aims to confirm whether the trend calculated above is statistically significant, i.e., to exclude the possibility of random fluctuations leading to false trends.

[0085] First, the system performs an s-statistic matrix construction step. The system reuses the difference matrix generated in the Sen slope calculation phase (i.e., the difference between the later phase matrix and the earlier phase matrix). The system applies a sign function to the difference matrix, mapping all numerical values in the matrix to a ternary set . Specifically, positive difference values are mapped to 1, negative difference values are mapped to -1, and zero difference values are mapped to 0. Subsequently, the system performs a summation operation on the sign matrix along the combination number dimension, generating an s-statistic matrix in two-dimensional space​​ Each element in this matrix This represents the difference in the number of positive and negative change pairs in the time series corresponding to the geographical location.

[0086] Next, the system performs variance correction and standardization. To improve the rigor of the test, the algorithm incorporates correction logic for tie data. A tie refers to a set of repeating data points in the time series where the values ​​are exactly the same. The system scans the time series vector of each pixel and counts the number of sets with repeating values. and the breadth of equal value groups .

[0087] The system calculates the variance according to the following formula. variance matrix of statistics ; in, This represents the effective time series length. For cases with no duplicate values, the summation term... It is 0.

[0088] based on Matrix and Matrix, system constructs standardized test statistic matrix To avoid the reduced computational efficiency caused by pixel-by-pixel conditional branches, the system uses mask matrix operations to implement piecewise function logic: Construct a positive value mask Calculate the corresponding position ; Construct a negative mask Calculate the corresponding position ; Constructing a zero-value mask The corresponding position is directly assigned the value 0. The three temporary matrices are superimposed to form the final matrix. matrix.

[0089] Finally, the system performs a multidimensional feature coupling classification step. The system introduces a saliency level threshold (typically set to...). Corresponding critical value The system initializes a classification result matrix with the same dimensions as the input image. The default value is set to 3 (representing stability).

[0090] The significance level threshold is the threshold set in statistical hypothesis testing. = 0.05 (i.e. 5% significance level). According to the standard normal distribution, the critical value (Z-score) corresponding to this level is 1.96. In the classification step of the system, this threshold is used to distinguish between significant features and random noise: features whose absolute value of the statistic exceeds 1.96 are considered significant and thus retained for further processing; otherwise, they are filtered out as noise.

[0091] The system utilizes the Sen slope matrix and the test statistic matrix to generate five mutually exclusive spatial masks through Boolean logical operations, and performs hierarchical assignment updates on Significant degradation mask: positions satisfying the condition are assigned a value of 5; Minor degradation mask: positions satisfying the condition are assigned a value of 4; Minor improvement mask: positions satisfying the condition are assigned a value of 2; Significant improvement mask: positions satisfying the condition are assigned a value of 1. Regions not covered by the above masks remain the default value of 3 (stable).

[0092] Finally, the statistical inference unit outputs the numerical classification result matrix To enhance the readability of human-computer interaction, the system binds a pre-set pseudo-color lookup table for this matrix, for example, rendering level 1 as dark green and level 5 as dark red, thus generating a spatial distribution map that intuitively reflects the regional ecological change pattern.

[0093] The visualization interaction module 150 is developed based on the MATLAB AppDesigner framework, and its front-end interface is built as a UI graphical container containing multiple functional partitions. The UI graphical container is mainly divided into: data input and configuration panel, execution control and state monitoring panel, and geographic spatial result preview panel.

[0094] The data input and configuration panel is deployed with a path selection control (Button) and a text display field. When the user triggers the path selection control, the system calls the file directory selection dialog box at the operating system bottom layer. Once the user confirms the path selection, the system obtains the path string through the callback function and immediately triggers the background pre-scanning thread. This pre-scanning thread quickly traverses the target folder, counts the number of raster files that meet the format requirements, and feeds back the statistical results to the status label on the panel in real time, so that the user can confirm the effectiveness of the data source. In addition, this panel also contains parameter setting components that allow users to input numerical limits for the start year and end year, as well as custom output result file storage paths. ​

[0095] The execution control and status monitoring panel contains a start processing instruction button and a progress visualization component. In order to solve the problem of interface freezing caused by large-scale matrix operation, the embodiment adopts an asynchronous interaction mode of separation of main thread and worker thread. When the user clicks the start processing button, the system does not directly perform calculation in the current UI drawing thread, but encapsulates the aforementioned heterogeneous data adaptive mapping module 120, multi-dimensional feature matrix construction module 130 and time series trend mining module 140 into an independent background task object through the parallel computing toolbox.

[0096] At the same time, the system starts a timing listener in the UI thread. The listener polls the execution state attribute of the background task object at a preset frequency (for example, every 0.5 seconds). When the background task completes the data processing of each year or completes a step (such as RSEI calculation completion, Sen slope calculation completion), it writes the current progress percentage value and log text information to the public data bus. The timing listener captures these information and drives the progress bar component (Gauge) on the interface to be dynamically filled, and at the same time, the log text is appended to the scroll log text box, realizing the visualization monitoring of the black box calculation process.

[0097] The core component of the geospatial result preview panel is a drawing coordinate axis with geospatial coordinate system support. When the background task completes the trend analysis and generates the classification result matrix After that, the system automatically triggers the visualization rendering subroutine. The subroutine first reads the matrix and loads the preset discretization color mapping table. The mapping table defines the correspondence between the values 1 to 5 and specific RGB color values (for example: value 1 corresponds to dark green, value 5 corresponds to dark red, and value 3 corresponds to gray).

[0098] The system uses imagesc or similar drawing instructions to map the classification matrix to a pseudo-color image and draw it on the coordinate axis. In order to preserve the geospatial context, the system reads the spatial reference object of the first input file, converts the pixel coordinates of the image to projected coordinates, and sets the X-axis and Y-axis scale labels of the coordinate axis to display the true longitude and latitude or meter coordinates. In addition, the panel also integrates an interactive toolbar, including zoom (Zoom), pan (Pan) and data probe tools. When the user activates the data probe and clicks on any point on the image, the system captures the row and column numbers of the point, indexes the original data cube, extracts the long-time series RSEl value change curve corresponding to the point, and pops up a floating window to draw the curve, thereby realizing interactive exploration from macro spatial distribution to micro time series change.

[0099] Refer to the attached Figure 3Specific application example: Long-Triangle Delta region long-term ecological quality trend monitoring (2015-2024).

[0100] This example aims to use the system provided by the present application to conduct a 10-year continuous monitoring of the ecological quality of a certain core area (covering an area of about 100 km x 100 km) in the Yangtze River Delta. The data sources mix Landsat-8 OLI / TIRS (2015-2020) and Sentinel-2 MSL (2021-2024) images.

[0101] User operation: The user selects the root directory D:_Delta in the system interface, which stores 50 original image files.

[0102] File scanning: The system traverses the directory and identifies 30.tif files.

[0103] Metadata parsing: The file name LC08_L1TP_118038_20150812_… is parsed by regular expression as: sensor = Landsat-8, year = 2015.

[0104] The file name S2A_MSIL2A_20220720_… is parsed by regular expression as: sensor = Sentinel-2, year = 2022.

[0105] Index construction: The system generates an index table, marking the data from 2015 to 2020 as TYPE_L8 and the data from 2021 to 2024 as TYPE_S2.

[0106] Processing Landsat-8 data (2015 scene): The system queries the mapping table vL8 and extracts the 2nd-7th and 10th bands. Action: Fill the physical 2nd band (blue light) into logical channel 1 and the physical 10th band (thermal infrared, 100m resolution) into logical channel 7. The system detects the resolution difference and automatically performs bicubic interpolation on the thermal infrared band, resampling it to 30m. Processing Sentinel-2 data (2022 scene): The system queries the mapping table and extracts the corresponding blue, green, red, near-infrared and short-wave infrared bands. Action: Since Sentinel-2 has no thermal infrared sensor, logical channel 7 is filled with a NaN mask. The red band (10m) and the short-wave infrared band (20m) are both resampled to the reference resolution of 30m to match the historical Landsat data.

[0107] The system performs batch processing calculations on the normalized logical vectors: Dynamic graph construction: 2015 data (including thermal infrared): Detect that logical channel 7 is valid, activate full feature mode.

[0108] Greenness: ; LST using the inverse Planck function.

[0109] Principal Component Analysis: Input Construct 4x4 Covariance Matrix.

[0110] 2022 data (missing thermal infrared): Detected logical channel 7 as , switched to default mode.

[0111] Feature Calculation: Skip LST calculation node.

[0112] Principal Component Analysis: Input only Construct 3x3 Covariance Matrix.

[0113] Feature Polar Correction: When calculating 2015 RSEI, the system solved for the first principal component eigenvector .

[0114] Judgment: Greenness Load , determined direction reversed.

[0115] Execution: , ensured high values represent good ecology, where represents the final RSEI value, represents the original RSEI value.

[0116] The system recombines the 10 RSEI images from 2015-2024 into a data cube of .

[0117] Vectorized Sen Slope Estimation: For each pixel , extract the time series .

[0118] Construct all pairs of time series .

[0119] Broadcast Difference: Calculate the slope for all pairs.

[0120] Result: Taking the median gives (0.023 annual improvement) with Figure 3 a black dashed line showing the fitted line based on this trend.

[0121] Mann-Kendall Significance Test: Calculate statistic and variance .

[0122] Normalized Z-Value: ; Assume the calculation of .

[0123] Decision logic: and (confidence 95%). Conclusion: the pixel is marked as Level 1 (significant improvement), as Figure 3 indicated by the top-left textbox, the system outputs the statistical conclusion.

[0124] Result rendering: refer to the attached Figure 3 , the urban center area (Shanghai urban area) appears dark gray or black (corresponding to significant degradation or stable low value); the surrounding wetland park and newly afforested area appear light gray or white (corresponding to significant improvement).

[0125] Interactive query: the user clicks on the Chongming East Beach area (light gray pixel) on the map, and the system pop-up window displays the "RSEI time series evolution trend chart" of the area (as Figure 3 shown). The chart uses a composite visualization form: The light gray column chart directly displays the RSEI annual observation values from 2015 to 2024, showing a fluctuating upward trend of the values over the years; The black dotted line identifies the Sen trend fitting line, clearly reflecting the overall improvement direction; The textbox in the upper left corner of the chart directly outputs the statistical results of the Mann-Kendall significance test (Z statistic is 2.15, greater than the critical value 1.96), and clearly determines the trend of the area as significant improvement (Level 1), thereby quantitatively verifying the effectiveness of the ecological protection project.

Claims

1. A system for batch index calculation and trend analysis of remote sensing images based on MATLAB GUI, characterized in that, include: The data acquisition and metadata parsing module (110) is used to traverse the target folder, parse the timestamp and sensor type identifier of the source data file, and generate a file index table containing the physical storage path; The heterogeneous data adaptive mapping module (120) is used to query the preset band mapping lookup table based on the sensor type identifier in the file index table, and remap source data files with different physical structures into standard logical spectral vectors of a unified format. A multidimensional feature matrix construction module (130) is used to calculate ecological feature components in parallel using matrix operators based on the standard logical spectral vector, and to perform principal component analysis transformation to synthesize a remote sensing ecological index matrix. The time-series trend mining module (140) is used to reorganize the remote sensing ecological index matrix of multiple time phases into a three-dimensional spatiotemporal data cube, perform vectorized Sen slope estimation and Mann-Kendall significance test along the time dimension, and generate a trend classification coding matrix. The visualization interaction module (150) is used to receive the user's start command, obtain the trend classification coding matrix through an asynchronous listening mechanism, and render it as a raster image with geographic coordinate information for display.

2. The remote sensing image batch index calculation and trend analysis system based on MATLABGUI according to claim 1, characterized in that, The heterogeneous data adaptive mapping module (120) generates a standard logical spectral vector that defines a three-dimensional matrix container with a predetermined depth. The logical channels corresponding to the predetermined depth are, in order, blue light channel, green light channel, red light channel, near-infrared channel, short-wave infrared first channel, short-wave infrared second channel and thermal infrared channel. The heterogeneous data adaptive mapping module (120) is used to extract the index mapping vector according to the sensor type identifier and use matrix slicing operation to fill the data in the source data file into the corresponding logical channel of the standard logical spectral vector; If the spatial resolution of the source data file is inconsistent with the system reference resolution, a bicubic interpolation resampling operation is performed before padding.

3. The remote sensing image batch index calculation and trend analysis system based on MATLABGUI according to claim 2, characterized in that, The multidimensional feature matrix construction module (130) includes a dynamic computation graph construction unit, which is used to scan the thermal infrared channels in the standard logical spectral vector before computation: If the thermal infrared channel data is valid, construct a full-band calculation map, activate the surface temperature calculation node, and set the principal component analysis transformation to receive four feature matrices of greenness, humidity, dryness, and heat as inputs; If the thermal infrared channel data is an invalid mask, a dependency path pruning operation is triggered, the surface temperature calculation node is marked as inactive, and the parameters of the principal component analysis transformation are reconstructed to only accept three feature matrices: greenness, humidity, and dryness.

4. The remote sensing image batch index calculation and trend analysis system based on MATLAB GUI according to claim 1, characterized in that, The multidimensional feature matrix construction module (130) is used to execute the automatic feature polarity correction logic when performing principal component analysis transformation: Calculate the covariance matrix of the input feature matrix set and solve for the eigenvectors of the first principal component; Extract the loading coefficients corresponding to the greenness component from the feature vector; The sign of the loading coefficient is determined. If the loading coefficient is less than zero, the direction of the first principal component is determined to be negatively correlated with the ecological quality. A numerical inversion operation is then performed on the initially calculated index matrix to generate the final remote sensing ecological index matrix.

5. The remote sensing image batch index calculation and trend analysis system based on MATLAB GUI according to claim 2, characterized in that, The multidimensional feature matrix construction module (130) performs the following full matrix vectorization operation when calculating ecological feature components: Extract the near-infrared channel matrix and red light channel matrix of the standard logical spectral vector, and perform matrix division to generate the greenness component; Based on the sensor type identifier, load the corresponding physical coefficient vector, and perform a linear weighted superposition operation on the relevant logical channel matrix to generate humidity components; Parallel calculation of the bare soil index matrix and the building index matrix, and taking the arithmetic mean to generate the dryness component; The digital quantization matrix of the thermal infrared channel is converted into a radiance matrix. The radiance temperature matrix is ​​calculated using the inverse Planck function transformation, and the thermal component is generated by correction using the emissivity matrix.

6. The remote sensing image batch index calculation and trend analysis system based on MATLAB GUI according to claim 1, characterized in that, The time-series trend mining module (140) is used to construct the three-dimensional spatiotemporal data cube through the following steps: Analyze the year information corresponding to all input remote sensing ecological index matrices and establish an ascending time index vector; Read the row and column dimensions of the first remote sensing ecological index matrix and the length of the time index vector, and pre-allocate three-dimensional continuous storage space in memory; Traverse the input remote sensing ecological index matrix sequence, map each two-dimensional matrix to the slice level corresponding to the three-dimensional continuous storage space, and bind the geographic reference information of the first matrix to the metadata field of the three-dimensional spatiotemporal data cube.

7. The remote sensing image batch index calculation and trend analysis system based on MATLABGUI according to claim 6, characterized in that, The time-series trend mining module (140) is used to perform vectorized Sen slope estimation on the three-dimensional spatiotemporal data cube: The three-dimensional spatiotemporal data cube is flattened in the spatial dimension and converted into a pixel time two-dimensional matrix; a full permutation time pair index is constructed, and a subsequent phase matrix and a previous phase matrix are generated based on the pixel time two-dimensional matrix; The slope matrix of all time pairs is calculated using broadcast difference operations, and row-level median reduction is performed along the combination dimension to obtain the Sen slope image.

8. The remote sensing image batch index calculation and trend analysis system based on MATLAB GUI according to claim 7, characterized in that, The time-series trend mining module (140) is also used to perform the Mann-Kendall significance test based on the Sen slope image: Based on the sign of the difference between the later phase matrix and the earlier phase matrix, calculate the S-statistic matrix and the variance matrix; A standardized test statistic matrix is ​​constructed, and combined with the Sen slope image and a preset significance level threshold, each pixel is classified and assigned a value through Boolean logic masking operation to generate a trend classification coding matrix containing five levels: significant improvement, slight improvement, stable, slight degradation, and significant degradation.

9. The remote sensing image batch index calculation and trend analysis system based on MATLABGUI according to claim 1, characterized in that, The visualization interaction module (150) adopts an architecture that separates the main thread and the background task thread. It is used to encapsulate and execute the calculation logic of the data acquisition and metadata parsing module (110), the heterogeneous data adaptive mapping module (120), the multidimensional feature matrix construction module (130), and the time series trend mining module (140) in the background task thread. The visualization interaction module (150) includes a timed listener, which is used to poll the execution status attributes of the background task thread, obtain the progress percentage value and log text information in real time, and update the progress bar component and log display field on the graphical user interface.

10. A method for batch index calculation and trend analysis of remote sensing images based on MATLAB GUI, characterized in that, The system for batch index calculation and trend analysis of remote sensing images based on MATLAB GUI as described in any one of claims 1-9 includes the following steps: The system receives the target folder path specified by the user through a graphical user interface and parses the timestamps and sensor types of the source data files under the path to generate a file index table. Based on the file index table and the preset band mapping lookup table, source data with different physical storage structures are remapped into standard logical spectral vectors of a unified dimension. Based on the standard logical spectral vector, the full feature input mode or the default feature input mode is adaptively selected, and the remote sensing ecological index matrix at each time node is calculated through principal component analysis transformation. The generated multiple remote sensing ecological index matrices are reassembled in chronological order into a three-dimensional spatiotemporal data cube in memory; The vectorized Sen slope estimation and Mann-Kendall significance test are performed on the three-dimensional spatiotemporal data cube to generate a trend classification coding matrix; The trend classification coding matrix is ​​converted into a raster image with geographic information and displayed in a graphical user interface.

Images