Image processing method and device for aerial remote sensing image and storage medium

Abstract

The invention discloses an image processing method and device for an aerial remote sensing image and a storage medium, and relates to the field of image analysis, and the method comprises the steps: obtaining an aerial image shot by an unmanned plane in real time, and carrying out the image preprocessing of the aerial image, and obtaining a preprocessed aerial image; performing rigid framework extraction on the preprocessed aerial image to obtain an original incomplete framework graph after gravel interference is removed; carrying out fracture framework repair and completion on the original incomplete framework diagram to obtain an unmanned aerial vehicle visual angle topological structure diagram; performing elastic topological mapping matching on the unmanned aerial vehicle visual angle topological structure diagram and a pre-stored satellite reference diagram to obtain a matching result; and carrying out pose calculation according to a matching result, and outputting accurate coordinate information of the unmanned aerial vehicle in the world coordinate system. According to the method and the device, high-precision positioning can be realized in an extreme environment in which the GPS is rejected and the surface texture is subjected to unstructured deformation.

Description

Technical Field

[0001] This application relates to the field of image analysis, and in particular to an image processing method, apparatus and storage medium for aerial remote sensing images. Background Technology

[0002] With the rapid development of drone technology, its application in emergency rescue, disaster assessment, and other fields is becoming increasingly widespread. Especially after geological disasters such as earthquakes and landslides, drones can quickly enter disaster areas for reconnaissance and material delivery, becoming an important tool in rescue efforts. However, in these extreme environments, drones face severe positioning challenges. On the one hand, GPS signals often malfunction due to electromagnetic interference from collapsed buildings, terrain obstruction, or human interference, rendering traditional satellite positioning systems ineffective. On the other hand, the drastic changes in the surface environment caused by disasters also pose unprecedented difficulties for visual positioning technology.

[0003] In existing technologies, visual positioning mainly relies on SLAM (Simultaneous Localization and Mapping) technology and feature-point-based image matching methods. SLAM technology constructs an environmental map in real time and locates its own position within the map, but it has high requirements for the continuity and stability of environmental features. In post-disaster scenarios, original road markings, building textures, and other features are completely covered by rubble, making it impossible for SLAM systems to establish effective feature associations. Feature-point matching methods based on algorithms such as SIFT and SURF achieve positioning by extracting key points in the image and matching them with a pre-stored reference map. These methods perform well in normal environments, but their core assumption is the relative stability of scene textures. When unstructured deformation occurs on the ground, there are huge texture differences between real-time images captured by drones and pre-disaster satellite reference maps. Originally clear roads are covered by rubble, building surfaces are damaged, and vegetation is buried, causing a fundamental change in the description of feature points, resulting in a sharp drop in the success rate of feature matching or even complete failure. Therefore, in extreme environments where GPS is denied and the surface texture undergoes unstructured deformation, existing visual positioning technologies rely too heavily on the stability of texture features, causing the matching between the UAV and the pre-stored reference map to fail, thus making it impossible to achieve high-precision positioning. Summary of the Invention

[0004] This application provides an image processing method, apparatus, and storage medium for aerial remote sensing images, used to achieve high-precision positioning in extreme environments where GPS is denied and the surface texture undergoes unstructured deformation.

[0005] To achieve the above objectives, the embodiments of this application adopt the following technical solutions: Firstly, an image processing method for aerial remote sensing images is provided, the method comprising: The system acquires aerial images captured in real time by drones and performs image preprocessing on the aerial images to obtain preprocessed aerial images. Rigid structure extraction is performed on the preprocessed aerial images to obtain the original incomplete structure diagram after removing debris interference; The original incomplete architecture diagram is repaired and completed to obtain the topology diagram from the perspective of the UAV. The topology map from the UAV's perspective is matched with a pre-stored satellite reference map using flexible topology mapping to obtain the matching result; Based on the matching results, pose calculation is performed, and the precise coordinate information of the UAV in the world coordinate system is output.

[0006] In one possible implementation of the first aspect, rigid architecture extraction is performed on the preprocessed aerial imagery to obtain the original fragmented architecture diagram after removing debris interference, including: Obtain the pixel grayscale distribution of the preprocessed aerial imagery, and construct a three-dimensional terrain matrix based on the pixel grayscale distribution; Determine the virtual fluid parameters and initialize the full-map uniform flow field of the 3D terrain matrix; Based on virtual fluid parameters, virtual fluid is iteratively flowed on a three-dimensional terrain matrix, where the virtual fluid parameters include the velocity decay coefficient and the viscosity coefficient. During the iterative flow process, regions where the flow velocity diverges or decays rapidly are identified in the three-dimensional terrain matrix, as well as streamline regions where the flow velocity converges and remains stable. Set the pixel gradient value to zero in areas where the flow velocity is diverging or rapidly decaying, and enhance the pixel value to a preset value in streamline areas where the flow velocity converges and remains stable, in order to remove gravel interference. After a preset number of iterations, the streamline convergence center is determined, and a binarized original incomplete architecture diagram is generated based on the streamline convergence center.

[0007] In another possible implementation of the first aspect, virtual fluid parameters are determined, and a uniform flow field across the entire three-dimensional terrain matrix is ​​initialized, including: To obtain the resolution and shooting altitude of aerial images; The initial value of the flow velocity attenuation coefficient is determined based on the resolution and shooting altitude of the aerial imagery; The initial value of the viscosity coefficient is determined based on the pre-stored characteristics of the gravel density in the post-disaster environment; An initial velocity vector is set at each pixel location of the 3D terrain matrix, and the direction of the initial velocity vector points in the gradient descent direction. The initial velocity vector is normalized to initialize a uniform flow field across the entire 3D terrain matrix.

[0008] In another possible implementation of the first aspect, setting the pixel gradient values ​​of regions with divergent or rapidly decaying flow rates to zero includes: Calculate the velocity divergence and velocity decay rate at each pixel location in the 3D terrain matrix; Determine whether the velocity divergence is greater than a preset velocity threshold, or whether the velocity decay rate is greater than a preset gradient threshold. If the velocity divergence is greater than the velocity threshold or the velocity decay rate is greater than the gradient threshold, the gradient value at the corresponding pixel position is set to zero to simulate the effect of sediment being washed away.

[0009] In another possible implementation of the first aspect, the original incomplete architecture diagram is repaired and completed to obtain a topology diagram from the perspective of the UAV, including: Identify all line break points in the original incomplete architecture diagram and calculate the tangent direction vector for each line break point; A virtual stress field is constructed between each pair of adjacent line breaks, and a probability score for connection between each pair of adjacent line breaks is calculated based on the virtual stress field. A growth threshold is set based on the probability score, and connecting lines are generated for line breakpoints with probability scores higher than the growth threshold. The generated connection lines are merged with the original incomplete architecture diagram to output the repaired topology diagram from the drone's perspective.

[0010] In another possible implementation of the first aspect, a virtual stress field is constructed between each pair of adjacent line breaks, and a probability score for connection between each pair of adjacent line breaks is calculated based on the virtual stress field, including: Calculate the Euclidean distance between each pair of adjacent breakpoints; Calculate the cosine of the angle between the tangent direction vectors of each pair of adjacent breakpoints; Extract the original texture residue of pixels between each pair of adjacent breakpoints, where the Euclidean distance, the cosine of the included angle, and the original texture residue constitute a virtual stress field; Based on the Euclidean distance, the cosine of the included angle, and the original texture residual, a weighted summation is used to calculate the probability score of connection between each pair of adjacent line breakpoints.

[0011] In another possible implementation of the first aspect, connecting lines are generated for line breakpoint pairs with probability scores higher than a growth threshold, including: For line breakpoint pairs with probability scores higher than the growth threshold, determine the starting and ending coordinates of the connecting line. Calculate the curvature parameter of the connecting line based on the tangent direction vectors at the two line break points; Connector lines are generated based on curvature parameters.

[0012] In another possible implementation of the first aspect, a flexible topology mapping matching is performed between the UAV-view topology map and a pre-stored satellite reference map to obtain the matching result, including: Based on satellite baseline images, standard roads and standard building structures are extracted, and a set of baseline gravity sources is constructed based on the standard roads and standard building structures; Define the nodes in the topology graph from the drone's perspective as the set of mass points to be matched; The set of particles to be matched is projected onto the coordinate system of the reference gravitational source set to obtain the initial rough position; Iteratively calculate the virtual gravitational resultant force of all reference gravitational sources on each particle in the set of particles to be matched, and move the position of the particle according to the virtual gravitational resultant force; During the movement of the particle's position, the total potential energy change rate of the topology map from the UAV's perspective is calculated in real time. When the total potential energy change rate is lower than the preset value, the iteration stops and the final position of the particle to be matched is output as the matching result.

[0013] Secondly, this application provides an image processing apparatus for aerial remote sensing images, comprising: The image preprocessing module is used to acquire aerial images captured in real time by the drone and to perform image preprocessing on the aerial images to obtain preprocessed aerial images. The rigid structure extraction module is used to extract the rigid structure from the preprocessed aerial images to obtain the original incomplete structure diagram after removing debris interference. The architecture repair module is used to repair and complete the broken architecture diagram of the original incomplete architecture diagram to obtain the topology diagram from the perspective of the UAV. The topology matching module is used to perform flexible topology mapping matching between the topology map from the perspective of the UAV and the pre-stored satellite reference map to obtain the matching result; The pose calculation module is used to calculate the pose based on the matching results and output the precise coordinate information of the UAV in the world coordinate system.

[0014] Thirdly, this application provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the above-described image processing method for aerial remote sensing images.

[0015] By introducing virtual fluid dynamics and virtual gravitational field models, the above technical solutions effectively address the challenge of UAV positioning in environments with GPS rejection and unstructured surface deformation. Virtual fluid simulation technology cleverly utilizes the fundamental difference in geometric continuity between rigid structures and debris accumulations to successfully extract stable road and building skeletons from chaotic scenes, overcoming the dependence of traditional feature extraction methods on texture stability. The fractured structure repair technology, based on a virtual stress field model, intelligently restores structural connections broken by disasters, reconstructing complete topological relationships. The elastic topological mapping matching method allows for local structural deformation and displacement, exhibiting stronger environmental adaptability compared to traditional rigid matching. Even under drastic surface changes, it can establish reliable matching relationships, ultimately achieving high-precision positioning in extreme environments. These technical solutions do not rely on GPS signals or require stable surface textures; they only require pre-stored satellite reference images, making them highly practical and providing reliable technical support for UAV navigation in post-disaster emergency rescue, significantly improving rescue efficiency and safety.

[0016] Other features and advantages of the embodiments of this application will be described in detail in the following detailed description section. Attached Figure Description

[0017] Figure 1 A flowchart illustrating an image processing method for aerial remote sensing images provided in an embodiment of this application; Figure 2 A schematic diagram of a steady-state determination process for a topological structure from the perspective of an unmanned aerial vehicle (UAV) based on the rate of change of potential energy, provided for an embodiment of this application; Figure 3 This is a schematic diagram of the structure of an image processing device for aerial remote sensing images provided in an embodiment of this application. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only for illustration and explanation of the embodiments of this application and are not intended to limit the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0019] It should be noted that if the embodiments of this application involve directional indicators (such as up, down, left, right, front, back, etc.), the directional indicators are only used to explain the relative positional relationship and movement of each component in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicators will also change accordingly.

[0020] Furthermore, if the embodiments of this application involve descriptions such as "first" or "second," these descriptions are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, features defined with "first" or "second" may explicitly or implicitly include at least one of those features. Additionally, the technical solutions of various embodiments can be combined with each other, but this must be based on the ability of those skilled in the art to implement them. If the combination of technical solutions is contradictory or impossible to implement, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed in this application.

[0021] Figure 1 The illustration schematically shows a flowchart of an image processing method for aerial remote sensing images according to an embodiment of this application. For example... Figure 1 As shown in the figure, this application provides an image processing method for aerial remote sensing images, which may include the following steps.

[0022] S110. Acquire aerial images captured in real time by the drone, and perform image preprocessing on the aerial images to obtain preprocessed aerial images. S120. Extract the rigid structure from the preprocessed aerial image to obtain the original incomplete structure diagram after removing debris interference. S130. Repair and complete the broken architecture diagram of the original incomplete architecture diagram to obtain the topology diagram from the perspective of the UAV. S140. Perform flexible topology mapping matching between the topology map from the UAV's perspective and the pre-stored satellite reference map to obtain the matching result; S150. Based on the matching results, perform pose calculation and output the precise coordinate information of the UAV in the world coordinate system.

[0023] In this embodiment, aerial images are preprocessed to obtain preprocessed aerial images, including: The aerial images are processed to obtain grayscale aerial images; Aerial grayscale images are enhanced by filtering at specific frequencies to suppress interference from dust and uneven lighting. Edge detection is performed on the filtered and enhanced image to extract the geometric edges of the objects; Morphological processing is performed on the geometric edges to obtain the preprocessed aerial image.

[0024] Specifically, the aerial grayscale image undergoes specific frequency filtering enhancement processing, including: The cutoff frequency of the filter is determined based on the dust particle size characteristics of the post-disaster environment. A bandpass filter is used to perform frequency domain filtering on aerial grayscale images to retain mid-frequency components and suppress high-frequency noise and low-frequency background. The filtered image is enhanced for contrast to highlight edge features.

