A multi-source surveying and mapping data intelligent fusion processing method and system

Abstract

The application discloses a kind of multi-source surveying and mapping data intelligent fusion processing method and system, it is related to digital image processing technical field;The method includes the following steps: obtaining the multi-temporal orthophoto of flexible membrane surface to be processed area, multi-temporal photogrammetry elevation and multi-temporal laser elevation, constructs joint elevation field and extracts the candidate valley line set carrying cross direction information, further determines the geometric low center and mirror flash migration amount, constructs valley line locking weight, extracts global lock valley center line by particle swarm algorithm, and according to this, multi-temporal orthophoto, photogrammetry elevation and laser elevation are fused, and fusion image and fusion elevation are obtained;The application can improve the stability of real valley line positioning of flexible membrane surface, weaken mirror flash migration interference, improve the time sequence consistency and geometric reliability of fusion result.

Description

Technical Field

[0001] This invention relates to the field of digital image processing technology, and in particular to a method and system for intelligent fusion processing of multi-source surveying and mapping data. Background Technology

[0002] During the inspection of wastewater treatment floating cover tanks, membrane-covered storage tanks, and other flexible membrane facilities, it is usually necessary to use multi-source mapping data such as UAV orthophotos, photogrammetric elevation, and laser elevation to identify and analyze wrinkles, local depressions, boundary shifts, and abnormal areas on the membrane surface. Because flexible membrane surfaces continuously generate morphological undulations and local valley structures under the influence of wind loads, hydrodynamic forces, surface attachment loads, and changes in internal support conditions, and because shallow water accumulation on the membrane surface can easily form specular flare bands under different lighting conditions, the same location exhibits both actual geometric changes and significant optical brightness drift at different time phases. Therefore, how to stably extract the true valley line positions from multi-source mapping data under multi-temporal conditions, reduce the interference of specular flare migration on structural interpretation, and obtain fused images and fused elevations that accurately reflect the long-term stable geometric state of the membrane surface has become a key technical problem in the intelligent inspection of flexible membrane facilities.

[0003] Existing technologies, when addressing the aforementioned issues, typically focus on brightness analysis of single image data or surface reconstruction of single elevation data, lacking a dedicated mechanism to handle the differences between multi-temporal optical drift and the true valley line position. When specular flare bands, localized high-brightness anomalies, or short-term deformations exist on the membrane surface, optically anomalous areas are easily misidentified as true valley line positions, or instantaneous elevation disturbances are directly incorporated into the fusion result, leading to unstable valley line positioning, blurred structural boundaries, and distorted fused elevation. Furthermore, existing technologies do not adequately utilize the combined stability of geometric centers and brightness drift fluctuations across multiple temporal phases, making it difficult to establish a stable and consistent valley line reference across the entire area. Consequently, in complex membrane surface scenarios, there are still shortcomings such as insufficient reliability of fusion results, weak anti-interference capabilities, and poor temporal consistency. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies, such as unstable extraction of the true valley line position within the processing area of ​​the flexible membrane surface, difficulty in removing mirror flash migration interference, and insufficient reliability of multi-source mapping data fusion results. Therefore, this invention proposes an intelligent fusion processing method and system for multi-source mapping data.

[0005] To address the problems existing in the prior art, the present invention adopts the following technical solution: A method for intelligent fusion processing of multi-source surveying and mapping data includes: S1. Acquire multi-temporal orthophotos, multi-temporal photogrammetric elevations, and multi-temporal laser elevations of the area to be processed on the flexible membrane surface; S2. Construct a joint elevation field based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and extract a set of candidate valley lines carrying cross-sectional direction information based on the joint elevation field, wherein the cross-sectional direction is the direction perpendicular to the local valley line extension direction; S3. Along the cross-sectional direction corresponding to the candidate valley line set, determine the geometric valley center and mirror flash migration on the joint elevation field and the multi-temporal orthophoto image, respectively. S4. Based on the fluctuation of the geometric valley center and the fluctuation of the mirror flash migration under multiple time phases, construct valley line locking weights; S5. Extract the global valley centerline based on the valley line locking weight using the particle swarm optimization algorithm; S6. Using the global valley centerline as a reference, and combining the valley line locking weight, perform data fusion on the multi-temporal orthophoto, the multi-temporal photogrammetric elevation, and the multi-temporal laser elevation to obtain the fused image and fused elevation.

[0006] Preferably, a joint elevation field is constructed based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and a candidate valley line set carrying cross-sectional direction information is extracted based on the joint elevation field, including: Calculate the absolute values ​​of the Laplace operator for the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, respectively; The first weighting coefficient and the second weighting coefficient are determined based on the reciprocal of the absolute value of the Laplacian operator of the multi-temporal photogrammetric elevation and the reciprocal of the absolute value of the Laplacian operator of the multi-temporal laser elevation, respectively. The first weighting coefficient is multiplied by the multi-temporal photogrammetric elevation, the second weighting coefficient is multiplied by the multi-temporal laser elevation, and the weighted sum is performed according to the normalized results of the first weighting coefficient and the second weighting coefficient to construct a joint elevation field. Determine the local valley line extension direction of the joint elevation field, and define the direction perpendicular to the local valley line extension direction as the transverse direction; Calculate the second derivative of the joint elevation field along the transverse direction, and determine the valley line intensity based on the negative value of the second derivative; Based on the continuous spatial extension of the valley line intensity, the location of the region carrying the cross-sectional direction information is extracted to form a candidate valley line set.

[0007] Preferably, determining the geometric valley center on the combined elevation field includes: For each position in the candidate valley line set, a one-dimensional profile is established along the corresponding transverse direction; In the joint elevation field corresponding to the one-dimensional profile, the location corresponding to the minimum elevation is searched and used as the geometric valley center.

[0008] Preferably, measuring mirror flicker migration on the multi-temporal orthophoto image includes: Based on the cross-sectional direction, a brightness profile is extracted from the multi-temporal orthophoto image, and the local mean of the brightness profile is calculated. The portion of the brightness profile that is higher than the local mean is identified as a significant portion; Using the geometric valley center as a reference, the mirror flash saliency center is calculated by performing a weighted integral on the cross-sectional distance based on the salient portion; The positional difference between the salient center of the mirror flash and the center of the geometric trough is calculated as the mirror flash migration amount.

[0009] Preferably, based on the fluctuation of the geometric trough center and the fluctuation of the mirror flash migration over multiple time phases, a trough locking weight is constructed, including: At the same spatial location, the first variance of the geometric trough centers in multiple time phases is statistically analyzed and used as the fluctuation amount of the geometric trough centers; At the same spatial location, the second variance of the mirror flash migration amount across multiple time phases is statistically analyzed and used as the fluctuation amount of the mirror flash migration amount; The first variance is transformed by its reciprocal to obtain the first reciprocal type quantity; Divide the first reciprocal type by the sum of the first reciprocal type and the second variance to obtain the valley line locking weight.

[0010] Preferably, the global valley centerline is extracted using a particle swarm optimization algorithm based on the valley centerline locking weight, including: The average valley line intensity is calculated based on the valley line intensity under multiple defined time phases; Define a continuous curve based on a set of control points; Calculate the integral of the product of the valley locking weight and the average valley intensity at each point on the continuous curve, and combine the integral with the smoothness penalty term of the continuous curve to construct the objective function; The control point set is iteratively optimized using the particle swarm optimization algorithm. The continuous curve corresponding to the control point set that makes the objective function optimal is used as the global valley centerline.

[0011] Preferably, the control point set is iteratively optimized using a particle swarm optimization algorithm, including: In each iteration of the particle swarm optimization algorithm, for each particle, the objective function value at the current position, the deviation between the objective function value at the current position and the objective function value at the individual's historical best position, and the deviation between the objective function value at the current position and the objective function value at the group's historical best position are calculated. Based on the normalized exponential function, the deviation between the objective function value at the current position and the objective function value at the individual's historical best position, and the deviation between the objective function value at the current position and the objective function value at the group's historical best position are processed to adaptively generate the inertia term weight, individual cognitive term weight, and group social term weight in the velocity update equation. The velocity and position of each particle are updated based on the velocity update equation, which includes the inertial term weight, the individual cognitive term weight, and the group social term weight.

[0012] Preferably, data fusion of the multi-temporal orthophotos includes: Obtain the reference transverse direction corresponding to the global valley centerline at the current spatial location; For each time phase and each spatial location, the orthophoto value of the current time phase is multiplied by the valley line locking weight to obtain the first component; In the reference cross-sectional direction, with the current position as the center, two symmetrical positions are determined according to the specular flash migration amount, and the mean value of the image values ​​corresponding to the two symmetrical positions is calculated. The difference obtained by subtracting the valley line locking weight is multiplied by the mean value to obtain the second component. The first component and the second component are added together to generate a corrected image for the corresponding time phase; The average value of the corrected images across all time phases is used to obtain the fused image.

[0013] Preferably, the data fusion of the multi-temporal photogrammetric elevation and the multi-temporal laser elevation includes: For each temporal phase and each spatial location, the valley line locking weight is used as the first multiplier and multiplied by the multi-temporal photogrammetric elevation to obtain the first product; The difference obtained by subtracting the valley line locking weight is used as the second multiplier and multiplied by the multi-temporal laser elevation to obtain the second product; Add the first product to the second product to obtain the weighted elevation for the corresponding time phase; The weighted elevations for all time phases are averaged to obtain the merged elevation.

