Remote sensing image change detection and damage assessment rapid method and system for disaster emergency

Abstract

The application provides a disaster emergency-oriented remote sensing image change detection and damage assessment rapid method and system, relates to the disaster emergency and remote sensing technical field, first obtains a post-disaster remote sensing image data set covering a disaster occurrence area, and calls a pre-disaster building shadow area topological description set and a surrounding ground object layout description set of a corresponding spatial position from a pre-disaster ground object structure knowledge base, then carries out spatial position geometry calculation of the shadow area according to an angle parameter and building geometry information, extracts abnormal offset information of the shadow area, then constructs an indirect evidence chain to mark a suspected damaged building object, generates a distribution map, then carries out damage mode classification processing on the suspected damaged building object, determines a damage type mark and calculates an initial damage degree score, finally generates a damage assessment result data set containing a damaged building spatial position coordinate and a damage level code, and the building damage situation in disaster emergency can be rapidly and accurately evaluated.

Description

Technical Field

[0001] This invention relates to the fields of disaster emergency response and remote sensing technology, and more specifically, to a rapid method and system for detecting changes and assessing damage in remote sensing images for disaster emergency response. Background Technology

[0002] In the field of disaster emergency response, timely and accurate information on the damage caused to buildings is crucial, as it directly affects the rational allocation of rescue resources, the effective implementation of rescue operations, and subsequent reconstruction planning. Remote sensing imagery technology, with its advantages of large-scale and rapid information acquisition, has become an important means of obtaining information on building damage in disaster emergencies.

[0003] Currently, methods for assessing building damage based on remote sensing imagery mainly rely on directly observing changes in the building's appearance in the images, such as obvious features like collapse or breakage. However, these methods have several limitations. Firstly, after a disaster, accurately identifying building damage directly from remote sensing images is often difficult due to weather conditions (such as cloud cover) and image resolution. For example, in severe weather conditions like heavy rain or dense fog, remote sensing images may be blurry and unclear, making it difficult to clearly present the building's detailed features, leading to misjudgments or omissions. Secondly, relying solely on direct observation cannot comprehensively and deeply analyze the underlying causes and extent of building damage. For instance, some buildings may not appear to have collapsed externally, but may have internal problems such as foundation damage or structural deformation, which are difficult to detect through direct image observation alone. Summary of the Invention

[0004] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a rapid method for remote sensing image change detection and damage assessment for disaster emergency response, the method comprising: Acquire a set of post-disaster remote sensing image data covering the disaster-stricken area, and retrieve a set of topological descriptions of the shadow area of ​​pre-disaster buildings and a set of descriptions of the layout of pre-disaster buildings from a pre-constructed knowledge base of pre-disaster ground features, which correspond to the spatial location of the post-disaster remote sensing image data set. Based on the solar incidence angle parameters and satellite observation angle parameters in the post-disaster remote sensing image data set, and combined with the building geometric information in the pre-disaster building shadow area topology description set, the spatial location of the shadow area is geometrically calculated to generate a post-disaster expected shadow area spatial distribution map. The post-disaster expected shadow area spatial distribution map is then overlaid with the actual shadow area in the post-disaster remote sensing image data set to extract abnormal shadow area offset information. Based on the abnormal offset information of the shadow area and the set of pre-disaster descriptions of the surrounding land features of the building, an indirect evidence chain of building structural damage is constructed. The spatial location of suspected damaged building objects in the post-disaster remote sensing image data set is determined by the indirect evidence chain, and a distribution map of suspected damaged buildings carrying indirect damage evidence markers is generated. Each suspected damaged building in the suspected damaged building distribution map is classified according to its damage pattern to determine the damage type identifier of each suspected damaged building and calculate its initial damage score. A damage assessment result dataset containing the spatial location coordinates of the damaged buildings and the damage level code is generated based on the damage type identifier and the initial damage degree score.

[0005] Furthermore, embodiments of the present invention also provide a rapid system for detecting changes and assessing damage in remote sensing images for disaster emergency response, characterized in that it includes: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to execute the aforementioned rapid method for remote sensing image change detection and damage assessment for disaster emergency response by executing the machine-executable instructions.

[0006] In another aspect, embodiments of the present invention also provide a computer program product, the computer program product including machine-executable instructions, the machine-executable instructions being stored in a computer-readable storage medium, a processor of a computer device reading the machine-executable instructions from the computer-readable storage medium, the processor executing the machine-executable instructions, causing the computer device to execute the aforementioned rapid method for detecting changes in remote sensing images and assessing damage for disaster emergency response.

[0007] Based on the above, by acquiring a set of post-disaster remote sensing imagery and combining it with a pre-constructed knowledge base of pre-disaster land cover structures, this method fully utilizes multi-source information from both before and after the disaster. It uses angular parameters from post-disaster remote sensing images and a set of pre-disaster building shadow area topological descriptions to perform geometric calculations of the spatial location of shadow areas, generating a spatial distribution map of the expected post-disaster shadow areas. This map is then overlaid with the actual shadow areas to extract anomalous offset information. Changes in shadow areas are used to indirectly reflect the damage to buildings, overcoming the limitations of direct observation images due to weather and resolution constraints, and improving the accuracy and reliability of damage detection. Based on the anomalous offset information in the shadow areas and the pre-disaster description set of surrounding land cover layouts, an indirect evidence chain is constructed. The spatial location of suspected damaged buildings is then determined, generating a distribution map of suspected damaged buildings carrying indirect damage evidence markers, further clarifying the location and extent of the damaged buildings. The system classifies suspected damaged buildings by damage pattern, identifies damage type identifiers, and calculates initial damage scores. Finally, it generates a damage assessment dataset containing the spatial coordinates of the damaged buildings and damage level codes, enabling a comprehensive and detailed assessment of building damage. This allows for rapid and accurate assessment of building damage during disaster emergencies. Attached Figure Description

[0008] Figure 1 This is a schematic diagram of the execution flow of a rapid method for detecting changes and assessing damage in remote sensing images for disaster emergency response, provided in an embodiment of the present invention.

[0009] Figure 2 This is a schematic diagram of exemplary hardware and software components of a rapid remote sensing image change detection and damage assessment system for disaster emergency response provided in an embodiment of the present invention. Detailed Implementation

[0010] Figure 1 This is a flowchart illustrating a rapid method for detecting changes and assessing damage in remote sensing images for disaster emergency response, provided by an embodiment of the present invention. A detailed description follows.

[0011] Step S110: Obtain a set of post-disaster remote sensing image data covering the disaster-stricken area, and retrieve a set of topological descriptions of the shadow area of ​​pre-disaster buildings and a set of descriptions of the layout of pre-disaster buildings from a pre-constructed knowledge base of pre-disaster land cover structures, which correspond to the spatial location of the post-disaster remote sensing image data set.

[0012] In this embodiment, after a strong typhoon passes, the emergency response processing platform first acquires a high-resolution post-disaster remote sensing image dataset covering the entire affected coastal city. This dataset contains multiple satellite image units, such as image unit ImgOpt_Post_01 taken by optical remote sensing satellite A and image unit ImgSar_Post_02 taken by synthetic aperture radar satellite B. Each image unit includes metadata containing information such as spatial coverage boundary coordinates, image acquisition timestamp, and sensor type. For example, the metadata of image unit ImgOpt_Post_01 records the four boundary coordinates of its coverage area: longitude Lon_minA to Lon_maxA, latitude Lat_minA to Lat_maxA, and image acquisition timestamp T_cap_01. Simultaneously, the processing accesses a pre-constructed pre-disaster ground feature structure knowledge base, which stores detailed ground feature information of the city before the disaster, organized using geographic grid cells as the basic storage unit. The processing method uses the spatial coverage boundary coordinates of each post-disaster remote sensing image unit as the query basis to accurately retrieve two core data corresponding to its spatial location from the pre-disaster land cover structure knowledge base: one is the set of topological descriptions of the shadow area of ​​pre-disaster buildings, describing the geometric relationship between pre-disaster buildings and their shadows; the other is the set of descriptions of the layout of land features around pre-disaster buildings, describing the spatial layout relationship between pre-disaster buildings and surrounding roads, vegetation, neighboring buildings, and other elements. The details of this retrieval process are implemented through the following sub-steps.

[0013] Step S111: Analyze the spatial coverage boundary coordinates and image acquisition timestamp of each post-disaster remote sensing image unit in the post-disaster remote sensing image data set; generate image spatial index code based on the spatial coverage boundary coordinates; and determine solar position parameters and satellite position parameters based on the image acquisition timestamp.

[0014] Specifically, for each post-disaster remote sensing image unit in the post-disaster remote sensing image dataset, such as an image file named ImgOpt_Post_01, metadata parsing is performed. The geographic coordinates of the four corner points of its spatial coverage area are extracted: minimum longitude Lon_min, minimum latitude Lat_min, maximum longitude Lon_max, and maximum latitude Lat_max, forming a boundary coordinate rectangle. Based on this boundary coordinate rectangle, the four double-precision floating-point numbers Lon_min, Lat_min, Lon_max, and Lat_max are hierarchically encoded using a geohashing algorithm to generate a string-based image spatial index code, such as "GeoHash_Alpha". Simultaneously, based on the image acquisition timestamp T_cap_01 and combined with a celestial mechanics model, the solar azimuth angle parameter Sun_Azi and solar altitude angle parameter Sun_Ele, as well as the satellite azimuth angle parameter Sat_Azi and altitude angle parameter Sat_Ele, are calculated over the urban area at time T_cap_01.

[0015] Step S112: Using the image spatial index code as the retrieval keyword, perform a spatial range matching query in the spatial index partition table of the pre-disaster ground feature structure knowledge base. The pre-disaster ground feature structure knowledge base stores parametric description data partitions of pre-disaster ground feature structures organized by geographic grid units.

[0016] Specifically, the generated image spatial index code "GeoHash_Alpha" is used as the search keyword to perform exact matching and range-expanded queries in the spatial index partition table of the pre-disaster ground feature structure knowledge base. This pre-disaster ground feature structure knowledge base is constructed by dividing the geographic space into multiple geographic grid units of predetermined sizes, for example, each unit is 1 km. A 1km square area is defined, with each geographic grid cell corresponding to a unique grid code, such as "Grid_Code_001" to "Grid_Code_N". Each geographic grid cell stores parametric description data of all land features within its range, forming an independent data partition. The spatial index partition table records the mapping relationship between each grid code and its corresponding data partition storage address. The query process aims to find all geographic grid cells represented by grid codes whose spatial extent intersects with the spatial extent represented by the image spatial index code "GeoHash_Alpha".

[0017] Step S113: Obtain the pre-disaster land cover structure parameterized description data partitions corresponding to all geographic grid units that have spatial overlap with the image spatial index code in the spatial range matching query results as the initial screening data partition set.

[0018] Specifically, after the spatial extent matching query is completed, a list of geographic grid cell codes is obtained, for example, containing four codes: "Grid_Code_125", "Grid_Code_126", "Grid_Code_145", and "Grid_Code_146". The spatial extent of these four geographic grid cells overlaps with the coverage area of ​​the ImgOpt_Post_01 image. Based on these four grid codes, their corresponding pre-disaster ground feature structure parameterized description data partitions, namely Partition_125, Partition_126, Partition_145, and Partition_146, are retrieved from the storage system of the pre-disaster ground feature structure knowledge base. These four data partitions are then combined to form the initial screening data partition set for the current post-disaster remote sensing image unit ImgOpt_Post_01.

[0019] Step S114: For each preliminary screening data partition unit in the preliminary screening data partition set, extract the building shadow area topology description sub-block stored in the preliminary screening data partition unit. The building shadow area topology description sub-block contains a sequence of shadow area outline polygons of multiple building objects under standard solar incidence angle conditions and the encoding of the topological connection relationship between each shadow area outline polygon and the corresponding building base polygon.

[0020] Specifically, each initial data partition unit in the initial data partition set is traversed; for example, partition data Partition_125 is processed first. The building shadow area topology description sub-block ShadowDesc_125 stored within Partition_125 is extracted. This building shadow area topology description sub-block ShadowDesc_125 is a complex data structure containing a list, where each element corresponds to a building object. For a building object with the ID Bldg_ID_001, its record includes: a sequence of polygons representing the shadow area outline, PolySeq_Shd_Std_B001, calculated by projection from a 3D model under preset standard solar incidence angle conditions (e.g., solar altitude angle Std_Ele, azimuth angle Std_Azi). This sequence may contain multiple polygons to describe the shadows cast by the building's main structure and its attached structures. It also includes a topological connection code, Link_B001, which is a tuple (base polygon ID_Base_B001, shadow polygon ID_Shd_B001) that explicitly indicates which building base polygon, Polygon_Base_B001, each polygon in the shadow area outline polygon PolySeq_Shd_Std_B001 is connected to.

[0021] Step S115: For each preliminary screening data partition unit in the preliminary screening data partition set, extract the building surrounding feature layout description sub-block stored in the preliminary screening data partition unit. The building surrounding feature layout description sub-block includes the road centerline vector line segment within a preset distance range around each building object, the base polygon of the adjacent building, the boundary polygon of the vegetation cover area, and the spatial distance parameters and azimuth parameters between each corresponding feature element and the central building object.

[0022] Specifically, continuing with the initial screening data partitioning unit Partition_125 as an example, the layout description sub-block LayoutDesc_125 surrounding the building is extracted from it. This layout description sub-block LayoutDesc_125 is also organized using the building object as the index. For the same building object Bldg_ID_001, its surrounding layout description information includes: the set of vector line segments of all road centerlines Set_Road_Seg_B001 within a buffer range with the geometric center of its base polygon Polygon_Base_B001 as the center and a preset distance threshold R_surround (e.g., 200M) as the radius; the set of base polygons of all adjacent buildings Set_Nearby_Base_B001; and the set of boundary polygons of all vegetation-covered areas Set_Veg_Boundary_B001. For each extracted feature, such as a road centerline segment Road_Segment_01, the system also stores the shortest distance value Dist_Road_B001 from the base polygon of the building Bldg_ID_001, as well as the azimuth parameter Angle_Road_B001 of the segment relative to the geometric center of the building Bldg_ID_001. These distance and azimuth parameters together constitute a refined description of spatial relationships.