[0025] In disaster relief scenarios, drones first collect real-time ground image data using their onboard high-definition cameras. The raw aerial images often contain a lot of noise interference, including dust particles stirred up after the disaster, uneven lighting, and sensor noise. The image preprocessing stage first converts the color RGB images into grayscale images, and calculates the grayscale value of each pixel using a weighted average method. This reduces the computational complexity of subsequent processing while preserving key structural information in the image.

[0026] To address the prevalent dust interference problem in post-disaster environments, the preprocessing module needs to determine filter parameters based on the statistical characteristics of dust particle size in the actual post-disaster environment. Specifically, by analyzing the dust particle size distribution in historical disaster data, it can be determined that dust is mainly concentrated in the high-frequency region of the image frequency domain, while useful building and road edge information is mainly distributed in the mid-frequency region. Therefore, a bandpass filter is used to perform frequency domain filtering on the grayscale image, and an appropriate cutoff frequency range is set to effectively suppress high-frequency noise and filter out low-frequency background changes, thereby preserving the edge structure information in the mid-frequency band.

[0027] After filtering, to further highlight edge features in the image, contrast enhancement is required. Histogram equalization or adaptive contrast enhancement algorithms are used to stretch the image's grayscale dynamic range, making the contrast between edges and background more pronounced. Subsequently, edge detection operators, such as the Canny or Sobel operators, are applied to extract the geometric edges of objects in the image. Edge detection results typically contain discontinuous edge points. To obtain continuous and complete edge contours, morphological processing is also necessary, including dilation and erosion operations, connecting broken edge fragments, and removing isolated noise points, ultimately yielding the preprocessed aerial image. This preprocessing effectively reduces the difficulty of subsequent processing, laying a solid foundation for rigid structure extraction.

[0028] The purpose of rigid structure extraction is to extract stable road and building skeletons from chaotic scenes filled with rubble and debris. In this embodiment, a three-dimensional terrain matrix is ​​first constructed based on the pixel grayscale distribution of the preprocessed aerial imagery, mapping the grayscale values ​​of the two-dimensional image to the third-dimensional height information; a larger grayscale value indicates a higher terrain height. This three-dimensional representation provides the terrain foundation for subsequent virtual fluid simulation. Next, a virtual fluid dynamics model is introduced for rubble interference removal. The velocity attenuation coefficient needs to be determined based on the resolution and shooting altitude of the aerial imagery. The higher the resolution and the lower the shooting altitude, the smaller the velocity attenuation coefficient should be to ensure that the fluid can flow fully on the detailed terrain. The viscosity coefficient is determined based on the pre-stored rubble density characteristics of the post-disaster environment; the greater the rubble density, the larger the viscosity coefficient should be to simulate the physical phenomenon of fluid flow being obstructed in rubble accumulation areas.

[0029] In the initialization phase, an initial velocity vector is set at each pixel location in the 3D terrain matrix, pointing towards the gradient descent direction at that location, i.e., flowing from high to low. All initial velocity vectors are then normalized to form a uniform flow field across the entire map. The iterative flow process of the virtual fluid simulates the flow behavior of real fluids on the terrain surface. During the flow, the fluid will experience velocity divergence or rapid attenuation in areas with accumulated gravel, because the gravel surface is rough and irregular, and the fluid cannot form stable streamlines. However, at rigid structures such as roads and building edges, where the surface is relatively flat and has a clear geometric orientation, the fluid will converge along these structures and form stable streamlines.

[0030] By calculating the velocity divergence and velocity decay rate at each pixel location, gravel areas and rigid structure areas can be identified. Areas with velocity divergence greater than a preset threshold or velocity decay rate greater than a preset gradient threshold are identified as gravel interference areas. The pixel gradient values ​​in these areas are then set to zero to simulate the effect of sediment being washed away.

[0031] Conversely, streamline regions where flow velocities converge and remain stable are identified as rigid structures. The pixel values ​​of these regions are then amplified to preset values ​​to make them more prominent in the image. After multiple iterations of flow, the streamlines gradually converge to the center of the rigid structure, forming a clear streamline convergence center. Based on these streamline convergence centers, a binarized original incomplete structure image is generated, where white pixels represent the extracted rigid structure and black pixels represent the background. This virtual fluid dynamics-based method cleverly utilizes the fundamental difference in geometric continuity between rigid structures and debris accumulation, effectively removing debris interference.

[0032] Although the original fragmented architecture diagram has removed debris interference, the damage caused by the disaster resulted in broken road and building structures, presenting discontinuous line fragments. The fractured architecture repair and completion steps aim to restore the fractured structure and reconstruct the complete topological connections. First, it is necessary to identify all line breaks in the original fragmented architecture diagram; these breaks are typically located at the endpoints of the lines. For each identified break, its tangent direction vector is calculated. This vector represents the direction of the line's extension at the break point, and the tangent direction can be estimated by fitting several pixels near the break point.

[0033] Next, a virtual stress field is constructed between each pair of adjacent line breakpoints. The virtual stress field is a physical field model that simulates the tendency of structural connections. The virtual stress field consists of three key parameters: First, the Euclidean distance, which represents the spatial distance between the two breakpoints; the closer the distance, the greater the possibility of connection. Second, the cosine value of the angle between the tangent direction vectors, which reflects whether the extension directions of the two breakpoints tend to be consistent; the closer the cosine value is to 1, the more consistent the directions, and the greater the possibility of connection. Third, the original texture residue, which assesses whether there are residual traces of roads or buildings in the area between the two breakpoints by analyzing the original image texture features. The higher the texture residue, the greater the possibility that a connecting structure originally existed.

[0034] Based on the three parameters mentioned above, a weighted summation is used to calculate the probability score of connecting each pair of adjacent breakpoints. Specifically, the distance factor uses a negative exponential decay function, with lower weights for greater distances; the orientation consistency factor directly uses the cosine of the included angle; and the texture persistence factor is quantified using local texture variance. The three factors are then weighted and summed according to preset weights to obtain the final probability score. A growth threshold is set based on the probability score; only breakpoint pairs with scores higher than this threshold are considered to be connected. For these high-scoring breakpoint pairs, connection lines need to be generated to complete the broken architecture.

[0035] It's important to note that the generation of connectors is not a simple straight-line connection. Instead, the curvature parameters of the connectors are calculated based on the tangent direction vectors at the two break points, ensuring a smooth transition between the lines at both ends and conforming to the natural direction of the road or building edges. Smooth connector curves are generated using Bézier curves or spline interpolation methods, based on the start and end coordinates and the curvature parameters. All generated connectors are then fused with the original incomplete architecture diagram to obtain a restored topology map from the drone's perspective. This restoration process effectively recovers the structural connections broken by the disaster, providing complete structural information for subsequent topology matching.

[0036] Elastic topology mapping matching is used to solve the localization problem after unstructured deformation of surface texture. It matches the topology map from the UAV's perspective with a pre-stored satellite reference map, but unlike traditional rigid registration methods, this embodiment employs an elastic matching strategy based on a virtual gravitational field.

[0037] First, standard road networks and standard building structures are extracted from satellite baseline images. These standard structures represent the pre-disaster surface state and have clearly defined geometric locations and topological relationships. Key nodes such as extracted standard road intersections, road turning points, and building corners are defined as baseline gravity sources, constructing a set of baseline gravity sources. Each baseline gravity source is assigned a virtual mass, which is related to the importance of the node in the topological structure; for example, the mass of a main road intersection is greater than that of a secondary road intersection. Simultaneously, key nodes extracted from the topological structure map from the UAV perspective are defined as mass points to be matched; these mass points represent structural nodes actually observed after the disaster.

[0038] The matching process first projects the set of particles to be matched onto the world coordinate system of the set of reference gravitational sources, obtaining initial positions through coarse scale estimation and orientation alignment. Then, an iterative optimization phase is initiated, calculating the virtual resultant gravitational force exerted on each particle by all reference gravitational sources. The magnitude of the virtual gravity is proportional to the mass of the gravitational source and inversely proportional to the square of the distance from the particle to the source, with its direction pointing towards the source.

[0039] In this system, the net force acting on each particle determines its direction and velocity of movement. Particles move towards the direction of stronger gravity, gradually being attracted to the vicinity of the best-matching reference gravitational source. During the particle's movement, the total potential energy of the entire system is calculated in real time, defined as the weighted sum of the distances from all particles to their nearest gravitational source. As iteration progresses, the particles gradually find the optimal matching position, and the total potential energy continuously decreases. When the rate of change of the total potential energy falls below a preset threshold, it indicates that the system has reached a stable state, and iteration stops. The final positions of the particles to be matched are output as the matching results, establishing a correspondence between the topology from the UAV's perspective and the satellite reference image. This elastic matching method based on a virtual gravitational field allows for a certain degree of deformation and displacement of the local structure, enabling it to adapt to unstructured changes in the post-disaster surface and exhibiting stronger robustness compared to traditional rigid matching methods.

[0040] Pose calculation can convert topology matching results into precise UAV position and attitude information. Based on the matching results obtained in the previous step, a mapping relationship has been established between key nodes in the UAV's topology map and corresponding nodes in the satellite reference map. These correspondences constitute a set of two-dimensional to two-dimensional point-to-point matchings.

[0041] Since nodes in the satellite reference image have known world coordinates, while nodes in the UAV's topology map have coordinates in the image coordinate system, the UAV's pose can be recovered by solving the perspective transformation matrix. The perspective transformation matrix contains geometric transformation parameters such as rotation, translation, and scaling, describing the mapping from the world coordinate system to the image coordinate system. Using multiple pairs of matched points, the optimal perspective transformation matrix is ​​solved using either the least squares method or the RANSAC algorithm. The RANSAC algorithm is more robust to mismatched point pairs, automatically removing outliers and retaining interior points for matrix solving.

[0042] The obtained perspective transformation matrix can be decomposed into the camera's intrinsic and extrinsic parameter matrices. The extrinsic parameter matrix is ​​further decomposed into a rotation matrix and a translation vector. The rotation matrix describes the rotation relationship between the UAV camera coordinate system and the world coordinate system, and can be converted into Euler angles or quaternions to represent the UAV's attitude. The translation vector describes the position coordinates of the UAV camera's optical center in the world coordinate system. Considering that the UAV's flight altitude information can be obtained through a barometer or laser rangefinder, combining the altitude information with the horizontal component of the translation vector yields the complete three-dimensional coordinates of the UAV in the world coordinate system.

[0043] Furthermore, by analyzing pose changes between consecutive frames and combining data from the inertial measurement unit, Kalman filtering or extended Kalman filtering can be performed to further improve the accuracy and stability of pose estimation. The final output precise coordinate information includes the UAV's three-dimensional position coordinates and three-axis attitude angles, which can be used for applications such as UAV navigation control, trajectory planning, and target localization.

[0044] This embodiment effectively solves the challenge of UAV positioning in environments with GPS rejection and unstructured surface deformation by introducing virtual fluid dynamics and virtual gravitational field models. Virtual fluid simulation technology cleverly utilizes the fundamental difference in geometric continuity between rigid structures and debris accumulation to successfully extract stable road and building skeletons from chaotic scenes, overcoming the dependence of traditional feature extraction methods on texture stability. The fractured structure repair technology, based on a virtual stress field model, intelligently restores structural connections broken by disaster damage, reconstructing complete topological relationships. The elastic topological mapping matching method allows for local structural deformation and displacement, exhibiting stronger environmental adaptability compared to traditional rigid matching. It can establish reliable matching relationships even under drastic surface changes, ultimately achieving high-precision positioning in extreme environments. This technical solution does not rely on GPS signals or require stable surface textures; it only requires pre-stored satellite reference images to function, making it highly practical and providing reliable technical support for UAV navigation in post-disaster emergency rescue, significantly improving rescue efficiency and safety.

[0045] In one embodiment of this invention, rigid structure extraction is performed on the preprocessed aerial image to obtain the original incomplete structure diagram after removing debris interference, including the following steps: S210. Obtain the pixel grayscale distribution of the preprocessed aerial imagery and construct a three-dimensional terrain matrix based on the pixel grayscale distribution. S220. Determine the virtual fluid parameters and initialize the uniform flow field of the entire three-dimensional terrain matrix; S230. Based on virtual fluid parameters, virtual fluid is iteratively flowed on a three-dimensional terrain matrix, where the virtual fluid parameters include velocity attenuation coefficient and viscosity coefficient. S240. During the iterative flow process, determine the regions where the flow velocity diverges or decays rapidly in the three-dimensional terrain matrix, as well as the streamline regions where the flow velocity converges and remains stable. S250: Set the pixel gradient value of the region where the flow velocity is diverging or rapidly decaying to zero, and enhance the pixel value of the streamline region where the flow velocity converges and remains stable to a preset value in order to remove gravel interference. S260. After a preset number of iterations, determine the streamline convergence center and generate a binarized original incomplete architecture diagram based on the streamline convergence center.

[0046] The preprocessed aerial imagery is stored in a two-dimensional matrix, with each pixel location corresponding to a grayscale value. To achieve rigid structure extraction based on virtual hydrodynamics, the two-dimensional grayscale image needs to be converted into a three-dimensional terrain representation. Each pixel location in the preprocessed aerial imagery is traversed, and the grayscale value at that location is read; this grayscale value reflects the brightness information at that location.