[0014] To address the aforementioned problems, the present invention also provides an intelligent fusion processing system for multi-source surveying and mapping data, comprising: The data acquisition module is used to acquire multi-temporal orthophotos, multi-temporal photogrammetric elevations, and multi-temporal laser elevations of the area to be processed on the flexible membrane surface. The valley line extraction module is used to construct a joint elevation field based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and to extract a set of candidate valley lines carrying cross-sectional direction information based on the joint elevation field, wherein the cross-sectional direction is the direction perpendicular to the local valley line extension direction. The offset measurement module is used to measure the geometric valley center and mirror flash migration along the cross direction corresponding to the candidate valley line set, respectively, on the joint elevation field and the multi-temporal orthophoto. The weight construction module is used to construct valley line locking weights based on the fluctuation of the geometric valley center and the fluctuation of the mirror flash migration in multiple time phases. The global solution module is used to extract the global valley centerline based on the valley line locking weight using the particle swarm optimization algorithm. The multi-source fusion module is used to fuse the multi-temporal orthophoto, multi-temporal photogrammetric elevation, and multi-temporal laser elevation data based on the global valley centerline and combined with the valley line locking weight, to obtain fused images and fused elevations.

[0015] Compared with the prior art, the beneficial effects of the present invention are: 1. This invention introduces multi-temporal orthophotos, multi-temporal photogrammetric elevation, and multi-temporal laser elevation to first construct a joint elevation field. Then, it extracts a set of candidate valley lines carrying cross-sectional information and measures the geometric valley center and mirror flash migration in the corresponding cross-sectional directions, thereby distinguishing the true concave structure of the flexible membrane surface from the brightness anomaly formed by shallow water reflection. Furthermore, by jointly quantifying the fluctuations of the geometric valley center and mirror flash migration across multiple temporal phases, a valley line locking weight is constructed, and a global valley center line is extracted accordingly. This ensures that the spatial position of the true valley line on the membrane surface can be stably constrained under multi-temporal conditions. This improves the accuracy and continuity of valley line positioning and reduces the interference of mirror flash drift on structural interpretation.

[0016] 2. This invention uses the global valley centerline as a unified structural benchmark, performs mirror flash stripping fusion on multi-temporal orthophotos, and performs adaptive weighted fusion on multi-temporal photogrammetric elevations and multi-temporal laser elevations. This allows high-confidence areas to retain the true texture and detail information of the membrane surface, and allows areas that are strongly affected by mirror flash or have large geometric fluctuations to be compensated by stable elevation information. This not only produces fused images with clear boundaries and stable temporal sequence, but also fused elevations with strong consistency and reliable geometric accuracy, thereby improving the ability to identify wrinkle extension, local depressions and abnormal areas during the inspection of flexible membrane facilities. Attached Figure Description

[0017] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this application, illustrate exemplary embodiments of the invention and, together with their description, serve to explain the invention and do not constitute an undue limitation thereof. In the drawings: Figure 1 This is a flowchart illustrating an intelligent fusion processing method for multi-source surveying and mapping data according to an embodiment of the present invention. Figure 2 This is a functional block diagram of a multi-source mapping data intelligent fusion processing system provided in an embodiment of the present invention. Detailed Implementation

[0018] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0019] This embodiment provides a method for intelligent fusion processing of multi-source surveying and mapping data. (See also...) Figure 1 Specifically, it includes the following steps: S1. Acquire multi-temporal orthophotos, multi-temporal photogrammetric elevations, and multi-temporal laser elevations of the area to be processed on the flexible membrane surface; S2. Construct a joint elevation field based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and extract a set of candidate valley lines carrying cross-sectional direction information based on the joint elevation field, wherein the cross-sectional direction is the direction perpendicular to the local valley line extension direction; S3. Along the cross-sectional direction corresponding to the candidate valley line set, determine the geometric valley center and mirror flash migration on the joint elevation field and the multi-temporal orthophoto image, respectively. S4. Based on the fluctuation of the geometric valley center and the fluctuation of the mirror flash migration under multiple time phases, construct valley line locking weights; S5. Extract the global valley centerline based on the valley line locking weight using the particle swarm optimization algorithm; S6. Using the global valley centerline as a reference, and combining the valley line locking weight, perform data fusion on the multi-temporal orthophoto, the multi-temporal photogrammetric elevation, and the multi-temporal laser elevation to obtain the fused image and fused elevation.

[0020] In one embodiment of the present invention, acquiring multi-temporal orthophotos, multi-temporal photogrammetric elevations, and multi-temporal laser elevations of the area to be processed on the flexible membrane surface includes: The flexible membrane facility to be treated is delineated into a treatment area. Multiple ground control points are evenly distributed around the treatment area, covering the perimeter and diagonal positions of the treatment area. The control points are set with forced centering markers. The plane and elevation coordinates of the control points are obtained through static observation of the global satellite navigation system to ensure that the spatial reference of all collected data is consistent. To address the inspection needs of flexible membrane structures, aerial survey data is collected in multiple consecutive time phases according to a pre-set time phase plan. This time phase plan is determined based on the membrane structure's operating cycle, historical inspection frequency, and local weather stability window. Preferably, the time interval between two adjacent collections is no less than one calendar day and no more than one month to ensure that the phased changes in the membrane structure are reflected while avoiding incomparable structural states due to excessively large time intervals. Each time phase is collected under meteorological conditions of low wind speed, no rainfall, and no newly added large-area flowing water accumulation on the membrane surface. Low wind speed refers to an average wind speed within ten minutes before and during the collection period that does not exceed 50% of the maximum allowable wind speed for similar UAV aerial survey operations in the local area. The collection time is uniformly limited to a period of time where the solar altitude angle changes gradually and there is no significant new projection obstruction in the area to be processed. Preferably, the start time of the first time phase is used as the reference time, and the deviation of subsequent time phases from this reference time does not exceed two hours to reduce specular flash background changes caused by differences in illumination geometry. The same set of standardized flight parameters is used for each time phase of the aerial survey flight. These parameters are determined based on the trial data collection results of the first time phase and remain consistent in subsequent time phases. Specifically, the flight altitude relative to the membrane surface is determined jointly based on the camera focal length, sensor size, and target ground sampling distance, ensuring that the ground sampling distance of the orthophoto is no greater than half the minimum wrinkle width to be identified. The forward overlap and lateral overlap are set to no less than 70% and no less than 60%, respectively, according to the requirements of multi-view stereo modeling. If subsequent time phases cannot fully maintain the flight altitude of the first time phase due to temporary airspace restrictions or changes in equipment status, the ground sampling distance change rate is kept to no more than 10%, and the data is uniformly resampled to the spatial resolution of the first time phase after generating the orthophoto and elevation data. This ensures that the basic data of each time phase remains comparable at the acquisition level, facilitating the unified calculation of mirror flash migration and valley line locking weights in subsequent phases. During each time-phase acquisition process, a sequence of aerial photographs of the area to be processed is acquired using an area-array aerial survey camera mounted on an UAV platform. After the aerial photographs are acquired, distortion correction and aerial triangulation adjustment are performed on the sequence of aerial photographs. Absolute orientation is then achieved by combining the aerial photographs with ground control points, generating a digital orthophoto of the area to be processed as the orthophoto of the corresponding time phase. The orthophotos of all time phases have a unified resolution and spatial coordinate range, forming a multi-time-phase orthophoto. For the sequence of aerial photographs after aerial triangulation adjustment for each time phase, a digital surface model of the area to be processed is generated using a multi-view stereo matching algorithm. Outlier removal and filtering are performed on the digital surface model to remove flight noise and elevation interference from non-membrane ground features, obtaining the photogrammetric elevation of the corresponding time phase. The photogrammetric elevations of all time phases maintain strict spatial registration with the orthophoto of the same time phase, forming a multi-time-phase photogrammetric elevation. Each phase of the aerial survey flight simultaneously acquires laser point cloud data of the area to be processed through an airborne lidar system mounted on the same UAV platform. The point cloud density of the lidar is determined according to the accuracy requirements of membrane surface elevation detection. After acquisition, the laser point cloud data undergoes trajectory calculation, point cloud positioning and adjustment processing. Absolute orientation is completed in conjunction with ground control points. The oriented point cloud is then subjected to noise reduction filtering and ground point classification processing to extract the effective point cloud corresponding to the membrane surface. A digital elevation model with the same resolution and coordinate range as the orthophoto is generated through spatial interpolation to obtain the laser elevation of the corresponding phase. The laser elevations of all phases maintain strict spatial registration with the orthophoto and photogrammetric elevations of the same phase, forming multi-phase laser elevations.

