A habitat area extraction method and system based on remote sensing images

Abstract

The application discloses a kind of based on remote sensing image's habitat area extraction method and system, belong to image extraction technical field.The method specifically includes: based on real-time remote sensing data is preprocessed and feature enhancement;Different features are classified and extracted using multi-dimensional extraction method to form a multi-dimensional data set;Based on data set, using lightweight integrated learning, habitat is classified to form habitat classification atlas;Based on multiscale space, combined with context mechanism, context precision optimization and boundary optimization are carried out on habitat classification atlas to form the final habitat result map.The application proposes a kind of lightweight method specially for habitat extraction, by the combination of multi-feature fusion and lightweight machine learning model, to realize efficient, accurate habitat area identification, provide important foundation.

Description

Technical Field

[0001] This invention relates to the field of image extraction technology, and in particular to a method and system for extracting habitat areas based on remote sensing images. Background Technology

[0002] The accurate extraction of habitat areas, as a core foundation for biodiversity conservation, ecological environment assessment, and territorial spatial planning, directly impacts the scientific validity of ecological protection strategies and sustainable ecosystem management. In the development of remote sensing technology applied to habitat area extraction, identification methods based on single features have played a crucial role. Among these, spectral features have been the most widely used due to their ease of extraction from remote sensing images and their direct indicative nature in distinguishing land cover types. By capturing the differences in reflection, absorption, and emission of electromagnetic waves of different wavelengths by different land cover types, land cover types can be identified. For example, the strong reflectivity of vegetation in the near-infrared band can be used to identify vegetated habitats such as forests and grasslands, while the strong absorption characteristics of water bodies in the infrared band can be used to distinguish aquatic habitats. This has accumulated rich practical experience for habitat area extraction.

[0003] However, as ecological protection work becomes more refined, the limitations of traditional extraction methods have become increasingly apparent. Because they rely heavily on single features and do not fully consider the unique multi-scale characteristics of habitats, they are less effective at handling habitats with complex structures and cannot meet actual extraction needs. Summary of the Invention

[0004] The purpose of this invention is to provide a method and system for habitat area extraction based on remote sensing images. The method uses original multispectral images as raw data and processes and enhances image features using multiple algorithm techniques. This provides stable data input for subsequent extraction of habitat images using multi-scale feature algorithms and techniques. This solves the problem that the traditional methods mentioned above do not fully consider the unique multi-scale features of habitats, resulting in poor processing capabilities when dealing with complex habitat structures, which makes it difficult to meet actual extraction needs.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] Radiometric calibration and atmospheric correction were performed on the original multispectral images. The Retinex algorithm was used to enhance the image features and calculate the multi-source vegetation index. Redundant spectral information was removed and the main feature components were retained. Finally, multi-feature stacked data was obtained.

[0007] Based on multi-feature stacked data, a multi-dimensional collaborative extraction method is used to classify and extract different features and establish a multi-dimensional dataset.

[0008] Based on a multidimensional dataset, a pre-trained random forest model is used as a classifier. The model parameters are fine-tuned using a transfer learning strategy. The fine-tuned model is used as a multidimensional dataset processor. The model identifies and classifies habitats into multiple categories, and finally combines them to obtain a preliminary habitat classification map.

[0009] Based on the preliminary habitat classification map, multi-scale morphological filtering is used to smooth the internal regions while preserving edge details. The boundaries are then refined using a fast boundary optimization method to obtain the final habitat extraction result map.

[0010] Furthermore, radiometric calibration and atmospheric correction of the original multispectral imagery include the following steps:

[0011] S1: Suppose the original multispectral image has K bands. Radiometric calibration converts the DN values ​​of the image into reflectance. For each band K, the apparent reflectance formula for each pixel (i,j) is: P k (i,j)=Gain k *DN k (i,j)+Bias k , where P k Indicates reflectivity; Gain k and Bias k DN represents the scaling factor; k This represents the DN value of the k-th band;

[0012] S2: Atmospheric correction is performed on the reflectance using the dark pixel method, resulting in an atmospherically corrected reflectance image [R1, R2...R]. k ].

[0013] Furthermore, the process of enhancing image features and calculating multi-source vegetation indices using the Retinex algorithm after correction includes the following steps:

[0014] S1: For reflectance image [R1, R2...R... k Perform individual band processing. For each band k of each reflectance image, construct a multi-scale Gaussian convolution kernel. The construction formula is as follows: G n (x,y) represents the Gaussian convolution kernel function at the nth scale; (x,y) represents the relative coordinate position within the convolution kernel; K n σ represents the normalization coefficient for the nth scale; n The standard deviation of the Gaussian kernel controls the convolution scale; x 2 +y 2 Represents the square of the radial distance; This represents the calculation of the squared negative distance after standardization; three convolution kernels [σ1, σ2, σ3] of different scales are obtained through multi-scale Gaussian convolution kernels, corresponding to small size → medium size → large size respectively;

[0015] S2: Perform multi-scale Retinex calculations for each pixel in each band. Calculate the Retinex output of the k-th band at each scale based on three scales [σ1, σ2, σ3]. Finally, perform a weighted summation to obtain the enhanced image of each band [E1, E2...E...]. k ];

[0016] S3: Image based on each enhanced band [E1, E2...E k ] Calculate the NDVI and LSWI indices as vegetation indices;