[0047] In post-disaster scenarios, areas with accumulated rubble and debris exhibit irregular, high-frequency fluctuations in grayscale values ​​due to their rough surfaces and complex shadows. In contrast, rigid structures such as roads and building edges, with their relatively smooth surfaces and well-defined geometric orientations, show relatively gentle and directional changes in grayscale values. When constructing the 3D terrain matrix, the row and column coordinates of the image are used as the X and Y axes, and the normalized grayscale values ​​are used as the Z-axis height. Normalization maps the original grayscale values ​​to a specific range, such as 0 to 1, ensuring that the height values ​​have a uniform dimension.

[0048] Locations with higher grayscale values ​​appear as higher peaks in 3D terrain, while locations with lower grayscale values ​​appear as lower valleys. Areas with accumulated gravel appear as rugged, undulating terrain in 3D, while the edges of roads and buildings appear as clearly defined gullies or ridges. This 3D terrain representation provides the physical basis for subsequent virtual fluid simulations, enabling the simulation of fluid flow behavior on the terrain surface, thus distinguishing between rigid structures and gravel disturbances.

[0049] The determination of virtual fluid parameters directly affects the effectiveness of rigid structure extraction and needs to be adaptively set according to the characteristics of the actual scene. The velocity attenuation coefficient describes the rate at which the fluid velocity decreases during flow, and this parameter is closely related to the resolution and shooting altitude of the aerial imagery. When the aerial imagery resolution is high, the image contains more detail and the frequency of terrain undulations is higher, requiring the fluid to flow at a finer scale. Therefore, the velocity attenuation coefficient should be set smaller to allow the fluid to fully explore terrain details. Conversely, when the resolution is low, the terrain undulations are relatively gentle, and the velocity attenuation coefficient can be set larger. The shooting altitude also affects the setting of the velocity attenuation coefficient; the lower the shooting altitude, the clearer the ground details, and the velocity attenuation coefficient should be reduced accordingly.

[0050] The viscosity coefficient describes the degree of viscosity of a fluid, reflecting the magnitude of resistance encountered during its flow. In post-disaster environments, debris density is a key factor influencing the setting of the viscosity coefficient. By analyzing pre-stored historical disaster data, the distribution characteristics of debris density in different types of disaster scenarios can be statistically analyzed. Higher debris density means a rougher surface, greater resistance to fluid flow, and therefore a higher viscosity coefficient should be set.

[0051] When initializing a uniform flow field across the entire map, an initial velocity vector is set at each pixel location in the 3D terrain matrix. The direction of the initial velocity vector is determined by calculating the gradient direction at that location. The gradient direction points in the direction of the fastest increase in grayscale value, while the velocity direction should point in the direction of gradient descent, i.e., from high to low, simulating the natural flow of fluid under gravity. To calculate the gradient, the Sobel operator or other gradient operators are used to calculate the gradient components in the X and Y directions at that location, and then synthesized to obtain the gradient vector. The opposite direction of this gradient is taken as the initial velocity direction. The magnitude of the initial velocity can be set to a uniform value, and then all initial velocity vectors are normalized to a magnitude of 1, forming a uniform flow field across the entire map with different directions but consistent velocity magnitudes.

[0052] The iterative flow simulation of virtual fluids gradually distinguishes between rigid structures and debris interference by simulating the dynamic flow behavior of fluids on a 3D terrain surface. In each iteration, every pixel position in the 3D terrain matrix is ​​traversed, updating the velocity vector at that location. The velocity update follows fundamental principles of fluid mechanics, comprehensively considering factors such as gravity, viscous drag, and terrain constraints. Gravity causes the fluid to accelerate along the gradient descent direction, with the magnitude of acceleration proportional to the gradient amplitude at that location. Viscous drag causes the velocity to gradually decay, with the decay rate determined by both the viscosity coefficient and the velocity decay coefficient.

[0053] In areas with gravel accumulation, the undulating and irregular terrain causes the fluid to constantly change direction during flow, resulting in significant differences in velocity direction between adjacent locations and causing velocity divergence. Simultaneously, the rough surfaces generate substantial frictional resistance, causing the velocity to decay rapidly. However, at rigid structures such as roads and building edges, the terrain has a clear geometric orientation, allowing the fluid to form stable streamlines along these structures. The velocity direction at adjacent locations tends to align, leading to velocity convergence and relative stability.

[0054] During the iteration process, the continuity constraint of the fluid also needs to be considered, that is, the flow rate into a certain location should be equal to the flow rate out of that location. By calculating the velocity divergence at each location, it can be determined whether that location is a convergence point or a divergence point of the fluid. A negative velocity divergence indicates fluid convergence, while a positive one indicates fluid divergence. After multiple iterations, the flow field gradually evolves, forming a clear streamline pattern at rigid structures, while the flow field remains chaotic in the gravel region.

[0055] During the virtual fluid iterative flow process, flow field characteristics are monitored and analyzed in real time to identify gravel regions where velocity diverges or decays rapidly, as well as rigid structure regions where velocity converges and remains stable. For each pixel location, its velocity divergence is calculated. The velocity divergence is obtained by calculating the divergence operator of the velocity vector at that location, specifically the sum of the partial derivatives of the velocity at that location in the X and Y directions. The physical meaning of velocity divergence is the net outflow from that location per unit time; a positive value indicates that the fluid is diverging outward from that location, while a negative value indicates that the fluid is converging towards that location.

[0056] In areas with gravel accumulation, due to irregular terrain, fluid cannot form a stable flow direction, resulting in large absolute values ​​of velocity divergence and frequent sign changes. A velocity divergence threshold is set; when the absolute value of the velocity divergence at a certain location exceeds this threshold, that location is determined to belong to a velocity divergence region. Simultaneously, the velocity decay rate is calculated, defined as the difference between the velocity magnitude of the current iteration and the velocity magnitude of the previous iteration, divided by the velocity magnitude of the previous iteration. The velocity decay rate reflects the rate of velocity change; in gravel areas, due to the greater resistance caused by surface roughness, the velocity decay rate is higher.

[0057] A velocity decay rate threshold is set. When the velocity decay rate at a certain location exceeds this threshold, that location is determined to be a rapidly decaying region. Regions where velocity converges and remains stable are characterized by a negative velocity divergence with a large absolute value, indicating fluid convergence at that location; simultaneously, a small velocity decay rate indicates relatively stable velocity. These regions correspond to rigid structures such as roads and building edges, along which fluids form stable streamlines. By setting combined judgment conditions, locations that simultaneously satisfy a velocity divergence less than a negative threshold and a velocity decay rate less than a positive threshold are identified as streamlined stable regions.

[0058] Based on the identified different region types, the 3D terrain matrix is ​​differentiated to remove debris interference and enhance the rigid structure. For debris regions identified as having divergent or rapidly decaying flow velocities, the pixel gradient values ​​at these locations are set to zero. Setting the gradient values ​​to zero simulates the physical process of sediment being washed away by fluid, eliminating the edge features of these regions in the image and treating them as a flat background in subsequent processing.

[0059] In practice, the corresponding values ​​of these locations in the gradient map are set to zero, or the grayscale values ​​of these locations in the original grayscale image are smoothed to the average value of the surrounding background. After processing, the irregular textures and edges caused by gravel accumulation are suppressed and no longer interfere with the extraction of rigid structures. For streamline regions identified as areas of convergent flow velocity and stable flow, the pixel values ​​at these locations are enhanced to preset values. The pixel value enhancement operation simulates the effect of rigid structures becoming more prominent after fluid erosion, making these areas form obvious highlight features in the image.

[0060] The preset value is typically set to the maximum value of the image's grayscale range, such as 255, to ensure that these locations are preserved as foreground in subsequent binarization processing. The enhancement operation not only increases the pixel values ​​at these locations but also appropriately expands their neighborhood, making the streamlined regions more continuous and complete through morphological dilation. After differential processing, the aerial image, originally mixed with debris interference, is transformed into an enhanced image that highlights rigid structures and suppresses debris noise.

[0061] After a predetermined number of iterations and differentiation processes, the flow field in the 3D terrain matrix has fully evolved, and streamlines form a clear convergence pattern at rigid structures. The streamline convergence center refers to the location where the velocity divergence is negative and its absolute value reaches a local maximum. These locations are the target points where the fluids ultimately converge, corresponding to the skeletal lines of road centerlines or building edges.

[0062] By traversing the enhanced image, pixel locations satisfying the velocity divergence condition are identified, and non-maximum suppression is applied. The location with the largest absolute value of velocity divergence within a local area is retained as the streamline convergence center. These convergence centers constitute the skeleton nodes of the rigid architecture. To generate continuous architectural lines, the convergence centers are connected. Based on the direction of the streamlines and the spatial relationship between the convergence centers, nearest neighbor connection or streamline-direction-based directional connection algorithms are used to connect adjacent convergence centers into line segments.

[0063] During the connection process, distance and direction consistency thresholds can be set to ensure the rationality of the connections. The generated set of line segments constitutes the complete skeleton network of the rigid architecture. Finally, binarization is performed, setting the pixel positions of the skeleton lines to white and the remaining positions to black, generating a binarized original incomplete architecture diagram. This architecture diagram clearly shows the roads and building structures preserved in the post-disaster scene, removing the interference of rubble and debris. However, due to disaster damage, some structures are in a broken state and need to be repaired and completed in subsequent steps.

[0064] This implementation method achieves automatic extraction of rigid structures by constructing a 3D terrain matrix and introducing a virtual fluid dynamics model. Adaptive setting of virtual fluid parameters ensures the method's adaptability to different resolutions and scene features. The iterative flow process utilizes the difference in geometric continuity between rigid structures and gravel deposits, identifying gravel interference regions and rigid structure regions through calculations of physical quantities such as velocity divergence and velocity decay rate. A differentiated processing strategy suppresses gravel noise and enhances rigid structure characteristics, improving the accuracy and robustness of structure extraction. The determination of streamline convergence centers and the generation of a binarized structure diagram provide high-quality input data for subsequent fracture structure repair.

[0065] In one embodiment of this invention, determining virtual fluid parameters and initializing a uniform flow field across the entire three-dimensional terrain matrix includes the following steps: S310: Acquire the resolution and shooting altitude of aerial images; S320. Determine the initial value of the flow velocity attenuation coefficient based on the resolution and shooting height of the aerial image. S330. Determine the initial value of the viscosity coefficient based on the pre-stored characteristics of the gravel density in the post-disaster environment; S340. Set an initial velocity vector at each pixel position of the three-dimensional terrain matrix. The direction of the initial velocity vector points in the gradient descent direction. S350. Normalize the initial velocity vector to initialize a uniform flow field across the entire three-dimensional terrain matrix.

[0066] The resolution and shooting altitude of aerial imagery are key factors affecting the setting of virtual fluid parameters and need to be accurately acquired before initializing the flow field. The resolution of aerial imagery is usually expressed in pixels, representing the width and height of the image, such as 3840×2160 pixels. It is also necessary to obtain the ground sampling distance, i.e., the actual ground size corresponding to each pixel.

[0067] Ground sampling distance is closely related to shooting altitude, camera focal length, and sensor size. Higher shooting altitudes result in greater ground sampling distances, meaning each pixel covers a larger area of ​​the ground, but contains relatively less detail in the image. Shooting altitude information can be directly read from the drone's flight control system, typically measured by a barometer, GPS altimeter, or laser rangefinder. In disaster relief scenarios, the typical flight altitude range for drones is 50 to 200 meters, corresponding to ground sampling distances of approximately 2 to 10 centimeters.

[0068] Resolution and shooting altitude together determine the spatial scale characteristics of an image. High resolution and low shooting altitude mean that the image contains richer terrain details and a higher frequency of terrain undulations. Virtual fluid simulation needs to be performed at a finer scale to accurately distinguish between rigid structures and debris interference. Conversely, low resolution and high shooting altitude make the terrain in the image relatively flat, and fluid simulation can be performed at a coarser scale. By obtaining these two key parameters, necessary input information is provided for the subsequent adaptive determination of the flow velocity attenuation coefficient, ensuring that the virtual fluid parameters can match the actual image features.

[0069] The velocity decay coefficient describes how quickly the velocity of a virtual fluid decays during flow. This parameter needs to be adaptively determined based on the resolution and shooting altitude of the aerial image. Physically, the velocity decay coefficient represents the proportion of velocity lost by the fluid per unit time or unit distance. A larger decay coefficient indicates faster velocity decay and a shorter flow distance.

[0070] At high resolution and low shooting altitude, the terrain details in the image are rich, and the geometric features of rigid structures are small in scale. For example, the width of a road may only be a few dozen pixels, and the details of building edges are also more refined. In order for the virtual fluid to fully explore these detailed features, the fluid needs to flow over a longer distance. Therefore, the flow rate attenuation coefficient should be set to a small value to slow down the attenuation process of the flow rate.

[0071] In specific calculations, an inverse relationship can be established between the velocity attenuation coefficient and the ground sampling distance; the smaller the ground sampling distance, the smaller the velocity attenuation coefficient. The influence of image resolution should also be considered; higher resolution means more pixels, resulting in a larger spatial range that the fluid needs to traverse, and consequently, a smaller velocity attenuation coefficient. An empirical formula can be used to calculate the initial value of the velocity attenuation coefficient. This formula comprehensively considers the effects of ground sampling distance and image resolution, balancing their weights by introducing adjustment parameters.