[0023] Step S116: Based on the solar position parameters and the satellite position parameters, perform projection transformation processing on the shadow area outline polygon sequence in the topological descriptor block of the building shadow area, transform the shadow area outline polygon sequence from the reference coordinate system under the standard solar incidence angle condition to the projection coordinate system under the actual observation geometry condition corresponding to the image acquisition timestamp, and generate the projection-transformed shadow area outline polygon sequence.

[0024] Specifically, at this point, the solar position parameters Sun_Azi and Sun_Ele and the satellite position parameters Sat_Azi and Sat_Ele at the actual imaging time are obtained. A geometric transformation needs to be performed on the sequence of shadow area contour polygons stored in the building shadow area topology descriptor block, based on the standard solar incidence angle condition. Taking the standard shadow polygon Polygon_Shd_Std_B001 of building Bldg_ID_001 as an example, its coordinates are defined in a local reference coordinate system O_std with the building base as the origin. First, a rotation matrix R_sun is constructed from the standard conditions to the actual illumination conditions based on the actual solar azimuth angle Sun_Azi and elevation angle Sun_Ele. Using R_sun, the coordinates of all vertices of the standard shadow polygon Polygon_Shd_Std_B001 are rotated to obtain the intermediate shadow polygon Polygon_Shd_Int_B001. Then, based on the actual satellite observation directions Sat_Azi and Sat_Ele, a perspective projection transformation matrix M_sat is constructed to project the intermediate shadow polygon Polygon_Shd_Int_B001 from the three-dimensional object space to the two-dimensional image space, generating a projection-transformed shadow area outline polygon Polygon_Shd_Proj_B001 that is geometrically consistent with the actual image. This transformation is performed on all buildings to obtain the projection-transformed shadow area outline polygon sequence Set_Polygon_Shd_Proj.

[0025] Step S117: Associate and combine the projection-transformed shadow area outline polygon sequence with the corresponding building base polygon to construct a pre-disaster building shadow area topological description set with the building base polygon as the index key.

[0026] Specifically, the sequence of shadow area contour polygons after projection transformation generated in step S116, Set_Polygon_Shd_Proj, is associated with the untransformed building base polygons directly obtained from the pre-disaster ground feature structure knowledge base. Taking building Bldg_ID_001 as an example, its base polygon is Polygon_Base_B001, and its corresponding projection-transformed shadow polygon is Polygon_Shd_Proj_B001. These two polygons are stored as a key-value pair (Key: Polygon_Base_B001, Value: Polygon_Shd_Proj_B001) in a hash map table HashMap_Shadow. At the same time, the topological connection relationship code Link_B001 extracted in step S114 is also attached as an attribute to this key-value pair. After traversing all building objects, HashMap_Shadow constitutes a complete set of pre-disaster building shadow area topological descriptions indexed by building base polygons. This set accurately describes the shadow area shape that each building should have produced under actual imaging conditions.

[0027] Step S118: Reconstruct the spatial topology of the road centerline vector segment, the base polygon of the adjacent building, and the boundary polygon of the vegetation cover area in the sub-block describing the layout of the building's surrounding features, and establish a topology map of the surrounding features layout with the central building object as the node, which serves as the pre-disaster building surrounding features layout description set.

[0028] Specifically, the road centerline vector segments, adjacent building base polygons, and vegetation cover area boundary polygons extracted from the LayoutDesc_125 sub-block of the building's surrounding features are initially stored in a list format centered on the buildings. To facilitate subsequent change analysis, their spatial topology needs to be reconstructed. The building Bldg_ID_001 is designated as the central node Node_B001, whose attributes include the base polygon Polygon_Base_B001. For each road centerline segment Road_Segment_01, it is abstracted into a road node Node_Road_01, whose attributes include its vector segment coordinates. Then, the topological relationship between building node Node_B001 and road node Node_Road_01 is calculated, including connectivity (i.e., if the distance is less than the threshold T_dist, they are considered connected) and spatial relationship parameters (distance Dist_Road_B001, azimuth angle Angle_Road_B001). This relationship is then added to the graph as an edge Edge_B001_Road01, with the edge's attributes including the aforementioned distance and azimuth angle. Similarly, perform similar operations on the base polygons of adjacent buildings and the boundary polygons of vegetation-covered areas to establish topological edges between building nodes and adjacent building nodes and vegetation nodes. Ultimately, this results in a topological graph Graph_Layout_B001, with the central building object as the core node and surrounding features as subordinate nodes, connected by edges with spatial relationship attributes. Performing this operation on all buildings yields a series of topological graphs, collectively forming a set of pre-disaster descriptions of the surrounding feature layout of buildings.

[0029] Step S120: Based on the solar incidence angle parameters and satellite observation angle parameters in the post-disaster remote sensing image data set, and combined with the building geometric information in the pre-disaster building shadow area topology description set, perform geometric calculation of the spatial location of the shadow area to generate a post-disaster expected shadow area spatial distribution map. Then, overlay the post-disaster expected shadow area spatial distribution map with the actual shadow area in the post-disaster remote sensing image data set to extract the shadow area abnormal offset information.

[0030] In this embodiment, after obtaining the topological description set of building shadow areas before the disaster, it is necessary to predict the location where building shadows should appear on the post-disaster image based on this, and compare it with the actual shadows on the image to discover shadow anomalies caused by changes in building height or shape (possibly due to damage). The specific calculation and analysis process is implemented by the following sub-steps.

[0031] Step S121: Extract the solar azimuth angle, solar altitude angle, satellite azimuth angle, and satellite altitude angle values ​​corresponding to each post-disaster remote sensing image unit from the metadata of the post-disaster remote sensing image data set, and use them as solar incidence angle parameters and satellite observation angle parameters.

[0032] Specifically, for each post-disaster remote sensing image unit in the post-disaster remote sensing image dataset, such as ImgOpt_Post_01, its metadata file is parsed again. The precise solar azimuth angle (Sun_Azi_Post), solar altitude angle (Sun_Ele_Post), and satellite azimuth angle (Sat_Azi_Post) and altitude angle (Sat_Ele_Post) for the image unit's imaging time are obtained directly from specific fields of the metadata or calculated using auxiliary ephemeris data. These four values ​​serve as unique solar incidence angle and satellite observation angle parameters for this image unit, and are used for subsequent shadow calculations for this image.

[0033] Step S122: Obtain the base polygon coordinate sequence of each building object in the pre-disaster building shadow area topology description set and the height attribute parameter of the building object. The height attribute parameter is stored in the pre-disaster ground feature structure knowledge base and associated with the building object identifier.

[0034] Specifically, from the pre-disaster building shadow area topological description set HashMap_Shadow constructed in step S117, each index key, i.e., the building base polygon, is traversed. For each base polygon, such as Polygon_Base_B001, a secondary query is performed in the core data table of the pre-disaster ground feature structure knowledge base, using its building object identifier Bldg_ID_001 as the key. This core data table stores the detailed attributes of each building object, including the height attribute parameter H_B001, which is usually in meters and represents the relative height or altitude of the building. The coordinate sequence (a set of ordered two-dimensional points) of the base polygon Polygon_Base_B001 is associated with the height attribute parameter H_B001 as input for subsequent individual building shadow calculations.

[0035] Step S123: Calculate the three-dimensional unit vector of the solar incidence direction based on the solar azimuth angle value and the solar altitude angle value, and calculate the three-dimensional unit vector of the satellite observation direction based on the satellite azimuth angle value and the satellite altitude angle value.

[0036] Specifically, using the conversion relationship between spherical and Cartesian coordinates, the three-dimensional unit vector Vec_Sun representing the direction of solar incidence is calculated based on the solar azimuth angle Sun_Azi_Post and the solar altitude angle Sun_Ele_Post. The calculation process is as follows: first, the angle is converted to radians; then, the x-component of Vec_Sun is equal to the cosine function value (solar altitude angle) multiplied by the sine function value (solar azimuth angle); the y-component is equal to the cosine function value (solar altitude angle) multiplied by the cosine function value (solar azimuth angle); and the z-component is equal to the sine function value (solar altitude angle). Similarly, the three-dimensional unit vector Vec_Sat representing the satellite observation direction is calculated based on the satellite azimuth angle Sat_Azi_Post and the satellite altitude angle Sat_Ele_Post. These two unit vectors, Vec_Sun and Vec_Sat, accurately describe the illumination direction and observation direction at the moment of imaging.

[0037] Step S124: For each building object, based on its height attribute parameter and the direction of sunlight, calculate the coordinates of the shadow projection points of the vertices of its base polygon boundary on the horizontal plane, and connect the coordinates of the shadow projection endpoints corresponding to all vertices to form the initial shadow area polygon generated by the building object under the illumination of sunlight.

[0038] Specifically, taking the building object Bldg_ID_001 as an example, its base polygon Polygon_Base_B001 consists of a series of vertices, assuming there are K vertices arranged counterclockwise as P1, P2, ..., Pk. For each vertex Pi (coordinates Xi, Yi, Zi_base, where Zi_base is the base elevation), considering the building height H_B001, the position of the top of the building corresponding to point Pi is Qi, with coordinates Xi, Yi, Zi_base + H_B001. The calculation of the shadow projection point Pi_shd essentially involves finding the intersection of the straight line passing through point Qi from the opposite direction of the sun's incident direction (i.e., the direction of the ray) with the horizontal plane (Z = Zi_base). Since the unit vector Vec_Sun of the sun's incident direction has already been obtained, the parametric equation of this line is Qi + t. Vec_Sun. Set its Z-coordinate to Zi_base, solve for the parameter t, and then calculate the X and Y coordinates of Pi_shd. Perform the above operation on all vertices to obtain a set of shadow projection points Pi_shd_1 to Pi_shd_k. Connect these k points in the same order as the vertices of the original base polygon to form a closed polygon, which is the initial shadow area polygon Polygon_Shd_Prelim_B001 produced by the building Bldg_ID_001 under sunlight.

[0039] Step S125: Perform visibility clipping processing on the preliminary shadow area polygon. Based on the three-dimensional unit vector of the satellite observation direction, determine the part of the preliminary shadow area polygon that is not visible due to being obscured by the building itself. Remove the invisible part from the preliminary shadow area polygon to generate the post-disaster expected shadow area polygon that meets the satellite observation geometry conditions.

[0040] Specifically, the initial shadow area polygon Polygon_Shd_Prelim_B001 is theoretically the complete shadow projected onto the ground. However, since the satellite is observing from a specific direction, some shadow areas may be obscured by the buildings themselves and cannot be seen by the sensors. These shadows should not appear in the final image, so visibility clipping is required. Using the 3D model of the building Bldg_ID_001 (which can be obtained by simply stretching the base polygon and height H_B001) and the satellite observation direction unit vector Vec_Sat as input, a hidden surface removal algorithm is executed. Specifically, starting from the satellite observation direction, the occlusion relationship between the building model and the initial shadow area polygon is determined. Methods based on ray projection or depth buffers can be used. For each point within the initial shadow polygon, it is determined whether the ray emitted from that point in the opposite direction of Vec_Sat intersects with the building model. If they intersect, the point is obscured and invisible. All visible shadow points are refitted into one or more polygons, ultimately generating the post-disaster expected shadow area polygon Polygon_Shd_Expect_B001 that conforms to the satellite observation geometry. The polygon eliminates the portion obscured by the building itself, and is consistent with the shadow shape that can be observed in the actual image.

[0041] Step S126: Merge the polygons of the expected shadow areas of all building objects after the disaster, eliminate the overlapping areas between adjacent shadow area polygons, and generate a spatial distribution map of the expected shadow areas of the disaster covering the entire image range.

[0042] Specifically, all building objects within the coverage area of ​​the current post-disaster remote sensing image unit are traversed. Steps S124 and S125 are executed for each building object to obtain its corresponding post-disaster expected shadow area polygons, such as Polygon_Shd_Expect_B001, Polygon_Shd_Expect_B002, etc. All of these polygons are then input into a spatial merging operator. This operator first performs intersection calculations on the boundaries of all polygons to identify overlapping areas. Then, through Boolean union operations, the polygons are merged into one or more non-overlapping complex polygons. The final output is a rasterized or vectorized layer, namely, a spatial distribution map of the post-disaster expected shadow area, Map_Shd_Expect, covering the entire image area. Each pixel or region in this map is labeled as either "expected shadow" or "unexpected shadow."

[0043] Step S127: Automatically extract shadow areas from the post-disaster remote sensing image data set. By analyzing the spectral and texture features of each pixel in the image, identify the pixel areas belonging to shadows, and aggregate and connect the identified shadow pixel areas to generate a spatial distribution map of the actual shadow areas after the disaster.

[0044] Specifically, shadow detection is performed on the original post-disaster remote sensing image unit ImgOpt_Post_01. First, the image is converted from the RGB color space to the HSV or HSI color space because shadow areas typically have higher hue and lower value / intensity. For each pixel P(i, j) in the image, its spectral feature vector F_spec is extracted, containing the reflectance values ​​of each band. Simultaneously, its texture features F_tex are extracted, such as contrast and entropy calculated through local binary mode or gray-level co-occurrence matrix. F_spec and F_tex are combined into a feature vector and input into a pre-trained pixel-level shadow classifier (such as random forest or support vector machine). The classifier outputs the probability P_shadow that the pixel belongs to the shadow. Pixels with a probability greater than a preset threshold Th_shadow are marked as shadow pixels. After classifying all pixels, a binary shadow labeling map is obtained. Then, connected component analysis is performed on this binary map, aggregating adjacent shadow pixels together to form independent connected regions, and assigning a unique identifier ID to each region. Finally, the boundaries of the above connected regions are vectorized and attributed to generate a spatial distribution map of the actual shadow area after the disaster, Map_Shd_Actual, which accurately describes the shadow location and shape actually observed on the image.

[0045] Step S128: Compare the expected spatial distribution map of the post-disaster shadow area with the actual spatial distribution map of the post-disaster shadow area pixel by pixel, mark the pixel positions where shadows are expected to exist but do not actually exist and the pixel positions where shadows are expected not to exist but actually exist, and generate a shadow area difference marker binary map.