[0017] S4: Image based on each enhanced band [E1, E2...E k Principal component analysis (PCA) dimensionality reduction of the image includes data centralization, covariance matrix calculation, eigenvalue decomposition, and feature selection reconstruction, ultimately yielding the principal component image (PC).

[0018] S5: Stack the results of steps S1-S4 above to form multi-feature stacked data F∈R. M×N×D .

[0019] Furthermore, the step of classifying and extracting different features includes the following:

[0020] S1: Based on multi-feature stacked data F∈R M×N×D Where M represents the image height; N represents the image width; D represents the number of features, with dimensions [B, 4, K], where B is the original number of bands; 4 is the number of exponential features; and K is the number of principal components, dividing the data into three types of features: spectral features, texture features, and geometric features.

[0021] S2: Perform statistical calculations on the D original feature channels of each pixel. The statistical data includes the mean, standard deviation, maximum value, and minimum value. Extract the NDVI and LSWI indices mentioned above, and add a simple ratio vegetation index, the formula of which is: SR=ρ NIR / ρ Red SR represents the simple ratio vegetation index; the higher the value, the more lush the vegetation and the higher the biomass. ρ NIR ρ represents the reflectance value of ground features in the near-infrared band. Red This represents the reflectance value of ground features in the red visible light band; a division operation is performed to highlight vegetation signals; finally, the above calculated parameters are combined to form the spectral feature matrix F1∈R. M×N×9 .

[0022] Furthermore, the method of using multi-dimensional collaborative extraction to classify and extract different features also includes the following steps:

[0023] S1: Using the GLCM algorithm, stacked multi-feature data F∈R M×N×D The image is quantized into 8 gray levels, and the quantization formula is:

[0024]

[0025] Where Q(x,y) represents the gray level; floor() represents the floor function; I(x,y) represents the pixel value of the original input image at position (x,y); I min Indicates the minimum pixel value of the current processing window or the entire image; I max This represents the maximum pixel value of the current processing window or the entire image; 8 represents the quantization level; this formula reduces the traditional 16 gray levels to 8, thus reducing its computational complexity.

[0026] S2: Perform multi-directional GLCM calculations, with the directions being 0°→45°→90°→135°. Set the window size to 7×7 and the step size to 1. Average the results from the four directions to obtain the comprehensive GLCM matrix P. Extract core features from the GLCM matrix P, including contrast, homogeneity, and entropy. Repeat the above calculations using different windows, including 3×3, 7×7, and 15×15, with scales ranging from small to medium to large. Finally, combine them to form the texture feature matrix F2∈R. M×N×9 ;

[0027] S3: Use the OTSU algorithm to binarize the mean image in the spectral features to obtain image B. Then, perform a distance transform on image B to calculate the Euclidean distance. The formula is as follows: Where D(i,j) represents the distance change result matrix; (i,j) represents the current foreground pixel coordinates; p, q represent the background pixel coordinates; min represents the minimum value function, which finds the nearest background pixel; √ represents the square root operation, which gives the target the true distance; the distance value obtained by the formula reflects the relative position of the pixel within the habitat patch;

[0028] S4: Calculate the gradient magnitude using the Sobel operator. Within a 7×7 sliding window, count the proportion of pixels whose gradient magnitude exceeds a threshold, normalize them to the range [0, 1], and obtain the edge density map ED(i,j). Finally, combine the above calculation parameters to form the geometric feature matrix F3∈R. M×N×3 ;

[0029] S5: Assume the above spectral feature matrix F1∈R M×N×9 Texture feature matrix F2∈R M×N×9 and geometric characteristic matrix F3∈R M×N×3Merge into F 综 ∈R M×N×21 Data set and output.

[0030] Furthermore, using a random forest model as a classifier and combining it with a transfer learning strategy to construct a preliminary habitat classification map includes the following steps:

[0031] S1: Based on F 综 ∈R M×N×21 The dataset was simplified using a multi-dimensional filtering strategy, including importance assessment, correlation analysis, and filtering strategies, resulting in the selection of the core feature dataset F. 核 ∈R M×N×K Where K=10 indicates that 10 core features are retained, and finally, the standardization process is performed and the output is completed.

[0032] S2: Set the initial parameters of the pre-trained random forest model. This model includes 100 decision trees, a similar habitat category system, and adjustable model parameters. Fine-tune the model parameters using a transfer learning mechanism. The transfer learning formula is: P 概 (c|X)=α×P 预 (c|X)+(1-α)×P 微 (c|X), where P 概 (c|X) represents the final category prediction probability; c represents the habitat index; X represents the input feature vector; α represents the weight coefficients of the pre-trained random forest model; P 预 (c|X) represents the predicted probability of the pre-trained model; P 微 (c|X) represents the predicted probability of the model after fine-tuning; the formula is used to calculate the fine-tuned pre-trained random forest model.

[0033] Furthermore, constructing a preliminary habitat classification map using a random forest model as a classifier and combining it with a transfer learning strategy also includes the following steps:

[0034] S1: Based on the fine-tuned pre-trained random forest model, a lightweight model with better computational efficiency is constructed. The training process includes parameter setting → parallel training → setting the prediction mechanism → setting the probability fusion mechanism. After training, a lightweight random forest classifier model RF is obtained.

[0035] S2: Using the lightweight random forest classifier model RF, class inference is performed on the entire image. The class inference process includes feature subset extraction → block prediction → class decision. Finally, the inference results are combined to form a preliminary habitat classification map.

[0036] Furthermore, the method of constructing the habitat result map using multi-scale morphological filtering and fast boundary optimization includes the following steps:

[0037] S1: Based on the preliminary habitat classification map, a morphological filtering algorithm is used to divide the map into a dual-scale structure of [3×3, 7×7] windows. Opening operation is performed on the 3×3 window scale, and closing operation is performed on the 7×7 window scale. The results of the two scales are fused using the majority voting principle to finally obtain a smoothed classification map.

[0038] S2: Boundary pixel identification is performed based on the smoothed classification map. Category changes are detected through a 3×3 window. When there are at least two different categories within the window, the center pixel is marked as a boundary pixel.

[0039] S3: For each boundary pixel, examine the category distribution within its 5×5 neighborhood. The strategy is to count the frequency of each category within the neighborhood. If the frequency of a neighboring category is greater than the frequency of the current category plus a threshold T, then the current pixel is reclassified as that neighboring category. The threshold T = 5.

[0040] S4: Apply 3×3 mean filtering to the boundary region, which is performed only when the boundary pixels are limited, and finally output the classification map after boundary optimization.

[0041] Furthermore, the method of constructing the habitat result map using multi-scale morphological filtering and fast boundary optimization also includes the following steps:

[0042] S1: Based on the optimized classification map, the flooding fill algorithm is used to calculate the number of pixels in each connected region. When the region area is less than 10 pixels, it is replaced with the surrounding main category to remove small patches.

[0043] S2: Set ecological rules to correct the category. The rule is to change the category of small water bodies completely surrounded by vegetation to wetlands and single pixels to replace the main category of the 3×3 neighborhood and check the consistency of river width.

[0044] S3: Using a 3×3 morphological closing operation, fill the remaining small holes to obtain the final habitat extraction result image.

[0045] Furthermore, the habitat area extraction system based on remote sensing imagery specifically includes:

[0046] The first module is used to perform radiometric calibration and atmospheric correction on the original multispectral images, enhance image features and calculate multi-source vegetation indices, remove redundant spectral information and retain the main feature components, and finally stack multi-feature data.

[0047] The second module is used to classify and extract different features and build a multi-dimensional dataset;

[0048] The third module is used to build a pre-trained random forest model as a classifier. The model is fine-tuned by combining transfer learning strategies. The fine-tuned model is used as a multi-dimensional dataset processor. The model identifies and classifies habitats into multiple categories, and finally combines them to form a preliminary habitat classification map.

[0049] The fourth module is used to smooth the internal regions while preserving edge details, refine the boundaries, and construct the final habitat extraction result map.

[0050] Compared with the prior art, the present invention has at least one of the following technical effects:

[0051] 1. This invention utilizes the Retinex algorithm to enhance multi-scale image features to eliminate the influence of terrain shadows on vegetation reflectance, and integrates spectral reflectance, vegetation index and principal component features to achieve multi-scale feature data stacking. Compared with traditional methods, it has the technical advantage of multi-scale feature fusion, and solves the problem that traditional methods rely on single features, resulting in poor processing capability when dealing with complex habitat structures.

[0052] 2. This invention utilizes a multi-scale feature extraction method to extract features from three layers: spectrum, texture, and geometry. It extracts only the core features of the features, simplifying the extraction complexity and introducing a collaborative mechanism. Compared with traditional methods, it has the technical advantages of being more lightweight and having clear ecological significance, thus solving the problems of insufficient lightweightness and unclear ecological significance of traditional methods.

[0053] 3. This invention is based on lightweight ensemble learning for habitat classification. It introduces a pre-trained random forest model as the foundation and incorporates a transfer learning mechanism to fine-tune the model's internal parameters, making its processing capabilities more accurate. Compared with traditional methods, it has the technical advantage of more accurate image classification processing models and solves the problem of high model complexity in traditional methods, which leads to wasted computing resources and poor transferability.

[0054] 4. This invention utilizes dual-scale morphological filtering as a method for smoothing habitat images and employs a fast boundary optimization method to optimize the boundaries of habitat classification maps. Compared with traditional methods, it has the technical advantage of simplifying the optimization process and solves the problem that traditional methods suffer from poor accuracy due to large computational loads caused by complex optimization processes.

[0055] 5. This invention forms a complete data link through data flow processing and feature management strategies, ensuring that data is uniformly managed at each processing stage. Compared with traditional methods, it has the technical advantage of data normalization and solves the problem of data silos and fragmentation in traditional systems. Attached Figure Description

[0056] To more clearly illustrate the technical solutions in the embodiments or examples of the present invention, the drawings used in the embodiments or examples will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained according to these drawings without creative effort.

[0057] Figure 1 This is a schematic diagram of the process for a habitat area extraction method based on remote sensing imagery according to the present invention.

[0058] Figure 2 This is a system architecture diagram for habitat area extraction based on remote sensing imagery, according to the present invention. Detailed Implementation

[0059] In the following description, specific details such as particular system architectures and techniques are set forth for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details. In other instances, detailed descriptions of well-known systems, apparatuses, circuits, and methods have been omitted so as not to obscure the description of this application with unnecessary detail.

[0060] It should be understood that, when used in this application specification and the appended claims, the term "comprising" indicates the presence of the described features, integrals, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components and / or a collection thereof.

[0061] It should also be understood that the term “and / or” as used in this application specification and the appended claims means any combination of one or more of the associated listed items and all possible combinations, and includes such combinations.

[0062] As used in this application specification and the appended claims, the term "if" may be interpreted, depending on the context, as "when," "once," "in response to determination," or "in response to detection." Similarly, the phrase "if determined" or "if detected [the described condition or event]" may be interpreted, depending on the context, as meaning "once determined," "in response to determination," "once detected [the described condition or event]," or "in response to detection [the described condition or event]."

[0063] Furthermore, in the description of this application and the appended claims, the terms "first," "second," "third," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0064] References to "one embodiment" or "some embodiments" as described in this specification mean that one or more embodiments of this application include a specific feature, structure, or characteristic described in connection with that embodiment. Therefore, the phrases "in one embodiment," "in some embodiments," "in other embodiments," "in still other embodiments," etc., appearing in different parts of this specification do not necessarily refer to the same embodiment, but rather mean "one or more, but not all, embodiments," unless otherwise specifically emphasized. The terms "comprising," "including," "having," and variations thereof mean "including but not limited to," unless otherwise specifically emphasized.