[0072] The adjustment parameters can be calibrated based on historical data and experimental results to ensure a suitable velocity attenuation coefficient is obtained under different shooting conditions. In practical applications, for aerial images with a shooting height of 100 meters, a ground sampling distance of 5 centimeters, and a resolution of 3840×2160, the typical initial value of the velocity attenuation coefficient is approximately 0.05 to 0.1, indicating that the fluid loses 5% to 10% of its velocity in each iteration. Through this adaptive determination method, the velocity attenuation coefficient can be matched with the spatial scale characteristics of the image, enabling virtual fluid simulation to achieve good rigidity extraction results under different shooting conditions.

[0073] The viscosity coefficient describes the viscosity of a virtual fluid, reflecting the magnitude of internal frictional resistance experienced by the fluid during flow. In post-disaster environments, the density of rubble and debris covering the surface is a key factor influencing the viscosity coefficient setting. The higher the rubble density and the rougher the surface, the greater the resistance experienced by the virtual fluid during flow; a larger viscosity coefficient should be set to simulate this hindering effect.

[0074] The pre-stored post-disaster environmental debris density feature data comes from statistical analysis of historical disaster scenarios. Aerial images of different types of disasters, such as earthquakes, landslides, and building collapses, are analyzed to extract texture features and density distribution patterns in debris-covered areas. Debris density can be quantitatively assessed by analyzing features such as local variance, texture complexity, and edge density. High-density debris areas exhibit high-frequency textures and irregular edges, with larger local variance; low-density debris areas have relatively smooth textures and smaller local variance.

[0075] A debris density feature database is established to store typical debris density values ​​corresponding to different disaster types and severity. When processing new aerial imagery, the texture features of the current image are first analyzed, and the texture complexity index of the entire image or a local area is calculated. Then, it is matched with feature patterns in the database to determine the debris density level of the current scene. Based on the debris density level, the corresponding recommended viscosity coefficient value is retrieved, or an accurate initial viscosity coefficient value is calculated using interpolation methods.

[0076] The viscosity coefficient typically ranges from 0.1 to 0.5, with smaller values ​​for lower gravel density and larger values ​​for higher gravel density. A larger viscosity coefficient causes the virtual fluid to decay rapidly in gravel regions, making it difficult to form stable streamlines and thus distinguishing these regions from rigid structures. Through this adaptive determination method based on pre-stored feature data, the viscosity coefficient can reflect the surface roughness characteristics of the actual post-disaster environment, improving the realism of the virtual fluid simulation and the accuracy of rigid structure extraction.

[0077] Setting the initial velocity vector requires determining the direction and magnitude of the velocity at each pixel location in the 3D terrain matrix. The initial velocity vector should point in the gradient descent direction, simulating the natural flow of fluid from high to low altitudes under gravity. The gradient descent direction is obtained by calculating the gradient vector of the 3D terrain matrix at that location and taking its opposite direction.

[0078] The gradient vector is calculated using numerical differentiation. For a pixel at coordinates (x, y), its partial derivatives in the X and Y directions are calculated. The partial derivative in the X direction is obtained by subtracting the height of the left pixel from the height of the pixel to the right of that position and then dividing by 2. The partial derivative in the Y direction is obtained by subtracting the height of the upper pixel from the height of the pixel below and then dividing by 2. These partial derivatives in the X and Y directions are used as the two components of the gradient vector, forming a two-dimensional gradient vector.

[0079] The gradient vector points in the direction of the fastest height increase, while the flow velocity should point in the direction of the fastest height decrease. Therefore, the opposite direction of the gradient vector is taken as the initial flow velocity direction. For image edge locations, due to the lack of adjacent pixels, the gradient can be calculated using one-sided difference or inward offset. The magnitude of the initial flow velocity can be set to a uniform constant value, such as 1.0, indicating that the initial flow velocity magnitude is the same at all locations, only the direction differs.

[0080] It is important to note that this setup simplifies the initialization process, resulting in a flow field with diverse directions but uniform velocity distribution in the initial state. In some implementations, the initial flow velocity can also be set based on the gradient magnitude; a larger gradient magnitude indicates a steeper terrain slope, allowing for a larger initial flow velocity. However, this approach increases computational complexity. By setting an initial flow velocity vector pointing in the gradient descent direction for each pixel location, the initial conditions for the virtual fluid flow on the 3D terrain surface are established.

[0081] Normalization ensures a consistent flow intensity in the initial state by unifying the magnitudes of all initial velocity vectors. The normalization process iterates through each pixel of the 3D terrain matrix, calculating and scaling the magnitude of the initial velocity vector at that location. The magnitude of the velocity vector is obtained by taking the square root of the sum of the squares of its X and Y components, representing the magnitude of the velocity.

[0082] The normalization operation divides each component of the velocity vector by its magnitude, making the normalized velocity vector have a magnitude of 1, while maintaining the direction. After processing, the velocity magnitude at all pixel locations is uniformly set to unit velocity, with only the velocity direction varying depending on the terrain gradient. The physical significance of normalization is to eliminate differences in initial velocity magnitude, preventing the virtual fluid from exhibiting unreasonable flow patterns due to excessively high or low initial velocities at certain locations during the initial flow phase.

[0083] In subsequent iterative flow processes, the flow velocity is dynamically adjusted based on factors such as terrain features, velocity attenuation coefficient, and viscosity coefficient, forming stable streamlines at rigid structures and rapidly attenuating in gravel areas. Normalization also simplifies the numerical calculation of the flow field, avoiding numerical instability caused by excessively large or small velocity ranges.

[0084] After normalization, the uniform flow field of the entire 3D terrain matrix is ​​initialized, with each pixel having a clear flow direction and uniform flow magnitude, providing a good initial state for the iterative flow of the virtual fluid. This initialization method ensures the uniformity and stability of the flow field, enabling the subsequent flow evolution process to realistically reflect the influence of terrain features on fluid behavior, thereby accurately distinguishing between rigid structures and debris interference.

[0085] This implementation achieves adaptive setting of virtual fluid parameters and precise initialization of a uniform flow field across the entire image through a systematic parameter determination and flow field initialization process. By acquiring the resolution and shooting altitude of aerial imagery and combining it with ground sampling distance calculations, the velocity attenuation coefficient is adaptively determined to match the spatial scale characteristics of the image. Based on pre-stored post-disaster environment debris density characteristic data, the viscosity coefficient is determined to reflect the surface roughness of the actual scene. By calculating the gradient descent direction to set the initial velocity vector and performing normalization processing, a uniform flow field with diverse directions but uniform velocity across the entire image is established, laying the foundation for subsequent iterative flow of the virtual fluid and ensuring the adaptability and robustness of the rigid structure extraction method under different shooting conditions and disaster scenarios.

[0086] In one embodiment of this invention, setting the pixel gradient value of a region with divergent or rapidly decaying flow velocity to zero includes the following steps: S410. Calculate the velocity divergence and velocity decay rate at each pixel location of the three-dimensional terrain matrix; S420. Determine whether the velocity divergence is greater than a preset velocity threshold, or whether the velocity decay rate is greater than a preset gradient threshold. S430. When the velocity divergence is greater than the velocity threshold or the velocity decay rate is greater than the gradient threshold, the gradient value at the corresponding pixel position is set to zero to simulate the effect of sediment being washed away.

[0087] Velocity divergence and velocity decay rate are key physical quantities for determining areas affected by debris interference, requiring precise calculation at each pixel location in the 3D terrain matrix. Velocity divergence describes the degree of dispersion or convergence of fluid at that location and is obtained by calculating the divergence operator of the velocity vector field. For a pixel with location coordinates (x, y), its velocity vector contains X-direction and Y-direction components, and the velocity divergence is equal to the sum of the partial derivatives of the velocity X-component in the X-direction and the partial derivatives of the velocity Y-component in the Y-direction.

[0088] The partial derivatives are calculated using the central difference method. The partial derivative of the velocity X component in the X direction is obtained by subtracting the velocity X component of the left pixel from the velocity X component of the pixel to the right of that position and then dividing by 2. Similarly, the partial derivative of the velocity Y component in the Y direction is obtained by subtracting the velocity Y component of the upper pixel from the velocity Y component of the lower pixel and then dividing by 2. A positive velocity divergence value indicates that the fluid is diverging outward from that position, while a negative value indicates that the fluid is converging towards that position. The larger the absolute value, the stronger the divergence or convergence.

[0089] In areas with gravel accumulation, the irregular terrain prevents fluid from forming a stable flow direction, resulting in significant differences in velocity direction between adjacent locations. This leads to a large absolute value of the velocity divergence and frequent sign changes. The velocity decay rate describes how quickly the velocity changes over time and is calculated by comparing the velocity magnitudes of the current and previous iterations. For each pixel location, the difference between the velocity magnitude of the current iteration and the velocity magnitude of the previous iteration is calculated, and then divided by the velocity magnitude of the previous iteration to obtain the normalized velocity decay rate.

[0090] A positive velocity decay rate indicates a decrease in flow velocity, while a negative rate indicates an increase. A larger absolute value indicates a more drastic change. In areas with gravel, the frictional resistance caused by surface roughness is significant, leading to rapid velocity decay and a high velocity decay rate. By simultaneously calculating these two physical quantities, the spatial distribution and temporal evolution characteristics of the flow field can be comprehensively considered to determine the gravel interference area, improving the accuracy and robustness of identification.

[0091] Determining whether velocity divergence and velocity decay rate exceed preset thresholds is used to identify areas affected by debris interference. The velocity and gradient thresholds need to be determined comprehensively based on the actual scene characteristics and virtual fluid parameters. The velocity threshold is used to determine whether the velocity divergence is significant; its physical meaning is the critical value of the net outflow from a unit area per unit time.

[0092] Setting the velocity threshold requires considering the initial velocity and velocity decay coefficient of the virtual fluid. The higher the initial velocity or the lower the velocity decay coefficient, the higher the velocity threshold should be. In typical applications, the velocity threshold can be set to 10% to 30% of the initial velocity. For example, when the initial velocity is normalized to 1.0, the velocity threshold can be set to 0.1 to 0.3. The gradient threshold is used to determine whether the velocity decay rate is significant; its physical meaning is the critical proportion of velocity change.

[0093] Setting the gradient threshold requires consideration of the viscosity coefficient and the number of iterations. A higher viscosity coefficient results in faster velocity decay, necessitating a correspondingly larger gradient threshold. In typical applications, the gradient threshold can be set between 20% and 50%, indicating that if the velocity decays beyond this percentage in a single iteration, it is considered a rapid decay. The judgment process uses a logical OR operation: if the absolute value of the velocity divergence is greater than the velocity threshold, or the velocity decay rate is greater than the gradient threshold, the location is identified as a gravel interference area. This judgment strategy has strong fault tolerance and can capture different types of gravel interference patterns, including areas with significant velocity divergence and areas with rapid velocity decay.

[0094] Setting the pixel gradient values ​​of areas with gravel interference to zero is used to remove gravel interference, simulating the physical effect of the surface becoming flat after sediment is washed away by fluid. For pixel locations identified as gravel interference areas, the corresponding location in the gradient map is first located. The gradient map stores the gradient magnitude of each location in the 3D terrain matrix, reflecting the terrain slope at that location.

[0095] Setting the gradient magnitude to zero at the corresponding location in the gradient map indicates that the terrain at that location becomes flat and no longer has edge features. The gradient zeroing operation can be applied to both the X and Y components of the gradient, making the gradient zero in both directions at that location. In some implementations, to avoid abrupt boundaries caused by the gradient zeroing operation, the zeroing region can be smoothed using Gaussian filtering or mean filtering to smooth the zeroing region and its neighborhood, allowing the gradient value to gradually transition from zero to the surrounding non-zero values, forming a more natural boundary.

[0096] Setting the gradient value to zero eliminates the edge response of gravel areas in the image, causing these areas to be treated as background during subsequent edge detection or architecture extraction, and no longer participating in the construction of rigid structures. This operation suppresses irregular textures and false edges caused by gravel accumulation, preventing these interfering features from being misidentified as road or building structures.

[0097] Meanwhile, by preserving the gradient values ​​of rigid structural regions where flow velocity converges and is stable, and even enhancing the pixel values ​​of these regions, a dual effect of interference suppression and feature enhancement is achieved, improving the signal-to-noise ratio and accuracy of rigid architecture extraction, and laying a high-quality data foundation for subsequent architecture repair and topology matching.

[0098] This implementation method comprehensively identifies debris interference regions by calculating two key physical quantities: velocity divergence and velocity decay rate, considering both the spatial distribution and temporal evolution of the flow field. By setting reasonable velocity and gradient thresholds and employing a logical OR operation-based judgment strategy, accurate identification of debris interference regions is achieved. Setting the pixel gradient values ​​of the identified debris regions to zero simulates the physical effect of sediment being washed away, eliminating irregular textures and false edge interference caused by debris accumulation, and improving the accuracy and robustness of rigid structure extraction.

[0099] In one embodiment of this invention, the original incomplete architecture diagram is repaired and completed to obtain a topology diagram from the perspective of the UAV, including the following steps: S510. Identify all line break points in the original incomplete architecture diagram and calculate the tangent direction vector for each line break point; S520. Construct a virtual stress field between each pair of adjacent line breaks and calculate the probability score of connection between each pair of adjacent line breaks based on the virtual stress field. S530. Set a growth threshold based on the probability score, and generate connecting lines for line breakpoints with probability scores higher than the growth threshold. S540: Merge the generated connection lines with the original incomplete architecture diagram to output the repaired UAV view topology diagram.