[0046] Specifically, the spatial distribution map of the expected shadow area generated in step S126 (Map_Shd_Expect, usually in vector or raster form) is spatially registered and compared pixel by pixel with the spatial distribution map of the actual shadow area generated in step S127 (Map_Shd_Actual, raster form). For each pixel location Pix(x, y), its value Val_exp in Map_Shd_Expect (1 indicates expected shadow, 0 indicates unexpected shadow) and its value Val_act in Map_Shd_Actual (1 indicates actual shadow, 0 indicates unrealistic shadow) are queried. Val_exp and Val_act are compared. If Val_exp equals 1 and Val_act equals 0, the pixel is marked as a first-class difference pixel, representing "missing shadow". If Val_exp equals 0 and Val_act equals 1, the pixel is marked as a second-class difference pixel, representing "abnormal shadow". All pixels marked as differing are assigned a value of 1 in the output image, while the remaining pixels are assigned a value of 0, thus generating a binary map of shadow area differences, Map_Diff. This binary map visually shows the locations where the expected shadow differs from the actual shadow.

[0047] Step S129: Calculate the spatial offset direction and spatial offset distance between each connected difference region in the binary map of the shadow area difference marker and the base polygon of the adjacent building. Combine the spatial offset direction and the spatial offset distance into the shadow area abnormal offset vector corresponding to each difference region. The shadow area abnormal offset vectors corresponding to all difference regions constitute the shadow area abnormal offset information.

[0048] Specifically, a connected component analysis is performed on the binary map Map_Diff, which identifies all connected difference regions, such as region Diff_Region_01. For each connected difference region, its geometric center coordinates, Center_Diff, are first calculated. Then, the nearest building base polygon to this center point is searched in space. For example, the base polygon Polygon_Base_B001 of building Bldg_ID_001 is found, with its geometric center at Center_B001. The vector Vec_Offset pointing from Center_B001 to Center_Diff is calculated. The direction of this vector is the spatial offset direction, Dir_Offset, and the magnitude of this vector is the spatial offset distance, Dist_Offset. Dir_Offset and Dist_Offset are combined to form the shadow area anomaly offset vector Vec_Anomaly corresponding to the difference region. If the difference region is a shadow change caused by building damage, this offset vector reflects the direction and distance of the shadow's abnormal movement relative to the building base. Record the abnormal offset vectors of all connected difference regions to form a list List_Anomaly_Vec, which is the shadow area abnormal offset information extracted from the current image unit.

[0049] Step S130: Construct an indirect evidence chain for building structural damage based on the abnormal offset information of the shadow area and the pre-disaster description set of the surrounding land features of the building. Use the indirect evidence chain to spatially locate the suspected damaged building objects in the post-disaster remote sensing image data set and generate a distribution map of suspected damaged buildings carrying indirect damage evidence markers.

[0050] In this embodiment, shadow anomalies alone are insufficient to definitively determine building damage. It is necessary to combine changes in surrounding features to form a mutually corroborating chain of indirect evidence, thereby more reliably identifying suspected damaged buildings. The specific construction and calibration process is implemented through the following sub-steps.

[0051] Step S131: Spatial matching is performed between each shadow area abnormal offset vector in the shadow area abnormal offset information and the building base polygon in the pre-disaster building shadow area topology description set to determine the candidate building object associated with each shadow area abnormal offset vector.

[0052] Specifically, from the list of shadow area anomaly offset information, List_Anomaly_Vec, obtained in step S129, each anomaly offset vector, Vec_Anomaly_K, is retrieved. This vector is associated with the center point of a difference region, Center_Diff_K. A spatial connectivity query is performed, using Center_Diff_K as the center and setting a search radius R_match (e.g., 50M), to find all building base polygons falling within this range. If only one building base polygon exists within the search range, such as Polygon_Base_B001, then that building, Bldg_ID_001, is identified as a candidate building object associated with this anomaly offset vector. If multiple exist, further filtering is needed based on the characteristics of the offset vector itself, such as selecting the building whose offset direction, Dir_Offset_K, is closest to the direction from which the expected shadow of that building points to its base, as the most likely associated object. Ultimately, each anomaly offset vector is associated with one or more candidate building objects, establishing a preliminary connection between the shadow anomaly and a specific building.

[0053] Step S132: For each candidate building object, count the number of all associated shadow area anomaly offset vectors, calculate the average offset direction and average offset distance of the shadow area anomaly offset vectors, and generate a shadow anomaly statistical feature descriptor for the candidate building object.

[0054] Specifically, taking the candidate building object Bldg_ID_001 as an example, all shadow area anomaly offset vectors associated with it determined in step S131 are collected. Assume there are M vectors in total, namely Vec_A_1, Vec_A_2, ..., Vec_A_M. First, the number of vectors, Count_Anomaly, is counted, which equals M. Then, the average offset direction, Dir_Avg, of these vectors is calculated. When calculating the average direction, the arithmetic mean of the angle values ​​cannot be simply calculated. Instead, the direction angle of each offset vector, Dir_A_m, should be converted into a unit vector. The x and y components of the unit vector are summed and normalized, and then the angle value is calculated back to obtain Dir_Avg. The average offset distance, Dist_Avg, is calculated, which is the arithmetic mean of the magnitudes Dist_A_m of the M vectors. The three values, Count_Anomaly, Dir_Avg, and Dist_Avg, are combined into a three-dimensional feature vector, Vec_ShadowStat_B001, which is the shadow anomaly statistical feature descriptor of the candidate building object.

[0055] Step S133: Obtain the topology map of the surrounding land features corresponding to the candidate building object in the pre-disaster building surrounding land feature layout description set, and extract the vertical distance parameter between the candidate building object and the center line of the surrounding road, the interval distance parameter between the candidate building object and the base of the adjacent building, and the area ratio parameter of the vegetation coverage area around the candidate building object from the surrounding land feature layout topology map.

[0056] Specifically, taking the candidate building object Bldg_ID_001 as an example, the corresponding topology graph Graph_Layout_B001 is retrieved from the set of surrounding feature layout topology graphs constructed in step S118. From this graph, the edge Edge_Road of type "Road" connected to the building node Node_B001 is queried. For each road edge, the distance parameter Dist_Road in its attributes is taken, and then the minimum value of all Dist_Roads is calculated as the vertical distance parameter Dist_Road_min_B001 between the building and the centerline of the surrounding road. Next, the edge Edge_Nearby of type "nearby building" is queried, and the distance parameter Dist_Nearby (i.e., the interval between buildings) in its attributes is taken. The average value of all Dist_Nearbys is calculated as the average interval distance parameter Dist_Nearby_avg_B001 between the building and the base of the neighboring building. Finally, the node Node_Veg of type "vegetation" is queried, the sum of the areas of all vegetation nodes Area_Veg_total is calculated, and then divided by the total buffer area Area_Buffer calculated based on the preset distance R_surround to obtain the area ratio parameter Ratio_Veg_B001 of the vegetation coverage area around the building. The above Dist_Road_min_B001, Dist_Nearby_avg_B001, and Ratio_Veg_B001 together constitute the static background features around the building.

[0057] Step S134: Extract the current status of ground features from the local image area covering the candidate building object in the post-disaster remote sensing image data set, identify the actual road centerline location, the base outline of adjacent buildings and the boundary of vegetation coverage area after the disaster, and generate a post-disaster surrounding ground feature layout update map.

[0058] Specifically, centered on the candidate building object Bldg_ID_001, a local image patch Patch_Post_B001 is cropped from the post-disaster remote sensing image unit ImgOpt_Post_01. Its range should be slightly larger than the buffer range R_surround used in step S133. For this image patch Patch_Post_B001, a series of deep learning models are applied to extract ground features. First, a pre-trained road extraction network (based on the U-Net architecture) is used. The image patch is input, and the output is a road probability map. After binarization and vectorization, a set of vector line segments representing the actual post-disaster road centerlines, Set_Road_Post_B001, is obtained. Then, a building extraction network is used to extract the base contours of adjacent buildings, resulting in a polygon set Set_Base_Post_B001. Finally, a vegetation extraction network is used to extract the boundaries of vegetation-covered areas, resulting in a polygon set Set_Veg_Post_B001. The three newly extracted feature layers are overlaid together and spatially registered to generate a post-disaster surrounding feature layout update map Map_Layout_Post_B001 centered on the building Bldg_ID_001.

[0059] Step S135: Compare the road centerline position in the post-disaster surrounding feature layout update map with the road centerline vector line segment in the surrounding feature layout topology map, and calculate the road centerline position offset and road continuity change index.

[0060] Specifically, the set of road centerlines, Set_Road_Post_B001, in the post-disaster surrounding feature layout update map Map_Layout_Post_B001 generated in step S134 is compared with the set of road centerlines, Set_Road_Pre_B001, in the pre-disaster surrounding feature layout topology map Graph_Layout_B001 generated in step S118. First, road centerline matching is performed to identify the line segments corresponding to the same road before and after the disaster. For each matched pair of road line segments, the Hausdorff distance between them is calculated as the offset of that road segment, Offset_Road_K. The offsets of all matched road segments are averaged to obtain the average road centerline offset, Offset_Road_avg. Regarding changes in road continuity, the number of connected components in the pre-disaster road network, Conn_Pre, is calculated, along with the number of connected components in the post-disaster road network, Conn_Post. If Conn_Post is greater than Conn_Pre, it indicates that the road has been broken due to damage, resulting in decreased continuity. The road continuity change index Change_Conn is defined as the difference between Conn_Post and Conn_Pre, with a positive value indicating an increase in breaks.

[0061] Step S136: Compare the outline of the base of the neighboring building in the updated map of the surrounding land features after the disaster with the polygon of the base of the neighboring building in the topology map of the surrounding land features, and calculate the rate of change of the base area of ​​the neighboring building and the offset of the base position of the neighboring building.

[0062] Specifically, the set of post-disaster adjacent building base outlines (Set_Base_Post_B001) is compared with the set of pre-disaster adjacent building base polygons (Set_Base_Pre_B001). First, building objects are matched, typically based on positional overlap (Intersection over Union). For each successfully matched adjacent building, such as building B002, its pre-disaster area (Area_Pre_B002) and post-disaster area (Area_Post_B002) are calculated. The area change rate (Rate_Area_B002) is defined as the absolute value of (Area_Post_B002 minus Area_Pre_B002) divided by Area_Pre_B002. The area change rate is calculated for all matched buildings, and the maximum or average value is taken as the adjacent building base area change rate index (Change_Area_Nearby). Simultaneously, the Euclidean distance between the pre-disaster and post-disaster geometric centers of the building is calculated as the building's position offset (Offset_Pos_B002). The average position offset of all matched buildings is calculated to obtain the average position offset of the base of the neighboring buildings, Offset_Pos_Nearby_avg.

[0063] Step S137: Compare the vegetation cover area boundary in the post-disaster surrounding land feature layout update map with the vegetation cover area boundary polygon in the surrounding land feature layout topology map, and calculate the vegetation cover area boundary change rate and the texture fragmentation index inside the vegetation cover area.

[0064] Specifically, the set of post-disaster vegetation cover area boundaries, Set_Veg_Post_B001, is compared with the set of pre-disaster vegetation cover area boundary polygons, Set_Veg_Pre_B001. First, the change in total vegetation cover area is calculated. Let the total pre-disaster vegetation area be Area_Veg_Pre_total, and the post-disaster vegetation area be Area_Veg_Post_total. The rate of change of vegetation cover area boundaries, Rate_Veg_Change, is defined as the absolute value of (Area_Veg_Post_total minus Area_Veg_Pre_total) divided by Area_Veg_Pre_total. Then, the texture fragmentation within the vegetation cover area is calculated. The pre-disaster and post-disaster vegetation areas are cropped from the original images, converted to grayscale, and their texture features are calculated, such as contrast, dissimilarity, and entropy calculated using the grayscale co-occurrence matrix. The fragmentation index Frag_Veg can be defined as the ratio of the entropy value Entropy_Post of the post-disaster vegetation texture to the entropy value Entropy_Pre of the pre-disaster vegetation texture. If the ratio is much greater than 1, it indicates that the vegetation area has become fragmented and disordered, and may have been affected by a disaster (such as fallen trees).

[0065] Step S138: Combine the road centerline position offset, road continuity change index, adjacent building base area change rate, adjacent building base position offset, vegetation cover area boundary change rate, and vegetation cover area internal texture fragmentation index to generate the surrounding land feature change feature vector of the candidate building object.

[0066] Specifically, the average road centerline offset (Offset_Road_avg) and road continuity change index (Change_Conn) calculated in step S135, the adjacent building base area change rate (Change_Area_Nearby) and average adjacent building base offset (Offset_Pos_Nearby_avg) calculated in step S136, and the vegetation cover area boundary change rate (Rate_Veg_Change) and vegetation cover area internal texture fragmentation index (Frag_Veg) calculated in step S137 are concatenated into a multi-dimensional vector in a fixed order. For example, a six-dimensional feature vector of surrounding land cover change (Vec_Change_B001) is constructed, with dimensions as follows: Offset_Road_avg, Change_Conn, Change_Area_Nearby, Offset_Pos_Nearby_avg, Rate_Veg_Change, and Frag_Veg. This vector quantifies the degree of change in the surrounding environment of candidate buildings before and after a disaster.

[0067] Step S139: The shadow anomaly statistical feature descriptor is concatenated and fused with the surrounding land feature change feature vector to generate the indirect evidence chain feature vector of the candidate building object.

[0068] Specifically, the shadow anomaly statistical feature descriptor Vec_ShadowStat_B001 (a three-dimensional vector) generated in step S132 is concatenated with the surrounding terrain change feature vector Vec_Change_B001 (a six-dimensional vector) generated in step S138 to form a higher-dimensional composite vector. The concatenated vector Vec_Evidence_B001 has nine dimensions, with the first three components representing shadow anomaly statistical features and the last six components representing surrounding terrain change features. This Vec_Evidence_B001 is the constructed indirect evidence chain feature vector for whether building Bldg_ID_001 has suffered structural damage; it integrates anomaly information from the building's own shadow and change information from the surrounding environment.

[0069] Step S1310: Input the indirect evidence chain feature vector into a pre-trained suspected damage classifier. The suspected damage classifier outputs a confidence score indicating that the candidate building object belongs to a suspected damaged building. Select candidate building objects whose confidence scores exceed a preset confidence threshold as suspected damaged building objects.

