A method for assessing regional accessibility potential considering multiple stress factors

Abstract

The present invention provides a method for evaluating regional accessibility potential taking into account multi-factor coercion, which includes: 1) extracting terrain elements and inverting terrain complexity by combining optical and terrain data; 2) collecting existing open source data to obtain national road network information and inverting surface cover type information; 3) evaluating areas with extreme temperature and rainfall based on historical time series climate data, and inverting medium and high resolution soil moisture information based on surface cover types and open source kilometer-level soil moisture products using deep learning and multi-source remote sensing data; 4) inverting traffic safety factors; 5) conducting a comprehensive analysis of each factor and evaluating accessibility potential based on the analytic hierarchy process. The present invention effectively solves the problems of insufficient consideration of geological and environmental influencing factors in the process of regional accessibility potential evaluation, low recognition accuracy of relevant remote sensing elements at a wide scale, insufficient efficiency, and poor generalization ability of intelligent methods.

Description

Technical Field

[0001] The invention relates to the technical field of multi-source remote sensing data, and in particular to a method for evaluating regional accessibility potential taking into account multi-factor stress. Background Art

[0002] With the rapid development of my country's economy, road network construction has achieved remarkable results, and vehicle navigation technology has become increasingly mature. However, although the road network coverage in cities and some rural areas is relatively complete, looking at the world, most areas of the land surface are still roadless. According to statistics, of the 132 million square kilometers of land in the world, about 105 million square kilometers have no road coverage, accounting for 80% of the earth's land surface. This situation is particularly prominent in field application scenarios such as geological exploration, emergency mapping, disaster relief, off-road self-driving and even battlefield environments. In these complex roadless environments, the lack of road information will not only lead to a decrease in traffic efficiency, but also increase traffic risks, seriously affecting the execution efficiency and safety risk management of related work.

[0003] The emergence of new technologies such as artificial intelligence and cloud computing has provided new technical approaches for the efficient processing of multi-source remote sensing data and the intelligent assessment of the potential accessibility of large-scale areas. The application of these technologies can not only significantly improve the computational efficiency of data processing, but also optimize the information extraction process of remote sensing data through intelligent algorithms, and improve the accuracy and efficiency of remote sensing target recognition in complex environments. In addition, the massive remote sensing data and computing resources provided by cloud computing, combined with the powerful data processing capabilities of artificial intelligence, make real-time or near-real-time regional accessibility assessment possible, providing strong technical support for access decision-making in field environments.

[0004] In recent years, with the continuous advancement of satellite technology, the types of remote sensing data obtained have become increasingly rich, providing important technical support for the assessment of regional accessibility potential. However, when considering regional accessibility, current research and development methods often only focus on terrain and vehicle performance parameters, while the consideration of key factors such as surface type, geological conditions, and environmental safety is relatively insufficient. In addition, the processing process of massive remote sensing data is complex and the processing cycle is long, which makes it difficult to meet the application requirements of large-scale and high-timeliness. In terms of remote sensing information extraction, it currently mainly relies on analysis methods such as visual interpretation and remote sensing index, which has obvious limitations in the accuracy and efficiency of remote sensing target recognition in complex environments. In summary, although my country has made great progress in road network construction and vehicle navigation technology, it still faces many problems in conducting regional accessibility potential assessments that take into account multi-factor coercion in roadless areas.

[0005] To sum up, it is necessary to make further innovations on existing technologies. Summary of the invention

[0006] In response to the technical problems existing in the above-mentioned background technology, the present invention proposes a regional accessibility potential assessment method taking into account multi-factor coercion, which effectively solves the problems of insufficient consideration of geological and environmental influencing factors in the process of regional accessibility potential assessment, low recognition accuracy of relevant remote sensing elements at a wide scale, insufficient efficiency, and poor generalization ability of intelligent methods.

[0007] In order to solve the above technical problems, the present invention provides a method for evaluating regional accessibility potential taking into account multi-factor stress, which specifically includes the following steps:

[0008] 1) Extract terrain elements and calculate terrain complexity by combining optical and terrain data;

[0009] 2) Collect open source data of existing road network information to obtain national road network information, combine remote sensing big data, cloud computing platform and deep learning model to invert land cover type information, and pay attention to wetlands and land freeze-thaw areas affected by seasonal changes;

[0010] 3) Based on historical time-series climate data, evaluate current climate data in areas with extreme temperatures and rainfall, and use deep learning and multi-source remote sensing data to invert medium- and high-resolution soil moisture information based on land cover types and open-source kilometer-level soil moisture products;

[0011] 4) Inversion of traffic safety factors, that is, statistical analysis is performed based on the geological disaster points that have occurred to determine the high-incidence areas of geological disasters, and then deep learning and InSAR data are used to extract abnormal surface deformation areas associated with geological disaster risks;

[0012] 5) Comprehensive analysis of various factors is conducted and potential assessment of trafficability is conducted based on the analytic hierarchy process.

[0013] The method for assessing regional accessibility potential taking into account multi-factor stress, wherein the terrain complexity in step 1) is calculated by weighting terrain roughness, surface texture and terrain state index, and the calculation process involves DEM data and optical remote sensing data. The specific inversion calculation process is as follows:

[0014] TCF=w1×R+w2×T+w3×SI (1);

[0015] Among them, w1, w2 and w3 in the above formula (1) are weight coefficients. It is usually recommended that w1 and w2 take values ​​of 0.4 and w3 take a value of 0.2; R is the surface roughness; T is the surface texture coefficient; SI is the surface shape index;

[0016] The surface roughness R in the above formula (1) is used to characterize the undulation of the terrain. The specific calculation formula is as follows:

[0017]

[0018] In the above formula (2), z i is the elevation value of each sampling point, is the average elevation of all sampling points, and N is the total number of sampling points;

[0019] The surface texture coefficient T in the above formula (1) is used to characterize the complexity of the surface texture in the optical image. The specific calculation formula is as follows:

[0020]

[0021] Through the calculation of the gray-level co-occurrence matrix, P(i,j) represents the frequency of occurrence of gray levels i and j under the conditions specified by distance and direction, and G is the number of gray levels;

[0022] The surface shape index SI in the above formula (1) is used to calculate the curvature change of the terrain in different axes. The specific formula is as follows:

[0023]

[0024] In the above formula (4), k1 is the curvature of the terrain in the horizontal direction, and k2 is the curvature of the terrain in the y longitudinal direction.

[0025] The method for assessing regional accessibility potential taking into account multi-factor stress, wherein: the calculation result of the terrain complexity needs to be normalized so that its value is within the range of 0-1.

[0026] The method for assessing regional accessibility potential taking into account multi-factor coercion, wherein: the terrain elements in step 1) include slope and terrain undulation; step 1) is to calculate the slope through DEM data, and distinguish whether it is a passable area according to the slope value, wherein 15° is agreed to be the critical value, that is, the area with a slope greater than 15° is an impassable area, and the area with a slope less than or equal to 15° is a passable area; the weight type of the slope is an interval type, when the slope is greater than 15, it is assigned a value of 0, and when the slope is less than 15, it needs to be divided by 15 for normalization so that its value range is 0-1.

[0027] The method for assessing regional accessibility potential taking into account multi-factor coercion, wherein: the data source for collecting open source data of existing road network information in step 2) mainly comes from OpenStreetMap or open source road network information published by official agencies; the weight type of the open source road network information is interval type, in which wide roads such as national highways are assigned a value of 1.0, main roads of moderate width are assigned a value of 0.8, secondary roads are assigned a value of 0.5, residential roads are assigned a value of 0.4, and narrow roads such as non-urban dirt roads are assigned a value of 0.2.