[0100] After rigid structure extraction, the original incomplete architectural diagram, while removing debris interference, still exhibits discontinuous line fragments due to road and building structural breaks caused by the disaster. Identifying line breakpoints requires traversing all white pixels in the original incomplete architectural diagram and determining whether each pixel is an endpoint of a line.

[0101] The determination of line breakpoints is based on pixel connectivity analysis. For the binarized architecture graph, an 8-connected neighborhood detection method is used to check the pixel values ​​of the target pixel's eight neighboring locations. If the target pixel is white, it indicates that it belongs to an architecture line, and only one of its eight neighbors is white; therefore, this pixel is determined to be a line breakpoint. This situation indicates that the line terminates at that position and does not continue in other directions. By traversing the entire architecture graph, all line breakpoint locations can be identified, and their coordinate information can be recorded.

[0102] After identifying the breakpoints, it is necessary to calculate the tangent direction vector for each breakpoint. This vector represents the direction of line extension at the breakpoint and is an important basis for subsequently determining the possibility of connection between breakpoints. The tangent direction is calculated using a local line fitting method. Starting from the breakpoint, the line is traced inwards at a certain distance, such as 5 to 10 pixels, to extract the coordinates of all pixels on this line segment.

[0103] A least-squares line fitting is performed on these coordinate points to obtain the slope and intercept of the fitted line. The direction of the fitted line is the tangent direction. The tangent direction is represented as a unit vector, with the X and Y components calculated from the slope and normalized to a magnitude of 1. The tangent direction vector points in the direction the line extends, reflecting the geometric orientation of the structure at the break point.

[0104] It is important to note that for certain special cases, such as lines near the breakpoint being too short or sharply curved, weighted fitting or curve fitting methods can be used to improve the accuracy of tangent direction estimation. By calculating the tangent direction vector for each breakpoint, a geometric feature description of the breakpoint is established, providing necessary input information for subsequent virtual stress field construction and connection probability assessment.

[0105] Virtual stress fields are physical field models that simulate the tendency of structural connections. By comprehensively considering multiple factors such as spatial distance, directional consistency, and texture residue, they quantify the likelihood of connection between each pair of adjacent breakpoints. First, adjacent breakpoint pairs need to be identified. All identified breakpoints are traversed, and for each breakpoint, other breakpoints within a certain radius are searched. For example, the search radius can be set to 50 to 200 pixels, with the specific value determined based on scene features such as image resolution and road width. Breakpoints within the search range are considered adjacent breakpoints and form a breakpoint pair with the current breakpoint.

[0106] For each pair of adjacent breakpoints, a virtual stress field is constructed and a connection probability score is calculated. The virtual stress field consists of three key parameters. The first parameter is the Euclidean distance, representing the straight-line distance between the two breakpoints, which is obtained by calculating the square root of the sum of the squares of the coordinate differences between the two breakpoints. The closer the distance, the greater the probability that the two breakpoints originally belonged to the same structure, and the stronger the tendency to connect. The distance factor is normalized using a negative exponential decay function, so that breakpoint pairs that are farther apart receive lower weights.

[0107] The second parameter is the cosine of the angle between the tangent direction vectors, obtained by calculating the dot product of the tangent direction vectors at the two breakpoints. The cosine value ranges from -1 to 1. The closer the value is to 1, the more consistent the two tangent directions are, meaning the extension directions of the two breakpoints tend to be parallel or collinear, resulting in a smooth line when connected. A cosine value closer to -1 indicates that the two tangent directions are opposite, making them unsuitable for connection. The direction consistency factor can directly use the cosine value or undergo a nonlinear transformation to enhance the weight of highly consistent regions.

[0108] The third parameter is the original texture residue, which assesses whether there are residual traces of roads or buildings in the area between two breakpoints by analyzing the texture features of the original aerial image. Specifically, the pixel grayscale values ​​of the area connecting the two breakpoints are extracted, and the local texture variance or edge density of the area is calculated. A smaller texture variance and lower edge density indicate that the area is relatively flat, possibly a road surface, and has a higher texture residue.

[0109] The three parameters are weighted and summed according to preset weights to obtain the final probability score. The weight coefficients can be adjusted according to experimental data and scene characteristics. A typical weight allocation is 40% for distance factor, 40% for orientation consistency factor, and 20% for texture persistence factor. Through this multi-factor comprehensive evaluation method, the probability score can comprehensively reflect the rationality of the breakpoint connection, taking into account both geometric relationships and texture information of the original image.

[0110] A growth threshold is set based on the probability score to filter valid connections. The setting of the growth threshold needs to balance the integrity and accuracy of the connections. If the threshold is set too low, too many false connections will be made, connecting irrelevant breakpoints incorrectly and creating false architectural lines; if the threshold is set too high, real connections will be missed, leaving broken structures that should be repaired still broken.

[0111] The growth threshold can be determined by analyzing the probability score distribution of all breakpoint pairs. The mean and standard deviation of the scores are calculated, and the threshold is set to the mean plus a certain multiple of the standard deviation, such as the mean plus 0.5 times the standard deviation. This allows for the screening of breakpoint pairs with scores significantly higher than the average level. In some implementations, a fixed threshold can also be used, for example, setting the threshold to 0.6 to 0.8, indicating that only breakpoint pairs with probability scores exceeding this value are considered to be connected.

[0112] For breakpoint pairs with a probability score higher than the growth threshold, connectors need to be generated to complete the broken architecture. Connector generation is not a simple straight-line connection; instead, the curvature parameters of the connector are calculated based on the tangent direction vectors of the two breakpoints, ensuring that the generated connectors smoothly connect the lines at both ends, conforming to the natural direction of the road or building edges.

[0113] The specific generation method uses Bézier curves or cubic spline interpolation. First, the coordinates of the start and end points of the connecting line are determined, which are the positions of the two breakpoints. Then, the coordinates of the control points are calculated based on the tangent direction vectors of the two breakpoints. The control points are located between the start and end points, and their positions are determined by the tangent direction and the breakpoint spacing. For cubic Bézier curves, two control points are required. The first control point is located at the start point, extending a certain distance along the tangent direction, and the second control point is located at the end point, extending a certain distance in the opposite direction of the tangent direction. The extension distance is usually set to one-third of the breakpoint spacing.

[0114] Based on the start point, end point, and control points, a smooth connection curve is generated using the Bézier curve formula. The coordinates of points on the curve are calculated using parametric equations, with parameters varying from 0 to 1, corresponding to the curve's movement from the start point to the end point. The tangent directions of the generated connection lines at the start and end points are consistent with the tangent directions at the break points, ensuring the smoothness of the connection. The generated connection lines are rasterized into pixel coordinates, and the pixels along the connection path are set to white in the architecture diagram, completing the repair of the broken architecture.

[0115] The generated connectors are then merged with the original incomplete architecture diagram to output a complete topology diagram. The merging process requires adding all generated connectors to the original architecture diagram while maintaining the integrity of the original architecture lines. Specifically, a new image of the same size as the original incomplete architecture diagram is created. First, all pixels of the original architecture diagram are copied to the new image, preserving the original architecture lines.

[0116] Then, iterate through all the generated connectors. For each connector, set the pixels along its path to white in the new image. If some pixels along the connector path are already white, it means there are already structural lines at that location, so leave them unchanged; if pixels are black, it means that the location was originally background, so set them to white to add the connector. To ensure a smooth transition between the connector and the original lines, morphological processing can be performed on the area near the endpoints of the connector, using a dilation operation to fully blend the connector with the original lines at the endpoints, avoiding breaks or discontinuities.

[0117] After fusion, the resulting image contains all the lines of the original architecture diagram as well as the generated connecting lines, forming a repaired and complete architecture network. To further optimize the topology, the fused image can be refined by using a skeleton extraction algorithm to refine the lines into single-pixel-width skeleton lines, eliminating uneven line thickness and making the architecture network more regular.

[0118] The refined image is the topology map from the drone's perspective. This map clearly shows the topological connections of roads and building structures in the post-disaster scene. Broken structures have been reasonably repaired and completed, providing high-quality input data for subsequent resilient topology matching. The output topology map not only retains the reliable structures identified in the rigid architecture extraction stage, but also repairs the broken parts caused by disaster damage through intelligent connections, restoring the integrity and continuity of the structure and improving the success rate and accuracy of subsequent matching and positioning.

[0119] This implementation method achieves intelligent completion of the original incomplete architecture diagram through a systematic fracture architecture repair process. By identifying line breakpoints and calculating tangent direction vectors, a geometric feature description of the breakpoints is established. Based on a virtual stress field model, multiple factors such as spatial distance, directional consistency, and texture residue are comprehensively considered to quantitatively evaluate the connection probability of breakpoint pairs. Valid connections are screened by setting a growth threshold, and smooth connection lines are generated using Bézier curves to ensure the rationality and naturalness of the repair. The connection lines are fused with the original architecture diagram to output a complete topology diagram from a UAV perspective, laying the foundation for subsequent flexible topology matching.

[0120] In one embodiment of this example, a virtual stress field is constructed between each pair of adjacent line breaks, and a probability score for connection between each pair of adjacent line breaks is calculated based on the virtual stress field, including the following steps: S610. Calculate the Euclidean distance between each pair of adjacent breakpoints; S620. Calculate the cosine of the angle between the tangent direction vectors of each pair of adjacent breakpoints. S630. Extract the original texture residue of pixels between each pair of adjacent breakpoints, where the Euclidean distance, the cosine of the included angle, and the original texture residue constitute a virtual stress field. S640. Based on the Euclidean distance, the cosine of the included angle, and the original texture residual, calculate the probability score of connection between each pair of adjacent line breakpoints by weighted summation.

[0121] Euclidean distance is a fundamental geometric quantity for measuring the spatial proximity of two discontinuities, directly reflecting the spatial feasibility of connecting discontinuity pairs. For each pair of adjacent discontinuities, let the coordinates of the first discontinuity be... The coordinates of the second breakpoint are The Euclidean distance is obtained by taking the square root of the sum of the squares of the differences between the coordinates of the two points.

[0122] In actual calculations, the coordinate differences in the X and Y directions are first calculated, then squared and summed, and finally the square root of the sum is taken to obtain the Euclidean distance in pixels. The physical meaning of Euclidean distance is the straight-line distance between two breakpoints; the smaller the distance, the closer the two breakpoints are in space, and the greater the probability that they originally belonged to the same structure.

[0123] In post-disaster scenarios, road or building structural fractures are typically localized, with relatively small gaps, ranging from a few pixels to tens of pixels. Therefore, fracture pairs with smaller Euclidean distances are more likely to be the two ends of the same structure fractured due to disaster damage and should be prioritized for connection.

[0124] To convert the Euclidean distance into a normalized factor that can be used for weighted summation, a nonlinear transformation of the distance values ​​is required. A negative exponential decay function is used for this transformation; the distance factor is equal to the value of an exponential function with the natural constant e as the base and the negative distance divided by the characteristic distance as the exponent. The characteristic distance is an adjustment parameter that reflects the rate of distance decay; the larger the characteristic distance, the slower the distance factor decays with increasing distance, allowing for a higher probability of connection between more distant discontinuities.

[0125] The feature distance needs to be determined based on the image resolution and typical break gaps, and is usually set to one-third to one-half of the expected maximum connection distance. After a negative exponential transformation, the distance factor ranges from 0 to 1, with the factor being 1 when the distance is zero and approaching 0 when the distance is infinite, forming a smooth decay curve. This transformation method preserves distance information and makes it suitable for weighted combination with other factors, providing reasonable distance weights for the calculation of probability scores.

[0126] The cosine of the angle between the tangent direction vectors reflects the consistency of the extension directions of the two breakpoints and is a key geometric feature for judging the rationality of the connection. For each pair of adjacent breakpoints, their respective tangent direction vectors have been calculated in the preceding steps. Let the tangent direction vector of the first breakpoint be... The tangent direction vector at the second break point is Both vectors are unit vectors.

[0127] The cosine of the included angle is obtained by calculating the dot product of the two vectors. The dot product is equal to the sum of the products of their corresponding components. Multiply Plus Multiply Since both vectors have been normalized to unit vectors, the result of the dot product is directly equal to the cosine of the angle between the two vectors, which ranges from -1 to 1. A cosine of 1 indicates that the two vectors have exactly the same direction, meaning that the extension directions of the two breakpoints are parallel and point in the same direction, and connecting them can form a smooth straight line or a gentle curve.

[0128] A cosine value of 0 indicates that the two vectors are perpendicular, their extension directions differ significantly, and the connection is less likely to be valid. A cosine value of -1 indicates that the two vectors are in completely opposite directions, meaning the two breakpoints extend towards each other. This is reasonable in some scenarios, such as the two ends of a road breaking towards each other. To distinguish between these two cases, the relative positional relationship of the breakpoints can be further analyzed to determine whether the two breakpoints are located ahead of their respective tangent extension directions.

[0129] It's important to note that if the second breakpoint is located ahead of the tangent extension direction of the first breakpoint, and the first breakpoint is located ahead of the tangent extension direction of the second breakpoint, then even if the cosine of the included angle is negative, it may still be a valid connection. In practical applications, the directional consistency factor can directly use the cosine of the included angle, or it can be nonlinearly transformed to enhance the weight of highly consistent regions. For example, the square or cube of the cosine value can be used, giving higher weights to breakpoint pairs with highly consistent directions, while rapidly reducing the weights of breakpoint pairs with significantly different directions. By calculating the cosine of the included angle, the degree of matching between breakpoint pairs in geometric directions is quantified, providing crucial directional information for connection probability scoring.