[0070] Specifically, a suspected damage classifier is constructed, such as a binary classification model based on a gradient boosting decision tree (e.g., XGBoost). This classifier model has been pre-trained using a large amount of historical disaster data. The features of the training data are the aforementioned nine-dimensional indirect evidence chain feature vector, labeled "damaged" or "not damaged." For the current candidate building object Bldg_ID_001, its indirect evidence chain feature vector Vec_Evidence_B001 is input into the trained suspected damage classifier. After layers of judgment in its internal decision tree, the classifier outputs a value Conf_B001 between 0 and 1, representing the confidence score that the building belongs to the "damaged" category. A preset confidence threshold Th_conf is used, for example, 0.75. Conf_B001 is compared with Th_conf. If Conf_B001 is greater than or equal to Th_conf, then the building Bldg_ID_001 is classified as a suspected damaged building object.

[0071] Step S1311: Mark the base polygon position of each suspected damaged building object in the spatial coordinate system of the post-disaster remote sensing image data set, and attach an indirect evidence chain feature vector as an attribute label to the base polygon of each suspected damaged building object to generate a distribution map of suspected damaged buildings carrying indirect damage evidence labels.

[0072] Specifically, for all suspected damaged building objects selected through step S1310, such as Bldg_ID_001, Bldg_ID_015, and Bldg_ID_078, a new vector layer, Map_Suspect, is created. In this layer, a new vector feature is created for each suspected damaged building object, its geometry being the base polygon of the building, such as Polygon_Base_B001. Simultaneously, a new field named "Evidence_Vector" is added to the attribute table of this feature, and the indirect evidence chain feature vector Vec_Evidence_B001 (a nine-dimensional value) corresponding to the building is stored in this field in a serialized form (such as a JSON array). Thus, the generated Map_Suspect layer not only indicates the spatial location of suspected damaged buildings but also carries detailed evidence data for making this judgment; therefore, it is called a distribution map of suspected damaged buildings carrying indirect damage evidence markers.

[0073] Step S140: For each suspected damaged building object in the suspected damaged building distribution map, damage mode classification processing is performed. Based on the abnormal offset information of the shadow area corresponding to the suspected damaged building object and the layout change characteristics of the surrounding land features, the damage type identifier of each suspected damaged building object is determined and its initial damage degree score is calculated through classification and regression analysis.

[0074] In this embodiment, after identifying suspected damaged buildings, it is necessary to further distinguish their specific damage patterns (such as complete collapse, partial collapse, tilting, etc.) and severity. The specific classification and scoring process is implemented by the following sub-steps.

[0075] Step S141: Extract the base polygon of each suspected damaged building object from the distribution map of suspected damaged buildings, and based on the base polygon, cut out all shadow area abnormal offset vectors covering the base polygon and its surrounding preset buffer range from the shadow area abnormal offset information to form the local abnormal offset vector field of the suspected damaged building object.

[0076] Specifically, taking the suspected damaged building object Bldg_ID_001 as an example, its base polygon Polygon_Base_B001 is extracted from the suspected damaged building distribution map Map_Suspect. Centered on this polygon, a preset buffer distance Buf_local (e.g., 30M) is extended outwards, forming a larger rectangular or circular region Region_Local_B001. Then, returning to the list of abnormal offset information in the shaded area List_Anomaly_Vec generated in step S129, vectors whose center points fall within the region_Local_B001 are selected. Assuming P vectors are selected, they collectively constitute a vector set Field_Anomaly_Local_B001. Each vector in this set carries a spatial location (its center point coordinates) and a vector value (offset direction, offset distance), thus forming a vector field distributed within a local space, called the local abnormal offset vector field.

[0077] Step S142: Perform spatial gradient analysis on the local abnormal offset vector field, calculate the spatial change rate of the offset vector and the directional consistency coefficient of the offset vector in the local abnormal offset vector field, and generate the spatial distribution feature descriptor of the offset field of the suspected damaged building object.

[0078] Specifically, spatial gradient calculations are performed on the local anomaly offset vector field Field_Anomaly_Local_B001. First, the spatial rate of change of the offset vector Grad_Field is calculated. This can be done by meshing the region and interpolating the offset vector at each mesh point. Then, the curl and divergence of the vector field are calculated. Curl characterizes the degree of rotation of the vector field, and divergence characterizes the degree of convergence and divergence. These two indicators together reflect the drastic change of the offset in space. Second, the directional consistency coefficient Consist_Field of the offset vector is calculated. Treating the direction angles of all vectors as points on a circle, the average resultant length of these angles is calculated. The closer this value is to 1, the more consistent the directions; the closer it is to 0, the more chaotic the directions. The calculated average curl Curl_avg, average divergence Div_avg, and directional consistency coefficient Consist_Field are combined into a three-dimensional vector Vec_Grad_B001, which is the spatial distribution feature descriptor of the offset field of the suspected damaged building object.

[0079] Step S143: Obtain the original surrounding land feature layout topology map of the suspected damaged building object from the pre-disaster building surrounding land feature layout description set based on the base polygon, and extract the post-disaster surrounding land feature layout update map of the suspected damaged building object from the post-disaster remote sensing image data set.

[0080] Specifically, taking the suspected damaged building object Bldg_ID_001 as an example, the original topology map Graph_Layout_Pre_B001 corresponding to it is retrieved again from the pre-disaster surrounding feature layout topology map set constructed in step S118. Simultaneously, the post-disaster surrounding feature layout update map Map_Layout_Post_B001 for this building, generated in step S134, is directly used as input. At this point, graph structure data describing the building's surrounding environment in both pre-disaster and post-disaster states is available.

[0081] Step S144: Input the original surrounding land feature layout topology map and the post-disaster surrounding land feature layout update map into a pre-trained graph change detection network. The graph change detection network calculates the differences between the node features and edge features in the two topology maps through a graph neural network and outputs a surrounding land feature layout change pattern encoding vector.

[0082] Specifically, a graph change detection network based on graph neural networks is constructed, for example, using a graph isomorphic network as the basic encoder. The pre-disaster topology graph Graph_Layout_Pre_B001 and the post-disaster topology graph Graph_Layout_Post_B001 are respectively input into two graph isomorphic network encoders sharing weights. For the pre-disaster graph, the encoder aggregates the features of each node (buildings, roads, vegetation) and its neighboring nodes (node ​​features include feature type, area, etc., and edge features include distance, azimuth, etc.), outputting a graph-level feature vector Vec_Pre_B001. Similarly, for the post-disaster graph, Vec_Post_B001 is output. Then, a difference calculation module is designed within the network, for example, by subtracting the two vectors or concatenating them before inputting them into a multilayer perceptron to calculate the change features between the two graphs. Finally, the graph change detection network outputs a fixed-length vector Vec_Change_Graph_B001, which is the encoding vector for the surrounding land cover layout change pattern, encoding the structural change pattern of the surrounding environment from pre-disaster to post-disaster.

[0083] Step S145: The offset field spatial distribution feature descriptor and the surrounding land feature layout change pattern encoding vector are concatenated to generate the damage feature fusion vector of the suspected damaged building object.

[0084] Specifically, the offset field spatial distribution feature descriptor Vec_Grad_B001 (a 3D vector) generated in step S142 is concatenated with the surrounding land feature layout change pattern encoding vector Vec_Change_Graph_B001 (assumed to be a 128-dimensional vector) generated in step S144 to form a higher-dimensional damage feature fusion vector Vec_Fusion_B001, which has a dimension of 131. This fusion vector integrates the spatial distribution details of the building's own shadow anomalies and the structural change information of the surrounding environment.

[0085] Step S146: Input the damage feature fusion vector into a pre-trained damage pattern classifier. The damage pattern classifier contains multiple binary sub-networks, each of which corresponds to a preset damage type. Each binary sub-network outputs the probability value of the suspected damaged building object belonging to the corresponding damage type based on the damage feature fusion vector.

[0086] Specifically, a damage pattern classifier is constructed. This classifier is not a single multi-classifier, but an ensemble model composed of multiple independent binary classification sub-networks. There are Q preset damage types, such as type T1: complete collapse, T2: partial structural collapse, T3: building tilt, T4: roof damage, etc. An independent binary classifier SubNet_q is trained for each type Tq. Each sub-network can be a support vector machine or a shallow neural network. For the current building Bldg_ID_001, its damage feature fusion vector Vec_Fusion_B001 is simultaneously input into the Q binary classification sub-networks. Each sub-network SubNet_q independently performs binary classification based on the input vector, outputting a probability value Prob_q, representing the probability that the building's damage pattern belongs to type Tq. Finally, a set of probability values ​​(Prob_1, Prob_2, ..., Prob_Q) is obtained.

[0087] Step S147: Select the damage type corresponding to the maximum probability value from all binary sub-network outputs as the damage type identifier of the suspected damaged building object, and use the maximum probability value as the initial score of the damage degree of the suspected damaged building object.

[0088] Specifically, for building Bldg_ID_001, the Q probability values ​​Prob_1 to Prob_Q output in step S146 are compared, and the maximum value, Prob_max, is found. Assuming Prob_max is Prob_3, corresponding to damage type T3 (building tilt), then T3 is determined as the damage type identifier Type_ID_B001 for the suspected damaged building. Simultaneously, this maximum probability value, Prob_max, is directly used as the initial damage score Score_Init_B001 for the building. A higher initial score indicates a greater certainty from the classifier that the building belongs to that damage type, and to some extent reflects the typicality or severity of the damage.

[0089] Step S148: Based on the damage type identifier, query the spatial distribution template of the area of ​​concern corresponding to the damage type from the preset damage type parameter table. The spatial distribution template of the area of ​​concern is used to indicate the main manifestation area of ​​the corresponding type of damage within the building base area.

[0090] Specifically, a damage type parameter table, Table_Type_Param, is pre-established. This table defines a spatial distribution template for the area of ​​interest for each damage type, Tq. For example, for the "complete collapse" type T1, its template might be a polygon Mask_T1 covering the entire building base; for the "partial structural collapse" type T2, its template might be a half-area template Mask_T2 pointing to one side of the building (such as the south side); for the "building tilt" type T3, its template might be a strip-shaped area Mask_T3 with a height related to the expected shaded area. These templates can be vector polygons or fixed-size grid masks, used to indicate which area within the building base should be focused on when judging this type of damage. Based on the damage type identifier Type_ID_B001 (i.e., T3) of the building Bldg_ID_001, the corresponding area of ​​interest template Mask_Focus_B001 is retrieved from Table_Type_Param.

[0091] Step S149: Spatial superposition of the spatial distribution template of the area of ​​interest and the local anomaly offset vector field; calculation of the average amplitude and average directional dispersion of the offset vector within the coverage area of ​​the spatial distribution template of the area of ​​interest; weighted correction of the initial damage score based on the average amplitude and average directional dispersion; and generation of the corrected damage score.

[0092] Specifically, the focus area template Mask_Focus_B001 (e.g., a mask representing the south side of a building) is spatially superimposed with the local anomaly offset vector field Field_Anomaly_Local_B001 generated in step S141. Anomaly offset vectors located within the area covered by Mask_Focus_B001 are selected. For the selected vectors, their average amplitude value Mag_Avg is calculated, which is the average magnitude of these vectors. Simultaneously, the average directional dispersion Disp_Angle of these vectors is calculated, which can be represented by the standard deviation of the vector direction angle; the smaller the dispersion, the more concentrated the offset direction. Based on the characteristics of damage type T3, weighting rules are pre-defined. For example, for the "tilted" damage mode, the anomaly offset vectors should be concentrated on the corresponding side of the tilt direction, and the amplitude should be positively correlated with the tilt angle. Therefore, a correction function F_corr can be designed, which takes Mag_Avg and Disp_Angle as input and outputs a correction coefficient K_corr between 0.5 and 1.5. For example, if Mag_Avg is large and Disp_Angle is small, then K_corr is greater than 1, indicating more severe damage; conversely, K_corr is less than 1. Ultimately, the corrected damage score Score_Final_B001 for the building Bldg_ID_001 is equal to the initial score Score_Init_B001 multiplied by the correction coefficient K_corr. This correction makes the score more consistent with the physical characteristics of the specific damage pattern.

[0093] Step S1410: The damage type identifier and the corrected damage level score are used as the final damage mode classification result for each suspected damaged building object.

[0094] Specifically, for each suspected damaged building object, its damage type identifier (such as Type_ID_B001) and the damage severity score (Score_Final_B001) after correction in step S149 are combined to form a tuple (Type, Score), which serves as the final damage mode classification result for the object.

[0095] Step S150: Generate a damage assessment result dataset containing the spatial location coordinates of the damaged building and the damage level code based on the damage type identifier and the initial damage degree score.

[0096] In this embodiment, after obtaining the detailed classification results for each suspected damaged building, it is necessary to convert them into a standardized damage assessment result dataset that is easy for the emergency command system to read and display. The specific generation process is implemented by the following sub-steps.

[0097] Step S151: Obtain the base polygon vertex coordinate sequence for each suspected damaged building object, and calculate the geometric center point coordinates of the base polygon vertex coordinate sequence as the spatial position coordinates of the suspected damaged building object.

[0098] Specifically, for each suspected damaged building object, such as Bldg_ID_001, its base polygon Polygon_Base_B001 consists of a set of vertex coordinates (X1, Y1), (X2, Y2), ..., (Xn, Yn). The arithmetic mean of these vertex coordinates is calculated. That is, the longitude coordinate Lon_Cen_B001 of the center point is equal to the sum of the X coordinates of all vertices divided by the number of vertices n, and the latitude coordinate Lat_Cen_B001 of the center point is equal to the sum of the Y coordinates of all vertices divided by the number of vertices n. The calculated (Lon_Cen_B001, Lat_Cen_B001) are used as the spatial coordinates of the building object for point marking on the map.

[0099] Step S152: Based on the damage type identifier of each suspected damaged building object, query the basic damage level code corresponding to the damage type identifier from the preset damage level mapping table.

[0100] Specifically, a damage level mapping table, Table_Level_Map, is pre-constructed. This table defines the base damage level code corresponding to each damage type identifier. For example, damage type T1 (complete collapse) maps to base level code L5 (highest level), T2 (partial structural collapse) maps to L4, T3 (building tilt) maps to L3, and T4 (roof damage) maps to L2. For building Bldg_ID_001, whose damage type identifier is T3, the corresponding base damage level code Base_Level_B001 obtained by querying Table_Level_Map is L3.