[0065] The specific implementation method is as follows:

[0066] See appendix Figure 1 -Appendix Figure 2 This invention discloses a method for extracting habitat areas based on remote sensing imagery, the method specifically comprising:

[0067] Radiometric calibration and atmospheric correction were performed on the original multispectral images. The Retinex algorithm was used to enhance the image features and calculate the multi-source vegetation index. Redundant spectral information was removed and the main feature components were retained. Finally, multi-feature stacked data was obtained.

[0068] Based on multi-feature stacked data, a multi-dimensional collaborative extraction method is used to classify and extract different features and establish a multi-dimensional dataset.

[0069] Based on a multidimensional dataset, a pre-trained random forest model is used as a classifier. The model parameters are fine-tuned using a transfer learning strategy. The fine-tuned model is used as a multidimensional dataset processor. The model identifies and classifies habitats into multiple categories, and finally combines them to obtain a preliminary habitat classification map.

[0070] Based on the preliminary habitat classification map, multi-scale morphological filtering is used to smooth the internal regions while preserving edge details. The boundaries are then refined using a fast boundary optimization method to obtain the final habitat extraction result map.

[0071] The above-mentioned method for extracting habitat areas based on remote sensing images consists of sequentially connected steps, each with specific features. The implementation of the above method can be carried out according to the steps, and after implementation, a complete image extraction and optimization chain can be formed based on original multispectral images - multi-feature data stacking processing - multi-scale feature extraction - classifier model construction - habitat identification and mapping - map optimization.

[0072] In practice, radiometric calibration and atmospheric correction of the original multispectral imagery include the following steps:

[0073] S1: Suppose the original multispectral image has K bands. Radiometric calibration converts the DN values ​​of the image into reflectance. For each band K, the apparent reflectance formula for each pixel (i,j) is: P k (i,j)=Gain k *DN k (i,j)+Bias k , where P k Indicates reflectivity; Gain k and Bias k DN represents the scaling factor; k This represents the DN value of the k-th band;

[0074] S2: Atmospheric correction is performed on the reflectance using the dark pixel method, resulting in an atmospherically corrected reflectance image [R1, R2...R]. k ].

[0075] The atmospheric correction of reflectance using the dark pixel method described above includes the following processing steps:

[0076] Step 1) Select the pixel with the lowest reflectance of 0.1% in the shortwave infrared band as the dark pixel;

[0077] Step 2) Perform path radiation calculation, the formula is: L b =Percentiles(ρ) 像 ,95%), of which L b ρ represents the atmospheric path radiative reflectance contribution of band b; 像 This represents a dark pixel; 95% represents the 95th percentile of the reflectance of a dark pixel; the purpose of the formula is to estimate the contribution of atmospheric path radiation to apparent reflectance and avoid the influence of outliers.

[0078] Step 3) Calculate the true surface reflectance. For each pixel (i,j) in each band b, the true reflectance formula is: ρ 反 (i,j)=ρ 像 (i,j)-L p , where ρ 反 (i,j) represents the true surface reflectance of band b; L b The atmospheric path radiative reflectance contribution of band b (obtained from step 2 above); ρ 像 (i,j) represents each pixel (i,j) of each band b; the formula subtracts the path radiation effect from the apparent reflectance to obtain the true surface reflectance;