[0130] Original texture residue reflects the texture features of the region between two breakpoints in the original aerial image, and is used to assess whether there are residual traces of roads or buildings in that area. Extracting texture residue requires going back to the original aerial image before image preprocessing and analyzing the pixel grayscale distribution characteristics of the line region between the two breakpoints.

[0131] First, the path connecting the two breakpoints is determined. A pixel sequence from the first breakpoint to the second breakpoint is generated using linear interpolation; this sequence constitutes the region to be analyzed. For each pixel on the path, its grayscale value is read from the original aerial image, forming a one-dimensional grayscale sequence. Texture persistence is calculated based on the statistical characteristics of this grayscale sequence, mainly including two indicators: local variance and edge density.

[0132] Local variance reflects the degree of fluctuation in grayscale values. It is calculated by taking the mean of the grayscale sequence, then squared the difference between each grayscale value and the mean, and averaging all squared values ​​to obtain the variance. A smaller variance indicates a more uniform grayscale distribution, suggesting the area is more likely to be a flat surface like a road, with higher texture retention. A larger variance indicates more drastic grayscale fluctuations, suggesting the area may be covered by gravel or other debris, with lower texture retention.

[0133] Edge density reflects the intensity of edge features in a region. It is calculated by performing gradient calculations on the grayscale sequence and counting the proportion of pixels with gradient magnitudes exceeding a certain threshold to the total number of pixels. Low edge density indicates that the region lacks obvious edge features, possibly a flat road surface with high texture retention. High edge density indicates that the region has a large number of edges, possibly gravel accumulation or complex textures with low texture retention.

[0134] The reciprocals of the local variance and the edge density are normalized, and then averaged or weighted to obtain the final texture residual factor, ranging from 0 to 1. A texture residual factor closer to 1 indicates that the area is more likely to be the original road or building surface, and the higher the plausibility of the connection. By extracting the original texture residual, the texture information of the original image is introduced, compensating for potential judgment biases that may arise from relying solely on geometric features, making the assessment of connection probability more comprehensive and accurate.

[0135] Based on three key parameters of the virtual stress field, a connection probability score is calculated through weighted summation to comprehensively evaluate the rationality of the breakpoints for connection. The probability score is calculated by multiplying the distance factor by the distance weight, adding the orientation consistency factor by the orientation weight, and adding the texture residual factor by the texture weight. The sum of the three weight coefficients is 1 to ensure score normalization.

[0136] The weighting coefficients need to be set based on a balance between the actual scene characteristics and the importance of each factor. In typical applications, the distance factor and orientation consistency factor are relatively important because they directly reflect the geometric rationality of the connection, and can be set to 0.4 and 0.4 respectively. The texture persistence factor provides auxiliary information, and its weight can be set to 0.2. In some scenarios, if the texture information of the original image is well preserved, the texture weight can be appropriately increased; if the original image is severely damaged or the texture information is unreliable, the texture weight can be decreased and the weight of the geometric factor can be increased.

[0137] By weighted summation, the contributions of the three factors are combined into a unified score, ranging from 0 to 1, with higher scores indicating a greater likelihood of connection. High-scoring breakpoint pairs demonstrate strong connection tendencies across multiple dimensions, including spatial distance, extension direction, and texture features, and should be prioritized for connection. Low-scoring breakpoint pairs indicate significant mismatches in at least one dimension, making the connection less plausible, and should be excluded or given lower priority.

[0138] This multi-factor weighted evaluation method can comprehensively reflect the rationality of the breakpoints in connection, taking into account both the constraints of geometric relationships and the texture information of the original image. This improves the accuracy and robustness of connection judgment and provides a reliable decision-making basis for subsequent growth threshold setting and connection line generation.

[0139] This implementation systematically evaluates the connectivity probability of breakpoint pairs by constructing a virtual stress field model. The impact of spatial proximity on connectivity is quantified by calculating Euclidean distance and performing a negative exponential transformation. The constraint of consistent extension directions on connectivity is quantified by calculating the cosine of the angle between tangent direction vectors. Image texture information is introduced as an auxiliary judgment criterion by extracting original texture residuals. The three factors are weighted and summed to form a unified probability score, comprehensively reflecting the rationality of connectivity, improving the accuracy of broken structure repair, and laying the foundation for generating high-quality UAV-view topology maps.

[0140] In one embodiment of this example, generating connecting lines for line breakpoint pairs with probability scores higher than the growth threshold includes the following steps: S710. For line breakpoints with a probability score higher than the growth threshold, determine the starting and ending coordinates of the connecting line. S720. Calculate the curvature parameter of the connecting line based on the tangent direction vectors of the two line break points; S730, generates connection lines based on curvature parameters.

[0141] For line breakpoint pairs with a probability score higher than the growth threshold, connecting lines need to be generated to complete the broken structure. First, the start and end positions of the connecting lines must be determined. The start and end coordinates correspond to the pixel coordinates of the two breakpoints in the breakpoint pair, respectively.

[0142] In practice, all breakpoint pairs that have passed the scoring and filtering are traversed. For each breakpoint pair, the coordinates of the first breakpoint are read as the starting coordinates of the connecting line, denoted as . Read the coordinates of the second breakpoint as the coordinates of the endpoint of the connecting line, denoted as . These two coordinate points correspond to the endpoints of the lines in the original incomplete structural diagram, representing the two ends of the broken structure.

[0143] After determining the start and end points, the directionality of the connection also needs to be considered, i.e., from which breakpoint to which breakpoint the connecting line will be generated. Although the connection is geometrically bidirectional, a clear order of start and end points needs to be specified during the actual generation process for subsequent curve parameterization calculations. The order can be determined based on the coordinate relationship of the breakpoints. For example, the breakpoint with the smaller X-coordinate can be chosen as the start point first, and if the X-coordinates are the same, the breakpoint with the smaller Y-coordinate can be chosen as the start point. This regularized order setting helps maintain the consistency of the algorithm.

[0144] In some implementations, the natural direction of the connection can be determined based on the tangent direction vector at the break point, and the order in which the connecting lines extend along the tangent direction can be selected. After determining the coordinates of the start and end points, these two coordinate points become the boundary constraints for subsequent curve generation. The connecting lines must start from the start point, follow a smooth path to the end point, thereby achieving the repair of the broken structure. The accurate determination of the start and end points is the basis for the generation of connecting lines, ensuring that the repaired structure can accurately connect to the original line endpoints, maintaining the continuity and integrity of the structure.

[0145] The curvature parameter determines the degree of curvature and shape characteristics of the connecting line. It needs to be calculated based on the tangent direction vectors of the two breakpoints to ensure that the generated connecting line smoothly connects the lines at both ends, conforming to the natural direction of the road or building edge. The tangent direction vector of each breakpoint has already been calculated in the preceding steps. Let the tangent direction vector of the starting point be... The tangent direction vector at the endpoint is These two vectors represent the direction of the line's extension at the starting and ending points, respectively.

[0146] The calculation of curvature parameters is based on the control point determination method for Bézier curves or spline curves. For cubic Bézier curves, two control points are needed to define the shape of the curve. The first control point is located at a certain distance extending from the starting point along its tangent direction; this distance is called the control distance, and the magnitude of the control distance directly affects the curvature of the curve. The control distance is usually set as a certain proportion of the Euclidean distance between the starting and ending points, with a typical proportionality coefficient of one-third to one-half.

[0147] The coordinates of the first control point are calculated by adding the starting coordinates to the tangent direction vector multiplied by the control distance, i.e. Where d is the control distance. The second control point is located at the position of the same control distance extending from the endpoint along its tangent in the opposite direction. The opposite direction of the tangent is obtained by taking the negative of the tangent direction vector. The coordinates of the second control point are calculated as follows: .

[0148] The two control points determined in this way ensure that the tangent direction of the Bézier curve at the starting point is consistent with the tangent direction vector at the starting point, and the tangent direction at the ending point is consistent with the tangent direction vector at the ending point, thus ensuring a smooth connection between the connecting line and the original line. The scaling factor of the control distance can be adaptively adjusted according to the distance between the breakpoints and the angle between the tangent directions. When the two tangent directions are at the same height, a larger scaling factor can be used to generate a smoother curve; when the tangent directions differ significantly, a smaller scaling factor can be used to generate a more curved curve to adapt to changes in direction.

[0149] Accurate calculation of curvature parameters is key to generating high-quality connectors. By setting the control point positions appropriately, the connectors can satisfy the tangent constraints at the endpoints and form a natural and smooth curve shape, avoiding sharp turns or unnatural bends.

[0150] Generating connecting lines based on curvature parameters is the process of converting mathematical curves into image pixel representations. This requires calculating the coordinates of discrete points on the curve using curve parameterization methods and then drawing the connecting lines in the architecture diagram. For a cubic Bézier curve, its parametric equation is that the coordinates of points on the curve are functions of the parameter t, where t ranges from 0 to 1. t = 0 corresponds to the starting point, and t = 1 corresponds to the ending point.

[0151] The parametric equation of a Bézier curve is expressed as follows: the x-coordinate of a point on the curve is equal to the x-coordinate of the starting point multiplied by... Add the X coordinate of the first control point multiplied by Add the X coordinate of the second control point multiplied by Add the endpoint X coordinate multiplied by The calculation method for the Y coordinate is the same; simply replace the X coordinate with the Y coordinate.

[0152] By gradually increasing the parameter t from 0 to 1 and sampling at a certain step size, such as 0.01, the coordinates of a series of discrete points on the curve can be calculated. A smaller step size results in denser sampling points and a smoother curve, but also increases the computational load. In practical applications, the step size can be adaptively determined based on the distance between the starting and ending points; the longer the distance, the larger the step size can be to control the total number of sampling points.

[0153] After calculating the coordinates of discrete points on the curve, these floating-point coordinates need to be converted into integer pixel coordinates through rounding or integer operations. Then, a linear interpolation method is used to connect adjacent discrete points. Using the Bresenham line algorithm or the DDA algorithm, a continuous pixel sequence is generated between every two adjacent points, ensuring that the connecting line has no breaks in the image.

[0154] All pixel sequences are merged to form a complete connecting line path. In the architecture diagram, all pixels on the path are set to white to complete the drawing of the connecting line. The generated connecting line smoothly connects with the original line at the start and end points, and the middle part presents a natural curved shape, which not only conforms to geometric constraints but also has a good visual effect.

[0155] It is important to note that in some implementations, to further optimize the quality of the connectors, post-processing can be performed on the generated connectors. For example, morphological closing operations can be used to fill any possible micro-gaps, or curve smoothing algorithms can be used to further smooth the connectors and remove any jagged edges or irregularities. By using a connector generation method based on curvature parameters, high-quality repair of fractured architectures is achieved. The generated connectors are not only geometrically sound but also visually natural, providing complete and accurate structural information for subsequent topology matching.

[0156] This implementation achieves precise connection of breakpoint pairs through a systematic connector generation process. By determining the coordinates of the start and end points, the boundary constraints of the connector lines are clarified. By calculating the curvature parameters based on the tangent direction vector, the control point positions of the Bézier curves are determined, ensuring a smooth transition between the connector lines and the original lines at the endpoints. Precise drawing of the connector lines is achieved by generating Bézier curves using parametric methods and rasterizing them into pixel representations. The generated connector lines satisfy geometric constraints and possess a natural, smooth shape, repairing broken architectural structures, improving the quality of architectural repair, and providing technical support for constructing a complete UAV-view topology map.

[0157] In one embodiment of this invention, a flexible topology mapping matching is performed between the UAV-view topology map and a pre-stored satellite reference map to obtain the matching result, including the following steps: S810: Based on satellite baseline images, extract standard roads and standard building structures, and construct a set of baseline gravity sources based on standard roads and standard building structures; S820. Define the nodes in the topology diagram from the perspective of the UAV as a set of mass points to be matched. S830. Project the set of particles to be matched onto the coordinate system where the reference gravitational source set is located to obtain the initial rough position; S840. Iteratively calculate the virtual gravitational resultant force of all reference gravitational sources on each particle in the set of particles to be matched, and move the position of the particle according to the virtual gravitational resultant force. S850: During the movement of the mass point, the total potential energy change rate of the topology map from the perspective of the UAV is calculated in real time. When the total potential energy change rate is lower than the preset value, the iteration is stopped and the final position of the mass point to be matched is output as the matching result.

[0158] In this embodiment, iteratively calculating the virtual gravitational resultant force exerted on each particle in the set of particles to be matched by all reference gravitational sources, and moving the particle's position according to the virtual gravitational resultant force, includes: Mass parameters are set for each reference gravity source, where the mass parameters of intersection nodes are greater than those of ordinary road nodes; Calculate the distance vector between each particle to be matched and all reference gravitational sources; Based on the distance vector and mass parameters, calculate the magnitude and direction of the virtual gravitational force on each particle to be matched; The resultant force vector is obtained by vector superposition of all virtual gravitational forces based on their magnitude and direction. The position coordinates of the particle to be matched are updated based on the resultant force vector and the preset movement step size, so as to move the position of the particle.