[0028] The method for assessing regional accessibility potential taking into account multi-factor stress, wherein: the data source for inverting the surface cover type information in step 2) involves multi-source remote sensing data. In order to meet the needs of real-time classification of multi-type remote sensing in a wide area, it is necessary to construct a multi-source multi-classification remote sensing classification model, and to perform lightweight transformation on the network based on the knowledge distillation method. The specific steps are as follows:

[0029] 2.1) Construction of multi-source remote sensing sample dataset for land cover classification

[0030] The multi-source remote sensing land cover classification sample dataset consists of multi-source remote sensing data and true value data. Across the country, typical areas were selected to collect multi-source remote sensing data. After the collected data were downloaded, they were manually interpreted to obtain the corresponding true value labels of land cover classification, which were then merged with the normalized multi-source remote sensing data to form original samples with 16 bands. The original samples were sequentially cropped horizontally and vertically with a sample frame size of 256×256 to obtain the final multi-source remote sensing land cover classification sample dataset.

[0031] 2.2) Land cover classification network training and distillation

[0032] The improved UNet network, namely OUNet network, is used as the surface cover classification network; the OUNet network makes three improvements based on Unet: first, the downsampling operation is reduced to retain more shallow spatial information; second, the number of basic convolution units in the encoder is increased, and the depth-separable convolution is used to replace the ordinary convolution; third, the attention mechanism is introduced in the decoder part;

[0033] Through the knowledge distillation technology, a small "student" model is trained to imitate the large "teacher" model; the student network DSNet chooses to use three groups of CBR convolution units, each of which consists of a convolution layer, a batch layer and a ReLU activation layer. No downsampling operation is used, and 1×1 convolution and Softmax activation are used at the output to obtain a binary mask result of the same size as the input;

[0034] Let U be the original U-Net model, i.e. the teacher model, S be the student model, and the input data be X; the output of the teacher model can be expressed as U(X), and the output of the student model can be expressed as S(X); and the loss function in the knowledge distillation process consists of the cross entropy loss function and the distillation loss function, which can be defined as:

[0035] L total =L CE +L distill (5);

[0036] Among them, L in the above formula (5) totalis the cross entropy loss function, and its calculation formula is:

[0037]

[0038] In the above formula (6), C is the summary of the category, y o,c is the one-hot encoding of the true label, p c Predict the probability that the sample belongs to category C for the model;

[0039] In the above formula (5), L distill is the distillation loss function, and its calculation formula is:

[0040]

[0041] In the above formula (7), T(X) c and S(X) c Respectively represent the probability that the teacher model and the student model predict that the input X belongs to category C;

[0042] The surface cover classification sample dataset containing multi-source remote sensing data and true value data constructed in the above step 2.1) is respectively fed into the teacher network and the student network for training to generate U(X) and S(X) respectively. The cross entropy loss function L total During the training process, it will gradually converge, making S(X) gradually approach U(X) and the true value;

[0043] 2.3) Real-time extraction of land cover types

[0044] The student network is reconstructed using JAVA language in the GEE platform, and the parameters of the trained student network are imported into the reconstructed network for initialization; then, by inputting the normalized multi-source remote sensing data with 15 bands, the labels of the corresponding land cover types can be obtained in real time;

[0045] 2.4) Classification of accessible / inaccessible surfaces

[0046] Based on the surface cover classification results, it is determined whether the area is drivable or not.

[0047] The method for assessing regional accessibility potential taking into account multi-factor stress, wherein: in the step 2.4), water bodies and ice and snow building areas in the surface cover classification results are classified as inaccessible areas, and the classification results are all assigned a value of 0; bare land, cultivated land, grassland, and woodland in the surface cover classification results are classified as accessible areas, and are assigned values ​​of 0.8, 0.5, 0.3, and 0.1, respectively; and in the identification process, it is necessary to pay attention to the fact that some areas have seasonal wetlands or land freeze-thaw conditions due to the influence of rainfall and temperature, and the weight coefficient can be reduced for such areas.

[0048] The method for assessing regional accessibility potential taking into account multi-factor stress, wherein: the time series meteorological data in step 3) is obtained from meteorological stations in the study area or through relevant open source remote sensing products; in terms of temperature, <-40° is assigned a value of 0.1, -40°-25° is assigned a value of 0.4, -25°-32° is assigned a value of 0.8, 32°-40° is assigned a value of 0.4, and >40° is assigned a value of 0.1; in terms of rainfall, >100 mm is assigned a value of 0.1, 50-100 mm is assigned a value of 0.2, 20-50 mm is assigned a value of 0.4, 5-20 mm is assigned a value of 0.6, and 0-5 mm is assigned a value of 0.8;

[0049] The data sources for inverting high-resolution soil moisture information based on surface cover types and open source kilometer-level soil moisture products in step 3) include Landsat8 visible and thermal infrared data, Sentinel-1 radar data and terrain data; Landsat8 data is used for calculating water body index and extracting emissivity information, and Sentinel-1 and terrain data are used to obtain backscatter coefficient data after terrain correction.

[0050] The method for assessing regional accessibility potential taking into account multi-factor stress, wherein the specific process of inverting high-resolution soil moisture information based on the surface cover type and the open source kilometer-level soil moisture product in step 3) is as follows:

[0051] 3.1) Data preparation and preprocessing:

[0052] Collect SMAP soil moisture content products, Landsat8 reflectance and thermal infrared images, Sentinel-1 backscatter coefficient data and terrain data; use Landsat data to calculate the water body index NDWI, the calculation formula is as follows:

[0053]

[0054] In the above formula (8), GREEN represents the green light band, and NIR represents the near infrared band.

[0055] The spatial resolution of the water body index NDWI, Landsat8 thermal infrared data, Sentinel-1 backscatter coefficient data, DEM data and the surface cover type data obtained in step 2) above was unified to 30 meters, and superimposed into a multi-source soil moisture sample dataset with 6 bands; based on the surface cover type data in the dataset, J areas with an area greater than 1km were selected. 2 The bare land is taken, and the coordinates of its center point are taken, and the set of NDWI index, thermal infrared emissivity data, backscattering coefficient data and terrain data corresponding to the coordinates are extracted as Xt, and the corresponding SMAP soil moisture value is taken as Yt; the sample data set of soil moisture inversion composed of J groups (Xt, Yt) of data;

[0056] 3.2) Model training and soil moisture inversion

[0057] The feature vector Xt and target variable Yt are trained using the random forest algorithm. The performance of the model is evaluated. If the requirements are met, the model can be applied to new data sets to achieve the inversion of medium- and high-resolution soil moisture content. The inverted results will be normalized to the range of 0-1.

[0058] The method for assessing regional accessibility potential taking into account multi-factor coercion comprises the following steps: the data source for extracting abnormal surface deformation areas associated with geological hazards using deep learning and InSAR data in step 4) comes from Sentinel-1 radar data, and surface deformation phase data is obtained through preprocessing. In combination with fine-tuning technology based on a large visual model, intelligent extraction of abnormal surface deformation areas associated with wide-area geological hazards is achieved.

[0059] The method for assessing regional traversability potential taking into account multi-factor stress, wherein the specific process of inverting the traversability safety factor in step 4) is as follows:

[0060] 4.1) Data preprocessing

[0061] The Sentinel-1 radar data were processed by the interferogram stacking time-series InSAR method to obtain the surface deformation phase data associated with geological hazards, namely the annual average surface deformation phase map. The linear phase rate p in the annual average surface deformation phase map is calculated as follows:

[0062]

[0063] In the above formula (9), Δt j is the time interval of the jth group of differential interference patterns; is the unwrapped differential interference phase of the jth group;