[0101] Step S153: Based on the initial damage score of each suspected damaged building object, perform level fine-tuning processing on the basis of the basic damage level code, map the initial damage score to the level fine-tuning increment value, and add the basic damage level code and the level fine-tuning increment value to generate the final damage level code of the suspected damaged building object.

[0102] Specifically, the initial damage rating for building Bldg_ID_001 is Score_Init_B001 (obtained in step S147, a value between 0 and 1). A linear mapping function is defined to map the initial rating interval [0, 1] to the fine-tuning increment interval [0, 2]. For example, the fine-tuning increment value Delta_Level_B001 is equal to the product of Score_Init_B001 and 2. Then, the base level code Base_Level_B001 is treated as a single value, for example, L3 is treated as the value 3. The final damage level code Final_Level_B001 is equal to the base level value plus the fine-tuning increment value. If Base_Level_B001 is 3 and Score_Init_B001 is 0.8, then Delta_Level is 1.6 and Final_Level is 4.6. This encoding method allows for more detailed differentiation based on the degree of damage under the same basic damage type. The final damage level code can be a floating-point number, or it can be rounded and combined with a subdivision suffix, such as L4-6.

[0103] Step S154: Extract the indirect evidence chain feature vector of each suspected damaged building object from the distribution map of suspected damaged buildings, perform hash encoding on the indirect evidence chain feature vector, and generate the evidence tracing index code associated with the suspected damaged building object.

[0104] Specifically, for building Bldg_ID_001, the indirect evidence chain feature vector Vec_Evidence_B001 (a nine-dimensional floating-point array) is read from the attribute table of the suspected damaged building distribution map Map_Suspect. This array is serialized into a string. Then, a cryptographic hash function, such as SHA-256, is applied to this string to generate a fixed-length hash value, such as "3f7e2c8a...". This hash value is used as the evidence tracing index code Hash_Evidence_B001 for the building object. This index code can be used to subsequently trace the original evidence features that generated the evaluation result without storing the complete feature vector, ensuring data verifiability and simplicity.

[0105] Step S155: Combine the spatial location coordinates, the final damage level code, and the evidence tracing index code to generate a damage assessment record unit for the suspected damaged building object.

[0106] Specifically, the spatial coordinates (Lon_Cen_B001, Lat_Cen_B001) of the building Bldg_ID_001, the final damage level code Final_Level_B001 (e.g., L4-6), and the evidence tracing index code Hash_Evidence_B001 are combined into a structured data record unit, Record_B001. This record unit can be represented in JSON format, for example: {"BuildingID": "B001", "Location": {"Lon": Lon_Cen_B001, "Lat": Lat_Cen_B001}, "DamageLevel": "L4-6", "EvidenceHash": "3f7e2c8a..."}. This record unit contains the core information of the evaluation results and is linked to the original chain of evidence through the hash code.

[0107] Step S156: Sort all damage assessment record units of suspected damaged buildings in space according to their spatial location coordinates to generate a damage assessment record sequence organized by spatial grid.

[0108] Specifically, damage assessment record units of all suspected damaged building objects are collected, for example, there are N such units. To improve the efficiency of subsequent spatial query and analysis, the record units are spatially sorted according to their spatial location coordinates (Lon, Lat). For example, a Hilbert curve or geohashing algorithm is used to calculate a spatial sorting code based on the coordinates, and then all record units are sorted according to the sorting code. After sorting, an ordered sequence of damage assessment records, List_Record_Sorted, is obtained. This sequence can be easily stored and indexed in blocks according to a spatial grid.

[0109] Step S157: Perform redundancy elimination processing on the damage assessment record sequence, check whether there are multiple damage assessment record units corresponding to the same building object. If so, retain the damage assessment record unit with the highest damage degree score after correction and delete the remaining redundant records to generate a set of redundancy-free damage assessment records.

[0110] Specifically, in the initially generated damage assessment record sequence List_Record_Sorted, algorithm errors or data overlap may cause the same physical building to be identified as two slightly spatially offset suspected objects, resulting in two records. Therefore, redundancy removal is necessary. The sorted record sequence is traversed, and for each record Record_A, its spatial coordinates are checked to see if they are very close to the spatial coordinates of the subsequent record Record_B, for example, if the Euclidean distance is less than a threshold Th_dist_dup (e.g., 5 meters). If the distance is less than the threshold, they are considered to likely point to the same building. At this point, the implicit corrected damage severity scores in the two records are compared (which can be calculated from the final damage level code or stored separately). The record with the higher score (i.e., more severe damage) is retained, and the other redundant record is deleted. After processing all records, the deredundant damage assessment record set Set_Record_Deduped is obtained.

[0111] Step S158: The set of damage assessment records after redundancy removal is encapsulated into a standard geographic information data format, and a data production timestamp and data version identifier are added to each damage assessment record to generate a damage assessment result dataset containing the spatial location coordinates of the damaged buildings and the damage level code.

[0112] Specifically, the deduplicated damage assessment record set, `Set_Record_Deduped`, is encapsulated into a standard geographic information data format, such as GeoJSON or Shapefile. During encapsulation, metadata is added to the entire dataset, including a data production timestamp (`Time_Stamp_Prod`, e.g., the time the processing was completed) and a data version identifier (`Version_ID`, e.g., "V1.0"). Simultaneously, it is ensured that the core fields of each record—spatial coordinates and the final damage level code—are correctly mapped to the corresponding fields in the output format. The resulting file or data stream is a damage assessment result dataset containing the spatial coordinates and damage level codes of the damaged buildings. This dataset can be distributed to emergency command centers, rescue teams, and disaster reconstruction planning departments.

[0113] Step S210: After generating a distribution map of suspected damaged buildings carrying indirect damage evidence markers, the method performs spatial clustering analysis on the suspected damaged building objects in the distribution map to generate a distribution map of damaged building cluster areas, including: In this embodiment, after generating suspected damage information for individual buildings, the spatial distribution patterns of damaged buildings are further analyzed from a macroscopic perspective to identify severely damaged clusters. This is of great value for coordinating rescue efforts and planning the deployment of supplies. The specific cluster analysis process is implemented through the following sub-steps.

[0114] Step S211: Extract the geometric center point of the base polygon of each suspected damaged building object in the distribution map of suspected damaged buildings as the spatial point representation of the suspected damaged building object.

[0115] Specifically, iterate through all suspected damaged building objects in the suspected damaged building distribution map Map_Suspect. For each object, such as Bldg_ID_001, calculate the coordinates (Lon_Cen_B001, Lat_Cen_B001) of its base polygon Polygon_Base_B001, using the same method as in step S151. Use a point feature Point_B001 containing these coordinates to represent the entire building for subsequent cluster analysis. All representative points form a point set Set_Points_Suspect.

[0116] Step S212: Perform density-based spatial clustering on all spatial point representations, using a preset spatial search radius as the neighborhood range parameter and a preset minimum number of neighborhood points as the density threshold parameter, to identify densely distributed point clusters as preliminary damaged aggregation areas.

[0117] Specifically, a density-based spatial clustering algorithm, such as the DBSCAN algorithm, is applied to the point set Set_Points_Suspect. Two key parameters are set: the spatial search radius Eps (e.g., 100M) and the minimum number of neighboring points MinPts (e.g., 5). After the algorithm starts, for each point, such as Point_B001, the number of other points contained within a circle centered at that point and with a radius of Eps is counted. If this number is greater than or equal to MinPts, the point is marked as a core point. Then, based on density connectivity, all mutually reachable core points and their neighboring boundary points are aggregated into a cluster Cluster_K. Each identified cluster Cluster_K corresponds to a densely distributed area of ​​suspected damaged buildings and is marked as the preliminary damage cluster Prelim_Cluster_K.

[0118] Step S213: Construct the convex hull polygon for each initial damage cluster area, and use the spatial range covered by the convex hull polygon as the initial boundary of the damage cluster area.

[0119] Specifically, for each preliminary damage cluster region Prelim_Cluster_K, the coordinates of representative points of all suspected damaged building objects within that region are obtained. Using a convex hull algorithm in computational geometry (such as the Graham scan method), these points are calculated to obtain a minimum convex polygon Convex_Hull_K that can contain all points. This convex polygon Convex_Hull_K is defined as the initial boundary of the damage cluster region.

[0120] Step S214: Obtain the set of damage type identifiers for all suspected damaged building objects within each initial damage cluster area, count the frequency of each damage type identifier appearing within the cluster area, and take the damage type identifier with the highest frequency as the dominant damage type identifier for that damage cluster area.

[0121] Specifically, for each Prelim_Cluster_K of preliminary damage clusters, all suspected damaged building objects contained within it are traversed. The damage type identifier (e.g., T1, T2, T3, etc.) determined in step S147 is read from the attributes of each object. A type frequency statistics table, such as a dictionary Dict_Type_Count, is established. During the traversal, each time a damage type Tq is encountered, the count value corresponding to Tq in Dict_Type_Count is incremented by 1. After traversing all objects within the region, the key with the largest count value in Dict_Type_Count is found; its corresponding damage type, such as T2, is the dominant damage type identifier Dom_Type_K for that damage cluster.

[0122] Step S215: Obtain the initial damage score set of all suspected damaged buildings within each initial damage cluster area, and calculate the average value of the initial damage score set as the regional average damage score of the damage cluster area.

[0123] Specifically, for each Prelim_Cluster_K of preliminary damage clusters, all suspected damaged building objects contained within it are traversed again. The initial damage score Score_Init, determined in step S147, is read from the attributes of each object. All scores are collected to form a list of values ​​List_Score_K. The arithmetic mean Avg_Score_K of all scores in List_Score_K is calculated. This average Avg_Score_K is the regional average damage score for the damage cluster, reflecting the overall average damage severity of the region.

[0124] Step S216: Based on the dominant damage type identifier and the regional average damage score, smooth and expand the initial boundary of the preliminary damage cluster area to generate a distribution map of the damaged building cluster area containing the dominant damage type identifier attribute label and the regional average damage score attribute label.

[0125] Specifically, for each initial damage cluster region Prelim_Cluster_K, its initial boundary is a convex hull polygon Convex_Hull_K. To better reflect the actual geographical distribution and visual perception, smoothing and expansion are performed. Smoothing can be achieved using the Douglas-Puk algorithm or Bézier curve fitting to remove sharp corners on the convex hull boundary, making it smoother. Expansion involves extending the smoothed polygon outward by a small buffer distance (e.g., 20M) to include potential impact areas that may be located at the edge of the region but not fully covered by the convex hull. After smoothing and expansion, the final region boundary Polygon_Final_K is obtained. Then, on a new vector layer Map_Cluster, features are created for the polygon Polygon_Final_K, and attribute fields are added: "Dominant Damage Type" is assigned the value Dom_Type_K, and "Average Damage Score" is assigned the value Avg_Score_K. Perform the above operations on all initially damaged clusters, and the resulting Map_Cluster is a distribution map of damaged building clusters that includes the dominant damage type identifier attribute and the regional average damage score attribute.

[0126] Step S217: Perform spatial overlay analysis on the preliminary damage cluster area, and generate a damage impact area function assessment report by combining the regional functional attribute data in the pre-disaster ground structure knowledge base.

[0127] Specifically, the initial boundary of the preliminary damaged cluster area generated in step S213 (or the final boundary obtained in step S216) is spatially overlaid with the regional functional attribute data stored in the pre-disaster ground feature structure knowledge base to assess the impact of the disaster on different functional zones. The detailed process is implemented by the following sub-steps.

[0128] Step S2171: Retrieve the regional functional attribute data partition covering the spatial range of the initially damaged cluster area from the pre-disaster ground structure knowledge base. The regional functional attribute data partition includes land use classification code, population density distribution raster layer, transportation network importance level vector layer, and lifeline pipeline distribution vector layer.

[0129] Specifically, for each preliminarily damaged cluster area Prelim_Cluster_K, its initial boundary polygon Convex_Hull_K is obtained. Using this polygon as the spatial query range, a request is sent to the pre-disaster land cover knowledge base. The pre-disaster land cover knowledge base returns data partitions covering the area. These partitions contain various functional attribute data: a vector layer of land use classification codes, where each land parcel has a use code LU_Code (e.g., R2 represents Class II residential land, B1 represents commercial land); a population density distribution raster layer, where the value of each cell represents the population density Pop_Dens at that location; a transportation network importance level vector layer, where each road feature has an importance level Road_Level (e.g., highway, arterial road, secondary arterial road); and a lifeline pipeline distribution vector layer, including water supply, power supply, and gas pipelines, where each pipeline has attributes such as diameter Dia and material Mat.

[0130] Step S2172: Perform spatial overlay analysis between the initial boundary of the preliminary damaged cluster area and the land use classification code, calculate the area proportion corresponding to different land use classification codes within the preliminary damaged cluster area, and generate a statistical table of land use type composition of the damaged area.

[0131] Specifically, a spatial overlay (Intersect) analysis is performed between the polygon Convex_Hull_K and the land use classification code layer. The area of ​​the intersection between the polygon Convex_Hull_K and each land parcel polygon is calculated. The areas are grouped according to the land use code LU_Code, and the total intersection area Area_LU corresponding to each LU_Code is summarized. Then, the percentage of the area of ​​each LU_Code to the total area (area of ​​the polygon Convex_Hull_K) is calculated. Finally, a statistical table Table_LU_Composition is generated, listing the area and percentage of various land uses within the damaged cluster area, such as "Residential Land: 60%", "Commercial Land: 30%", and "Green Space: 10%".

[0132] Step S2173: Spatial overlay analysis is performed between the initial boundary of the initially damaged cluster area and the population density distribution raster layer, and the sum of the population density raster values ​​within the initially damaged cluster area is used as an estimate of the affected population size.

[0133] Specifically, the polygon Convex_Hull_K is used as a mask to clip and perform zonal statistics on the population density distribution raster layer. The values ​​of all raster cells within the polygon, Pop_Dens_i, are extracted. Since each cell represents the population per unit area (e.g., 100 square meters), the sum of all these cell values, Total_Pop_K, is the estimated size of the affected population within the damaged cluster area. This estimate, based on pre-disaster resident population data, can be used to estimate the number of people requiring relocation and resettlement.