[0159] Reference Figure 2 It illustrates a schematic diagram of the steady-state determination process for the topology structure from the UAV's perspective based on the rate of change of potential energy. During the movement of the particle's position, the total rate of change of potential energy in the topology structure from the UAV's perspective is calculated in real time, and the iteration stops when the total rate of change of potential energy falls below a preset value. This includes: Calculate the total potential energy of all unmatched particles in the topology graph from the UAV's perspective in the current iteration; Calculate the rate of change of total potential energy relative to the previous iteration step; Determine whether the rate of change of total potential energy is lower than the preset value; If the rate of change is lower than a preset small value, the topology diagram from the drone's perspective is determined to have reached a steady state, and the iteration stops. If the rate of change is not lower than the preset value, continue to execute the next iteration.

[0160] Based on the matching results, pose calculation is performed, and the precise coordinate information of the UAV in the world coordinate system is output, including: Based on the mapping relationship between the final coordinates of the particle to be matched in steady state and the coordinates of the original image in the matching results, the homography matrix is ​​calculated; The camera's rotation matrix and translation vector are solved by inversely solving the homography matrix; Based on the rotation matrix and translation vector, and combined with the UAV's flight altitude information, calculate the latitude and longitude coordinates of the UAV in the world coordinate system; Output latitude and longitude coordinates and altitude information as the precise coordinate information of the UAV.

[0161] Satellite baseline images are pre-stored pre-disaster surface images containing complete road networks and building distribution information, serving as a reference benchmark for resilient topology matching. Extracting standard roads and standard building structures is a prerequisite for constructing a baseline gravity source set, requiring systematic feature extraction and structured representation of the satellite baseline images.

[0162] Satellite reference images are typically derived from high-resolution satellite imagery or aerial mapping data, possessing high spatial accuracy and complete georeferenced information. Standard road extraction employs road detection algorithms. First, the satellite reference imagery is preprocessed, including grayscale conversion, contrast enhancement, and noise filtering. Then, edge detection algorithms are used to extract road edges, and Hough transform or road tracing algorithms are used to identify road centerlines, forming a vector representation of the road network. The road network includes road segments and intersection nodes. Road segments represent continuous roads between two intersections, and intersection nodes represent the junctions of multiple roads.

[0163] The extraction of standard building structures employs building detection algorithms. By analyzing the texture, shadow, and geometric features of satellite imagery, the outlines and locations of buildings are identified, and key feature points such as corners, edges, and center points are extracted. The extracted standard roads and building structures are stored as vector data, with each feature point or node having precise world coordinates, typically represented by latitude and longitude or projected coordinates.

[0164] Based on the extracted standard road and building structures, a set of benchmark gravity sources is constructed. These benchmark gravity sources are the gravitational emission points in the virtual gravitational field model, representing key nodes in the Earth's surface structure. The construction process first identifies all intersection nodes in the road network. Intersection nodes are the junctions of multiple roads, possessing significant topological meaning and playing a crucial role in navigation and positioning.

[0165] All intersection nodes are added to the baseline gravity source set, and a high mass parameter is assigned to each node, reflecting its importance in the topology. The mass parameter of arterial road intersections is greater than that of secondary arterial road intersections, and the mass parameter of complex intersections with multiple roads is greater than that of simple crossroads. In addition to intersection nodes, ordinary road nodes are also set up along the road segments. These nodes are evenly distributed along the road centerline, with the spacing determined based on the road length and curvature characteristics, typically ranging from 20 to 50 meters.

[0166] The mass parameters of ordinary road nodes are set relatively small, but they provide constraint information on road orientation, helping to guide the mass points to be matched along the correct road direction. For building structures, corner points and important edge points of the buildings are added to the baseline gravity source set. These points provide geometric constraints on the buildings, and the mass parameters are set according to the size and importance of the buildings. The baseline gravity source set constructed in this way includes both topologically critical nodes and geometrically constrained points, forming a complete gravitational field distribution.

[0167] The topology map from the UAV's perspective is a structured representation obtained after extracting rigid structures and repairing fractured structures. It includes the topological connections of roads and building structures in the post-disaster scenario. Defining the nodes in the UAV's topology map as a set of mass points to be matched is a key step in establishing a virtual gravitational field matching model.

[0168] Node extraction employs topological analysis, traversing all lines in the topological structure map from the UAV's perspective to identify line intersections, endpoints, and turning points. Intersections are the points where multiple lines meet, corresponding to road intersections or building edges, and have significant topological meaning. Endpoints are the terminating points of lines, potentially corresponding to the ends of roads or corners of building edges. Turning points are locations where the direction of a line changes significantly, reflecting road curvature or building corners. Through topological analysis, all these key nodes can be identified, and the image coordinates of each node can be recorded.

[0169] All identified nodes are defined as points to be matched, each point having position coordinates and topological connectivity. Initially, the position coordinates are the pixel coordinates of the node in the topological structure map from the UAV's perspective, represented as... , where i is the index of the mass point. Topological connectivity describes the adjacency relationships between mass points. By analyzing the connection of lines, it can be determined which mass points are directly connected by lines. These connections constitute the skeleton of the topological structure graph.

[0170] The set of particles to be matched not only includes the spatial location information of the nodes but also preserves the topological connections between them. This structured representation allows for consideration of not only the positional matching of individual nodes but also the consistency of the overall topology during the matching process. To facilitate subsequent gravitational calculations and position updates, a unit mass is assigned to each particle to be matched, representing its response capability in the virtual gravitational field. By defining the nodes of the topological structure map from the UAV's perspective as the set of particles to be matched, a mapping from image features to the physical model is established.

[0171] The set of particles to be matched is initially located in the image coordinate system of the topological structure map from the UAV's perspective, while the set of reference gravitational sources is located in the world coordinate system of the satellite reference map; the two are in different coordinate systems. In order to perform matching, the set of particles to be matched needs to be projected onto the world coordinate system where the set of reference gravitational sources is located to obtain an initial coarse position.

[0172] The projection process requires establishing an initial mapping between the image coordinate system and the world coordinate system. This mapping is estimated based on the UAV's flight parameters and camera parameters. The UAV's flight parameters include flight altitude, flight attitude, and approximate geographical location, which can be obtained from the UAV's flight control system or inertial measurement unit. Camera parameters include focal length, principal point position, and distortion coefficients, which are predetermined through camera calibration.

[0173] The initial projection uses a simplified perspective projection model, assuming the ground is flat and the drone camera shoots vertically downwards. Based on the flight altitude and camera focal length, the ground coverage area and ground sampling distance of the image can be calculated. The ground sampling distance is equal to the flight altitude multiplied by the sensor pixel size divided by the focal length, representing the actual ground size corresponding to each pixel.

[0174] Based on the ground sampling distance, the image coordinates can be converted to a ground distance relative to the image center. Then, combined with the approximate geographical location of the UAV, the ground distance can be converted to latitude and longitude or projected coordinates in the world coordinate system. In the specific calculation, let the world coordinates corresponding to the image center be... For image coordinates The particle, whose offset relative to the image center is . Multiply by the ground sampling distance to get the ground offset, add the world coordinates of the image center to get the initial world coordinates of the mass point.

[0175] The purpose of the initial coarse position is to provide a reasonable starting location for the particles to be matched, placing them in the vicinity of the reference gravitational source set, thus providing a good initial value for subsequent iterative optimization. If the initial position is too deviated, it may lead to slow convergence of the iterative process or getting trapped in local optima.

[0176] It should be noted that in some implementations, more precise initial projection methods can be used, such as using attitude information provided by an inertial measurement unit for attitude compensation, or using a digital elevation model for terrain correction, to further improve the accuracy of the initial position.

[0177] The core optimization process of elastic topology matching involves iteratively calculating the virtual resultant gravitational force acting on the particles to be matched and moving the particles' positions. This process simulates the motion of particles under the influence of gravity in a physical gravitational field. Through multiple iterations, the positions of the particles are gradually adjusted to move them closer to the best-matching reference gravitational source, eventually reaching a stable state. Each iteration includes three main steps: gravity calculation, resultant force superposition, and position update.

[0178] First, a mass parameter is assigned to each baseline gravitational source. This mass parameter reflects the importance and attractive force of the source within the topology. Intersection nodes, as key topological nodes in the road network, have a high mass parameter, typically set to 5 to 10 units. The mass parameter of arterial road intersections is greater than that of secondary arterial road intersections, and the mass parameter of complex intersections with multiple roads is greater than that of simple crossroads. The specific value of the mass parameter can be set according to the connectivity of the intersection; the higher the connectivity, the greater the mass parameter.

[0179] The mass parameters of ordinary road nodes are set relatively small, typically 1 to 2 units of mass. They provide constraints on road orientation but do not dominate the matching process. The mass parameters of building corner nodes are set according to the size of the building; larger buildings have larger mass parameters at their corners, while smaller buildings have smaller mass parameters. By setting differentiated mass parameters, important topological nodes exert a stronger attraction on the mass points to be matched, guiding the mass points to preferentially approach these key nodes.

[0180] Next, the distance vector between each particle to be matched and all reference gravitational sources is calculated. For the i-th particle to be matched, its current position coordinates are... For the j-th reference gravitational source, its position coordinates are The distance vector is the vector pointing from the point mass to the gravitational source, calculated as follows: The magnitude of the distance vector represents the Euclidean distance between the point mass and the gravitational source, calculated as the square root of the sum of the squares of the X and Y components of the distance vector. The direction of the distance vector indicates the direction of gravity, from the point mass towards the gravitational source.

[0181] Based on the distance vector and mass parameters, the magnitude and direction of the virtual gravitational force experienced by each particle to be matched are calculated. The virtual gravity follows laws similar to universal gravitation: the magnitude is directly proportional to the mass of the gravitational source and inversely proportional to the square of the distance. Specifically, the formula is that the magnitude of the gravitational force experienced by the i-th particle from the j-th gravitational source is equal to the gravitational constant multiplied by the mass of the gravitational source divided by the square of the distance.

[0182] The gravitational constant is an adjustment parameter that controls the overall strength of gravity. It needs to be calibrated according to the actual scenario, with a typical value between 0.1 and 1.0. To avoid excessive gravity and numerical instability when a point mass is very close to the gravitational source, a small regularization term, such as 0.01, can be added to the square of the distance. The corrected magnitude of gravity is equal to the gravitational constant multiplied by the mass of the gravitational source, divided by the square of the distance, plus the regularization term. The direction of gravity is the unit direction vector of the distance vector, obtained by dividing the distance vector by its magnitude. The gravitational vector is equal to the magnitude of gravity multiplied by the gravitational direction vector, represented as a two-dimensional vector containing X-axis and Y-axis components.

[0183] The resultant force vector is obtained by vector superimposing all virtual gravitational forces acting on each particle to be matched. The particle is acted upon by gravitational forces from all reference gravitational sources, and these gravitational vectors need to be vector-summed. The X-component of the resultant force vector equals the sum of the X-components of all gravitational vectors, and the Y-component equals the sum of the Y-components of all gravitational vectors. The direction of the resultant force vector indicates the direction in which the particle should move, and the magnitude of the resultant force vector indicates the strength of the driving force for movement.

[0184] In some implementations, to account for the constraints of topological connections, an elastic constraint term can be introduced into the resultant force calculation. For pairs of particles connected by lines, a virtual spring force is applied between them. The magnitude of the spring force is proportional to the change in the distance between the particles relative to the original distance, and the direction is the direction for restoring the original distance. This elastic constraint ensures that the relative positional relationship between the particles does not change drastically during the matching process, maintaining the overall consistency of the topology.

[0185] The position coordinates of the particle to be matched are updated based on the resultant force vector and a preset step size. The step size controls the adjustment range of the particle's position in each iteration. An excessively large step size may cause the particle's position to oscillate or exceed the optimal position, while an excessively small step size will result in slow convergence. The step size can be set to a fixed value, such as 0.5 to 2.0 pixels, or an adaptive strategy can be used to dynamically adjust the step size according to the magnitude of the resultant force. A larger step size is used when the resultant force is large to accelerate convergence, and a smaller step size is used when the resultant force is small to improve accuracy.

[0186] The formula for updating the particle's position is: the new position coordinates equal to the current position coordinates plus the unit direction vector of the resultant force vector multiplied by the step size. The unit direction vector of the resultant force vector is obtained by dividing the resultant force vector by its magnitude, ensuring the correct direction of movement and that the movement distance is controlled by the step size. The updated particle position serves as the input for the next iteration. Through multiple iterations, the particle gradually moves closer to the best-matching reference gravitational source, and the spatial distribution of the entire particle set gradually becomes consistent with the distribution of the reference gravitational source set.

[0187] During the movement of the particle's position, it is necessary to monitor the convergence state of the matching process in real time to determine whether a stable match has been reached, thus determining the termination condition for the iteration. The total potential energy is a key indicator for evaluating the matching quality and convergence state. The total potential energy is defined as the weighted sum of the distances of all particles to be matched to their nearest reference gravitational source.

[0188] To calculate the total potential energy of all unmatched particles in the topology diagram from the UAV's perspective in the current iteration, firstly, for each unmatched particle, calculate its distance to all reference gravitational sources and find the minimum distance, which represents the deviation between the particle and the best-matching gravitational source. Summate the minimum distances of all particles to obtain the total potential energy. In some implementations, the distances can be weighted, with weights set according to the mass parameters of the gravitational sources or the importance of the particles, so that the matching deviation of key nodes contributes more to the total potential energy. The physical meaning of the total potential energy is the potential energy of the entire particle system in the gravitational field; the lower the potential energy, the closer the particle's position is to the gravitational source, and the higher the matching mass.