[0064] Then, the InSAR deformation data of different causes are identified and made into true value labels for the identification of InSAR surface deformation anomaly areas, which are then merged with the annual average surface deformation phase map into a sample data set for intelligent identification of surface deformation phase, providing basic data for subsequent hidden danger identification applications.

[0065] 4.2) Intelligent identification of abnormal surface deformation areas

[0066] In the Pytorch environment, an intelligent recognition model InSARNet for abnormal surface deformation areas is built. The intelligent recognition model InSARNet for abnormal surface deformation areas consists of three parts: an image encoder, a keyword decoder and an image decoder. Then, based on the surface deformation phase intelligent recognition sample data set obtained in the above step 4.1), the intelligent recognition model InSARNet for abnormal surface deformation areas is trained. During the training process, the parameters of the image encoder are frozen, and only the embeddable fine-tuning module and the image decoder are trained. After the intelligent recognition model InSARNet for abnormal surface deformation areas is trained, the abnormal surface deformation areas associated with geological disaster risks in unknown areas can be identified by inputting InSAR deformation phase data. The identification result is expressed as a probability value of 0-1. The higher the probability value, the greater the possibility that it is an abnormal surface deformation area.

[0067] The method for assessing regional accessibility potential taking into account multi-factor stress, wherein the main steps of processing Sentinel-1 radar data by interferogram stacking time-series InSAR method in step 4.1) include:

[0068] 4.1.1) Splice the three sub-strips of Sentinel-1 TOPS data into a complete SLC image;

[0069] 4.1.2) Interference image pair selection;

[0070] 4.1.3) Differential interferometry calculation;

[0071] 4.1.4) Phase unwrapping;

[0072] 4.1.5) Phase superposition processing;

[0073] 4.1.6) Trend stripe removal;

[0074] 4.1.7) Geocoding, finally obtaining the annual average surface deformation phase map;

[0075] The image encoder is composed of a SAM pre-trained ViT model and an embeddable fine-tuning module; the VIT model adopts the ViT-H / 16 version, which consists of four isometric global attention modules; the embeddable fine-tuning module consists of two sets of multi-layer perceptrons and two layers of ReLU activation layers; the input of the image encoder is an image of size 1024×1024×3; the output of the image encoder is a 16×16 image feature embedding processed by the patch embedding operation and the SAM pre-trained ViT module;

[0076] The image decoder consists of a ViT module and two layers of multi-layer perceptrons, which are used to restore the image feature embedding information into label data consistent with the input image size; the image encoder of InSARNet is initialized with parameters from a large visual model trained on massive natural images.

[0077] The method for assessing regional traversability potential taking into account multi-factor stress, wherein the step 5) involves using the analytic hierarchy process to assess regional traversability potential, and the specific steps are as follows:

[0078] 5.1) Determine the hierarchical relationship

[0079] Construct a hierarchical relationship based on "target layer-factor layer-element layer";

[0080] 5.2) Constructing the comparison matrix

[0081] Construct a matrix M, the elements in it are marked as a ij ; Compare each element pairwise and determine the appropriate scale, then fill the determined scale into the corresponding position of the matrix M. The final matrix M is the pairwise comparison matrix; where a ij It represents the scale of i compared to j; a ij The corresponding position is the i-th row and j-th column in the matrix M; if i=j, it means that i and j are equally important, that is, a ij =1, which means that the scales of the main diagonal elements of the matrix M are all 1; a ij >0 and a ij ×a ji =1;

[0082] 5.3) Calculation of weight coefficient

[0083] The weights are calculated using the arithmetic mean method, geometric mean method and eigenvalue method respectively, and the arithmetic mean of the calculation results is calculated;

[0084] For the arithmetic mean method to calculate the weight, the calculation formula is as follows:

[0085] Assume that the pairwise comparison matrix M is:

[0086]

[0087] The weight vector w is obtained by arithmetic mean method i for:

[0088]

[0089] For the arithmetic mean method to calculate the weight, the calculation formula is as follows:

[0090] Assume that the pairwise comparison matrix M is:

[0091]

[0092] The weight vector w is obtained by geometric mean method i for:

[0093]

[0094] For the eigenvalue method to find the weight, the calculation formula is as follows:

[0095] First assume that the maximum eigenvalue of the pairwise comparison matrix M is λ max and its corresponding eigenvector is ξ, then:

[0096]

[0097] The weight vector w is obtained by eigenvalue method i for:

[0098]

[0099] 5.4) Consistency Check

[0100] Calculate the largest eigenvalue λ of the eigenvector of the comparison matrix max :

[0101] CW = C*w (16);

[0102]

[0103] In the above formulas (16)-(17), C is the comparison matrix, w is the indicator weight, and n represents the comparison matrix order;

[0104] Calculate the consistency index CI of the comparison matrix:

[0105]

[0106] When CI is close to 0, it can be preliminarily judged that the constructed pairwise comparison matrix has satisfactory consistency; then calculate the consistency ratio CR:

[0107]

[0108] In the above formula (19), RI is the average consistency index; usually when CR < 0.1, the consistency of the pairwise comparison matrix can be considered acceptable;

[0109] 5.5) Assessment of traversability potential

[0110] The factor layer and the element layer are calculated according to the above steps 5.2) to 5.4) respectively to obtain the weight coefficients of the factor layer and the element layer, and the initial values ​​are weighted and summed to obtain the traversable potential assessment score; according to the traversable potential assessment score, the k-means clustering algorithm is used to classify the traversable potential assessment score.

[0111] The method for assessing regional accessibility potential taking into account multi-factor coercion, wherein: the target layer in step 5.1) is the regional accessibility potential score; the factor layer includes topography, surface cover, climate environment, and traffic safety; the element layer includes slope, terrain complexity, roads, surface cover type, temperature, rainfall, soil moisture content, geological disaster risk information, and existing geological disaster points.

[0112] The method for evaluating regional accessibility potential taking into account multi-factor stress, wherein the process of using the k-means clustering algorithm to classify the evaluation results in step 5.5) is as follows: predetermine the number of classification categories K, randomly select K objects {a1, a2, ..., a K} as given n data {x1,x2,…,x n}; calculate the distance between each data point and each cluster center, assign it to each cluster center according to the principle of closest distance, and recalculate each category center, so that the sum of the squares of the distances from each data point to the nearest cluster center is minimized;

[0113] And the sum of the squares of the distances of each data point to the nearest cluster center is the objective function W n :

[0114]

[0115] Among them, W in the above formula (20) n is the distance objective function of clustering, a j is the initial cluster center, K is the number of classification categories, X i is the data to be clustered, and n is the number of data.

[0116] By adopting the above technical solution, the present invention has the following beneficial effects:

[0117] The present invention aims to invert the terrain factors, surface state factors, climate and environmental factors, and traffic safety factors that affect regional accessibility based on multi-phase optical, radar, terrain and other multi-source remote sensing data, taking into account geological and meteorological environmental factors, quantify the spatiotemporal variation characteristics of each factor and its impact on regional accessibility, construct a potential accessibility assessment model under a complex environment, and combine spatial analysis technology to realize the mapping of regional scale accessibility potential areas, provide basic data support for route planning in complex unstructured environments, improve traffic efficiency and safety, and effectively serve the needs of geological exploration, emergency mapping, disaster relief, and war zone navigation.

[0118] The present invention realizes for the first time an accurate assessment of the regional accessibility potential through a comprehensive inversion of multiple factors such as terrain factors, surface cover factors, climate environment and traffic safety.