[0021] It should be noted that flexible membrane surfaces refer to continuous covering structures formed by flexible polymer membrane materials or composite membrane materials. These structures can undergo elastic deformation or slow displacement under their own weight, surface loads, internal support conditions, wind loads, and the support of liquids or gases, exhibiting geometric changes such as local undulations, depressions, wrinkle extension, and boundary displacement. Unlike rigid panels, the surface morphology of flexible membrane surfaces dynamically adjusts with changes in stress and operating conditions, thus easily presenting both actual geometric changes and optical reflection anomalies in the surveying and acquisition results. It should be noted that multi-temporal orthophotos refer to a collection of images acquired at different times for the same area to be processed, which are then processed through distortion correction, aerial triangulation adjustment, georegistration, and orthorectification to form an image set with a unified spatial reference, unified projection relationship, and unified resolution. The temporal images in this image set can be mapped pixel-by-pixel or raster-by-raster in the same planar coordinate system, and are used to characterize the surface texture, brightness distribution, and reflection changes of the area to be processed at different times. It should be noted that multi-temporal photogrammetric elevation refers to a set of elevation data obtained for the same area to be processed at different acquisition times, based on a sequence of aerial photographs or multi-view image data, through aerial triangulation, multi-view stereo matching, surface reconstruction, and elevation interpolation. Each temporal result in this elevation set maintains a spatial registration relationship with the orthophoto of the corresponding temporal phase and reflects the geometric undulations, local concavity, and wrinkle morphology changes of the surface of the area to be processed at different times. Multi-temporal laser elevation refers to a set of elevation data obtained for the same area to be processed at different acquisition times using laser point cloud data acquired by airborne or other carrier-mounted laser ranging equipment, and processed through positioning calculation, point cloud registration, filtering and denoising, target surface extraction, and elevation rasterization. Each temporal result in this elevation set is within a unified spatial coordinate framework with the orthophoto and photogrammetric elevation of the corresponding temporal phase. It can be used to characterize the true elevation distribution and temporal changes of the surface of the area to be processed and provides independent elevation data to suppress the interference of optical imaging anomalies on geometric measurement results.

[0022] In one embodiment of the present invention, a joint elevation field is constructed based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and a candidate valley line set carrying cross-sectional direction information is extracted based on the joint elevation field, including: For each time phase, the multi-temporal photogrammetric elevation and multi-temporal laser elevation are processed in parallel across the entire region by traversing the grid one by one. First, the photogrammetric elevation and laser elevation of the current time phase are preprocessed separately. The overall trend surface of the elevation of the entire region is fitted by a quadratic polynomial. The original elevation data is subtracted from the fitted trend surface value to complete the detrending process and eliminate the interference of the overall slope of the membrane surface on the subsequent local curvature calculation. The absolute values ​​of the Laplacian operator were calculated for the preprocessed photogrammetric and laser elevation data. For the raster-format elevation data, a 3×3 window four-neighbor discrete Laplacian convolution kernel was used to perform the convolution operation. The center element of the four-neighbor discrete Laplacian convolution kernel was -4, the four adjacent elements (top, bottom, left, and right) were 1, and the other four diagonal elements were 0. After the convolution operation, the absolute value of the convolution result for each raster position was taken to obtain the absolute Laplacian values ​​for the photogrammetric and laser elevation data, respectively. The core purpose of using the Laplacian operator for calculation here is to quantify the degree of local curvature change in the elevation data through second-order differential operations. The smaller the absolute value of the Laplacian operator, the smoother the elevation surface, which is more in line with the characteristics of the real concave structure of the flexible membrane surface. The larger the absolute value of the Laplacian operator, the more obvious the abrupt change in the elevation surface, which is more likely to be noise or non-structural abrupt change. Based on this, the weights can be constructed to prioritize the contribution of the real geometric structure. The first weighting coefficient and the second weighting coefficient are determined based on the reciprocals of the absolute values ​​of the Laplacian operator for photogrammetric elevation and laser elevation, respectively. The formula for calculating the first weighting coefficient is as follows: ; The formula for calculating the second weighting coefficient is: ; In the formula, The first weighting coefficient for the current grid position. The second weighting coefficient for the current grid position. The absolute value of the Laplacian operator for the photogrammetric elevation of the current grid location. This represents the absolute value of the Laplacian operator for the laser elevation at the current grid position. The value is a preset minimum positive number, which is one-thousandth of the elevation mean error of the entire region in the current time phase elevation data. If the elevation mean error is not predetermined, it is one-hundred-thousandth of the range of the entire region in the current time phase elevation data, in order to avoid the calculation anomaly of zero denominator. The purpose of using the reciprocal form to construct the weights here is to enhance the weights of the elevation smooth region and suppress the weights of the elevation abrupt region, so that the subsequent construction of the joint elevation field is more in line with the real membrane geometry. After completing the grid-by-grid calculation of the first and second weighting coefficients, the first weighting coefficient is multiplied by the photogrammetric elevation at the same grid location, and the second weighting coefficient is multiplied by the laser elevation at the same grid location. Then, the results are weighted and summed according to the normalized results of the first and second weighting coefficients to construct the joint elevation field for the current time phase. The formula for calculating the joint elevation field is as follows: ; In the formula, This is the joint elevation field value for the current grid location. The photogrammetric elevation of the current grid location after preprocessing. The laser elevation is the preprocessed elevation of the current grid position. The meanings of the other parameters are the same as those in the aforementioned formula. Here, a joint elevation field is constructed by normalization and weighted summation, which can realize the adaptive fusion of photogrammetric elevation and laser elevation. It retains the effective information of the two data sources in the real structure area with smooth elevation, and suppresses the interference of abnormal data in the noisy area with abrupt elevation changes, thereby improving the stability and reliability of elevation data.

[0023] It should be noted that the joint elevation field refers to the elevation distribution result formed by adaptively weighting the photogrammetric elevation and laser elevation acquired in the same area to be processed at the same time phase according to the local curvature changes at their respective corresponding locations. This elevation field is not a simple superposition of the two types of elevation data, but rather adjusts the contribution ratio of the two types of elevation data according to the elevation smoothness at different locations. This ensures that locations with smooth elevations that better match the actual concave shape of the flexible membrane surface are effectively preserved, while locations with significant elevation abrupt changes and more likely to be affected by noise are given a lower participation. This results in a unified elevation expression that takes into account both surface geometric continuity and local structural stability.

[0024] The local valley line extension direction of the joint elevation field is determined grid-by-grid, and the direction perpendicular to the local valley line extension direction is determined as the transverse direction. Specifically, for each grid position in the joint elevation field, the Hessian matrix at that position is calculated using a 3×3 window and the central difference method. The second-order partial derivative in the x-direction is calculated using the central difference of adjacent grids in the x-direction, the second-order partial derivative in the y-direction is calculated using the central difference of adjacent grids in the y-direction, and the mixed partial derivative in the xy-direction is calculated using the central difference of diagonally adjacent grids. The elements of the Hessian matrix are obtained from the above... The second-order partial derivatives and mixed partial derivatives are used to perform eigenvalue decomposition on the Hessian matrix, resulting in two eigenvalues ​​and two corresponding orthogonal eigenvectors. The direction of the eigenvector corresponding to the eigenvalue with the smaller absolute value is the local valley extension direction. This direction is the direction of minimum local curvature of the elevation surface and also the natural extension direction of the valley line of the membrane surface folds. The direction perpendicular to this local valley extension direction is the transverse direction. This direction is the direction of maximum local curvature of the elevation surface and also the normal section direction of the membrane surface concave structure. The corresponding transverse direction information is stored synchronously at each grid position.

[0025] It should be noted that the transverse direction refers to the direction perpendicular to the valley line extension direction at a local location in the joint elevation field. This direction is used to characterize the change direction of local wrinkles or depressions on the flexible membrane surface along the valley line extension path. It is also the reference direction for establishing one-dimensional profiles, locating the geometric valley center, and measuring mirror flash migration. When observing along the transverse direction, the degree of local depression of the membrane surface, the position of the valley line boundary, and the offset of optical anomalies relative to the true geometric center can be reflected more directly. Therefore, the transverse direction is a key reference direction connecting valley line geometric analysis and mirror flash interference analysis.

[0026] After completing the full-area calculation in the transverse direction, the second derivative of the joint elevation field along the transverse direction is calculated grid by grid, and the valley line intensity is determined based on the negative value of the second derivative. The formula for calculating the valley line intensity is as follows: ; In the formula, The valley line intensity at the current grid position. To combine the second derivative of the elevation field along the cross-sectional direction corresponding to the current grid, the negative value of the second derivative is used here and non-negative truncation is taken to calculate the valley line intensity. The core reason is that the second derivative of the concave structure corresponding to the valley line of the membrane wrinkle along the cross-sectional direction is negative. The larger the absolute value of the negative value, the more significant the concavity, and the more obvious the valley line characteristics. By using non-negative truncation, the positive value results corresponding to the convex structure can be removed, and only the contribution of the concave structure that conforms to the valley line characteristics can be retained. Based on the continuous spatial extension of valley line intensity, the location of regions carrying transverse direction information is extracted to form a candidate valley line set. The specific implementation process is as follows: the Otsu method is used to adaptively calculate the valley line intensity of the entire region to obtain the valley line intensity threshold. The valley line intensity map is then binarized and segmented to obtain the valley line candidate regions. Morphological refinement processing is then performed on the valley line candidate regions to extract a linear skeleton with a single pixel width. Subsequently, the linear skeleton is continuously tracked to determine the minimum length threshold. The minimum length threshold is one-thousandth of the diagonal length of the entire range of the region to be processed. Short line segments with a length less than the minimum length threshold are removed, and spatially continuous valley line segments are retained. Finally, each grid position on all retained valley line segments, together with the corresponding transverse direction information, is stored to form a candidate valley line set carrying transverse direction information.