[0079] Step 4) Set the reflectivity constraint range, and the constraint strategy is: when ρ 反 If (i,j)<0, then ρ 反 (i,j)=0.001; when ρ反 If (i,j)>1, then ρ 反 (i,j)=0.999; the purpose is to ensure that the reflectivity is within the physically reasonable range of [0,1] to complete atmospheric correction.

[0080] In specific implementation, the process of enhancing image features and calculating multi-source vegetation indices using the Retinex algorithm after correction includes the following steps:

[0081] S1: For reflectance image [R1, R2...R... k Perform individual band processing. For each band k of each reflectance image, construct a multi-scale Gaussian convolution kernel. The construction formula is as follows: G n (x,y) represents the Gaussian convolution kernel function at the nth scale; (x,y) represents the relative coordinate position within the convolution kernel; K n σ represents the normalization coefficient for the nth scale; n The standard deviation of the Gaussian kernel controls the convolution scale; x 2 +y 2 Represents the square of the radial distance; This represents the calculation of the squared negative distance after standardization; three convolution kernels [σ1, σ2, σ3] of different scales are obtained through multi-scale Gaussian convolution kernels, corresponding to small size → medium size → large size respectively;

[0082] S2: Perform multi-scale Retinex calculations for each pixel in each band. Calculate the Retinex output of the k-th band at each scale based on three scales [σ1, σ2, σ3]. Finally, perform a weighted summation to obtain the enhanced image of each band [E1, E2...E...]. k ];

[0083] S3: Image based on each enhanced band [E1, E2...E k ] Calculate the NDVI and LSWI indices as vegetation indices;

[0084] S4: Image based on each enhanced band [E1, E2...E k Principal component analysis (PCA) dimensionality reduction of the image includes data centralization, covariance matrix calculation, eigenvalue decomposition, and feature selection reconstruction, ultimately yielding the principal component image (PC).

[0085] S5: Stack the results of steps S1-S4 above to form multi-feature stacked data F∈R. M×N×D .

[0086] The above-mentioned multi-scale Retinex calculation and weighted fusion for each pixel in each band includes the following steps:

[0087] Step 1) Calculate multi-scale Retinex components based on three scales [σ1, σ2, σ3]. The formula for calculating the Retinex components is: S k σ(x,y)=log(R k (x,y))-log(Gσ(x,y)*R k (x,y)), where S k σ(x,y) represents the Retinex output component of band k at scale σ; R k (x,y) represents the reflectance value of band k at position (x,y) (i.e., the kth band of the atmospherically corrected reflectance image); (Gσ(x,y)*R k (x,y)) represents the Gaussian convolution operation, which represents the illumination estimation at scale σ; log represents the natural logarithm function (to enhance contrast); the Retinex components at three scales (i.e., [Sσ1,Sσ2,Sσ3]) are calculated using the formula.

[0088] Step 2) Perform a weighted summation of the Retinex components [Sσ1, Sσ2, Sσ3] obtained in Step 1. The weighting formula is: MSR k (x,y)=w1×Sσ1(x,y)+w2×Sσ2(x,y)+w3×Sσ3(x,y), where MSR k (x,y) represents the multi-scale Retinex fusion result of band k; Sσ1, Sσ2 and Sσ3 represent the Retinex components of the three scales; w1, w2 and w3 represent the weights of the three scales (the weights are distributed as follows: 0.4→0.35→0.25); by fusing multi-scale information through formulas, detail enhancement and overall brightness balance are achieved.

[0089] Step 3) Based on the Retinex component fusion result and the atmospherically corrected reflectance image, perform color restoration and contrast stretching, using the formula: p k (x,y)=α×MSR k (x,y)×C k (x,y)+β, where p k (x,y) represents the enhanced reflectivity of band k at position (x,y); α represents the contrast enhancement factor (α=1.2 to avoid excessive enhancement and noise); β represents the brightness shift coefficient (β=0.1 to prevent negative reflectivity); MSR k (x,y) represents the multi-scale Retinex fusion result of band k (obtained from step 2 above); C k (x,y) represents the color retention factor; the formula restores color accuracy and adjusts contrast and brightness to a suitable range.

[0090] The formula for calculating the color retention factor is as follows:

[0091]

[0092] Where C k (x,y) represents the color retention factor; k represents the current band index; (x,y) represents the pixel spatial coordinates; b represents the band cycle index; B represents the total number of bands; ∑ represents summation (i.e., summation over all bands); Let represent the surface reflectance of the k-th band at location (x,y); let represent the surface reflectance of the b-th band at location (x,y). The table shows the reflectance (each value in the summation operation represents the reflectance contribution of each band); the significance of division is to calculate the proportion of the current band's reflectance in the total reflectance.

[0093] The above calculations use NDVI and LSWI as vegetation indices. The formula for calculating the NDVI index is as follows:

[0094]

[0095] Where NDVI(i,j) represents the NDVI value at position (i,j); ρ NIR (i,j) represents the reflectance value of the near-infrared band at position (i,j); ρ Red (i,j) represents the reflectance value of the red light band at position (i,j); (i,j) represents the pixel spatial coordinates; this value quantifies vegetation cover and growth status, and the value range is between [-1, 1] (high value 0.3-0.8: healthy and dense vegetation, medium value 0.1-0.3: sparse or stressed vegetation, low value <0.1: non-vegetated area); the LSWI index calculation formula is:

[0096]

[0097] Where LSWI(i,j) represents the LSWI value at position (i,j); ρ NIR (i,j) represents the reflectance value of the near-infrared band at position (i,j); ρ SWIR (i,j) represents the reflectance value of the shortwave infrared band at position (i,j); (i,j) represents the pixel spatial coordinates; this value can detect surface moisture content and vegetation water stress, identify water bodies and flooded areas, and monitor soil moisture and vegetation moisture status. The value range is between [-1, 1] (high positive value > 0.2: water body or extremely wet surface, medium positive value 0-0.2: wet vegetation or soil, low value or negative value: dry surface or water-stressed vegetation).

[0098] The above-mentioned principal component analysis (PCA) dimensionality reduction for each enhanced band image includes data centralization → covariance matrix calculation → eigenvalue decomposition → feature selection and reconstruction. The specific implementation steps are as follows:

[0099] Step 1) For each enhanced band image, reconstruct the 3D image data into a 2D matrix (P×B), calculate the mean of each band, and subtract the corresponding band mean from each band value to obtain the centered data matrix X∈R. （P×B） ;

[0100] Step 2): Based on the result of Step 1), calculate its transpose matrix, perform matrix multiplication, and then divide by the total number of pixels P to normalize and obtain the covariance matrix. The formula is: C = (1 / P) × X T ×X∈R （P×B） Where C represents the covariance matrix; X T Represents the transpose of the centered matrix; P represents the total number of pixels; B represents the number of bands; X∈R （P×B） This represents the centralized data matrix (obtained from step 1 above).