[0119] The present invention proposes the concept and calculation method of terrain complexity; uses convolutional neural networks and knowledge organization technology to realize the inversion of surface cover factors based on the GEE platform; realizes the intelligent extraction of surface state factors; combines multimodal remote sensing data and machine learning methods to realize the inversion of medium and high resolution soil moisture content; realizes the identification of abnormal surface deformation areas associated with geological disaster risks through visual large model fine-tuning technology; for the first time, through the comprehensive inversion of terrain factors, surface cover factors, climate environment and traffic safety and other factors, an accurate assessment of the regional accessibility potential is realized; effectively solves the problems of insufficient consideration of geological and environmental influencing factors in the process of regional accessibility potential assessment, low recognition accuracy of relevant remote sensing elements at a wide scale, insufficient efficiency, and poor generalization ability of intelligent methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0120] In order to more clearly illustrate the specific implementation methods of the present invention or the technical solutions in the prior art, the drawings required for use in the specific implementation methods or the description of the prior art will be briefly introduced below. Obviously, the drawings in the following description are some implementation methods of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0121] Figure 1 A schematic diagram of the OUNet teacher network structure involved in the regional accessibility potential assessment method taking into account multi-factor stress of the present invention;

[0122] Figure 2 A schematic diagram of the DSNet student network structure involved in the regional accessibility potential assessment method taking into account multi-factor stress of the present invention;

[0123] Figure 3A schematic diagram of the InSARNet network structure involved in the regional accessibility potential assessment method taking into account multi-factor stress of the present invention;

[0124] Figure 4 It is a hierarchical relationship diagram of traversability potential assessment involved in the regional traversability potential assessment method taking into account multi-factor stress of the present invention. DETAILED DESCRIPTION

[0125] The technical solution of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0126] The present invention is further explained below in conjunction with specific implementation modes.

[0127] This embodiment provides a method for evaluating regional accessibility potential taking into account multi-factor stress, which specifically includes the following steps:

[0128] S100, extracting terrain elements, mainly including slope and terrain undulation, and calculating terrain complexity by combining optical and terrain data;

[0129] S200. Collect open source data of existing road network information to obtain national road network information; combine remote sensing big data, cloud computing platform and deep learning model to invert the surface cover type information of bare land, grassland, forest land, cultivated land, water body, building and ice and snow, and pay attention to wetlands and land freeze-thaw areas affected by seasonal changes;

[0130] S300, based on historical time-series climate data, evaluates current climate data for areas with extreme temperatures and rainfall, and inverts medium- and high-resolution soil moisture information using deep learning and multi-source remote sensing data based on land cover types and open-source kilometer-level soil moisture products;

[0131] S400, inversion of traffic safety factors, on the one hand, statistical analysis is performed based on the geological disaster points that have occurred to determine the high-incidence areas of geological disasters, and on the other hand, deep learning and InSAR data are used to extract abnormal surface deformation areas associated with geological disaster risks;

[0132] S500, comprehensive analysis of various factors, and evaluation of traversability potential based on analytic hierarchy process.

[0133] The above step S100 involves the calculation of terrain complexity (TCF), which is obtained by weighting terrain roughness, surface texture and terrain state index. The calculation process involves DEM data and Sentinel-2 optical remote sensing data. The specific calculation process is as follows:

[0134] TCF=w1×R+w2×T+w3×SI (1);

[0135] Among them, w1, w2 and w3 in the above formula (1) are weight coefficients. It is usually recommended that w1 and w2 take values ​​of 0.4 and w3 take a value of 0.2; R is the surface roughness; T is the surface texture coefficient; SI is the surface shape index;

[0136] The surface roughness R in the above formula (1) is used to characterize the undulation of the terrain. The specific calculation formula is as follows:

[0137]

[0138] In the above formula (2), z i is the elevation value of each sampling point, is the average elevation of all sampling points, and N is the total number of sampling points;

[0139] The surface texture coefficient T in the above formula (1) is used to characterize the complexity of the surface texture in the optical image. The specific calculation formula is as follows:

[0140]

[0141] Through the calculation of the gray level co-occurrence matrix (GLCM), P(i,j) represents the frequency of gray levels i and j under the conditions of specified distance and direction, and G is the number of gray levels;

[0142] The surface shape index SI in the above formula (1) is used to calculate the curvature change of the terrain in different axes. The specific formula is as follows:

[0143]

[0144] In the above formula (4), k1 is the curvature of the terrain in the horizontal direction, and k2 is the curvature of the terrain in the y longitudinal direction.

[0145] At the same time, the calculation results of the terrain complexity (TCF) mentioned above need to be normalized so that their values ​​are within the range of 0-1.

[0146] In the above step S100, the slope is calculated by DEM data, and whether it is a passable area is determined according to the slope value, wherein 15° is agreed to be the critical value, that is, the area with a slope greater than 15° is an impassable area, and the area with a slope less than or equal to 15° is a passable area; the weight type of the slope is interval type, when the slope is greater than 15, it is assigned a value of 0, and when the slope is less than 15, it needs to be divided by 15 for normalization so that its value range is 0-1.

[0147] The above step S200 involves the collection of open source data of existing road network information, and the data source mainly comes from the open source road network information published by OpenStreetMap or official institutions. The weight type of the open source road network (Roadmap) is interval type, in which the value of wide roads such as national roads is 1.0, the value of main roads with moderate width is 0.8, the value of secondary roads is 0.5, the value of residential roads is 0.4, and the value of narrow roads such as non-urban dirt roads is 0.2.

[0148] The above step S200 involves the inversion of surface cover information. The data sources include multi-source remote sensing data such as Sentinel-2 optical, Sentinel-1 radar and DEM. In order to meet the needs of wide-area multi-type remote sensing real-time classification, it is necessary to build a multi-source multi-classification remote sensing classification model and perform lightweight transformation on the network based on knowledge distillation technology. The specific steps are as follows:

[0149] S210, Construction of multi-source remote sensing sample dataset for land cover classification

[0150] The multi-source remote sensing land cover classification sample dataset consists of two parts: multi-source remote sensing data and true value data. The multi-source remote sensing data is a combination of Sentinel-2 reflectance data (13 bands), Sentinel-1 intensity data (1 band), and terrain data (1 band), containing a total of 15 bands. The true value data is obtained by manual interpretation and contains a total of 8 land feature classes. The land class names and values ​​are as follows: bare land 1, grassland 2, forest land 3, cultivated land 4, water body 5, building 6, ice and snow 7, and others 8.

[0151] Typical areas across the country were selected to collect multi-source remote sensing data, and this process was carried out on a cloud platform. After the collected data were downloaded, they were manually interpreted to obtain the corresponding true value labels of the surface cover classification, and merged with the normalized multi-source remote sensing data to form original samples with 16 bands. With a sample frame of 256×256 size, the original samples were sequentially cropped horizontally and vertically to obtain the final multi-source remote sensing surface cover classification sample dataset.

[0152] S220, land cover classification network training and distillation

[0153] The improved UNet network (OUNet) is used for the land cover classification network; Figure 1 As shown in the figure, the OUNet network has made three improvements on the basis of Unet: first, it reduces the downsampling operation and retains more shallow spatial information to improve the edge accuracy of remote sensing object classification results; second, it increases the number of basic convolution units in the encoder and uses depthwise separable convolution instead of ordinary convolution; third, it introduces the attention mechanism in the decoder part to enable better fusion of features of different scales.

[0154] Knowledge distillation is a model compression technique that reduces the size of the model by training a small "student" model to mimic a large "teacher" model. Figure 2 The design principle of the proposed method is to maintain sufficient learning ability and have a small number of parameters. Therefore, three groups of CBR convolutional units are selected. Each group of CBR convolutional units consists of a convolutional layer, a batch processing layer, and a ReLU activation layer. No downsampling operation is used. 1×1 convolution and Softmax activation are used at the output to obtain a binary mask result that is the same size as the input.