[0027] Furthermore, to avoid instability in the candidate valley line extraction results caused by differences in the scale of different regions to be processed, after obtaining the valley line intensity threshold using the Otsu method, the length, average valley line intensity, and line segment direction consistency of each linear skeleton in the binarized result are calculated simultaneously. Only when the line segment length is not less than the minimum length threshold, the average valley line intensity on the line segment is not less than the average valley line intensity of the entire region, and the mean square value of the direction change of adjacent skeleton points is not greater than the preset upper limit of direction fluctuation, is the line segment retained as a valid candidate valley line segment. The upper limit of direction fluctuation is determined by the statistical median of the mean square values ​​of the direction changes of all candidate line segments in the current time phase. By simultaneously constraining length, intensity, and direction consistency, it is possible to avoid relying solely on the length threshold to retain pseudo valley line segments formed by noise, local reflection disturbances, or boundary anomalies, thereby improving the stability of the candidate valley line set in representing real folded valley lines. Furthermore, when a grid location is not extracted as a candidate valley line in two adjacent time phases, but is extracted as a single candidate valley line point in the current time phase, the grid location is marked as a temporary outlier and is not directly included in the candidate valley line set; only when the location meets the candidate valley line conditions in at least two consecutive time phases is it retained as a valid candidate valley line location; this processing method is used to reduce the impact of single-time phase random noise on the subsequent determination of geometric valley center.

[0028] It should be noted that the candidate valley line set refers to the set of linear positions that may correspond to the actual valley lines of the flexible membrane surface, extracted from the joint elevation field based on the strength of the concave response along the transverse direction and its spatial continuity. Each position in this set includes not only spatial coordinate information but also transverse direction information corresponding to that position, which is used to support the subsequent determination of the geometric valley center, calculation of mirror flash migration, and solution of the global valley center line. The candidate valley line set is essentially a preliminary screening result of the actual valley line structure in the area to be processed. Its role is to preferentially delineate areas with stable concave characteristics and continuous extension characteristics from the elevation data of the entire region, providing a reliable structural starting point for subsequent refined fusion processing.

[0029] In one embodiment of the present invention, determining the geometric valley center and mirror flash migration along the cross-sectional direction corresponding to the candidate valley line set on the joint elevation field and the multi-temporal orthophoto image respectively includes: For each location in the candidate valley line set, a one-dimensional local coordinate system for the profile is established, with that location as the origin and the corresponding transverse direction as the positive direction of the local coordinate axis. The half-width of the profile is adaptively determined based on the average valley width of the local region where the current candidate valley line is located. Specifically, near the current candidate valley line location, the joint elevation field is pre-scanned along the transverse direction to determine the two boundary locations where the second derivative of the elevation changes from negative to zero or from negative to positive. Half the distance between the two boundaries is taken as the estimated local valley width. If this estimated value is less than three times the ground sampling distance, the half-width of the profile is set to three times the ground sampling distance; if this estimated value is greater than ten times the ground sampling distance, the half-width of the profile is set to ten times the ground sampling distance; otherwise, this estimated value is directly used as the half-width of the profile. The sampling step size is one ground sampling distance, and the number of sampling points is automatically determined based on the half-width of the profile and the sampling step size, and adjusted to the nearest odd number to ensure that the origin location is the center location of the sampling sequence. Sampling is performed along both the positive and negative directions of the transverse direction on the joint elevation field and the simultaneous multi-temporal orthophoto image at equally spaced sampling steps. During the sampling process, bilinear interpolation is used to obtain values ​​at non-raster center locations to ensure the spatial continuity and registration consistency of the sampled data, generating one-dimensional elevation profiles for the corresponding joint elevation field and one-dimensional brightness profiles for the corresponding multi-temporal orthophoto image. The use of equally spaced sampling centered on candidate valley lines ensures that the subsequent calculations of geometric valley centers and specular flash saliency centers are based on the same spatial reference, avoiding calculation errors caused by spatial misalignment. After constructing the one-dimensional elevation profile, the sampling point corresponding to the minimum elevation value is searched in the joint elevation field sampling sequence corresponding to the one-dimensional profile. If there are multiple elevation minima in the profile, the minima closest to the origin is selected as the effective minimum point. The position of the minimum point in the local coordinate system is determined as the geometric valley center. By determining the geometric valley center through the minimum search of the one-dimensional profile, the true geometric concavity center of the membrane wrinkles along the transverse direction can be accurately located, providing a stable geometric reference for the subsequent calculation of mirror flash migration.

[0030] It should be noted that the geometric valley center refers to the position of the minimum elevation determined in the joint elevation field corresponding to a one-dimensional profile established along the corresponding cross-sectional direction at a certain location in the candidate valley line set. This position is used to characterize the true geometric center of the local concave structure of the flexible membrane surface in the cross-sectional direction. It is the core benchmark reflecting the actual spatial position of the fold valley line. Since the geometric valley center comes from the minimum elevation search result of the joint elevation field, it mainly reflects the true concave distribution of the membrane surface morphology, rather than the apparent bright center formed by the influence of local brightness changes.

[0031] After extracting the one-dimensional brightness profile, the local mean of the full sampling sequence of the brightness profile is calculated. The formula for calculating the local mean is as follows: ; In the formula, This represents the local mean of the current brightness profile. This represents the total number of sampling points in the brightness profile. Let be the brightness value corresponding to the i-th sampling point. The purpose of using the full-profile sampling sequence to calculate the mean is to accurately distinguish between the brightness abnormality area caused by mirror flash and the normal background area of ​​the film surface, based on the background brightness of the entire cross-section, thus avoiding the reference deviation caused by local window calculation. The difference between the brightness value of each sampling point in the brightness profile and the local mean is calculated. The difference result is truncated to non-negative, retaining only the portion where the brightness value is higher than the local mean. This portion is determined as the significant part, and the formula for calculating the significant part is: ; In the formula, is the salient value corresponding to the i-th sampling point. The meanings of the other parameters are the same as those in the aforementioned formula. Here, a non-negative truncation method is used to extract the salient part, which can remove areas with brightness lower than the background mean and retain only the high-brightness abnormal areas corresponding to the specular flare band, thus eliminating the interference of irrelevant background on the calculation of the specular salient center. If the sum of the significant portions is less than a preset lower limit of brightness, it is considered that there is no stable specular flare band in the profile. In this case, the specular flare salient center is directly set as the geometric trough center, and the specular flare migration at that location in that phase is set to zero. The preset lower limit of brightness is one-thousandth of the dynamic range of the brightness of the orthophoto image in the current phase. This processing method is used to solve the problem that the specular flare salient center cannot be reliably calculated when the brightness abnormality is not significant due to cloudy weather, weak reflection, or local film contamination. After extracting the salient parts, the cross-sectional distance is calculated using the geometric trough center as a reference and the values ​​of the salient parts as a weighted integral to obtain the mirror saliency center. The formula for calculating the mirror saliency center is as follows: ; In the formula, Let be the location of the salient center of the mirror flash in the local coordinate system. Let be the x-intercept distance of the i-th sampling point in the local coordinate system. The value is a preset minimum positive number, which is one ten-thousandth of the current dynamic range of the orthophoto's brightness. This is used to avoid calculation anomalies where the denominator is zero. The meanings of the other parameters are the same as those in the aforementioned formula. The core purpose of using the significant part values ​​as weights to perform weighted summation on the cross-sectional distance is to calculate the brightness centroid position of the specular flare band region. This position is the significant center of the specular flare optical characteristics, which can accurately characterize the actual distribution center of the specular flare band in the cross-sectional direction and form a clear positional correspondence with the center of the true geometric valley.

[0032] It should be noted that the saliency center of mirror flare refers to the center location obtained by extracting the brightness profile from multi-temporal orthophotos along the corresponding cross-sectional direction at a certain position in the candidate valley line concentration, taking the significant portion of the brightness profile above the local mean as the effective brightness anomalous region, and performing a weighted integral of the cross-sectional distance based on this significant portion. This location is used to characterize the concentrated distribution position of mirror flare bands or high brightness anomalous regions in the cross-sectional direction, and is a quantitative expression of the spatial landing point of mirror flare anomalies. The saliency center of mirror flare mainly reflects the distribution center of high brightness regions under illumination reflection conditions, so its position may not be consistent with the true geometric concavity center of the film surface.

[0033] The positional difference between the salient center of the specular flash and the center of the geometric valley in the same local coordinate system is determined by subtracting the positional value of the salient center of the specular flash from the center of the geometric valley in the same local coordinate system. This difference is defined as the specular flash migration amount. The sign of the specular flash migration amount represents the migration direction of the salient center of the specular flash relative to the center of the geometric valley in the cross-sectional direction, while the absolute value represents the migration distance. This amount is used to characterize the degree of offset of the specular flash optical anomaly relative to the true geometric valley center in the cross-sectional direction. Its sign represents the offset direction, and its absolute value represents the offset distance. The larger the specular flash migration amount, the more obvious the deviation of the optical highlight anomaly from the true concave structure, indicating a higher degree of specular flash interference at that location. The smaller the specular flash migration amount, the closer the specular flash anomaly is to the true geometric center, and the stronger the consistency between the optical performance and the geometric structure at that location.