[0101] Step 3) Based on the results of Step 2), perform eigenvalue decomposition on the covariance matrix C, arrange the eigenvalues ​​and eigenvectors in descending order of eigenvalue size, ensuring that the eigenvectors are unit vectors, and finally obtain the eigenvalue vector λ=[λ1,λ2...λB] and the eigenvector matrix V=[v1,v2...vB]∈R. （B×B） ;

[0102] Step 4) Based on the results of Step 3), calculate the cumulative contribution rate to determine the number of principal components K. The formula for calculating the cumulative contribution rate is: , where η K λ represents the cumulative contribution rate of the first K principal components; K represents the number of principal components selected; λ k ∑ represents the k-th eigenvalue; ∑ represents summation (the formula describes the summation of the numerator from the 1st to the Kth component; the denominator sums over all principal components); the purpose of division is to describe the proportional relationship between the partial components and the whole components;

[0103] Further, a projection matrix composed of the first K feature vectors is constructed to project the centered data onto the principal component space. Finally, the result is reconstructed into a three-dimensional image format to obtain the principal component image PC.

[0104] In practice, the steps for classifying and extracting different features include the following:

[0105] S1: Based on multi-feature stacked data F∈R M×N×DWhere M represents the image height; N represents the image width; D represents the number of features, with dimensions [B, 4, K], where B is the original number of bands; 4 is the number of exponential features; and K is the number of principal components, dividing the data into three types of features: spectral features, texture features, and geometric features.

[0106] S2: Perform statistical calculations on the D original feature channels of each pixel. The statistical data includes the mean, standard deviation, maximum value, and minimum value. Extract the NDVI and LSWI indices mentioned above, and add a simple ratio vegetation index, the formula of which is: SR=ρ NIR / ρ Red SR represents the simple ratio vegetation index; the higher the value, the more lush the vegetation and the higher the biomass. ρ NIR ρ represents the reflectance value of ground features in the near-infrared band. Red This represents the reflectance value of ground features in the red visible light band; a division operation is performed to highlight vegetation signals; finally, the above calculated parameters are combined to form the spectral feature matrix F1∈R. M×N×9 .

[0107] In practice, the method of using multi-dimensional collaborative extraction to classify and extract different features also includes the following steps:

[0108] S1: Using the GLCM algorithm, stacked multi-feature data F∈R M×N×D The image is quantized into 8 gray levels, and the quantization formula is:

[0109]

[0110] Where Q(x,y) represents the gray level; floor() represents the floor function; I(x,y) represents the pixel value of the original input image at position (x,y); I min Indicates the minimum pixel value of the current processing window or the entire image; I max This represents the maximum pixel value of the current processing window or the entire image; 8 represents the quantization level; this formula reduces the traditional 16 gray levels to 8, thus reducing its computational complexity.

[0111] S2: Perform multi-directional GLCM calculations, with the directions being 0°→45°→90°→135°. Set the window size to 7×7 and the step size to 1. Average the results from the four directions to obtain the comprehensive GLCM matrix P. Extract core features from the GLCM matrix P, including contrast, homogeneity, and entropy. Repeat the above calculations using different windows, including 3×3, 7×7, and 15×15, with scales ranging from small to medium to large. Finally, combine them to form the texture feature matrix F2∈R. M×N×9 ;

[0112] S3: Use the OTSU algorithm to binarize the mean image in the spectral features to obtain image B. Then, perform a distance transform on image B to calculate the Euclidean distance. The formula is as follows: Where D(i,j) represents the distance change result matrix; (i,j) represents the current foreground pixel coordinates; p, q represent the background pixel coordinates; min represents the minimum value function, which finds the nearest background pixel; √ represents the square root operation, which gives the target the true distance; the distance value obtained by the formula reflects the relative position of the pixel within the habitat patch;

[0113] S4: Calculate the gradient magnitude using the Sobel operator. Within a 7×7 sliding window, count the proportion of pixels whose gradient magnitude exceeds a threshold, normalize them to the range [0, 1], and obtain the edge density map ED(i,j). Finally, combine the above calculation parameters to form the geometric feature matrix F3∈R. M×N×3 ;

[0114] S5: Assume the above spectral feature matrix F1∈R M×N×9 Texture feature matrix F2∈R M×N×9 and geometric characteristic matrix F3∈R M×N×3 Merge into F 综 ∈R M×N×21 Data set and output.

[0115] The above multi-directional GLCM calculations, averaging the results from the four directions to obtain the comprehensive GLCM matrix P, and extracting the core features of contrast, homogeneity, and entropy, are implemented as follows:

[0116] Step 1) For the four directions, calculate and average using the GLCM general formula, which is:

[0117]

[0118] Where Pavg(i,j) represents the average co-occurrence probability; (i,j) represents the gray level; ∑ represents summation (i.e., summing the results over the four directions); θ represents the set of directions; I(x,y) represents the gray value at the current coordinate (x,y); I(x+dx θ ,y+dy θ ) represents the gray level of the neighborhood location; dx,dy represent the directional offset; δ represents the condition count (the count is 1 when the gray level pairing condition is met, and 0 otherwise); 1 / 4 represents the average coefficient of the four directions; the average GLCM of the four directions can be calculated using this general formula;

[0119] Step 2): Calculate the contrast based on the results of Step 1), using the following formula: C 对The formula represents the contrast result; (i,j) represents the gray level index (i.e., the quantized gray level, with a value range of 0-7, a total of 8 levels); P(i,j) represents the GLCM probability matrix element (i.e., the probability that gray levels i and j appear simultaneously in a specific spatial relationship); ∑∑ represents double summation (the purpose is to traverse all gray level combinations, i.e., 8×8=64 combinations); the result calculated by this formula is used to measure the severity of local gray level changes in the image;

[0120] Step 3) Calculate homogeneity based on the results of Step 1), using the following formula: H 同 This represents the homogeneity calculation result; (i,j) 2 P(i,j) represents the squared difference in gray levels; P(i,j) represents the element of the GLCM probability matrix (the specific meaning is consistent with the description in step 2 above); when i=j, the denominator is 1 and has the largest weight; when |ij| increases, the denominator increases rapidly and the weight decreases; the result calculated by this formula is used to measure the uniformity and smoothness of local areas of the image.

[0121] Step 4) Calculate the entropy based on the result of Step 1), using the following formula:

[0122]

[0123] Where E 熵 The result represents the entropy calculation result; P(i,j) represents the element of the GLCM probability matrix (the specific meaning is consistent with the description in steps 2 and 3 above); log represents the natural logarithm function (used to measure the amount of information); ε represents a very small positive number (to prevent mathematical overflow of log(0), take 1e-8); the purpose of using the negative sign is to ensure that the entropy value is positive; the result calculated by this formula is used to measure the complexity and randomness of the image texture.