[0134] Step S2174: Perform spatial overlay analysis on the initial boundary of the preliminary damage cluster area and the traffic network importance level vector layer to identify the traffic network elements at all levels contained in the preliminary damage cluster area, and extract the importance level identifier and element length of each included traffic network element to generate a list of affected traffic networks.

[0135] Specifically, the polygon Convex_Hull_K is spatially selected (Select by Location) or spatially connected to the traffic network importance level vector layer to identify all traffic network elements (road segments) that intersect with or are entirely within the polygon. For each selected road segment, its importance level identifier Road_Level and the length of the portion of the segment inside the polygon Length_In_Region are extracted (if the segment is only partially inside the polygon, its clipped length is calculated). All information is then compiled to generate an Affected Traffic Network List List_Affected_Roads, containing information such as road name (if any), level, and affected length, used to assess the accessibility of emergency lanes and repair priorities.

[0136] Step S2175: Perform spatial overlay analysis on the initial boundary of the preliminary damage cluster area and the lifeline pipeline distribution vector map layer to identify various lifeline pipeline elements contained in the preliminary damage cluster area, and extract the pipe diameter parameters, material parameters and start and end point coordinates of each included lifeline pipeline element to generate a list of affected lifeline pipelines.

[0137] Specifically, a spatial overlay analysis is performed between the polygon Convex_Hull_K and the lifeline pipeline distribution vector layer. All water, electricity, and gas pipeline elements intersecting with the polygon are identified. For each selected pipeline, its type (water / electricity / gas), diameter parameter Dia, material parameter Mat, and the coordinates of the start and end points (Start_Point, End_Point) within the polygon are extracted. All information is then compiled to generate an affected lifeline pipeline list List_Affected_Utilities. This list is crucial for organizing emergency repairs of water, electricity, and gas infrastructure.

[0138] Step S2176: The estimated affected population size is weighted and corrected based on the dominant damage type identifier and the regional average damage score to generate a corrected assessment value of the affected population size.

[0139] Specifically, considering the varying substantial impacts of different damage types and degrees on the population—for example, all residents in a completely collapsed area may require resettlement, while residents in an area with only roof damage may only require partial resettlement—a correction coefficient function F_pop_adj is defined. This function takes the dominant damage type identifier Dom_Type_K and the area's average damage score Avg_Score_K as input and outputs a correction coefficient K_pop_adj between 0 and 1. For example, if the dominant type is complete collapse, K_pop_adj might be 0.95; if it's roof damage, K_pop_adj might be 0.3. Then, the total population estimated in step S2173, Total_Pop_K, is multiplied by K_pop_adj to obtain the corrected assessment value of the affected population size: Adj_Pop_K = Total_Pop_K. K_pop_adj. This revision makes the population assessment results closer to the actual disaster situation.

[0140] Step S2177: Combine the statistical table of land use types in the damaged area, the revised assessment value of the affected population size, the list of affected transportation networks, and the list of affected lifeline pipelines to generate a functional assessment report of the damaged area.

[0141] Specifically, all the above analysis results are integrated into a structured document or data package, namely the Damage Impact Area Function Assessment Report_Function_K. This report includes the following parts: area identifier (corresponding to Prelim_Cluster_K), land use type composition statistics table Table_LU_Composition, revised affected population size assessment value Adj_Pop_K, list of affected transportation networks List_Affected_Roads, and list of affected lifeline utilities List_Affected_Utilities.

[0142] Step S310, before classifying the suspected damaged building objects according to their damage patterns, the method performs image quality enhancement processing on the post-disaster remote sensing image dataset to generate an enhanced post-disaster remote sensing image dataset, including: In this embodiment, the original post-disaster remote sensing images may be of poor quality due to weather (such as clouds and fog), sensor noise, or resolution limitations, and directly using them for analysis would affect the accuracy of subsequent steps. Therefore, before classifying damage patterns, the images are first subjected to a series of enhancement processes. The specific enhancement process is implemented by the following sub-steps.

[0143] For example, in step S311, atmospheric scattering effect suppression processing is performed on each post-disaster remote sensing image unit in the post-disaster remote sensing image data set. The atmospheric light value and atmospheric transmittance distribution map of each post-disaster remote sensing image unit are estimated by the dark channel prior algorithm. Atmospheric correction is performed on each post-disaster remote sensing image unit according to the atmospheric light value and the atmospheric transmittance distribution map to generate a defogging post-disaster remote sensing image unit.