[0034] In one embodiment of the present invention, a valley line locking weight is constructed based on the fluctuation of the geometric valley center and the fluctuation of the mirror flash migration over multiple time phases, including: First, the spatial reference of all time-phase calculation results is unified. The geometric valley centers and mirror flash migrations calculated for each time phase are mapped to a unified spatial coordinate system, ensuring consistent spatial resolution and coordinate correspondence across all time phases, with each grid position corresponding to a unique fixed spatial coordinate. For each grid position in the unified spatial coordinate system, the calculation results of all time phases are iterated, extracting the effective geometric valley center position value and effective mirror flash migration value for that position in each time phase. The effective geometric valley center position value is the coordinate value of that position along the corresponding transverse direction in the unified spatial coordinate system for that time phase. An effective time phase satisfies the following conditions: that position in that time phase has an effective candidate valley line, an effective one-dimensional profile, an effective geometric valley center, and an effective mirror flash migration. To avoid statistical distortion of variance due to insufficient effective samples, the minimum number of effective time phases is set to 60% of the larger of the total number of time phases and three, rounded up. If the number of effective time phases at the current position is lower than this minimum number of effective time phases, the valley line locking weight is not directly reset to zero. Instead, it is handled according to the following rules: if the spatial distance between the current position and the nearest effective valley line position is no more than two grid cells, the valley line locking weight of the nearest effective valley line position is used as the initial value, and then multiplied by the distance decay coefficient to obtain the valley line locking weight of the current position; if the distance exceeds the above range, the valley line locking weight of the current position is reset to zero. The distance decay coefficient is the normalized result of the reciprocal of the Euclidean distance between the current position and the nearest effective valley line position. This processing method can avoid large-area weight holes in boundary areas or locally occluded areas due to insufficient effective samples, while also preventing unsupported data from being forcibly included in high-confidence areas. After filtering the valid data, calculate the first variance of the geometric trough center position values ​​for all valid time phases at that location. Use this first variance as the fluctuation of the geometric trough center. The formula for calculating the first variance is: ; In the formula, Let M be the fluctuation amount at the center of the first variance, i.e., the geometric trough, and M be the number of valid time phases at the current position. Let be the geometric trough center position value of the current location in the t-th valid time phase. This is the arithmetic mean of the geometric trough center position values ​​of all valid time phases at the current location. Here, variance is used to statistically measure the fluctuation of the geometric trough center. The core purpose is to quantify the dispersion of the geometric trough center across multiple time phases through variance. The smaller the variance value, the more stable the position of the geometric trough center at this location is across different time phases, and the more it conforms to the fixed geometric structure characteristics of the true membrane surface fold valley line. The larger the variance value, the more obvious the fluctuation of the geometric position at this location, and the more likely it is a false result caused by noise or unstable structure. For the same grid location, extract the mirror flash migration values ​​for all valid time phases, and calculate the second variance of the mirror flash migration values ​​for all valid time phases at that location. Use this second variance as the fluctuation of the mirror flash migration. The formula for calculating the second variance is: ; In the formula, This represents the fluctuation of the second variance, i.e., the mirror flash migration. Let be the mirror flash migration value at the current position in the t-th valid time phase. This is the arithmetic mean of all valid temporal mirror flash migration values ​​at the current location. The meanings of the other parameters are the same as those in the aforementioned formula. Here, variance is used to statistically measure the fluctuation of mirror flash migration. The core purpose is to quantify the degree of migration dispersion of mirror flash bands under multiple temporal phases through variance. The larger the variance value, the more obvious the mirror flash migration phenomenon at this location is, the more it conforms to the characteristics of time-varying optical anomalies, and the higher the degree of interference from mirror flash. The smaller the variance value, the more stable the mirror flash characteristics at this location are, and the lower the degree of interference from optical anomalies. Taking the reciprocal of the first variance yields the first reciprocal type quantity. The formula for calculating the first reciprocal type quantity is: ; In the formula, For the first reciprocal type quantity, The value is a preset minimum positive number, which is one-thousandth of the statistical mean of the first variance of all valid locations in the entire region. This is used to avoid the calculation abnormality of the denominator being zero. The meanings of the other parameters are the same as those in the aforementioned formula. The core purpose of using the reciprocal transformation here is to reverse the numerical relationship of the geometric fluctuation, so that the more stable the geometric trough center is, the larger the corresponding first reciprocal value is, which can obtain a higher contribution ratio in the subsequent weight calculation, and achieve priority reinforcement of stable geometric structures. Dividing the first reciprocal variable by the sum of the first reciprocal variable and the second variance yields the valley locking weight at the current position. The formula for calculating the valley locking weight is as follows: ; In the formula, The valley lock weight is the valley line locking weight at the current position. The meanings of the other parameters are the same as those in the aforementioned formula. The core purpose of using this fractional structure to calculate the valley lock weight is to jointly quantify the stability of the geometric structure and the volatility of mirror flash migration. This makes the value of the valley lock weight positively correlated with geometric stability and negatively correlated with the degree of mirror flash volatility. When the value of the valley lock weight is larger, it indicates that the geometric structure at this position is stable and minimally affected by mirror flash, making it a highly reliable true valley line position. When the value of the valley lock weight is smaller, it indicates that the geometric structure at this position is unstable or greatly affected by mirror flash migration, making it unsuitable as a reliable benchmark for subsequent fusion processing. All grid positions in the entire area are calculated according to the above process to obtain the valley lock weight covering the entire area to be processed.

[0035] It should be noted that the valley line locking weight refers to a weight parameter calculated based on the fluctuation of the geometric valley center and the fluctuation of mirror flash migration at the same spatial location across multiple temporal phases. This weight is used to characterize the reliability of the spatial location as a stable location of the true valley line. The smaller the fluctuation of the geometric valley center, the more stable the concave structure center corresponding to this location across multiple temporal phases, and the stronger its ability to represent the true valley line. The larger the fluctuation of the mirror flash migration, the more significant the influence of optical highlight anomalies on this location, and the stronger the interference on the true valley line location. By jointly quantifying the stability of the geometric valley center and the degree of change in mirror flash migration, the valley line locking weight can simultaneously reflect the geometric stability and the strength of mirror flash interference. The larger the valley line locking weight, the closer the location is to the true, stable, and a high-reliability valley line location that can be used for subsequent global valley center line solving and multi-source data fusion. The smaller the valley line locking weight, the more likely the location is to be affected by mirror flash migration interference or has unstable geometric location, and it is not suitable to be directly used as the core benchmark for fusion processing.