[0124] In practice, using a random forest model as a classifier and combining it with a transfer learning strategy to construct a preliminary habitat classification map includes the following steps:

[0125] S1: Based on F 综 ∈R M×N×21 The dataset was simplified using a multi-dimensional filtering strategy, including importance assessment, correlation analysis, and filtering strategies, resulting in the selection of the core feature dataset F. 核 ∈R M×N×K Where K=10 indicates that 10 core features are retained, and finally, the standardization process is performed and the output is completed.

[0126] S2: Set the initial parameters of the pre-trained random forest model. This model includes 100 decision trees, a similar habitat category system, and adjustable model parameters. Fine-tune the model parameters using a transfer learning mechanism. The transfer learning formula is: P概 (c|X)=α×P 预 (c|X)+(1-α)×P 微 (c|X), where P 概 (c|X) represents the final category prediction probability; c represents the habitat index; X represents the input feature vector; α represents the weight coefficients of the pre-trained random forest model; P 预 (c|X) represents the predicted probability of the pre-trained model; P 微 (c|X) represents the predicted probability of the model after fine-tuning; the formula is used to calculate the fine-tuned pre-trained random forest model.

[0127] The above-mentioned multi-dimensional screening strategy for simplifying the dataset includes importance assessment, correlation analysis, and screening strategy. The specific implementation steps are as follows:

[0128] Step 1) Construct a temporary random forest (using 30 trees) and train it on the training samples. For each tree, calculate the Gini impurity reduction of each feature at all split nodes. The formula for calculating the Gini impurity reduction is as follows:

[0129]

[0130] Gini 基 Indicates the reduction in impurity of the gin; N 左样 N represents the number of samples in the left child node. 右样 Indicates the number of samples in the right child node; Gini 左子 Indicates the Gini impurity of the left child node; 右子 This represents the Gini impurity of the right child node; the purpose of the division operation is to measure the proportion of the total samples of the parent node shared by the left and right child nodes. and This represents the proportion of the left and right child nodes to the total number of samples of the parent node, respectively. Then, the average of all trees is calculated to obtain the final importance score of each feature, and normalization is performed so that the sum of the importance of all features is 1.

[0131] Step 2) Calculate the 20×20 feature correlation coefficient matrix. The formula for the feature correlation coefficient (Pearson correlation coefficient) is:

[0132]

[0133] Where Corr(d1,d2) represents the correlation coefficient between features d1 and d2; Xd1 and Xd2 represent the vectors of all sample values ​​of features d1 and d2; μd1 and μd2 represent the sample means of features d1 and d2; ∑ represents the summation operator; P represents the total number of samples; σd1 and σd2 represent the sample standard deviations of features d1 and d2; the numerator and denominator represent the covariance of the two features multiplied by the total number of samples minus 1; the purpose of the division operation is to standardize the correlation measure to the range of [-1,1]; further, a correlation threshold |Corr|>0.85 is set as highly correlated, and for highly correlated feature pairs, the feature with higher importance is retained;

[0134] Step 3) For the most important features, retain the top 15 features (75% of the total 20 features). For each highly relevant feature group, if all features in the group are in the candidate features, only retain the highest-ranking feature in the group and remove the other features. If the number of features after filtering is greater than 12, remove features from low to high importance until 12 features remain. If the number of features after filtering is less than 12, supplement features from high to low importance until 12 features are reached.

[0135] In practice, using a random forest model as a classifier and combining it with a transfer learning strategy to construct a preliminary habitat classification map also includes the following steps:

[0136] S1: Based on the fine-tuned pre-trained random forest model, a lightweight model with better computational efficiency is constructed. The training process includes parameter setting → parallel training → setting the prediction mechanism → setting the probability fusion mechanism. After training, a lightweight random forest classifier model RF is obtained.

[0137] S2: Using the lightweight random forest classifier model RF, class inference is performed on the entire image. The class inference process includes feature subset extraction → block prediction → class decision. Finally, the inference results are combined to form a preliminary habitat classification map.

[0138] The above-mentioned lightweight model with better computational efficiency includes a training process that includes parameter setting, parallel training, setting a prediction mechanism, and setting a probability fusion mechanism. The specific implementation steps are as follows:

[0139] Step 1) The parameter settings are as follows: Decision trees = 50 trees, maximum tree depth limit = 10, minimum number of split samples per node = 20, number of features used per decision tree = 8, auto-sampling = TRUE and random seed = 42;

[0140] Step 2) Parallel training is specifically set up as follows: main process task allocation and result collection, worker process implementation of 4-8 parallel processes, each process processing 6-12 trees, and communication mechanism using shared memory to reduce data copying;

[0141] Step 3) The prediction mechanism is specifically set as follows: the large dataset is divided into 256×256 pixel blocks, the peak memory usage is controlled within <2GB, and pipeline processing is implemented (i.e., the next block of data is loaded when the current block is predicted). During prediction, the feature subset mapping during training is directly used to avoid repeated feature selection calculations and maintain the consistency between training and prediction.

[0142] Step 4) Set the probability fusion mechanism as follows: Use soft voting (i.e. probability averaging) and the final probability = 0.7 × pre-trained model probability + 0.3 × lightweight random forest probability.

[0143] In practice, the method of constructing the habitat result map using multi-scale morphological filtering and fast boundary optimization includes the following steps:

[0144] S1: Based on the preliminary habitat classification map, a morphological filtering algorithm is used to divide the map into a dual-scale structure of [3×3, 7×7] windows. Opening operation is performed on the 3×3 window scale, and closing operation is performed on the 7×7 window scale. The results of the two scales are fused using the majority voting principle to finally obtain a smoothed classification map.

[0145] S2: Boundary pixel identification is performed based on the smoothed classification map. Category changes are detected through a 3×3 window. When there are at least two different categories within the window, the center pixel is marked as a boundary pixel.

[0146] S3: For each boundary pixel, examine the category distribution within its 5×5 neighborhood. The strategy is to count the frequency of each category within the neighborhood. If the frequency of a neighboring category is greater than the frequency of the current category plus a threshold T, then the current pixel is reclassified as that neighboring category. The threshold T = 5.

[0147] S4: Apply 3×3 mean filtering to the boundary region, which is performed only when the boundary pixels are limited, and finally output the classification map after boundary optimization.

[0148] In practice, the method of constructing the habitat result map using multi-scale morphological filtering and fast boundary optimization also includes the following steps:

[0149] S1: Based on the optimized classification map, the flooding fill algorithm is used to calculate the number of pixels in each connected region. When the region area is less than 10 pixels, it is replaced with the surrounding main category to remove small patches.