[0155] Let U be the original U-Net model, i.e. the teacher model, S be the student model, and the input data be X; the output of the teacher model can be expressed as U(X), and the output of the student model can be expressed as S(X); the loss function in the knowledge distillation process consists of two parts, one is the cross entropy loss, and the other is the distillation loss, which can be defined as:

[0156] L total =L CE +L distill (5);

[0157] Among them, L in the above formula (5) total is the cross entropy loss function, and its calculation formula is:

[0158]

[0159] In the above formula (6), C is the summary of the category, y o,c is the one-hot encoding of the true label, p c Predict the probability that the sample belongs to category C for the model;

[0160] In the above formula (5), L distill is the distillation loss function, and its calculation formula is:

[0161]

[0162] In the above formula (7), T(X) c and S(X) c Respectively represent the probability that the teacher model and the student model predict that the input X belongs to category C;

[0163] The surface cover classification sample data set containing multi-source remote sensing data and true value data constructed in the above step S210 is respectively fed into the teacher model and the student model for training to generate U(X) and S(X) respectively. The loss function L total During the training process, it will gradually converge, making S(X) gradually approach U(X) and the true value data.

[0164] S230, real-time extraction of land cover types

[0165] Since the GEE platform does not provide Tensorflow or Pytorch programming interfaces, it is necessary to reconstruct the student network using JAVA in GEE and import the parameters of the trained student network into the reconstructed network for initialization; by inputting normalized multi-source remote sensing data with 15 bands composed of Sentinel-2, Sentinel-1 and terrain, the corresponding surface cover type labels can be obtained in real time.

[0166] S240. Classification of passable and impassable surfaces

[0167] According to the surface cover classification results, it is distinguished whether it is a passable area. Water bodies and ice and snow building areas are classified as impassable areas, and the classification results are all assigned a value of 0; bare land, cultivated land, grassland, and woodland are passable areas, and are assigned values ​​of 0.8, 0.5, 0.3, and 0.1 respectively; in the division process, it is necessary to pay attention to the fact that some areas have seasonal wetlands or land freeze-thaw due to the influence of rainfall and temperature. For such areas, the weight coefficient can be appropriately reduced.

[0168] The above step S300 involves using time series climate data to evaluate extreme temperature and rainfall areas; wherein, the time series meteorological data comes from the meteorological stations in the study area, and can also be obtained through relevant open source remote sensing products; in terms of temperature, <-40° is assigned a value of 0.1, -40°—25° is assigned a value of 0.4, -25°-32° is assigned a value of 0.8, 32°-40° is assigned a value of 0.4, and >40° is assigned a value of 0.1; in terms of rainfall, >100 mm is assigned a value of 0.1, 50-100 mm is assigned a value of 0.2, 20-50 mm is assigned a value of 0.4, 5-20 mm is assigned a value of 0.6, and 0-5 mm is assigned a value of 0.8.

[0169] The data sources for inverting high-resolution soil moisture information based on the surface cover type and the open-source kilometer-level soil moisture product in step S300 involve Landsat8 data, Sentinel-1 data and terrain data; among them, Landsat8 data mainly uses its optical and thermal infrared bands for calculating water index and extracting emissivity information, Sentinel-1 mainly uses backscatter coefficient data to better reflect water information, and terrain data is used for terrain correction; the specific steps are as follows:

[0170] S310, data preparation and preprocessing:

[0171] Collect SMAP soil water content products, Landsat8 reflectance and thermal infrared data from Landsat8 data, backscatter coefficient data and terrain data from Sentinel-1 data; use Landsat8 data to calculate the water body index NDWI, the calculation formula is as follows:

[0172]

[0173] In the above formula (8), GREEN represents the green light band, and NIR represents the near infrared band.

[0174] The spatial resolution of the data obtained after inverting the water body index NDWI, the thermal infrared data in the Landsat8 data, the backscatter coefficient data in the Sentinel-1 data, the DEM data and the surface cover type information in the above step S200 is unified to 30 meters, and superimposed into a 6-band multi-source soil moisture sample data set; based on the surface cover type data in the data set obtained after inverting the surface cover type information in the above step S200, J areas with an area greater than 1km are selected 2 For bare land, the NDWI index corresponding to the coordinates of its center point, the thermal infrared data in the Landsat8 data, the backscatter coefficient data in the Sentinel-1 data, and the terrain data are taken as the feature vector Xt, and the corresponding SMAP soil moisture content value is taken as the target variable Yt, and a sample data set of soil moisture inversion consisting of J groups of (Xt, Yt) data is obtained.

[0175] S320, model training and soil moisture inversion

[0176] The feature vector Xt and target variable Yt are trained using the random forest algorithm; the performance of the random forest model is evaluated, and if it meets the requirements (usually the accuracy reaches more than 90%), the model can be applied to new data sets to achieve the inversion of medium and high resolution soil moisture content; the inverted results will be normalized to the range of 0-1.

[0177] The above step S400 involves using deep learning and InSAR data to intelligently extract abnormal surface deformation areas associated with geological hazards; the data source comes from the radar data of Sentinel-1, and the surface deformation phase data is obtained by preprocessing the Sentinel-1 radar data, and combined with the fine-tuning technology based on the visual large model, the intelligent extraction of abnormal surface deformation areas associated with wide-area geological hazards is realized; specifically, the following steps are included:

[0178] S410, data preprocessing

[0179] The Sentinel-1 radar data is processed by the interferogram stacking time-series InSAR method to obtain the surface deformation phase data associated with geological hazards, that is, the annual average surface deformation phase map. The stacking InSAR technology controls the time and space baseline through multiple main images, obtains multiple unwrapped differential interferograms, and then accumulates the interference phase, effectively solving the temporal and spatial decorrelation problems and reducing the impact of atmospheric disturbances. It is one of the methods to improve the accuracy of deformation solution; the calculation formula of the linear phase rate p in the annual average surface deformation phase map is:

[0180]

[0181] In the above formula (9), Δt j is the time interval of the jth group of differential interference patterns; is the unwrapped differential interference phase of the jth group.

[0182] The main steps of processing Sentinel-1 data using the above-mentioned Stacking InSAR technology include:

[0183] S411, stitching the three sub-strips of Sentinel-1TOPS data into a complete SLC image;

[0184] S412, interference image pair selection;

[0185] S413, differential interference calculation;

[0186] S414, phase unwrapping;

[0187] S415, phase superposition processing;

[0188] S416, trend streak removal;

[0189] S417, geocoding, and finally obtain the annual average surface deformation phase map.

[0190] In order to accurately identify the InSAR surface deformation characteristics associated with landslide hazards, it is also necessary to identify the InSAR surface deformation phase data of different causes based on expert knowledge, and make true value labels for the identification of InSAR surface deformation anomaly areas, which can be merged with the surface deformation phase map data into a surface deformation phase intelligent identification sample data set to provide basic data for subsequent hazard identification applications.

[0191] S420, Intelligent identification of abnormal surface deformation areas