[0036] In one embodiment of the present invention, the global valley centerline is extracted using a particle swarm optimization algorithm based on the valley centerline locking weight, including: The valley line intensities calculated for each time phase are all mapped to a unified spatial coordinate system to ensure that the valley line intensities for all time phases have a consistent spatial resolution and grid coordinate correspondence. For each grid position in the unified spatial coordinate system, the valley line intensities of all valid time phases are traversed, and the arithmetic mean of the valley line intensities of all valid time phases at that position is calculated to obtain the average valley line intensity. The formula for calculating the average valley line intensity is as follows: ; In the formula, The average valley line intensity at spatial location x. Let x be the number of effective time phases at spatial location x. Let be the valley line intensity at spatial location x in the t-th effective time phase. The core purpose of using multi-time phase arithmetic average to calculate the average valley line intensity is to weaken the interference of single-time phase noise on the valley line features through time-series averaging, strengthen the real valley line features that exist stably throughout the entire time phase, and provide a stable feature basis for the subsequent extraction of the global center line. Based on the control point set, a continuous curve is defined. Principal component analysis is performed on all segments of the candidate valley line set. The direction corresponding to the first principal component is taken as the main extension direction of the candidate valley line set. Then, the total extension length of all segments of the candidate valley line set is calculated. Initial control points are evenly spaced along the main extension direction, with the spacing between control points being one percent of the total extension length. The number of control points is no less than five and no more than fifty. After the layout is completed, a control point set containing multiple control points is formed. Each control point corresponds to a two-dimensional plane coordinate in a unified spatial coordinate system. Based on the control point set, a uniform cubic B-squared curve is used. Spline curve fitting generates a continuous and smooth curve. The B-spline curve is of order four, and the node vector adopts a uniform node sequence from zero to one. The first and second derivatives of the curve are continuous, and the sampling step size of the curve is consistent with the ground sampling distance of the orthophoto in a unified spatial coordinate system. The core purpose of using a uniform cubic B-spline curve to define a continuous curve based on control points is to control the overall shape of the curve with a small number of control points, ensuring the spatial continuity and smoothness of the curve, avoiding curve breakage and shape distortion caused by point-by-point optimization, and reducing the solution dimension of subsequent optimization algorithms. The continuous curve is discretized into multiple sampling points with equal arc length intervals according to the sampling step size. The valley locking weight and average valley intensity of the corresponding spatial location of each sampling point are obtained. The trapezoidal numerical integration method is used to calculate the line integral of the product of the valley locking weight and the average valley intensity at each point on the continuous curve along the curve. Then, the integral result is combined with the smoothness penalty term of the continuous curve to obtain the objective function. The formula for calculating the objective function is as follows: ; In the formula, For the set of control points The corresponding objective function value, Let K be the set of control points to be optimized. Let K be the two-dimensional coordinates of the k-th control point. For the set of control points The fitted continuous curve, Weights are locked for the valley line at spatial location x. Let be the differential arc length of the curve. The smoothing coefficient, set to one percent of the maximum average valley line intensity across the entire region, is used to balance curve fitting accuracy and smoothness. The core logic of constructing the objective function here is that the first term is the fitting term, used to drive the curve to converge towards regions with high valley line locking weights and high average valley line intensity, making the curve closely resemble the spatial distribution of the actual valley lines. The second term is the smoothness penalty term, which constrains the curvature change of the curve through the second-order difference of control points, avoiding sharp inflection points and breaks in the curve, and ensuring the spatial continuity of the global valley locking centerline. In this embodiment, the objective function is solved using a maximization method; if minimization is used as the default solution method, the objective function is inversely calculated before solving. The control point set is iteratively optimized using a particle swarm optimization (PSO) algorithm. During PSO initialization, each particle corresponds to a set of control points, with the particle dimension being twice the number of control points. Each dimension corresponds to the x-coordinate or y-coordinate of a control point. The PSO size is adaptively determined based on the number of control points, specifically ten times the number of control points, and limited to between twenty and two hundred. The initial particle position is randomly generated within a ±5% neighborhood of the initial control point coordinates, and the initial velocity is set to zero. Simultaneously, the particle position boundary is set to the coordinate range of the region to be processed, and the velocity boundary is set to 1% of the position boundary range. When a particle's position or velocity exceeds its corresponding boundary, boundary truncation is used for constraint. The individual historical best position of each particle is initialized as its initial position, and the swarm's historical best position is initialized as the initial position of the particle with the largest objective function value in the PSO swarm. To prevent particles from being concentrated in a small neighborhood in the initial stage, which would lead to insufficient search, after the initial position is generated, the average Euclidean distance between all particles is calculated. If the average Euclidean distance is less than one percent of the diagonal length of the initial control point set, the initial positions of the particles are randomly generated again until the above dispersion requirement is met. The maximum number of iterations for the algorithm is set to twenty times the number of control points, with a minimum of one hundred and a maximum of one thousand iterations. Convergence is determined using a dual condition: termination occurs when the maximum number of iterations is reached; or, the algorithm is considered convergent and terminates when the relative change in the objective function of the group's historical best position is less than one-thousandth for ten consecutive iterations, and the average displacement of the set of historical best control points is less than one ground sampling distance. By simultaneously constraining both the change in the objective function and the displacement of the control points, the problem of local curve drift caused by relying solely on the change in the objective function for convergence determination can be avoided. In each iteration of the particle swarm optimization algorithm, for each particle in the swarm, the objective function value corresponding to the particle's current position is calculated. Then, the absolute deviation between the objective function value at the current position and the objective function value at the particle's historical best position is calculated. Simultaneously, the absolute deviation between the objective function value at the current position and the objective function value at the group's historical best position is also calculated. Subsequently, the above three calculation results are processed based on the softmax normalized exponential function to adaptively generate the inertia term weight, individual cognitive term weight, and group social term weight in the velocity update equation. The calculation formula for the adaptive generation of weights is as follows: ; In the formula, k is the current iteration number, and i is the current particle number. Let be the weight of the inertia term for the i-th particle in the k-th iteration. Let be the weight of the individual cognitive term of the i-th particle in the k-th iteration. The group social term weight of the i-th particle in the k-th iteration. Let i be the objective function value at the current position of the i-th particle in the k-th iteration. Let be the individual historical best position of the i-th particle in the k-th iteration. Let be the group's historical best position in the k-th iteration, and softmax be a normalized exponential function. Specifically, it's calculated by performing exponential operations on the three input values, then dividing each result by the sum of the three results to obtain three non-negative weights. The core purpose of using adaptive weight generation is to avoid the problems of poor convergence and easy getting trapped in local optima caused by manually pre-setting fixed weights. This allows the algorithm to automatically adjust its iteration strategy based on the particle's current fitness state. When the particle's current fitness is high, the inertia term weight is increased to maintain the particle's iteration direction. When the particle's current position deviates little from the individual best, the individual cognition term weight is increased to strengthen the guidance of individual experience. When the particle's current position deviates little from the group best, the group social term weight is increased to strengthen the guidance of group experience, thus improving the algorithm's convergence speed and global optimization ability. If the objective function value of a particle's current position does not improve for five consecutive iterations, a random perturbation is applied to the particle's velocity vector while maintaining the group's historical best position. The perturbation amplitude is taken as 0.5% of the current position's boundary range to improve the ability of locally stagnant particles to escape. After completing the adaptive calculation of the three weights, the velocity and position of each particle are updated based on the velocity update equation, which includes inertial term weights, individual cognitive term weights, and group social term weights. The velocity update equation and the position update equation are as follows: ; ; In the formula, Let be the velocity of the i-th particle in the k-th iteration. and The value is a uniformly distributed random number between zero and one, used to increase the global search capability of the algorithm. The meanings of the other parameters are the same as those in the aforementioned formula. The core purpose of using this velocity update equation is to combine the current motion state of the particle, the individual historical optimal experience and the global optimal experience of the group to realize the iterative update of the particle position and drive the particle swarm to converge towards the direction of the optimal objective function value. After each iteration updates the velocity and position of all particles, boundary constraints are applied to the updated particle positions and velocities. For each particle, if the objective function value of its current position is greater than the objective function value of its individual historical best position, then the individual historical best position is updated to its current position. For the entire particle swarm, if the objective function value of any particle in the current iteration is greater than the objective function value of the swarm's historical best position, then the swarm's historical best position is updated to that particle's current position. Furthermore, if the angle between the optimized global valley centerline and the main extension direction of the candidate valley line set exceeds 45 degrees, or if the proportion of the global valley centerline's length outside the candidate valley line set's envelope exceeds 20%, the optimization result is considered abnormal. In this case, the control point set corresponding to the swarm's second-best historical position is used as the replacement result, and five more local iterations are performed for fine-tuning. This abnormal rollback process prevents the particle swarm algorithm from generating centerline results that significantly deviate from the true valley line distribution due to interference from local abnormal data. Repeat the above iterative process until the preset maximum number of iterations is reached or the convergence threshold is met. After the iteration terminates, the set of control points corresponding to the final historical best position of the population is taken as the optimal control point set, and the uniform cubic B-spline curve generated by fitting the optimal control point set is taken as the global valley centerline.

[0037] It should be noted that the global valley centerline refers to a continuous centerline within the processing area, obtained by optimizing a continuous curve defined by a set of control points based on valley line intensity and valley line locking weights across multiple temporal phases. This centerline is not an instantaneous linear result directly determined by a single temporal phase or a single local location, but rather a globally optimal linear structure obtained by combining the stability of valley line responses across multiple temporal phases, the consistency of geometric structure, and the smoothness constraints of the continuous curve. The global valley centerline reflects the continuous extension path of the real wrinkled valley lines on the flexible membrane surface throughout the processing area. It serves as a unified structural benchmark for subsequently determining the reference cross-sectional direction, implementing mirror flash stripping image fusion, and conducting elevation fusion processing. The closer the global valley centerline is to the region with high valley line intensity and high valley line locking weights, the more accurate its representation of the real valley line morphology, and the higher its reliability as a reference benchmark for subsequent multi-source fusion.