[0150] S2: Set ecological rules to correct the category. The rule is to change the category of small water bodies completely surrounded by vegetation to wetlands and single pixels to replace the main category of the 3×3 neighborhood and check the consistency of river width.

[0151] S3: Using a 3×3 morphological closing operation, fill the remaining small holes to obtain the final habitat extraction result image.

[0152] The above steps involve opening operations on a 3×3 window and closing operations on a 7×7 window. The results from the two scales are then merged using a majority voting principle. The specific implementation steps are as follows:

[0153] Step 1), the opening operation formula is: C 3×3 =Dilation(Erosion(C,S 3×3 ),S 3×3 ), where C 3×3 Represents the result of the opening operation; Dilation represents the dilation operation; Erosion represents the erosion operation; C represents the preliminary habitat classification map; S 3×3 This represents a 3×3 window scale structure (i.e., a 3x3 row 3 column structure); the opening operation is dilation → erosion, therefore the Dilation function is used as the outer operation function.

[0154] Step 2), the closing operation formula is: C 7×7 =Erosion(Dilation(C,C 7×7 ),C 7×7 ), where C 7×7 Represents the result of the closing operation; Dilation represents the dilation operation; Erosion represents the erosion operation; C represents the preliminary habitat classification map; C 7×7 This represents a 7×7 window scale structure (i.e., a 7-row, 7-column structure); the opening operation is erosion → dilation, therefore the Erosion function is used as the outer operation function;

[0155] Step 3) Based on the results of Step 1) and Step 3), a vote is performed. The voting strategy is to count the occurrence frequency of each of the three input categories. If the occurrence frequency of a category is ≥2, the category is selected as the final result. Otherwise, it is considered that all three categories are different, and C(i,j) is selected as the final result (i.e., the original classification is maintained). The decision logic is summarized as follows: unanimous agreement (3 votes in the same category) is adopted directly; majority agreement (2 votes in the same category) is adopted; complete disagreement (3 votes in the same category) adopts a conservative strategy and maintains the original classification.

[0156] In specific implementation, a habitat area extraction system based on remote sensing imagery includes:

[0157] The first module is used to perform radiometric calibration and atmospheric correction on the original multispectral images, enhance image features and calculate multi-source vegetation indices, remove redundant spectral information and retain the main feature components, and finally stack multi-feature data.

[0158] The second module is used to classify and extract different features and build a multi-dimensional dataset;

[0159] The third module is used to build a pre-trained random forest model as a classifier. The model is fine-tuned by combining transfer learning strategies. The fine-tuned model is used as a multi-dimensional dataset processor. The model identifies and classifies habitats into multiple categories, and finally combines them to form a preliminary habitat classification map.

[0160] The fourth module is used to smooth the internal regions while preserving edge details, refine the boundaries, and construct the final habitat extraction result map.

[0161] Working principle: This method converts raw images into deep feature data through physical correction and mathematical enhancement, providing reliable input for subsequent analysis; it comprehensively describes habitat characteristics from three dimensions: spectral, texture, and geometry, and captures multiple attributes of habitats through feature fusion; it achieves classification with limited labeled data while maintaining computational efficiency by utilizing ensemble learning and transfer learning techniques; it optimizes the spatial continuity and ecological rationality of classification results through spatial context information and ecological knowledge; and it realizes deep extraction of habitat structure maps, providing a practical habitat image extraction technology for applications such as habitat monitoring, biodiversity conservation, and ecological environment assessment.

[0162] Those skilled in the art will recognize that the modules, units, and algorithm implementation steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, computer software, or a combination of electronic hardware and software. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

Claims