[0192] In the Pytorch environment, we build an intelligent recognition model for abnormal surface deformation areas, InSARNet (such as Figure 3 As shown in Figure 1, the model InSARNet mainly consists of three parts: image encoder, keyword decoder and image decoder. The image encoder is composed of the Vision Transformer (ViT) module pre-trained by SAM and the embeddable fine-tuning module. Specifically, VIT uses the ViT-H / 16 version, which consists of four equidistant global attention modules, which can better capture the long-distance dependencies between image slices. The embeddable fine-tuning module consists of two sets of multi-layer perceptrons and two layers of ReLU activation layers. The embeddable fine-tuning module allows the model InSARNet to fully utilize the extraction capabilities of the ViT model pre-trained by SAM while being more adaptable to downstream tasks, thereby achieving transfer learning at a lower computational cost. The input of the image encoder is an image of size 1024×1024×3, and the output of the image encoder is a 16×16 image feature embedding processed by the patch embedding operation and the ViT module pre-trained by SAM. The image decoder consists of a ViT module and two layers of multi-layer perceptrons, which are used to restore the image feature embedding information into label data of the same size as the input image. The image encoder of InSARNet is initialized with parameters from a large visual model trained on massive natural images, so it has strong general semantic segmentation performance. The embeddable fine-tuning module can improve the adaptability of the model to the task of surface deformation anomaly recognition.

[0193] The InSARNet model is trained based on the surface deformation phase intelligent recognition sample data set obtained in the above step S410. During the training process, the parameters of the encoder are frozen, and only the embeddable fine-tuning module and the image decoder are trained. After the InSARNet model is trained, the surface deformation anomaly area associated with geological disaster hazards in unknown areas can be identified by inputting the InSAR surface deformation phase data. Since the large visual model has good migration ability, combined with the embeddable fine-tuning module, the generalization of the surface deformation intelligent recognition ability under a wide-area scale is achieved. The recognition result is expressed as a probability value of 0-1, and the higher the probability value, the greater the possibility that it is a surface deformation anomaly area.

[0194] The above step S500 involves the use of the analytic hierarchy process to evaluate the regional traversability potential, and the specific steps are as follows:

[0195] S510: Determine hierarchical relationship

[0196] Construct a hierarchical relationship based on "target layer-factor layer-element layer", such as Figure 4 As shown. Among them, the target layer is the regional traversable potential score; the factor layer includes topography, surface cover, climate environment, and traffic safety; the element layer includes slope, terrain complexity, roads, surface cover type, temperature, rainfall, soil moisture content, geological disaster hidden danger information, and existing geological disaster points. On this basis, the present invention uses soil bearing capacity as an additional factor layer, allowing this factor to be included in the comprehensive evaluation of traversable potential when sufficient soil bearing capacity data is available. The initial numerical range of each factor selected by the present invention is shown in Table 1.

[0197] Table 1: Elements of roadability potential assessment

[0198]

[0199]

[0200]

[0201] S520, construct comparison matrix

[0202] Construct a matrix M, the elements in it are marked as a ij Combined with expert opinions and relevant information, compare each element pairwise, refer to Table 2 to determine the appropriate scale, and then fill the determined scale into the corresponding position of the matrix M. The final matrix M is the pairwise comparison matrix. ij The meaning of is the importance or scale of i compared with j. ij The corresponding position is the i-th row and j-th column in the matrix M. If i=j, it means that i and j are equally important, that is, a ij=1, which means that the scale of the main diagonal elements of the matrix M is 1. ij >0 and a ij ×a ji =1.

[0203] Table 2 Scale and meaning of the pairwise comparison matrix of the AHP

[0204]

[0205] S530, weight coefficient calculation

[0206] The weights are calculated by using the arithmetic mean method, geometric mean method and eigenvalue method respectively, and the arithmetic mean of the obtained weights is calculated to improve the accuracy of the calculation;

[0207] For the arithmetic mean method to calculate the weight, the specific calculation formula is as follows:

[0208] Assume that the pairwise comparison matrix M is:

[0209]

[0210] The weight vector w is obtained by arithmetic mean method i for:

[0211]

[0212] For the arithmetic mean method to calculate the weight, the specific calculation formula is as follows: Assume that the pairwise comparison matrix M is:

[0213]

[0214] The weight vector w is obtained by geometric mean method i for:

[0215]

[0216] For the eigenvalue method to calculate the weight, the specific calculation formula is as follows:

[0217] First assume that the maximum eigenvalue of the pairwise comparison matrix M is λ max and its corresponding eigenvector is ξ, then:

[0218]

[0219] The weight vector w is obtained by eigenvalue method i for:

[0220]

[0221] S540, consistency check

[0222] Calculate the maximum eigenvalue λ of the comparison matrix eigenvectormax :

[0223] CW = C*w (16);

[0224]

[0225] In the above formulas (16)-(17), C is the comparison matrix, w is the indicator weight, and n represents the comparison matrix order;

[0226] Calculate the comparison matrix consistency index CI:

[0227]

[0228] When CI is close to 0, it can be preliminarily judged that the constructed pairwise comparison matrix has satisfactory consistency; calculate the consistency ratio CR:

[0229]

[0230] In the above formula (19), RI is the average consistency index; usually when CR < 0.1, the consistency of the pairwise comparison matrix can be considered acceptable.

[0231] S550, Trafficability Potential Assessment

[0232] The above steps S520 to S540 are calculated for the factor layer and the element layer respectively to obtain the weight coefficients of the factor layer and the element layer, and the initial values ​​are weighted and summed to obtain the passable potential evaluation score; according to the passable potential evaluation score, the k-means clustering algorithm is used to classify the passable potential evaluation score to improve the hierarchical expression of the passable potential information. Among them, the calculation steps of the k-means clustering algorithm are:

[0233] Determine the number of categories K in advance, and randomly select K objects {a1, a2, …, a K} as given n data {x1,x2,…,x n}; calculate the distance between each data point and the initial cluster center, assign it to the initial cluster center according to the principle of closest distance, and recalculate the center of each category, so that the sum of the squares of the distances between each data point and the nearest cluster center is minimized;

[0234] And the sum of the squares of the distances of each data point to the nearest cluster center is the objective function W n :

[0235]

[0236] Among them, W in the above formula (20) n is the distance objective function of clustering, a j is the initial cluster center, K is the number of classification categories, Xi is the data to be clustered, and n is the number of data.

[0237] The present invention effectively solves the problems of insufficient consideration of geological and environmental influencing factors in the process of regional accessibility potential assessment, low recognition accuracy of relevant remote sensing elements at a wide scale, insufficient efficiency, and poor generalization ability of intelligent methods.

[0238] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for assessing regional accessibility potential taking into account multi-factor stress, comprising the following steps: 1) Extract terrain elements and calculate terrain complexity by combining optical and terrain data; 2) Collect open source data of existing road network information to obtain national road network information, combine remote sensing big data, cloud computing platform and deep learning model to invert land cover type information, and pay attention to wetlands and land freeze-thaw areas affected by seasonal changes; 3) Based on historical time-series climate data, evaluate current climate data in areas with extreme temperatures and rainfall, and use deep learning and multi-source remote sensing data to invert medium- and high-resolution soil moisture information based on land cover types and open-source kilometer-level soil moisture products; 4) Inversion of traffic safety factors, that is, statistical analysis is performed based on the geological disaster points that have occurred to determine the high-incidence areas of geological disasters, and then deep learning and InSAR data are used to extract abnormal surface deformation areas associated with geological disaster risks; 5) Comprehensive analysis of terrain factors, surface conditions factors, climate and environmental factors, and traffic safety factors, and evaluation of traffic potential based on the analytic hierarchy process; The data source for inverting the land cover type information in step 2) involves multi-source remote sensing data. In order to meet the needs of real-time classification of multi-type remote sensing in a wide area, it is necessary to build a multi-source multi-classification remote sensing classification model, and perform lightweight transformation on the network based on the knowledge distillation method. The specific steps are as follows: 2.1) Construction of multi-source remote sensing sample dataset for land cover classification The multi-source remote sensing land cover classification sample dataset consists of multi-source remote sensing data and true value data. Across the country, typical areas were selected to collect multi-source remote sensing data. After the collected data were downloaded, they were manually interpreted to obtain the corresponding true value labels of land cover classification, which were then merged with the normalized multi-source remote sensing data to form original samples with 16 bands. The original samples were sequentially cropped horizontally and vertically with a sample frame size of 256×256 to obtain the final multi-source remote sensing land cover classification sample dataset. 2.2) Land cover classification network training and distillation The improved UNet network, namely OUNet network, is used as the surface cover classification network; the OUNet network makes three improvements based on Unet: first, the downsampling operation is reduced to retain more shallow spatial information; second, the number of basic convolution units in the encoder is increased, and the depth-separable convolution is used to replace the ordinary convolution; third, the attention mechanism is introduced in the decoder part; Through the knowledge distillation technology, a small "student" model is trained to imitate the large "teacher" model; the student network DSNet chooses to use three groups of CBR convolution units, each of which consists of a convolution layer, a batch layer and a ReLU activation layer. No downsampling operation is used, and 1×1 convolution and Softmax activation are used at the output to obtain a binary mask result of the same size as the input; Let U be the original U-Net model, i.e. the teacher model, S be the student model, and the input data be X; the output of the teacher model can be expressed as U(X), and the output of the student model can be expressed as S(X); and the loss function in the knowledge distillation process consists of the cross entropy loss function and the distillation loss function, which can be defined as: L total =L CE +L distill (5); Among them, L in the above formula (5) total is the cross entropy loss function, and its calculation formula is: In the above formula (6), C is the summary of the category, y o,c is the one-hot encoding of the true label, p c Predict the probability that the sample belongs to category C for the model; In the above formula (5), L distill is the distillation loss function, and its calculation formula is: In the above formula (7), T(X) c and S(X) c Respectively represent the probability that the teacher model and the student model predict that the input X belongs to category C; The surface cover classification sample dataset containing multi-source remote sensing data and true value data constructed in the above step 2.1) is respectively fed into the teacher network and the student network for training to generate U(X) and S(X) respectively. The cross entropy loss function L total During the training process, it will gradually converge, making S(X) gradually approach U(X) and the true value; 2.3) Real-time extraction of land cover types The student network is reconstructed using JAVA language in the GEE platform, and the parameters of the trained student network are imported into the reconstructed network for initialization; then, by inputting the normalized multi-source remote sensing data with 15 bands, the labels of the corresponding land cover types can be obtained in real time; 2.4) Classification of accessible / inaccessible surfaces Based on the surface cover classification results, it is determined whether the area is drivable or not.