[0038] In one embodiment of the present invention, using the global valley centerline as a reference, and combining the valley line locking weight, data fusion is performed on the multi-temporal orthophoto, the multi-temporal photogrammetric elevation, and the multi-temporal laser elevation to obtain a fused image and a fused elevation, including: To achieve spatial benchmark unification for all data, the global valley centerline, valley locking weights, orthophotos of all time phases, and mirror flash migrations are all mapped to the same spatial coordinate system. This ensures that all data have consistent spatial resolution and raster coordinate correspondence, with each raster position corresponding to a unique fixed spatial coordinate. The global valley centerline is sampled at equal intervals, with the sampling step size matching the ground sampling distance of the orthophoto, resulting in a continuous sequence of sampling points covering the entire centerline. For each sampling point, the curve tangent direction at that location is calculated using the coordinate difference between that sampling point and its adjacent sampling points. This tangent direction represents the local extension direction of the global valley centerline at that location. The direction perpendicular to this local extension direction is determined as the reference cross-sectional direction. The unit vector corresponding to the reference cross-sectional direction is calculated and bound to the corresponding spatial position for storage. For each raster position in the unified spatial coordinate system, the reference cross-sectional direction unit vector corresponding to the global valley centerline sampling point is matched using the nearest neighbor principle to obtain the reference cross-sectional direction at that spatial position. This completes the assignment of reference cross-sectional directions for all raster positions across the entire region. For each valid temporal phase and each spatial grid location, a grid-by-grid calculation of the corrected image is performed. A valid temporal phase is defined as a phase where there is both valid mirror flash migration and valid orthophoto value, and the number of valid temporal phases is not less than two-thirds of the total number of temporal phases. The first component is obtained by multiplying the orthophoto value of the current spatial location at the current temporal phase by the valley locking weight at that location. The formula for calculating the first component is as follows: ; In the formula, Let x be the first component at the t-th temporal spatial location. Weights are locked for the valley line at spatial location x. Let x be the orthophoto value at the t-th temporal spatial location. The core purpose of constructing the first component here is to retain the original image information of high-confidence locations based on the valley line locking weight. The higher the valley line locking weight, the less the location is affected by mirror flash migration, the higher the confidence of the original image value, and the larger its proportion in the corrected image. After calculating the first component, two symmetrical positions are determined along the reference cross-section, centered on the current spatial location, based on the mirror flash migration amount at the current time phase. The positive symmetrical position is obtained by offsetting the absolute value of the mirror flash migration along the positive direction of the unit vector along the reference cross-section, using the current spatial location as the origin. The negative symmetrical position is obtained by offsetting the absolute value of the mirror flash migration along the negative direction of the unit vector along the reference cross-section. For both positive and negative symmetrical positions, if the position coordinates correspond to the raster center of the orthophoto image, the orthophoto image value of the corresponding raster is directly read. If the position coordinates do not correspond to the raster center, bilinear interpolation is used to obtain the orthophoto image value at that location. If the symmetrical position exceeds the boundary of the area to be processed, the nearest valid raster image value in the neighborhood of the current location is used to fill the gap. The system does not automatically assume that forward-symmetrical and reverse-symmetrical positions are immune to mirror flicker interference. Instead, it further performs a valid background check. The valid background check includes: determining whether the image brightness at the symmetrical position is higher than the upper limit of brightness determined by the sum of the mean and standard deviation of the background in the neighborhood at that position; if it is higher than the upper limit, the symmetrical position is determined to still be affected by mirror flicker interference; if it is not higher than the upper limit, it is determined to be a valid background position. When both the forward and reverse symmetrical positions are valid background positions, the arithmetic mean of their image values ​​is used to construct the second component. When only one symmetrical position is a valid background position, the image value of that valid background position is used directly to construct the second component. When both symmetrical positions are determined to still be affected by mirror flash, the search continues to expand to both sides along the reference transverse direction with a step size of one ground sampling distance, centered on the current position, until the first pair of valid background positions is found, or the maximum search radius is reached. The maximum search radius is the larger of three times the absolute value of the current mirror flash migration and five ground sampling distances. If no valid background position is found within the maximum search radius, the second component of the current position at the current time is set to zero, and only the first component is retained to participate in the correction image calculation for that time phase. The formula for calculating the second component is: ; In the formula, Let x be the second component at the t-th temporal spatial position. Let x be the mirror flash migration at the t-th temporal spatial location. Let x be the unit vector of the reference transverse direction at spatial location x. The absolute value of the mirror flash migration is given, and the meanings of the other parameters are the same as those in the formula above. The core purpose of constructing the second component is to compensate for the low confidence position that is severely affected by mirror flash migration by using the background image in the reference cross direction. Through the above effective background check and extended search mechanism, the assumption that there is no mirror flash at the symmetrical position is no longer relied upon, thereby ensuring that mirror flash removal and image correction can still be achieved in complex scenes. The first component and the second component are added grid by grid to generate the corrected image corresponding to the current time. The formula for calculating the corrected image is: ; In the formula, This is the corrected image value at the t-th temporal spatial location x, and the meanings of the other parameters are the same as those in the aforementioned formula. After completing the correction image calculation for all valid time phases according to the above process, the arithmetic mean of the correction images for all time phases is calculated grid-by-grid at the same spatial location to obtain the final fused image. The calculation formula for the fused image is as follows: ; In the formula, The fused image value at spatial location x. The number of valid temporal phases at the current position is given. The meanings of the other parameters are the same as those in the aforementioned formula. The core purpose of using the arithmetic mean of multi-temporal corrected images to calculate the fused image is to further weaken the random noise and optical interference of single-temporal residuals through the fusion of multi-temporal data, enhance the texture features of the true structure of the film surface, and obtain a fused image result that is temporally stable, has clear boundaries, and is free from mirror flicker interference. The calculated valley line locking weights for the entire region, the multi-temporal photogrammetric elevations and multi-temporal laser elevations for all time phases are all mapped to a unified spatial coordinate system to ensure that all data have a consistent spatial resolution and grid coordinate correspondence. Each grid position corresponds to a unique fixed spatial coordinate. At the same time, invalid values ​​of all data are preprocessed, and elevation values ​​that exceed the reasonable range and areas with no data are marked as invalid values. For each grid location in a unified spatial coordinate system, the effective time phases are first screened. An effective time phase is one where effective photogrammetric elevation, effective laser elevation, and effective valley line locking weights all exist simultaneously. The minimum number of effective time phases is determined using the same rules as described above. If the number of effective time phases at the current location meets the requirements, weighted elevation calculations are performed on all effective time phases. If the requirements are not met, a backtracking process is performed in the following order: the average valley line locking weight of the eight nearest effective grids in the neighborhood of the current location is calculated as the alternative valley line locking weight for the current location; this alternative valley line locking weight is applied to the photogrammetric elevation and laser elevation of each time phase at the current location to calculate the backtracked weighted elevation; the average of all backtracked weighted elevations is calculated as the fused elevation for the current location; if there are no fewer than three effective grids in the neighborhood of the current location, the arithmetic mean of the effective laser elevations of all time phases at the current location is directly used as the fused elevation. By using a hierarchical rollback method, the structural information of spatially adjacent locations can be used first to maintain the continuity of the fused elevation. Only when the neighborhood support is insufficient will it degenerate into the pure laser elevation mean result, thereby reducing the risk of elevation breakage caused by the lack of local samples. For each valid temporal phase and each spatial grid location, weighted elevation is calculated in parallel, grid-by-grid. The valley line of the current spatial location is locked as the weight, and multiplied by the multi-temporal photogrammetric elevation of the previous spatial location at the current time phase to obtain the first product. The formula for calculating the first product is as follows: ; In the formula, The first product at the t-th temporal spatial position x. Weights are locked for the valley line at spatial location x. Let x be the multi-temporal photogrammetric elevation at the t-th temporal spatial location. The core purpose of constructing the first product is to adaptively weight the photogrammetric elevation based on the valley line locking weight. The higher the valley line locking weight, the stronger the temporal stability of the geometric valley center, the less optical interference from the mirror flash migration zone, and the higher the degree of fit between the photogrammetric elevation and the real geometry of the membrane surface. Therefore, it is given a higher contribution ratio in the weighted elevation. After calculating the first product, the difference obtained by subtracting the valley line locking weight is used as the second multiplier. This second multiplier is then multiplied by the multi-temporal laser elevation of the current spatial position to obtain the second product. The formula for calculating the second product is as follows: ; In the formula, The second product at the t-th temporal spatial position x. Let x be the multi-temporal laser elevation at the t-th temporal spatial location. The meanings of the other parameters are consistent with the aforementioned formula. The core purpose of constructing the second product is to effectively supplement the location with low valley line locking weight using laser elevation. The lower the valley line locking weight, the more severe the interference from the mirror flare migration zone, and the more obvious the distortion caused by optical effects on photogrammetric elevation. Laser elevation is not affected by illumination conditions and mirror flare anomalies, and has stronger measurement stability for film surface geometry. Therefore, it is given a higher contribution ratio in the weighted elevation. Furthermore, when the absolute value of the difference between the photogrammetric elevation and the laser elevation at a certain time point exceeds twice the statistical median of the elevation differences of all valid time phases at that location, that time phase is marked as an elevation conflict phase. For elevation conflict phases, they are not directly included in the normal weighting; instead, the valley line locking weight of that phase is compressed to half of its original value before weighting calculation. This processing method can suppress the influence of outlier elevation values ​​caused by instantaneous wind loads, short-term surface water accumulation, or local mismatch anomalies on the final fused elevation.

[0039] After calculating the first and second products, the first and second products at the same time and spatial location are directly added together to obtain the weighted elevation corresponding to that time and location. The core logic of constructing the weighted elevation using this weighted summation method is to achieve adaptive fusion of photogrammetric elevation and laser elevation. It can automatically adjust the contribution ratio of the two elevation sources according to the degree of mirror flash interference at different locations on the membrane surface. In the stable geometric region of the real valley line, the high-resolution texture features of the photogrammetric elevation are preserved, while in the region with severe mirror flash interference, the stability and accuracy of the geometric measurement are ensured by relying on the laser elevation. After completing the weighted elevation calculation for all valid time phases across the entire region according to the above process, the arithmetic mean of the weighted elevations for all valid time phases is calculated grid-by-grid at the same spatial location to obtain the final merged elevation. The formula for calculating the merged elevation is as follows: ; In the formula, Let x be the fused elevation at spatial location x, and M be the number of valid temporal phases at the current spatial location. The core purpose of calculating the fused elevation using the arithmetic mean of multi-temporal weighted elevations is to further weaken the interference of wind load and hydrodynamic fluctuations on the membrane surface temporal deformation caused by single-temporal loads, as well as the influence of random noise in the elevation measurement process, through the fusion of multi-temporal data. The final result is a fused elevation that reflects the long-term stable geometry of the membrane surface, has strong temporal consistency, and reliable geometric accuracy. For the boundary locations of the area to be processed, if invalid data exceeding the boundary occurs during the weighted elevation calculation, the average elevation of the three nearest valid grid cells in the neighborhood of that location is used to fill the gap, ensuring the continuity and integrity of the fused elevation across the entire area.