1. A method for habitat area extraction based on remote sensing imagery, characterized in that, The method specifically includes: Radiometric calibration and atmospheric correction were performed on the original multispectral images. The Retinex algorithm was used to enhance the image features and calculate the multi-source vegetation index. Redundant spectral information was removed and the main feature components were retained. Finally, multi-feature stacked data was obtained. Based on multi-feature stacked data, a multi-dimensional collaborative extraction method is used to classify and extract different features and establish a multi-dimensional dataset. Based on a multidimensional dataset, a pre-trained random forest model is used as a classifier. The model parameters are fine-tuned using a transfer learning strategy. The fine-tuned model is used as a multidimensional dataset processor. The model identifies and classifies habitats into multiple categories, and finally combines them to obtain a preliminary habitat classification map. Based on the preliminary habitat classification map, multi-scale morphological filtering is used to smooth the internal region while preserving edge details. The boundary is then refined using a fast boundary optimization method to obtain the final habitat extraction result map. Radiometric calibration and atmospheric correction of raw multispectral images include the following steps: S1: Suppose the original multispectral image has K bands. Radiometric calibration converts the DN values ​​of the image into reflectance. For each band k, the apparent reflectance formula for each pixel (i,j) is: P k (i,j)=Gain k *DN k (i,j)+Bias k , where P k Indicates reflectivity; Gain k and Bias k DN represents the scaling factor; k This represents the DN value of the k-th band; S2: Atmospheric correction is performed on the reflectance using the dark pixel method, resulting in an atmospherically corrected reflectance image [R1, R2...R]. k ]; The process of enhancing image features and calculating multi-source vegetation indices using the Retinex algorithm after correction includes the following steps: S1: For reflectance image [R1, R2...R] k Perform individual band processing, and construct a multi-scale Gaussian convolution kernel for each band of each reflectance image. The construction formula is as follows: G n (x,y) represents the Gaussian convolution kernel function at the nth scale; (x,y) represents the relative coordinate position within the convolution kernel; K n σ represents the normalization coefficient for the nth scale; n The standard deviation of the Gaussian kernel controls the convolution scale; x 2 +y 2 Represents the square of the radial distance; This represents the calculation of the squared negative distance after standardization; three convolution kernels [σ1, σ2, σ3] of different scales are obtained through multi-scale Gaussian convolution kernels, corresponding to small size → medium size → large size respectively; S2: Perform multi-scale Retinex calculations for each pixel in each band. Calculate the Retinex output of the k-th band at each scale based on three scales [σ1, σ2, σ3]. Finally, perform a weighted summation to obtain the enhanced image of each band [E1, E2...E...]. k ]; S3: Image based on each enhanced band [E1, E2...E k ] Calculate the NDVI and LSWI indices as vegetation indices; S4: Image based on each enhanced band [E1, E2...E k Principal component analysis (PCA) dimensionality reduction of the image includes data centralization, covariance matrix calculation, eigenvalue decomposition, and feature selection reconstruction, ultimately yielding the principal component image (PC). S5: Stack the results of steps S1-S4 above to form multi-feature stacked data F∈R. M×N×D .

2. The method for extracting habitat areas based on remote sensing imagery according to claim 1, characterized in that, The steps for classifying and extracting different features include the following: S1: Based on multi-feature stacked data F∈R M×N×D Where M represents the image height; N represents the image width; D represents the number of features, with dimensions [B, 4, K], where B is the original number of bands; 4 is the number of exponential features; and K is the number of principal components, dividing the data into three types of features: spectral features, texture features, and geometric features. S2: Perform statistical calculations on the D original feature channels of each pixel. The statistical data includes the mean, standard deviation, maximum value, and minimum value. Extract the NDVI and LSWI indices mentioned above, and add a simple ratio vegetation index, the formula of which is: SR=ρ NIR / ρ Red SR represents the simple ratio vegetation index; the higher the value, the more lush the vegetation and the higher the biomass. ρ NIR ρ represents the reflectance value of ground features in the near-infrared band. Red This represents the reflectance value of ground features in the red visible light band; a division operation is performed to highlight vegetation signals; finally, the above calculated parameters are combined to form the spectral feature matrix F1∈R. M×N×9 .

3. The method for extracting habitat areas based on remote sensing imagery according to claim 2, characterized in that, The method of using multi-dimensional collaborative extraction to classify and extract different features also includes the following steps: S1: Using the GLCM algorithm, stacked multi-feature data F∈R M×N×D The image is quantized into 8 gray levels, and the quantization formula is: Where Q(x,y) represents the gray level; floor() represents the floor function; I(x,y) represents the pixel value of the original input image at position (x,y); I min Indicates the minimum pixel value of the current processing window or the entire image; I max This represents the maximum pixel value of the current processing window or the entire image; 8 represents the quantization level; this formula reduces the traditional 16 gray levels to 8, thus reducing its computational complexity. S2: Perform multi-directional GLCM calculations, with the directions being 0°→45°→90°→135°. Set the window size to 7×7 and the step size to 1. Average the results from the four directions to obtain the comprehensive GLCM matrix P. Extract core features from the GLCM matrix P, including contrast, homogeneity, and entropy. Repeat the above calculations using different windows, including 3×3, 7×7, and 15×15, with scales ranging from small to medium to large. Finally, combine them to form the texture feature matrix F2∈R. M×N×9 ; S3: Use the OTSU algorithm to binarize the mean image in the spectral features to obtain image B. Then, perform a distance transform on image B to calculate the Euclidean distance. The formula is as follows: Where D(i,j) represents the distance change result matrix; (i,j) represents the current foreground pixel coordinates; p, q represent the background pixel coordinates; min represents the minimum value function, which finds the nearest background pixel; √ represents the square root operation, which gives the target the true distance; the distance value obtained by the formula reflects the relative position of the pixel within the habitat patch; S4: Calculate the gradient magnitude using the Sobel operator. Within a 7×7 sliding window, count the proportion of pixels whose gradient magnitude exceeds a threshold, normalize them to the range [0, 1], and obtain the edge density map ED(i,j). Finally, combine the above calculation parameters to form the geometric feature matrix F3∈R. M×N×3 ; S5: Assume the above spectral feature matrix F1∈R M×N×9 Texture feature matrix F2∈R M×N×9 and geometric characteristic matrix F3∈R M×N×3 Merge into F 综 ∈R M×N×21 Data set and output.

4. The method for extracting habitat areas based on remote sensing imagery according to claim 3, characterized in that, Using a random forest model as a classifier and combining it with a transfer learning strategy to construct a preliminary habitat classification map includes the following steps: S1: Based on F 综 ∈R M×N×21 The dataset was simplified using a multi-dimensional filtering strategy, including importance assessment, correlation analysis, and filtering strategies, resulting in the selection of the core feature dataset F. 核 ∈R M×N×K Where K=10 indicates that 10 core features are retained, and finally, the standardization process is performed and the output is completed. S2: Set the initial parameters of the pre-trained random forest model. This model includes 100 decision trees, a similar habitat category system, and adjustable model parameters. Fine-tune the model parameters using a transfer learning mechanism. The transfer learning formula is: P 概 (c|X)=α×P 预 (c|X)+(1-α)×P 微 (c|X), where P 概 (c|X) represents the final category prediction probability; c represents the habitat index; X represents the input feature vector; α represents the weight coefficients of the pre-trained random forest model; P 预 (c|X) represents the predicted probability of the pre-trained model; P 微 (c|X) represents the predicted probability of the model after fine-tuning; the formula is used to calculate the fine-tuned pre-trained random forest model.

5. The method for extracting habitat areas based on remote sensing imagery according to claim 4, characterized in that, The method of constructing a preliminary habitat classification map using a random forest model as a classifier and combining it with a transfer learning strategy also includes the following steps: S1: Based on the fine-tuned pre-trained random forest model, a lightweight model with better computational efficiency is constructed. The training process includes parameter setting → parallel training → setting prediction mechanism → setting probability fusion mechanism. After training, a lightweight random forest classifier model RF is obtained. S2: Using the lightweight random forest classifier model RF, category inference is performed on the entire image. The inference process includes feature subset extraction → block prediction → category decision. Finally, the inference results are combined to form a preliminary habitat classification map.

6. The method for extracting habitat areas based on remote sensing imagery according to claim 5, characterized in that, The method of constructing the habitat map using multi-scale morphological filtering and fast boundary optimization includes the following steps: S1: Based on the preliminary habitat classification map, a morphological filtering algorithm is used to divide the map into a dual-scale structure of [3×3, 7×7] windows. Opening operation is performed on the 3×3 window scale, and closing operation is performed on the 7×7 window scale. The results of the two scales are fused using the majority voting principle to finally obtain a smoothed classification map. S2: Boundary pixel identification is performed based on the smoothed classification map. Category changes are detected through a 3×3 window. When there are at least two different categories within the window, the center pixel is marked as a boundary pixel. S3: For each boundary pixel, examine the category distribution within its 5×5 neighborhood. The strategy is to count the frequency of each category within the neighborhood. If the frequency of a neighboring category is greater than the frequency of the current category plus a threshold T, then the current pixel is reclassified as that neighboring category. The threshold T = 5. S4: Apply 3×3 mean filtering to the boundary region, which is performed only when the boundary pixels are limited, and finally output the classification map after boundary optimization.

7. The method for extracting habitat areas based on remote sensing imagery according to claim 6, characterized in that, The method of constructing the habitat map using multi-scale morphological filtering and fast boundary optimization also includes the following steps: S1: Based on the optimized classification map, the flooding fill algorithm is used to calculate the number of pixels in each connected region. When the region area is less than 10 pixels, it is replaced with the surrounding main category to remove small patches. S2: Set ecological rules to correct the category. The rule is to change the category of small water bodies completely surrounded by vegetation to wetlands and single pixels to replace the main category of the 3×3 neighborhood and check the consistency of river width. S3: Using a 3×3 morphological closing operation, fill the remaining small holes to obtain the final habitat extraction result image.

Images