2. The method for evaluating regional accessibility potential taking into account multi-factor stress according to claim 1, characterized in that: The calculation of terrain complexity in step 1) is obtained by weighting terrain roughness, surface texture and terrain state index. The calculation process involves DEM data and optical remote sensing data. The specific inversion calculation process is as follows: TCF=w1×R+w2×T+w3×SI (1); Wherein, in the above formula (1), w1, w2 and w3 are weight coefficients, w1 and w2 are 0.4, w3 is 0.2; R is the surface roughness; T is the surface texture coefficient; SI is the surface shape index; The surface roughness R in the above formula (1) is used to characterize the undulation of the terrain. The specific calculation formula is as follows: In the above formula (2), z i is the elevation value of each sampling point, is the average elevation of all sampling points, and N is the total number of sampling points; The surface texture coefficient T in the above formula (1) is used to characterize the complexity of the surface texture in the optical image. The specific calculation formula is as follows: Through the calculation of the gray-level co-occurrence matrix, P(i,j) represents the frequency of occurrence of gray levels i and j under the conditions specified by distance and direction, and G is the number of gray levels; The surface shape index SI in the above formula (1) is used to calculate the curvature change of the terrain in different axes. The specific formula is as follows: In the above formula (4), k1 is the curvature of the terrain in the horizontal direction, and k2 is the curvature of the terrain in the y longitudinal direction.

3. The method for evaluating regional accessibility potential taking into account multiple stress factors as claimed in claim 2, characterized in that: The calculation result of the terrain complexity needs to be normalized so that its value is within the range of 0-1.

4. The method for evaluating regional accessibility potential taking into account multiple stress factors as claimed in claim 1, characterized in that: The terrain elements in step 1) include slope and terrain relief; The step 1) is to calculate the slope through DEM data, and distinguish whether it is a passable area according to the slope value, wherein 15° is agreed to be the critical value, that is, the area with a slope greater than 15° is an impassable area, and the area with a slope less than or equal to 15° is a passable area; the weight type of the slope is interval type, when the slope is greater than 15, it is assigned a value of 0, and when the slope is less than 15, it needs to be divided by 15 for normalization so that its value range is 0-1.

5. The method for evaluating regional accessibility potential taking into account multiple stress factors as claimed in claim 1, characterized in that: The data source for the collection of open source data of existing road network information in the step 2) mainly comes from the open source road network information published by OpenStreetMap or official agencies; the weight type of the open source road network information is interval type, in which the wide roads such as national roads are assigned a value of 1.0, the main roads with moderate width are assigned a value of 0.8, the secondary roads are assigned a value of 0.5, the residential roads are assigned a value of 0.4, and the narrow roads such as non-urban dirt roads are assigned a value of 0.

2.

6. The method for evaluating regional accessibility potential taking into account multiple stress factors as claimed in claim 1, characterized in that: In the step 2.4), water bodies and ice and snow building areas in the surface cover classification results are classified as inaccessible areas, and the classification results are all assigned a value of 0; bare land, cultivated land, grassland, and woodland in the surface cover classification results are classified as accessible areas, and are assigned values ​​of 0.8, 0.5, 0.3, and 0.1, respectively; and in the identification process, it is necessary to pay attention to the fact that some areas have seasonal wetlands or land freeze-thaw conditions due to the influence of rainfall and temperature, and the weight coefficient can be reduced for such areas.

7. The method for evaluating regional accessibility potential taking into account multiple stress factors as claimed in claim 1, characterized in that: The time series meteorological data in step 3) are derived from meteorological stations in the study area or obtained through relevant open source remote sensing products; in terms of temperature, a value of 0.1 is assigned when the temperature is ≤-40°C, a value of 0.4 is assigned when the temperature is -40°C < temperature ≤-25°C, a value of 0.8 is assigned when the temperature is -25°C < temperature ≤32°C, a value of 0.4 is assigned when the temperature is 32°C < temperature ≤40°C, and a value of 0.1 is assigned when the temperature is >40°C; a value of 0.1 is assigned when the rainfall is >100 mm, a value of 0.2 is assigned when the rainfall is 50 mm < rainfall ≤100 mm, a value of 0.4 is assigned when the rainfall is 20 mm < rainfall ≤50 mm, a value of 0.6 is assigned when the rainfall is 5 mm < rainfall ≤20 mm, and a value of 0.8 is assigned when the rainfall is 0 mm ≤ rainfall ≤5 mm; The data sources for inverting high-resolution soil moisture information based on surface cover types and open source kilometer-level soil moisture products in step 3) include Landsat8 visible and thermal infrared data, Sentinel-1 radar data and terrain data; Landsat8 data is used for calculating water body index and extracting emissivity information, and Sentinel-1 and terrain data are used to obtain backscatter coefficient data after terrain correction.