[0040] like Figure 2 The diagram shown is a functional block diagram of a multi-source mapping data intelligent fusion processing system provided in an embodiment of the present invention. In this embodiment, the functions of each module / unit are as follows: The data acquisition module is used to acquire multi-temporal orthophotos, multi-temporal photogrammetric elevations, and multi-temporal laser elevations of the area to be processed on the flexible membrane surface. The valley line extraction module is used to construct a joint elevation field based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and to extract a set of candidate valley lines carrying cross-sectional direction information based on the joint elevation field, wherein the cross-sectional direction is the direction perpendicular to the local valley line extension direction. The offset measurement module is used to measure the geometric valley center and mirror flash migration along the cross direction corresponding to the candidate valley line set, respectively, on the joint elevation field and the multi-temporal orthophoto. The weight construction module is used to construct valley line locking weights based on the fluctuation of the geometric valley center and the fluctuation of the mirror flash migration in multiple time phases. The global solution module is used to extract the global valley centerline based on the valley line locking weight using the particle swarm optimization algorithm. The multi-source fusion module is used to fuse the multi-temporal orthophoto, multi-temporal photogrammetric elevation, and multi-temporal laser elevation data based on the global valley centerline and combined with the valley line locking weight, to obtain fused images and fused elevations.

[0041] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for intelligent fusion processing of multi-source surveying and mapping data, characterized in that, Includes the following steps: S1. Acquire multi-temporal orthophotos, multi-temporal photogrammetric elevations, and multi-temporal laser elevations of the area to be processed on the flexible membrane surface; S2. Construct a joint elevation field based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and extract a set of candidate valley lines carrying cross-sectional direction information based on the joint elevation field, wherein the cross-sectional direction is the direction perpendicular to the local valley line extension direction; S3. Along the cross-sectional direction corresponding to the candidate valley line set, determine the geometric valley center and mirror flash migration on the joint elevation field and the multi-temporal orthophoto image, respectively. S4. Based on the fluctuation of the geometric valley center and the fluctuation of the mirror flash migration under multiple time phases, construct valley line locking weights; S5. Extract the global valley centerline based on the valley line locking weight using the particle swarm optimization algorithm; S6. Using the global valley centerline as a reference, and combining the valley line locking weight, perform data fusion on the multi-temporal orthophoto, the multi-temporal photogrammetric elevation, and the multi-temporal laser elevation to obtain the fused image and fused elevation.

2. The intelligent fusion processing method for multi-source surveying and mapping data according to claim 1, characterized in that, A joint elevation field is constructed based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and a candidate valley line set carrying cross-sectional direction information is extracted from the joint elevation field, including: Calculate the absolute values ​​of the Laplace operator for the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, respectively; The first weighting coefficient and the second weighting coefficient are determined based on the reciprocal of the absolute value of the Laplacian operator of the multi-temporal photogrammetric elevation and the reciprocal of the absolute value of the Laplacian operator of the multi-temporal laser elevation, respectively. The first weighting coefficient is multiplied by the multi-temporal photogrammetric elevation, the second weighting coefficient is multiplied by the multi-temporal laser elevation, and the weighted sum is performed according to the normalized results of the first weighting coefficient and the second weighting coefficient to construct a joint elevation field. Determine the local valley line extension direction of the joint elevation field, and define the direction perpendicular to the local valley line extension direction as the transverse direction; Calculate the second derivative of the joint elevation field along the transverse direction, and determine the valley line intensity based on the negative value of the second derivative; Based on the continuous spatial extension of the valley line intensity, the location of the region carrying the cross-sectional direction information is extracted to form a candidate valley line set.

3. The intelligent fusion processing method for multi-source surveying and mapping data according to claim 2, characterized in that, Determining the center of the geometric valley on the joint elevation field includes: For each position in the candidate valley line set, a one-dimensional profile is established along the corresponding transverse direction; In the joint elevation field corresponding to the one-dimensional profile, the location corresponding to the minimum elevation is searched and used as the geometric valley center.

4. The intelligent fusion processing method for multi-source surveying and mapping data according to claim 3, characterized in that, Measuring mirror flicker migration on the multi-temporal orthophotos includes: Based on the cross-sectional direction, a brightness profile is extracted from the multi-temporal orthophoto image, and the local mean of the brightness profile is calculated. The portion of the brightness profile that is higher than the local mean is identified as a significant portion; Using the geometric valley center as a reference, the mirror flash saliency center is calculated by performing a weighted integral on the cross-sectional distance based on the salient portion; The positional difference between the salient center of the mirror flash and the center of the geometric trough is calculated as the mirror flash migration amount.

5. The intelligent fusion processing method for multi-source surveying and mapping data according to claim 1, characterized in that, Based on the fluctuations of the geometric trough center and the mirror flash migration over multiple time phases, a trough-locking weight is constructed, including: At the same spatial location, the first variance of the geometric trough centers in multiple time phases is statistically analyzed and used as the fluctuation of the geometric trough centers; At the same spatial location, the second variance of the mirror flash migration amount across multiple time phases is statistically analyzed and used as the fluctuation amount of the mirror flash migration amount; The first variance is transformed by its reciprocal to obtain the first reciprocal type quantity; Divide the first reciprocal type by the sum of the first reciprocal type and the second variance to obtain the valley line locking weight.

6. The intelligent fusion processing method for multi-source surveying and mapping data according to claim 5, characterized in that, Using the particle swarm optimization algorithm, the global valley centerline is extracted based on the valley locking weight, including: The average valley line intensity is calculated based on the valley line intensity under multiple defined time phases; Define a continuous curve based on a set of control points; Calculate the integral of the product of the valley locking weight and the average valley intensity at each point on the continuous curve, and combine the integral with the smoothness penalty term of the continuous curve to construct the objective function; The control point set is iteratively optimized using the particle swarm optimization algorithm. The continuous curve corresponding to the control point set that makes the objective function optimal is used as the global valley centerline.

7. The intelligent fusion processing method for multi-source mapping data according to claim 6, characterized in that, The control point set is iteratively optimized using the particle swarm optimization algorithm, including: In each iteration of the particle swarm optimization algorithm, for each particle, the objective function value at the current position, the deviation between the objective function value at the current position and the objective function value at the individual's historical best position, and the deviation between the objective function value at the current position and the objective function value at the group's historical best position are calculated. Based on the normalized exponential function, the deviation between the objective function value at the current position and the objective function value at the individual's historical best position, and the deviation between the objective function value at the current position and the objective function value at the group's historical best position are processed to adaptively generate the inertia term weight, individual cognitive term weight, and group social term weight in the velocity update equation. The velocity and position of each particle are updated based on the velocity update equation, which includes the inertial term weight, the individual cognitive term weight, and the group social term weight.

8. The intelligent fusion processing method for multi-source surveying and mapping data according to claim 1, characterized in that, Data fusion of the multi-temporal orthophotos includes: Obtain the reference transverse direction corresponding to the global valley centerline at the current spatial location; For each time phase and each spatial location, the orthophoto value of the current time phase is multiplied by the valley line locking weight to obtain the first component; In the reference cross-sectional direction, with the current position as the center, two symmetrical positions are determined according to the specular flash migration amount, and the mean value of the image values ​​corresponding to the two symmetrical positions is calculated. The difference obtained by subtracting the valley line locking weight is multiplied by the mean value to obtain the second component. The first component and the second component are added together to generate a corrected image for the corresponding time phase; The average value of the corrected images across all time phases is used to obtain the fused image.

9. The intelligent fusion processing method for multi-source surveying and mapping data according to claim 1, characterized in that, Data fusion of the multi-temporal photogrammetric elevation and the multi-temporal laser elevation includes: For each temporal phase and each spatial location, the valley line locking weight is used as the first multiplier and multiplied by the multi-temporal photogrammetric elevation to obtain the first product; The difference obtained by subtracting the valley line locking weight is used as the second multiplier and multiplied by the multi-temporal laser elevation to obtain the second product; Add the first product to the second product to obtain the weighted elevation for the corresponding time phase; The weighted elevations for all time phases are averaged to obtain the merged elevation.

10. A multi-source mapping data intelligent fusion processing system, used to execute the multi-source mapping data intelligent fusion processing method according to any one of claims 1-9, characterized in that, include: The data acquisition module is used to acquire multi-temporal orthophotos, multi-temporal photogrammetric elevations, and multi-temporal laser elevations of the area to be processed on the flexible membrane surface. The valley line extraction module is used to construct a joint elevation field based on the multi-temporal photogrammetric elevation and the multi-temporal laser elevation, and to extract a set of candidate valley lines carrying cross-sectional direction information based on the joint elevation field, wherein the cross-sectional direction is the direction perpendicular to the local valley line extension direction. The offset measurement module is used to measure the geometric valley center and mirror flash migration along the cross direction corresponding to the candidate valley line set, respectively, on the joint elevation field and the multi-temporal orthophoto. The weight construction module is used to construct valley line locking weights based on the fluctuation of the geometric valley center and the fluctuation of the mirror flash migration in multiple time phases. The global solution module is used to extract the global valley centerline based on the valley line locking weight using the particle swarm optimization algorithm. The multi-source fusion module is used to fuse the multi-temporal orthophoto, multi-temporal photogrammetric elevation, and multi-temporal laser elevation data based on the global valley centerline and combined with the valley line locking weight, to obtain fused images and fused elevations.

Images