[0144] Specifically, taking the post-disaster remote sensing image unit ImgOpt_Post_01 as an example, the dark channel prior dehazing algorithm is applied. First, for each pixel in the image, within its local neighborhood (e.g., 15... A 15-window image (Dark Channel) is obtained by calculating the minimum values ​​of three bands (R, G, B). Then, the top 0.1% of the brightest pixels in the Dark Channel image are selected, corresponding to the most opaque areas of the image (i.e., dense fog / haze areas). The maximum intensity value of these pixels in the original image is taken as the global atmospheric light value A_global. Next, the atmospheric transmittance t(x) of each pixel is estimated according to the atmospheric scattering model. The estimation formula is based on the Dark Channel prior, that is, t(x) equals 1 minus the quotient of omega multiplied by the Dark Channel value divided by A_global, where omega is a constant (usually taken as 0.95). By calculating pixel by pixel, an atmospheric transmittance distribution map Map_T is obtained. Finally, the original image is corrected using the inverse process of the atmospheric scattering model: for each pixel, its pixel value is subtracted from A_global, divided by the maximum value of t(x) (to prevent division by zero), and then added back to A_global to obtain the dehazed pixel value. After all pixels have been processed, the dehazed post-disaster remote sensing image unit Img_Dehazed_01 is generated.

[0145] Step S312: Perform sensor noise suppression processing on the defogging post-disaster remote sensing image unit, transform the defogging post-disaster remote sensing image unit to the frequency domain, identify and attenuate high-frequency noise components in the frequency domain, and then transform the attenuated frequency domain data back to the spatial domain to generate a denoised post-disaster remote sensing image unit.

[0146] Specifically, noise suppression is performed on the dehazed image unit Img_Dehazed_01. First, a two-dimensional Fast Fourier Transform (FFT) is performed on the image to transform it from the spatial domain to the frequency domain, resulting in a spectrogram Freq_Map. In the frequency domain, image information is mainly concentrated in the low- and mid-frequency regions, while random noise (such as salt-and-pepper noise and Gaussian noise) typically manifests as high-frequency components. Analyzing the Freq_Map reveals the high-frequency regions representing noise. A low-pass filter, such as a Butterworth low-pass filter, is designed, with its transfer function smoothly transitioning at the cutoff frequency D0. The Freq_Map is then multiplied element-wise with the filter function to attenuate high-frequency noise components, resulting in a filtered spectrogram Freq_Filtered. Next, an inverse Fourier transform is performed on Freq_Filtered to transform the data back from the frequency domain to the spatial domain, yielding the denoised post-disaster remote sensing image unit Img_Denoised_01. This step effectively suppresses graininess and random noise in the image.

[0147] Step S313: Perform spatial resolution enhancement processing on the denoised post-disaster remote sensing image unit, input the denoised post-disaster remote sensing image unit into a pre-trained super-resolution reconstruction network, and use the super-resolution reconstruction network to jointly model the low-frequency structural information and high-frequency detail information of the image, and output a high-resolution post-disaster remote sensing image unit with enhanced spatial resolution.

[0148] Specifically, to obtain clearer edges and details of ground features, super-resolution reconstruction is performed on the denoised image unit Img_Denoised_01. This is then fed into a pre-trained super-resolution reconstruction network, such as a model based on an enhanced deep residual network. The network takes a low-resolution image (e.g., the original resolution) as input and outputs a high-resolution image magnified by a certain factor (e.g., 2x or 4x). Internally, the network learns an end-to-end mapping from low to high resolution through multiple residual blocks. During training, the network is trained to minimize pixel-level differences (e.g., L1 loss) and perceptual loss between the reconstructed high-resolution image and the real high-resolution image. In this step, Img_Denoised_01 propagates forward layer by layer in the network. The bottom convolutional layers extract shallow features, the intermediate residual block groups learn high-frequency residual information, and the final upsampling module (e.g., sub-pixel convolutional layers) magnifies and reassembles the feature map to generate the final high-resolution output. This output is the high-resolution post-disaster remote sensing image unit Img_SR_01 with enhanced spatial resolution.

[0149] Step S314: Perform radiometric normalization processing on the high-resolution post-disaster remote sensing image unit, statistically analyze the pixel value histogram distribution of the high-resolution post-disaster remote sensing image unit, map the pixel value histogram distribution to the range of the standard radiometric reference distribution, and generate a radiometric normalized high-resolution post-disaster remote sensing image unit.

[0150] Specifically, to ensure consistent radiometric characteristics across images acquired at different times and from different sensors, facilitating subsequent quantitative analysis and comparison, radiometric normalization was performed on the super-resolution reconstructed image Img_SR_01. First, the grayscale value distribution of all pixels in Img_SR_01 across each band was statistically analyzed to obtain its histogram, Hist_SR. A standard radiometric reference distribution was pre-defined, such as the standard histogram Hist_Ref obtained from a large number of cloudless images. Then, a histogram matching algorithm was applied. For each grayscale level in Hist_SR, the algorithm found the grayscale level in Hist_Ref whose cumulative distribution function value was closest, constructing a grayscale mapping table, Map_Gray. Each pixel in the image was traversed, and based on its original grayscale value, it was searched through the mapping table Map_Gray and replaced with a new grayscale value. After processing, a new image was obtained, whose pixel value histogram distribution was basically consistent with the standard reference distribution Hist_Ref, i.e., the radiometrically normalized high-resolution post-disaster remote sensing image unit Img_Norm_01.

[0151] Step S315: The high-resolution post-disaster remote sensing image unit after radiometric normalization is subjected to edge sharpening processing. The high-frequency edge information of the high-resolution post-disaster remote sensing image unit after radiometric normalization is extracted by the Laplacian operator. The high-frequency edge information is superimposed back onto the original image to enhance the edge contrast, thereby generating the high-resolution post-disaster remote sensing image unit with enhanced edges.

[0152] Specifically, to highlight the edges of features such as buildings and roads, making them easier to identify and extract, edge sharpening is performed on the radiometrically normalized image Img_Norm_01. First, a 3D model is defined. A Laplacian kernel of 3, for example, a kernel with a center value of 8 and surrounding values ​​of -1, is used to convolve the image Img_Norm_01, resulting in a Laplacian image Img_Lap. This image highlights high-frequency edges and details in the original image. The Laplacian image is then superimposed back onto the original image with certain weights. The sharpened image Img_Sharp is equal to the original image Img_Norm_01 plus the Laplacian image Img_Lap multiplied by the sharpening factor Alpha (e.g., Alpha = 0.5). This operation enhances contrast at edges, making ground feature outlines clearer. The processed image is the high-resolution post-disaster remote sensing image unit Img_Sharp_01 with enhanced edges.

[0153] Step S316: All the processed edge-enhanced high-resolution post-disaster remote sensing image units are stitched together according to their original spatial locations to generate an enhanced post-disaster remote sensing image data set covering the entire disaster-affected area.

[0154] Specifically, for each image unit (such as ImgOpt_Post_01, ImgSar_Post_02, etc.) in the post-disaster remote sensing image dataset, steps S311 to S315 are sequentially performed to obtain their corresponding edge-enhanced high-resolution image units (such as Img_Sharp_01, Img_Sharp_02). Then, based on the original spatial coverage boundary coordinates of each image unit, image stitching technology is used to seamlessly stitch the processed image units together. The stitching process may involve image registration, color equalization, and feathering to eliminate stitching seams. Finally, an enhanced post-disaster remote sensing image dataset covering the entire disaster-affected area with comprehensively improved quality is generated. This dataset will serve as the basic input data for all subsequent analysis steps (such as shadow extraction, ground feature identification, etc.).

[0155] Step S410: After generating the damage assessment result dataset, the method performs emergency response decision support processing on the damage assessment result dataset to generate a dynamic deployment instruction sequence for different emergency response stages, including: In this embodiment, after generating the final damage assessment result dataset, it is further applied to emergency response decision support, transforming the assessment results into executable resource allocation and action instructions to guide rescue work at different stages. The specific decision support process is implemented through the following sub-steps.

[0156] For example, in step S411, the spatial location coordinates and damage level codes of the damaged buildings in the damage assessment result dataset are input into the emergency resource demand prediction model. The emergency resource demand prediction model outputs a list of estimated demand for various emergency resources based on the number of damaged buildings and the distribution of damage levels.

[0157] Specifically, from the damage assessment result dataset Set_Record_Deduped, the total number of damage records N_total and the distribution of the number of buildings (N_L1, N_L2, ...) for different damage levels (e.g., L1-L5) are calculated. These statistical data are then used as feature vectors and input into a pre-built emergency resource demand prediction model. This model can be a regression model or a neural network trained on historical disaster data. Internally, the model, through learned mapping relationships, outputs a multi-dimensional vector based on the input damage quantity and level distribution, namely, a list of estimated demand quantities for various emergency resources List_Resource_Demand. This list includes, for example: the number of rescue personnel Demand_Staff, the number of engineering machinery Demand_Mach, the number of first aid kits Demand_MedKit, the number of tents Demand_Tent, and the amount of drinking water Demand_Water, etc.

[0158] Step S412: Generate a heat map of the spatial distribution of damaged points based on the spatial coordinates of the damaged buildings in the damage assessment result dataset. Overlay the heat map of the spatial distribution of damaged points with a pre-stored spatial distribution layer of emergency resource reserve points to calculate the average spatial distance and accessibility path time between each emergency resource reserve point and the surrounding damaged points.

[0159] Specifically, the spatial coordinates of all damaged buildings are used as point features, and kernel density analysis is applied to generate a heatmap (Heatmap_Damage) showing the spatial distribution of damaged points. This map uses color depth to visually represent the clustering degree of damaged points. Then, a pre-stored emergency resource reserve point spatial distribution layer (Layer_Depot) is loaded, containing the locations and attributes of all material reserves, rescue team outposts, etc., throughout the city. The Heatmap_Damage (or the original damaged points) is overlaid with Layer_Depot for analysis. For each reserve point (Depot_M), all damaged points within a certain radius of it are searched. The average Euclidean distance from these damaged points to Depot_M is calculated as the average spatial distance (Dist_Avg_M). Simultaneously, the average time required for the optimal path from Depot_M to each damaged point is calculated using a route planning service or based on road network data, and this average time is used as the reachability path time (Time_Avg_M).

[0160] Step S413: Assess the spatial proximity and temporal accessibility of each emergency resource reserve point based on the average spatial distance and the accessibility path time, and sort the resource allocation priorities by combining the assessment results to generate an emergency resource allocation order list sorted by priority.

[0161] Specifically, for each reserve point Depot_M, its average spatial distance Dist_Avg_M and reachability path time Time_Avg_M are used as two indicators. First, these two indicators are normalized, for example, using min-max normalization to map them to the [0,1] interval, with smaller values ​​indicating better performance (closer distance, shorter time). Then, a comprehensive priority score Score_Prio_M is defined, which is a weighted sum of the normalized spatial proximity score and the normalized time reachability score, for example, with weights of 0.3 and 0.7 respectively, indicating a greater emphasis on time efficiency. After calculating the comprehensive score for all reserve points, all reserve points are sorted from highest to lowest score (higher score, higher priority, meaning faster arrival at the disaster area), generating a priority-sorted emergency resource allocation order list List_Depot_Prio. This list indicates which reserve points should be allocated resources from first.

[0162] Step S414: Extract the list of affected transportation networks from the functional assessment report of the damaged area, and calculate the capacity reduction coefficient of each transportation route based on the importance level identifier and length of the transportation network elements in the list of affected transportation networks.

[0163] Specifically, from the damage impact area functional assessment report (e.g., Report_Function_K) generated in step S217, the list of affected traffic networks, List_Affected_Roads, is extracted. This list contains the grade and affected length of each road within the damaged area. For each potential transport route Route_R (the path connecting the reserve point and the damaged point), the affected roads along its path are analyzed. Based on the road's importance grade (e.g., highway, arterial road), a base loss factor, Base_Factor, is preset (e.g., 0.2 for highways, 0.4 for arterial roads). For Route_R, its capacity loss coefficient, Coeff_Impair_R, is equal to 1 minus the sum of the products of the base loss factor of all affected road segments on the path multiplied by (the affected length of the road segment divided by the total length of the road segment). The closer the loss coefficient is to 0, the more severe the capacity loss.

[0164] Step S415: Dynamically adjust the allocation routes in the emergency resource allocation order list according to the traffic capacity loss coefficient, allocate transportation routes with lower loss coefficients to allocation tasks with higher priority, and generate an optimized emergency resource allocation route planning scheme.

[0165] Specifically, the route planning is dynamically adjusted by combining the priority list of reserve points (List_Depot_Prio) from step S413 and the route loss coefficient from step S414. For the highest priority allocation task (i.e., transporting supplies from the highest priority reserve point to the most severely damaged area), multiple possible routes from the reserve point to the target damaged area are first calculated, and the capacity loss coefficient of each route is calculated using the method in step S414. The route with the highest loss coefficient (i.e., the least capacity loss) is selected as the recommended route. For tasks with lower priority, when planning routes, it is necessary to consider whether the planned routes will cause congestion or resource conflicts, and to coordinate and optimize them if necessary. Finally, a complete emergency resource allocation route planning scheme (Plan_Routes) that considers priority and road conditions is generated. This scheme specifies the starting reserve point, target area, specific transportation route, and estimated time for each allocation task.

[0166] Step S416: Extract the affected population size assessment value from the functional assessment report of the damaged area, and calculate the site area requirement and daily consumption of living materials for the temporary resettlement site based on the affected population size assessment value.

[0167] Specifically, the revised assessment value of the affected population size, Adj_Pop_K, is extracted from the functional assessment report of the damaged area (e.g., Report_Function_K). Based on the per capita resettlement area standard (e.g., 5 square meters per person), the required site area for temporary resettlement sites, Req_Area_K, is calculated, which is Adj_Pop_K multiplied by 5. Simultaneously, based on the per capita consumption standard for daily necessities (e.g., 2 liters of drinking water and 1 kilogram of food per person per day), the daily consumption of various daily necessities is calculated. For example, the water requirement, Req_Water_Daily_K, equals Adj_Pop_K multiplied by 2 liters, and the food requirement, Req_Food_Daily_K, equals Adj_Pop_K multiplied by 1 kilogram. These requirements will guide the planning of resettlement sites and the continuous supply of materials.

[0168] Step S417: Input the site area requirement and daily consumption of living materials of the temporary resettlement site into the resettlement site site selection optimization model, and output the recommended temporary resettlement site location coordinate sequence by combining the heat map of the spatial distribution of damaged points and the spatial distribution map of available open space.

[0169] Specifically, a site selection optimization model for temporary resettlement sites is constructed. The model's inputs include: the site area requirement (Req_Area_K) and daily water consumption (Req_Water_Daily_K) calculated in step S416, a heatmap (Heatmap_Damage) showing the spatial distribution of damaged sites (indicating the demand center), and a pre-prepared map (Map_OpenSpace) showing the distribution of available open space (including parks, squares, stadiums, etc., suitable for temporary resettlement sites, each with an area attribute). The model's optimization objective is to select a set of locations that are as close as possible to the center of the damage heatmap (i.e., the demand center), while ensuring the total area requirement is met, and that the sites themselves have suitable conditions. The model can employ a multi-objective optimization algorithm or a clustering method to select the center point. Finally, the model outputs a set of recommended temporary resettlement site location coordinates (List_Shelter_Sites), for example, containing three recommended locations, each with coordinates (Lon, Lat) and an estimated capacity.

[0170] Step S418: Generate a set of rapid response instructions for the initial stage of emergency response based on the emergency resource allocation route planning scheme and the recommended temporary resettlement site location coordinate sequence. The set of rapid response instructions includes resource departure time instructions, transportation route instructions, and resettlement site activation instructions.

[0171] Specifically, the emergency resource allocation route planning scheme Plan_Routes generated in step S415 and the recommended temporary resettlement site location coordinate sequence List_Shelter_Sites generated in step S417 are converted into executable operation instructions. For each allocation task, a "resource departure time instruction" is generated, such as "Instruction: Depart from the Chengbei reserve warehouse at T0+1 hour to allocate 5 tons of drinking water"; at the same time, a corresponding "transportation route instruction" is generated, attaching the specific route coordinate sequence planned in step S415. For each recommended resettlement site, a "resettlement site activation instruction" is generated, such as "Instruction: Activate the Chengxi Stadium as a temporary resettlement site, which can accommodate 2000 people". All instructions are organized according to time sequence (such as "execute immediately", "execute after 1 hour") and type to form a rapid response instruction set Set_Cmd_Initial for the initial stage of emergency response.

[0172] Step S419: Based on the damage type identifiers of the damaged buildings in the damage assessment result dataset, calculate the proportion of each type of damage, and generate a list of post-disaster reconstruction material requirements by combining the building structural parameters in the pre-disaster ground structure knowledge base.

[0173] Specifically, during the mid-stage of the emergency response (reconstruction planning phase), the frequency and proportion of various damage type identifiers (T1, T2, T3...) are statistically analyzed from the damage assessment dataset Set_Record_Deduped. Then, for each damage type, typical structural parameters of that type of building are retrieved from the pre-disaster building structure knowledge base, such as average building area, main building materials (brick-concrete, frame), and floor height. Combining the proportion of damage types and typical structural parameters, a list of various materials needed for post-disaster reconstruction, List_Material_Demand, is estimated. For example, if the statistics show 100 T1 type (completely collapsed) buildings, with an average requirement of S_ton of steel and C_m3 of concrete per building, then the total steel requirement is 100. S_ton, total concrete requirement is 100 C_m3. This list is used to guide the procurement and stockpiling of materials for post-disaster reconstruction.

[0174] Step S4110: Generate a set of reconstruction material allocation instructions for the mid-term of emergency response based on the post-disaster reconstruction material demand list and the spatial distribution layer of emergency resource reserve points. The set of reconstruction material allocation instructions includes material category instructions, allocation quantity instructions, and receiving point instructions.

[0175] Specifically, the post-disaster reconstruction material demand list List_Material_Demand generated in step S419 is used as the target, combined with the spatial distribution layer Layer_Depot of emergency resource reserve points (at which point the reserve points may have already been replenished with building materials), to conduct mid-term material allocation planning. Unlike the initial rapid response, mid-term allocation places more emphasis on the systematic and balanced nature of the plan. The planning algorithm will formulate detailed allocation plans for each type of material (such as steel and cement) based on factors such as the inventory of each reserve point, its distance from the reconstruction area, and transportation capacity. For each allocation task, a "material category instruction" (such as "allocate rebar"), a "allocation quantity instruction" (such as "500 tons"), and a "receiving point instruction" (specifying one or more reconstruction material distribution centers or specific construction site locations) are generated. All instructions are summarized to form a set of reconstruction material allocation instructions Set_Cmd_MidTerm for the mid-term of the emergency response. This set of reconstruction material allocation instructions, together with the initial rapid response instruction set Set_Cmd_Initial, constitutes a dynamic allocation instruction sequence for different emergency response stages, providing comprehensive guidance for all aspects of work from emergency rescue to post-disaster reconstruction.

[0176] Based on the same inventive concept, please refer to Figure 2 The diagram shows a schematic block diagram of a rapid remote sensing image change detection and damage assessment system 100 for disaster emergency response provided in an embodiment of this application. It may include a communication unit 110, a machine-readable storage medium 120, and a processor 130.

[0177] In this embodiment, the machine-readable storage medium 120 can also be integrated into the processor 130 and can communicate and interact with external systems through the communication unit 110. The machine-readable storage medium 120 stores machine-executable instructions for executing the scheme of this application, and the processor 130 executes the machine-executable instructions stored in the machine-readable storage medium 120 to implement the inspection video stream processing method provided in the aforementioned method embodiments.

[0178] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A rapid method for detecting changes and assessing damage in remote sensing images for disaster emergency response, characterized in that, The method includes: Acquire a set of post-disaster remote sensing image data covering the disaster-stricken area, and retrieve a set of topological descriptions of the shadow area of ​​pre-disaster buildings and a set of descriptions of the layout of pre-disaster buildings from a pre-constructed knowledge base of pre-disaster ground features, which correspond to the spatial location of the post-disaster remote sensing image data set. Based on the solar incidence angle parameters and satellite observation angle parameters in the post-disaster remote sensing image data set, and combined with the building geometric information in the pre-disaster building shadow area topology description set, the spatial location of the shadow area is geometrically calculated to generate a post-disaster expected shadow area spatial distribution map. The post-disaster expected shadow area spatial distribution map is then overlaid with the actual shadow area in the post-disaster remote sensing image data set to extract abnormal shadow area offset information. Based on the abnormal offset information of the shadow area and the set of pre-disaster descriptions of the surrounding land features of the building, an indirect evidence chain of building structural damage is constructed. The spatial location of suspected damaged building objects in the post-disaster remote sensing image data set is determined by the indirect evidence chain, and a distribution map of suspected damaged buildings carrying indirect damage evidence markers is generated. Each suspected damaged building in the suspected damaged building distribution map is classified according to its damage pattern to determine the damage type identifier of each suspected damaged building and calculate its initial damage score. A damage assessment result dataset containing the spatial location coordinates of the damaged buildings and the damage level code is generated based on the damage type identifier and the initial damage degree score; After generating a distribution map of suspected damaged buildings carrying indirect damage evidence markers, the method performs spatial clustering analysis on the suspected damaged building objects in the distribution map to generate a distribution map of damaged building cluster areas, including: Extract the geometric center point of the base polygon of each suspected damaged building object in the distribution map of suspected damaged buildings as the spatial point representation of the suspected damaged building object; Density-based spatial clustering is performed on all spatial point representations. A preset spatial search radius is used as the neighborhood range parameter, and a preset minimum number of neighborhood points is used as the density threshold parameter. Densely distributed point clusters are identified as the initial damaged cluster areas. For each initial damage cluster, construct its convex hull polygon, and use the spatial range covered by the convex hull polygon as the initial boundary of the damage cluster; Obtain the set of damage type identifiers for all suspected damaged buildings within each initial damage cluster area, count the frequency of each damage type identifier within the cluster area, and take the damage type identifier with the highest frequency as the dominant damage type identifier for that damage cluster area; Obtain an initial set of damage severity scores for all suspected damaged buildings within each initial damage cluster area, and calculate the average value of the initial damage severity scores as the regional average damage severity score for that damage cluster area; Based on the dominant damage type identifier and the regional average damage score, the initial boundary of the preliminary damage cluster area is smoothed and expanded to generate a distribution map of the damaged building cluster area containing the dominant damage type identifier attribute label and the regional average damage score attribute label; A spatial overlay analysis was performed on the initially damaged clustered areas, and a functional assessment report of the damaged area was generated by combining the regional functional attribute data in the pre-disaster ground structure knowledge base.

2. The rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in claim 1, characterized in that, The step of retrieving the pre-constructed pre-disaster ground feature structure knowledge base, which corresponds to the spatial location of the post-disaster remote sensing image data set, from the pre-built pre-disaster ground feature structure knowledge base, includes: The spatial coverage boundary coordinates and image acquisition timestamp of each post-disaster remote sensing image unit in the post-disaster remote sensing image dataset are analyzed. An image spatial index code is generated based on the spatial coverage boundary coordinates. The solar position parameters and satellite position parameters are determined based on the image acquisition timestamp. Using the image spatial index code as the retrieval keyword, a spatial range matching query is performed in the spatial index partition table of the pre-disaster ground feature structure knowledge base. The pre-disaster ground feature structure knowledge base stores parametric description data partitions of pre-disaster ground feature structures organized by geographic grid units. The pre-disaster land cover structure parameterized description data partitions corresponding to all geographic grid units that have spatial overlap with the image spatial index code in the spatial range matching query results are used as the initial screening data partition set; For each preliminary screening data partition unit in the preliminary screening data partition set, extract the building shadow area topology description sub-block stored in the preliminary screening data partition unit. The building shadow area topology description sub-block contains a sequence of shadow area outline polygons of multiple building objects under standard solar incidence angle conditions and the encoding of the topological connection relationship between each shadow area outline polygon and the corresponding building base polygon. For each preliminary screening data partition unit in the preliminary screening data partition set, extract the building surrounding feature layout description sub-block stored in the preliminary screening data partition unit. The building surrounding feature layout description sub-block includes the road centerline vector line segment within a preset distance range around each building object, the base polygon of the adjacent building, the boundary polygon of the vegetation cover area, and the spatial distance parameters and azimuth parameters between each corresponding feature element and the central building object. Based on the solar position parameters and the satellite position parameters, the shadow area outline polygon sequence in the topological descriptor block of the building shadow area is subjected to projection transformation processing. The shadow area outline polygon sequence is transformed from the reference coordinate system under the standard solar incidence angle condition to the projection coordinate system under the actual observation geometry condition corresponding to the image acquisition timestamp, and the projection-transformed shadow area outline polygon sequence is generated. The sequence of shadow area outline polygons after projection transformation is associated and combined with the corresponding building base polygons to construct a set of pre-disaster building shadow area topological descriptions with the building base polygons as index keys. The spatial topological relationships of the road centerline vector segments, the base polygons of adjacent buildings, and the boundary polygons of vegetation cover areas in the sub-block describing the layout of features around the building are reconstructed to establish a topological map of the layout of features around the building with the central building object as the node, which serves as the pre-disaster description set of the layout of features around the building.

3. The rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in claim 1, characterized in that, The process involves using the solar incidence angle parameters and satellite observation angle parameters from the post-disaster remote sensing image dataset, combined with the building geometric information from the pre-disaster building shadow area topology description dataset, to perform geometric calculations of the spatial location of the shadow area, generating a post-disaster expected shadow area spatial distribution map. This post-disaster expected shadow area spatial distribution map is then overlaid with the actual shadow areas in the post-disaster remote sensing image dataset to extract abnormal shadow area offset information, including: The solar azimuth angle, solar altitude angle, satellite azimuth angle, and satellite altitude angle values ​​corresponding to each post-disaster remote sensing image unit are parsed from the metadata of the post-disaster remote sensing image data set and used as solar incidence angle parameters and satellite observation angle parameters. Obtain the base polygon coordinate sequence and height attribute parameter of each building object in the pre-disaster building shadow area topology description set. The height attribute parameter is stored in the pre-disaster ground feature structure knowledge base and associated with the building object identifier. Calculate the three-dimensional unit vector of the solar incidence direction based on the solar azimuth angle value and the solar altitude angle value; calculate the three-dimensional unit vector of the satellite observation direction based on the satellite azimuth angle value and the satellite altitude angle value. For each building object, based on its height attribute parameter and the direction of solar incidence, calculate the coordinates of the shadow projection points of the vertices of its base polygon boundary on the horizontal plane, and connect the coordinates of the shadow projection endpoints of all vertices to form the initial shadow area polygon generated by the building object under the illumination of the solar incident light. The initial shadow area polygon is subjected to visibility clipping. Based on the three-dimensional unit vector of the satellite observation direction, the invisible part of the initial shadow area polygon is determined by the building itself. The invisible part is removed from the initial shadow area polygon to generate the post-disaster expected shadow area polygon that meets the satellite observation geometry. Merge the polygons of the expected shadow area of ​​all building objects after the disaster, eliminate the overlapping areas between adjacent shadow area polygons, and generate a spatial distribution map of the expected shadow area of ​​the disaster covering the entire image range. The post-disaster remote sensing image dataset is subjected to automatic shadow region extraction processing. By analyzing the spectral and texture features of each pixel in the image, the pixel regions belonging to shadow are identified, and the identified shadow pixel regions are aggregated and connected to generate a spatial distribution map of the actual shadow areas after the disaster. The spatial distribution map of the expected shadow area after the disaster is compared pixel by pixel with the spatial distribution map of the actual shadow area after the disaster. The pixel positions where shadows were expected to exist but did not actually exist and the pixel positions where shadows were expected not to exist but actually existed are marked, and a binary map of shadow area difference marking is generated. For each connected difference region in the binary image of the shadow area difference marker, calculate its spatial offset direction and spatial offset distance with the base polygon of the adjacent building. Combine the spatial offset direction and spatial offset distance into the shadow area abnormal offset vector corresponding to each difference region. The shadow area abnormal offset vectors corresponding to all difference regions constitute the shadow area abnormal offset information.

4. The rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in claim 1, characterized in that, The process involves constructing an indirect evidence chain for building structural damage based on the abnormal offset information of the shadow area and the pre-disaster description set of surrounding features. This indirect evidence chain is then used to spatially locate suspected damaged building objects in the post-disaster remote sensing image dataset, generating a distribution map of suspected damaged buildings carrying indirect damage evidence markers. This includes: Each shadow area anomaly offset vector in the shadow area anomaly offset information is spatially matched with the building base polygon in the pre-disaster building shadow area topology description set to determine the candidate building object associated with each shadow area anomaly offset vector; For each candidate building object, count the number of all associated shadow area anomaly offset vectors, and calculate the average offset direction and average offset distance of the shadow area anomaly offset vectors to generate a shadow anomaly statistical feature descriptor for the candidate building object. Obtain the topology map of the surrounding land features corresponding to the candidate building object from the pre-disaster building surrounding land feature layout description set. Extract the vertical distance parameter between the candidate building object and the center line of the surrounding road, the interval distance parameter between the candidate building object and the base of the adjacent building, and the area ratio parameter of the vegetation coverage area around the candidate building object from the surrounding land feature layout topology map. The local image area covering the candidate building object in the post-disaster remote sensing image data set is processed by extracting the current status of ground features to identify the actual road centerline location, the base outline of adjacent buildings and the boundary of vegetation coverage area after the disaster, and generate a post-disaster surrounding ground feature layout update map. The road centerline position in the post-disaster surrounding feature layout update map is compared with the road centerline vector line segment in the surrounding feature layout topology map, and the road centerline position offset and road continuity change index are calculated. The outline of the base of the neighboring building in the updated map of the surrounding features after the disaster is compared with the polygon of the base of the neighboring building in the topology map of the surrounding features, and the change rate of the base area of ​​the neighboring building and the offset of the base position of the neighboring building are calculated. The boundary of the vegetation cover area in the post-disaster surrounding land feature layout update map is compared with the polygon of the vegetation cover area boundary in the surrounding land feature layout topology map, and the vegetation cover area boundary change rate and the texture fragmentation index inside the vegetation cover area are calculated. The road centerline offset, road continuity change index, adjacent building base area change rate, adjacent building base position offset, vegetation cover area boundary change rate, and vegetation cover area internal texture fragmentation index are combined to generate the surrounding land feature change vector of the candidate building object. The shadow anomaly statistical feature descriptor is concatenated and fused with the surrounding land feature change feature vector to generate the indirect evidence chain feature vector of the candidate building object; The indirect evidence chain feature vector is input into a pre-trained suspected damage classifier. The suspected damage classifier outputs a confidence score that the candidate building object belongs to a suspected damaged building. Candidate building objects whose confidence scores exceed a preset confidence threshold are selected as suspected damaged building objects. In the spatial coordinate system of the post-disaster remote sensing image data set, the base polygon position of each suspected damaged building object is marked, and an indirect evidence chain feature vector is attached to the base polygon of each suspected damaged building object as an attribute label to generate a distribution map of suspected damaged buildings carrying indirect damage evidence labels.

5. The rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in claim 1, characterized in that, The process of classifying the damage patterns of each suspected damaged building in the distribution map of suspected damaged buildings, determining the damage type identifier of each suspected damaged building, and calculating its initial damage score includes: Extract the base polygon of each suspected damaged building object from the distribution map of suspected damaged buildings, and based on the base polygon, cut out all shadow area abnormal offset vectors covering the base polygon and its surrounding preset buffer range from the shadow area abnormal offset information to form the local abnormal offset vector field of the suspected damaged building object; Spatial gradient analysis is performed on the local abnormal offset vector field to calculate the spatial change rate of the offset vector and the directional consistency coefficient of the offset vector in the local abnormal offset vector field, and an offset field spatial distribution feature descriptor of the suspected damaged building object is generated. Based on the base polygon, the original surrounding land feature layout topology map of the suspected damaged building object is obtained from the pre-disaster building surrounding land feature layout description set, and the post-disaster surrounding land feature layout update map of the suspected damaged building object is extracted from the post-disaster remote sensing image data set. The original surrounding land feature layout topology map and the post-disaster surrounding land feature layout updated map are input into a pre-trained graph change detection network. The graph change detection network calculates the differences between the node features and edge features in the two topology maps through a graph neural network and outputs a surrounding land feature layout change pattern encoding vector. The offset field spatial distribution feature descriptor and the surrounding land feature layout change pattern encoding vector are concatenated to generate the damage feature fusion vector of the suspected damaged building object. The damage feature fusion vector is input into a pre-trained damage pattern classifier. The damage pattern classifier contains multiple binary sub-networks, each of which corresponds to a preset damage type. Each binary sub-network outputs the probability value of the suspected damaged building object belonging to the corresponding damage type based on the damage feature fusion vector. Select the damage type corresponding to the maximum probability value output by all binary sub-networks as the damage type identifier of the suspected damaged building object, and use the maximum probability value as the initial score of the damage degree of the suspected damaged building object; Based on the damage type identifier, the spatial distribution template of the area of ​​concern corresponding to the damage type is queried from the preset damage type parameter table. The spatial distribution template of the area of ​​concern is used to indicate the main manifestation area of ​​the corresponding type of damage within the building base area. The spatial distribution template of the area of ​​interest is spatially superimposed with the local anomaly offset vector field. The average amplitude value and average directional dispersion of the offset vector within the coverage area of ​​the spatial distribution template of the area of ​​interest are calculated. The initial damage score is weighted and corrected based on the average amplitude value and the average directional dispersion to generate a corrected damage score. The damage type identifier and the revised damage severity score are used as the final damage pattern classification result for each suspected damaged building object.

6. The rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in claim 1, characterized in that, The damage assessment result dataset, which includes the spatial location coordinates of the damaged building and the damage level code, is generated based on the damage type identifier and the initial damage severity score. Obtain the base polygon vertex coordinate sequence for each suspected damaged building object, and calculate the coordinates of the geometric center point of the base polygon vertex coordinate sequence as the spatial position coordinates of the suspected damaged building object; Based on the damage type identifier of each suspected damaged building object, the basic damage level code corresponding to that damage type identifier is retrieved from the preset damage level mapping table; Based on the initial damage score of each suspected damaged building object, the level is fine-tuned on the basis of the basic damage level code. The initial damage score is mapped to the level fine-tuning increment value. The basic damage level code and the level fine-tuning increment value are added to generate the final damage level code of the suspected damaged building object. Extract the indirect evidence chain feature vector of each suspected damaged building object from the distribution map of suspected damaged buildings, perform hash encoding on the indirect evidence chain feature vector, and generate the evidence tracing index code associated with the suspected damaged building object; The spatial location coordinates, the final damage level code, and the evidence tracing index code are combined to generate a damage assessment record unit for the suspected damaged building object; All damage assessment record units of suspected damaged buildings are spatially sorted according to their spatial location coordinates to generate a damage assessment record sequence organized by a spatial grid. Redundancy elimination processing is performed on the damage assessment record sequence. It is checked whether there are multiple damage assessment record units corresponding to the same building object. If so, the damage assessment record unit with the highest damage degree score after correction is retained and the remaining redundant records are deleted to generate a set of redundancy-free damage assessment records. The deduplicated damage assessment record set is encapsulated into a standard geographic information data format, and a data production timestamp and data version identifier are added to each damage assessment record to generate a damage assessment result dataset containing the spatial location coordinates of the damaged buildings and the damage level code.

7. The rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in claim 1, characterized in that, The process involves performing spatial overlay analysis on the initially damaged clustered areas, and generating a functional assessment report of the damaged area by combining regional functional attribute data from the pre-disaster ground structure knowledge base. This includes: Retrieve regional functional attribute data partitions covering the spatial range of the initially damaged clustered area from the pre-disaster ground structure knowledge base. The regional functional attribute data partitions include land use classification codes, population density distribution raster layers, transportation network importance level vector layers, and lifeline pipeline distribution vector layers. Spatial overlay analysis is performed between the initial boundary of the preliminary damaged cluster area and the land use classification code, and the area proportion corresponding to different land use classification codes in the preliminary damaged cluster area is calculated to generate a statistical table of land use type composition of the damaged area. Spatial overlay analysis is performed between the initial boundary of the initially damaged cluster area and the population density distribution raster layer, and the sum of the population density raster values ​​within the initially damaged cluster area is used as an estimate of the affected population size. Spatial overlay analysis is performed between the initial boundary of the preliminary damage cluster area and the traffic network importance level vector layer to identify the traffic network elements at all levels contained in the preliminary damage cluster area, and the importance level identifier and element length of each included traffic network element are extracted to generate a list of affected traffic networks. Spatial overlay analysis is performed between the initial boundary of the preliminary damage cluster area and the lifeline pipeline distribution vector map layer to identify various lifeline pipeline elements contained in the preliminary damage cluster area, and the pipe diameter parameters, material parameters and start and end point coordinates of each included lifeline pipeline element are extracted to generate a list of affected lifeline pipelines. The estimated affected population size is weighted and corrected based on the dominant damage type identifier and the regional average damage score to generate a corrected assessment value of the affected population size. The damage impact area functional assessment report is generated by combining the statistical table of land use types in the damaged area, the revised assessment value of the affected population size, the list of affected transportation networks, and the list of affected lifeline pipelines.

8. A rapid system for detecting changes and assessing damage in remote sensing images for disaster emergency response, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in any one of claims 1 to 7 by executing the machine-executable instructions.

9. A computer program product, characterized in that, The computer program product includes machine-executable instructions stored in a computer-readable storage medium. A processor of a computer device reads the machine-executable instructions from the computer-readable storage medium and executes the machine-executable instructions, causing the computer device to perform the rapid method for remote sensing image change detection and damage assessment for disaster emergency response as described in any one of claims 1 to 7.

Images