8. The method for evaluating regional accessibility potential taking into account multi-factor stress according to claim 1 or 7, characterized in that: The specific process of inverting high-resolution soil moisture information based on the surface cover type and the open source kilometer-level soil moisture product in step 3) is as follows: 3.1) Data preparation and preprocessing: Collect SMAP soil moisture content products, Landsat8 reflectance and thermal infrared images, Sentinel-1 backscatter coefficient data and terrain data; use Landsat data to calculate the water body index NDWI, the calculation formula is as follows: In the above formula (8), GREEN represents the green light band, and NIR represents the near infrared band; The spatial resolution of the water body index NDWI, Landsat8 thermal infrared data, Sentinel-1 backscatter coefficient data, DEM data and the surface cover type data obtained in step 2) above was unified to 30 meters, and superimposed into a multi-source soil moisture sample dataset with 6 bands; based on the surface cover type data in the dataset, J areas with an area greater than 1km were selected. 2 The bare land is taken, and the coordinates of its center point are taken, and the set of NDWI index, thermal infrared emissivity data, backscattering coefficient data and terrain data corresponding to the coordinates are extracted as Xt, and the corresponding SMAP soil moisture value is taken as Yt; the sample data set of soil moisture inversion composed of J groups (Xt, Yt) of data; 3.2) Model training and soil moisture inversion The feature vector Xt and target variable Yt are trained using the random forest algorithm. The performance of the model is evaluated. If the requirements are met, the model can be applied to new data sets to achieve the inversion of medium- and high-resolution soil moisture content. The inverted results will be normalized to the range of 0-1.

9. The method for evaluating regional accessibility potential taking into account multiple stress factors as claimed in claim 1, characterized in that: The data source for extracting abnormal surface deformation areas associated with geological disaster risks using deep learning and InSAR data involved in step 4) comes from Sentinel-1 radar data. Surface deformation phase data is obtained through preprocessing, and combined with fine-tuning technology based on a large visual model, intelligent extraction of abnormal surface deformation areas associated with wide-area geological disaster risks is achieved.

10. The method for evaluating regional accessibility potential taking into account multiple stress factors according to claim 1, characterized in that: The specific process of inverting the traffic safety factor in step 4) is as follows: 4.1) Data preprocessing The Sentinel-1 radar data were processed by the interferogram stacking time-series InSAR method to obtain the surface deformation phase data associated with geological hazards, namely the annual average surface deformation phase map. The linear phase rate p in the annual average surface deformation phase map is calculated as follows: In the above formula (9), Δt j is the time interval of the jth group of differential interference patterns; is the unwrapped differential interference phase of the jth group; Then, the InSAR deformation data of different causes are identified and made into true value labels for the identification of InSAR surface deformation anomaly areas, which are then merged with the annual average surface deformation phase map into a sample data set for intelligent identification of surface deformation phase, providing basic data for subsequent hidden danger identification applications. 4.2) Intelligent identification of abnormal surface deformation areas In the Pytorch environment, an intelligent recognition model InSARNet for abnormal surface deformation areas is built. The intelligent recognition model InSARNet for abnormal surface deformation areas consists of three parts: an image encoder, a keyword decoder and an image decoder. Then, based on the surface deformation phase intelligent recognition sample data set obtained in the above step 4.1), the intelligent recognition model InSARNet for abnormal surface deformation areas is trained. During the training process, the parameters of the image encoder are frozen, and only the embeddable fine-tuning module and the image decoder are trained. After the intelligent recognition model InSARNet for abnormal surface deformation areas is trained, the abnormal surface deformation areas associated with geological disaster risks in unknown areas can be identified by inputting InSAR deformation phase data. The identification result is expressed as a probability value of 0-1. The higher the probability value, the greater the possibility that it is an abnormal surface deformation area.

11. The method for evaluating regional accessibility potential taking into account multiple stress factors according to claim 10, characterized in that: The main steps of processing the Sentinel-1 radar data by the interferogram stacking time series InSAR method in step 4.1) include: 4.1.1) Splice the three sub-strip data of Sentinel-1TOPS into a complete SLC image; 4.1.2) Interference image pair selection; 4.1.3) Differential interferometry calculation; 4.1.4) Phase unwrapping; 4.1.5) Phase superposition processing; 4.1.6) Trend stripe removal; 4.1.7) Geocoding, finally obtaining the annual average surface deformation phase map; The image encoder is composed of a ViT model pre-trained by SAM and an embeddable fine-tuning module; the ViT model adopts the ViT-H / 16 version, which consists of four isometric global attention modules; the embeddable fine-tuning module consists of two sets of multi-layer perceptrons and two layers of ReLU activation layers; the input of the image encoder is an image of size 1024×1024×3; the output of the image encoder is a 16×16 image feature embedding processed by the patch embedding operation and the ViT module pre-trained by SAM; The image decoder consists of a ViT module and two layers of multi-layer perceptrons, which are used to restore the image feature embedding information into label data consistent with the size of the input image; the image encoder of InSARNet is initialized with parameters from a large visual model trained on massive natural images.

12. The method for evaluating regional accessibility potential taking into account multiple stress factors according to claim 1, characterized in that: The step 5) involves using the analytic hierarchy process to evaluate the regional traversability potential, and the specific steps are as follows: 5.1) Determine the hierarchical relationship Construct a hierarchical relationship based on "target layer-factor layer-element layer"; 5.2) Constructing the comparison matrix Construct a matrix M, the elements in it are marked as a ij ; Compare each element pairwise and determine the appropriate scale, then fill the determined scale into the corresponding position of the matrix M. The final matrix M is the pairwise comparison matrix; where a ij It represents the scale of i compared to j; a ij The corresponding position is the i-th row and j-th column in the matrix M; if i=j, it means that i and j are equally important, that is, a ij =1, which means that the scales of the main diagonal elements of the matrix M are all 1; a ij >0 and a ij ×a ji =1; 5.3) Calculation of weight coefficient The weights are calculated using the arithmetic mean method, geometric mean method and eigenvalue method respectively, and the arithmetic mean of the calculation results is calculated; Assume that the pairwise comparison matrix M is: For the arithmetic mean method, we can get the weight vector w i for: Find the weight vector w for the geometric mean method i , the calculation formula is as follows: Find the weight vector w for the eigenvalue method i , the calculation formula is as follows: 5.4) Consistency Check Calculate the largest eigenvalue λ of the eigenvector of the comparison matrix max : CW = C*w (16); In the above formulas (16)-(17), C is the comparison matrix, w is the indicator weight, and n represents the comparison matrix order; Calculate the consistency index CI of the comparison matrix: When CI is close to 0, it can be preliminarily judged that the constructed pairwise comparison matrix has satisfactory consistency; then calculate the consistency ratio CR: In the above formula (19), RI is the average consistency index; usually when CR < 0.1, the consistency of the pairwise comparison matrix can be considered acceptable; 5.5) Assessment of traversability potential The factor layer and the element layer are calculated according to the above steps 5.2) to 5.4) respectively to obtain the weight coefficients of the factor layer and the element layer, and the initial values ​​are weighted and summed to obtain the traversable potential assessment score; according to the traversable potential assessment score, the k-means clustering algorithm is used to classify the traversable potential assessment score.

13. The method for evaluating regional accessibility potential taking into account multiple stress factors according to claim 12, characterized in that: The target layer in step 5.1) is the regional traversability potential score; the factor layer includes topography, surface cover, climate environment, and traffic safety; the element layer includes slope, terrain complexity, roads, surface cover type, temperature, rainfall, soil moisture content, geological disaster risk information, and existing geological disaster points.

14. The method for evaluating regional accessibility potential taking into account multiple stress factors according to claim 12, characterized in that: The process of using the k-means clustering algorithm to classify the evaluation results in step 5.5) is as follows: predetermine the number of classification categories K, randomly select K objects {a1, a2, ..., a K } as given n data {x1,x2,…,x n }; calculate the distance between each data point and each cluster center, assign it to each cluster center according to the principle of closest distance, and recalculate each category center, so that the sum of the squares of the distances from each data point to the nearest cluster center is minimized; And the sum of the squares of the distances of each data point to the nearest cluster center is the objective function W n : Among them, W in the above formula (20) n is the distance objective function of clustering, a j is the initial cluster center, K is the number of classification categories, X i is the data to be clustered, and n is the number of data.

Images