[0189] The rate of change of total potential energy relative to the previous iteration is calculated. This rate of change is defined as the current iteration's total potential energy minus the previous iteration's total potential energy, divided by the previous iteration's total potential energy, and expressed as a percentage. A negative rate of change indicates a decrease in total potential energy and improved matching quality; a positive rate of change indicates an increase in total potential energy and deteriorated matching quality. This latter situation should generally be avoided, as it may indicate an improper step size setting or local oscillations. The absolute value of the rate of change represents the drastic nature of the potential energy change; a larger absolute value indicates that the matching process is still rapidly adjusting, while a smaller absolute value indicates that the matching process is becoming more stable.

[0190] The process checks if the rate of change of total potential energy is below a preset value, which is a threshold for determining convergence. A typical value is 0.1% to 1%, indicating that the matching process is considered convergent when the change in total potential energy is less than this threshold. If the absolute value of the rate of change of total potential energy is below the preset value, the topology map from the UAV's perspective is considered to have reached a steady state, and the particle positions have converged to the optimal or near-optimal matching positions. At this point, iteration stops, and the final positions of the particles to be matched are output as the matching results. The final positions record the coordinates of each particle in the world coordinate system, establishing a correspondence between the UAV's topology map and the satellite reference image.

[0191] If the rate of change of total potential energy is not lower than a preset value, it indicates that the matching process has not yet converged, and the next iteration needs to be performed, repeating the process of gravitational calculation, resultant force superposition, and position update until the convergence condition is met. To avoid the iteration process from getting stuck in an infinite loop, a maximum number of iterations needs to be set, for example, 100 to 500 iterations. If convergence is still not achieved after reaching the maximum number of iterations, the iteration is forcibly stopped, the current particle position is output as the matching result, and it is marked that the matching may not have fully converged.

[0192] Based on the matching results, pose calculation is performed to convert the final position of the mass point into precise coordinate information of the UAV in the world coordinate system. The matching results establish a mapping relationship between nodes in the topology map from the UAV's perspective and corresponding nodes in the satellite reference map. These mapping relationships constitute a set of two-dimensional to two-dimensional point pairs. Specifically, for each mass point to be matched, its coordinates in the original image are... The final world coordinates in steady state are These two sets of coordinates form a point-to-point correspondence.

[0193] Based on these point-to-point correspondences, the homography matrix is ​​calculated. The homography matrix is ​​a 3×3 matrix describing the projection transformation relationship between two planes, mapping points in the image coordinate system to points in the world coordinate system. The homography matrix is ​​calculated using either a direct linear transformation method or the least squares method. For each pair of corresponding points, two linear equations can be established, and the eight unknown parameters of the homography matrix are obtained by solving the overdetermined system of equations.

[0194] To improve robustness, the RANSAC algorithm can be used to eliminate mismatched point pairs, retaining only inliers for homography matrix calculation. The RANSAC algorithm calculates candidate homography matrices by randomly sampling the minimum point set, then counts the number of inliers that match the matrix, and selects the matrix with the largest number of inliers as the optimal solution.

[0195] The camera's rotation matrix and translation vector can be deduced using the homography matrix. The homography matrix can be decomposed into a combination of the camera intrinsic parameter matrix, the rotation matrix, and the translation vector. A matrix decomposition algorithm can extract the rotation matrix and translation vector from the homography matrix. The rotation matrix describes the rotation relationship between the camera coordinate system and the world coordinate system; it is a 3×3 orthogonal matrix. The translation vector describes the position of the camera's optical center in the world coordinate system; it is a 3-dimensional vector. The decomposition process requires knowledge of the camera's intrinsic parameter matrix, which includes parameters such as focal length and principal point position, and is predetermined through camera calibration.

[0196] Based on the rotation matrix and translation vector, combined with the UAV's flight altitude information, the latitude and longitude coordinates of the UAV in the world coordinate system are calculated. The translation vector provides the horizontal position of the camera's optical center in the world coordinate system, i.e., the X and Y coordinates, which can be converted to latitude and longitude. If the world coordinate system uses a projected coordinate system, such as UTM coordinates, the projected coordinates can be converted to latitude and longitude through inverse projection transformation.

[0197] Flight altitude information is provided by a barometer or laser rangefinder. Combined with the Z component of the translation vector, the complete three-dimensional position of the UAV in the world coordinate system can be obtained. The rotation matrix can be converted into Euler angles or quaternions to represent the UAV's attitude, including pitch, roll, and yaw angles. The output latitude and longitude coordinates and altitude information serve as the UAV's precise coordinate information, which can be used for applications such as navigation control, trajectory planning, and target localization.

[0198] This implementation achieves flexible topology matching between the UAV-view topology map and the satellite baseline map by introducing a virtual gravitational field model. A baseline set of gravitational sources, including intersection nodes, road nodes, and building corner points, is constructed by extracting standard roads and building structures from the satellite baseline map. Nodes in the UAV-view topology map are defined as mass points to be matched, establishing a mapping between the physical model and image features. Initial projection places the mass points in reasonable starting positions, providing good initial values ​​for iterative optimization.

[0199] The iterative calculation process simulates the motion of a point mass in a gravitational field. By setting differentiated mass parameters, key topological nodes play a dominant role in the matching process. The calculation of virtual gravity comprehensively considers distance and mass factors, and the superposition of resultant forces and position updates enable the gradual adjustment of the point mass. By monitoring the rate of change of total potential energy in real time, the convergence state is accurately determined, ensuring that the matching process terminates in a timely manner when it reaches stability.

[0200] Based on the pose calculation of the matching results, the position and attitude information of the UAV were accurately recovered through homography matrix calculation and decomposition. The elastic topology matching process allows for deformation and displacement of local structures, which has stronger environmental adaptability compared with traditional rigid matching, and solves the positioning problem after unstructured deformation of the surface texture.

[0201] Reference Figure 3 This application also provides an image processing apparatus for aerial remote sensing images, comprising: The image preprocessing module is used to acquire aerial images captured in real time by the drone and to perform image preprocessing on the aerial images to obtain preprocessed aerial images. The rigid structure extraction module is used to extract the rigid structure from the preprocessed aerial images to obtain the original incomplete structure diagram after removing debris interference. The architecture repair module is used to repair and complete the broken architecture diagram of the original incomplete architecture diagram to obtain the topology diagram from the perspective of the UAV. The topology matching module is used to perform flexible topology mapping matching between the topology map from the perspective of the UAV and the pre-stored satellite reference map to obtain the matching result; The pose calculation module is used to calculate the pose based on the matching results and output the precise coordinate information of the UAV in the world coordinate system.

[0202] This application also provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the above-described image processing method for aerial remote sensing images.

[0203] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0204] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, as well as combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create a machine for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0205] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0206] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0207] In a typical configuration, a computing device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0208] Memory may include non-persistent memory in computer-readable media, such as random access memory (RAM) and / or non-volatile memory, such as read-only memory (ROM) or flash RAM. Memory is an example of computer-readable media.

[0209] Computer-readable media includes both permanent and non-permanent, removable and non-removable media that can store information using any method or technology. Information can be computer-readable instructions, data structures, modules of programs, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transferable medium that can be used to store information accessible by a computing device. As defined herein, computer-readable media does not include transient computer-readable media, such as modulated data signals and carrier waves.

[0210] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0211] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. An image processing method for aerial remote sensing images, characterized in that, include: The system acquires aerial images captured in real time by drones and performs image preprocessing on the aerial images to obtain preprocessed aerial images. Rigid structure extraction is performed on the preprocessed aerial images to obtain the original incomplete structure diagram after removing debris interference; The original incomplete architecture diagram is repaired and completed to obtain the topology diagram from the perspective of the UAV. The topology map from the UAV's perspective is matched with a pre-stored satellite reference map using flexible topology mapping to obtain the matching result; Based on the matching results, pose calculation is performed, and the precise coordinate information of the UAV in the world coordinate system is output.

2. The method according to claim 1, characterized in that, Rigid structure extraction is performed on the preprocessed aerial imagery to obtain the original incomplete structure diagram after removing debris interference, including: Obtain the pixel grayscale distribution of the preprocessed aerial imagery, and construct a three-dimensional terrain matrix based on the pixel grayscale distribution; Determine the virtual fluid parameters and initialize the full-map uniform flow field of the 3D terrain matrix; Based on virtual fluid parameters, virtual fluid is iteratively flowed on a three-dimensional terrain matrix, where the virtual fluid parameters include the velocity decay coefficient and the viscosity coefficient. During the iterative flow process, regions where the flow velocity diverges or decays rapidly are identified in the three-dimensional terrain matrix, as well as streamline regions where the flow velocity converges and remains stable. Set the pixel gradient value to zero in areas where the flow velocity is diverging or rapidly decaying, and enhance the pixel value to a preset value in streamline areas where the flow velocity converges and remains stable, in order to remove gravel interference. After a preset number of iterations, the streamline convergence center is determined, and a binarized original incomplete architecture diagram is generated based on the streamline convergence center.

3. The method according to claim 2, characterized in that, Determine the virtual fluid parameters and initialize the full-map uniform flow field of the 3D terrain matrix, including: Obtain the resolution and shooting altitude of aerial imagery; The initial value of the flow velocity attenuation coefficient is determined based on the resolution and shooting altitude of the aerial imagery; The initial value of the viscosity coefficient is determined based on the pre-stored characteristics of the gravel density in the post-disaster environment; An initial velocity vector is set at each pixel location of the 3D terrain matrix, and the direction of the initial velocity vector points in the gradient descent direction. The initial velocity vector is normalized to initialize a uniform flow field across the entire 3D terrain matrix.

4. The method according to claim 2, characterized in that, Set the pixel gradient values ​​to zero in regions where the flow rate diverges or decays rapidly, including: Calculate the velocity divergence and velocity decay rate at each pixel location in the 3D terrain matrix; Determine whether the velocity divergence is greater than a preset velocity threshold, or whether the velocity decay rate is greater than a preset gradient threshold. If the velocity divergence is greater than the velocity threshold or the velocity decay rate is greater than the gradient threshold, the gradient value at the corresponding pixel position is set to zero to simulate the effect of sediment being washed away.

5. The method according to claim 1, characterized in that, The original incomplete architecture diagram is repaired and completed to obtain a topology diagram from the perspective of the UAV, including: Identify all line breakpoints in the original incomplete architecture diagram and calculate the tangent direction vector for each line breakpoint; A virtual stress field is constructed between each pair of adjacent line breaks, and a probability score for connection between each pair of adjacent line breaks is calculated based on the virtual stress field. A growth threshold is set based on the probability score, and connecting lines are generated for line breakpoints with probability scores higher than the growth threshold. The generated connection lines are merged with the original incomplete architecture diagram to output the repaired topology diagram from the drone's perspective.

6. The method according to claim 5, characterized in that, A virtual stress field is constructed between each pair of adjacent line breaks, and a probability score for connection between each pair of adjacent line breaks is calculated based on the virtual stress field, including: Calculate the Euclidean distance between each pair of adjacent breakpoints; Calculate the cosine of the angle between the tangent direction vectors of each pair of adjacent breakpoints; Extract the original texture residue of pixels between each pair of adjacent breakpoints, where the Euclidean distance, the cosine of the included angle, and the original texture residue constitute a virtual stress field; Based on the Euclidean distance, the cosine of the included angle, and the original texture residual, a weighted summation is used to calculate the probability score of connection between each pair of adjacent line breakpoints.

7. The method according to claim 5, characterized in that, For line breakpoints with a probability score higher than the growth threshold, generate connecting lines, including: For line breakpoint pairs with probability scores higher than the growth threshold, determine the starting and ending coordinates of the connecting line. Calculate the curvature parameter of the connecting line based on the tangent direction vectors at the two line break points; Connecting lines are generated based on curvature parameters.

8. The method according to claim 1, characterized in that, The topology map from the UAV's perspective is matched with a pre-stored satellite reference map using flexible topology mapping to obtain the matching results, including: Based on satellite baseline images, standard roads and standard building structures are extracted, and a set of baseline gravity sources is constructed based on the standard roads and standard building structures; Define the nodes in the topology graph from the drone's perspective as the set of mass points to be matched; The set of particles to be matched is projected onto the coordinate system of the reference gravitational source set to obtain the initial rough position; Iteratively calculate the virtual gravitational resultant force of all reference gravitational sources on each particle in the set of particles to be matched, and move the position of the particle according to the virtual gravitational resultant force; During the movement of the particle's position, the total potential energy change rate of the topology map from the UAV's perspective is calculated in real time. When the total potential energy change rate is lower than the preset value, the iteration stops and the final position of the particle to be matched is output as the matching result.

9. An image processing apparatus for aerial remote sensing images, applied to the image processing method for aerial remote sensing images according to any one of claims 1-8, characterized in that, include: The image preprocessing module is used to acquire aerial images captured in real time by the drone and to perform image preprocessing on the aerial images to obtain preprocessed aerial images. The rigid structure extraction module is used to extract the rigid structure from the preprocessed aerial images to obtain the original incomplete structure diagram after removing debris interference. The architecture repair module is used to repair and complete the broken architecture diagram of the original incomplete architecture diagram to obtain the topology diagram from the perspective of the UAV. The topology matching module is used to perform flexible topology mapping matching between the topology map from the perspective of the UAV and the pre-stored satellite reference map to obtain the matching result; The pose calculation module is used to calculate the pose based on the matching results and output the precise coordinate information of the UAV in the world coordinate system.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the image processing method for aerial remote sensing images as described in any one of claims 1 to 8.

Images