Tower camera image high-precision positioning method based on sky-base cooperation

Abstract

The invention discloses a sky-base collaborative tower camera image high-precision positioning method, which comprises the following steps of: carrying out virtual squint projection on a satellite image based on a high-resolution satellite orthoimage of an area around a tower camera, and mapping the satellite image to a tower camera squint coordinate system; performing multi-scale phase fusion terrain enhancement feature extraction and cross-view matching on the satellite reference image and the tower camera image to obtain high-precision homonymy points and meter-scale control points; block adjustment and absolute positioning are carried out on a reference database image collected by a tower camera; and based on the positioned reference database, carrying out high-precision geographic positioning on the natural resource illegal behaviors. According to the method, high-precision positioning of the tower camera can be efficiently and automatically realized, the positioning efficiency of the tower camera image is remarkably improved, and protection of natural resources is facilitated.

Description

TECHNICAL FIELD

[0001] The application relates to the field of natural resource law enforcement, in particular to a space-based cooperative tower camera image high-precision positioning method. BACKGROUND

[0002] A natural resource stereoscopic monitoring system centered on multi-source visual perception technology is being rapidly constructed. Tower camera, as a wide-area monitoring node covering complex terrains such as forests and mountains, its meter-level spatial positioning capability is the key foundation to support natural resource "space-ground-land" collaborative observation and realize "global coverage, all-time response" intelligent monitoring.

[0003] The PTZ device relied on by the tower camera is often disturbed by environmental factors such as wind vibration and temperature change during long-term operation, resulting in an unavoidable drift of its initial positioning parameters. This drift often causes tens of meters or even hundreds of meters of error in spatial positioning directly relying on PTZ parameters in practice, thereby severely limiting its application feasibility in emergency response scenarios with extremely high positioning accuracy requirements. Although theoretically, the positioning accuracy can be improved to the meter level by means of visual positioning and other technical means, in actual operation, the fixed-view imaging characteristics of tower images make it difficult to directly apply the multi-view intersection strategy commonly used in unmanned aerial vehicle or satellite images, and therefore a series of technical obstacles still exist in high-precision positioning.

[0004] Satellite images, as high-precision geographic references, can provide wide-area and low-cost positioning benchmarks for tower images, but due to the fact that the former usually adopts top view orthographic projection while the latter is mostly near-ground oblique perspective imaging, there are essential differences in view angle and structure expression between the two, plus the differences in radiation characteristics and texture expression, making it difficult for traditional feature descriptors to work stably during matching, which seriously affects the matching accuracy and the generalization ability of the model.

[0005] The method proposes a space-based tower camera image high-precision positioning scheme: first, for the cross-view geometric-semantic matching problem of space-based heterogeneous images, a terrain-adaptive virtual oblique projection method is proposed, which reduces the geometric distortion caused by terrain undulation and view angle difference through digital elevation model fusion and dynamic mapping of the camera coordinate system; then, multi-scale phase moment principal component analysis and terrain feature enhancement are fused, combined with a hierarchical matching strategy, to effectively solve the mismatching problem caused by scale difference and affine deformation of heterogeneous images; next, a fixed optical center adjustment geometry model is constructed, the rotation parameters are optimized by measuring the residual error of the direction vector angle and Huber weighted loss, and the satellite image control points are combined to realize meter-level high-precision positioning of the tower image, thereby injecting technical innovation power into the construction of real scene three-dimensional China and empowering natural resource emergency monitoring and other needs. SUMMARY

[0006] The application provides a kind of tower camera image high-precision positioning method based on space-ground cooperation, which can automatically obtain meter-level control points of satellite orthographic image, and carry out high-precision positioning on tower image, efficiently monitor illegal behavior of natural resources, and empower emergency response monitoring.

[0007] In order to achieve the above-mentioned purpose, the application provides a kind of tower camera image high-precision positioning method based on space-ground cooperation, which includes the following steps: Step 1: Obtain high-resolution satellite orthographic image of the area around the tower as reference image, rotate 360° at the initial downward viewing angle to obtain the initial tower camera image, and combine the digital elevation model of the area with the initial downward viewing angle of the tower camera to perform virtual oblique projection on the satellite reference image, map the region of interest of the satellite image to the oblique coordinate system of the tower camera, and generate a virtual satellite image consistent with the initial downward viewing angle of the tower camera; Step 2: Perform multi-scale phase consistency analysis and occlusion-aware feature fusion on the satellite reference image and the tower camera image to generate an edge sketch of the satellite reference image and the tower camera image, and realize terrain-enhanced feature extraction; Step 3: Through a hierarchical matching framework, cross-view matching is performed on the satellite reference image and the tower camera image to obtain matching point pairs; Step 4: Control the tower camera to automatically adjust the downward viewing angle according to the preset sequence, and perform 360° rotation shooting at each downward viewing angle to collect a multi-angle image sequence covering the full-range scene around the tower, which together constitutes the original image of the reference dataset; Step 5: Based on the reference dataset, a tower image area network adjustment model with fixed optical center constraint is constructed to calculate the relative orientation parameters of adjacent images in the tower camera reference dataset; Among them, the relative orientation parameters of adjacent images specifically refer to the three-axis rotation angles of each image relative to the unified camera coordinate system under the condition of fixed optical center constraint; Step 6: Using the matching point pairs obtained in step 3, combining the reference dataset of step 4 and the relative orientation parameters of step 5, the meter-level control points are obtained and the absolute positioning parameters of the images in the reference dataset are calculated through space forward intersection and backward intersection calculation, the absolute positioning parameters refer to the position and attitude parameters required to convert the tower camera image from its own coordinate system to the global unified geographic coordinate system; Step 7: Match the tower camera image to be positioned with the images in the reference dataset, and realize high-precision positioning of the image to be positioned according to the known positioning parameters.

[0008] Further, the specific steps of the virtual oblique projection in step 1 include: A three-dimensional point set is constructed in the geographic coordinate system of satellite imagery, and the elevation value corresponding to each pixel is obtained by interpolation through a digital elevation model; Based on the camera's external parameters, points in the geographic coordinate system are transformed to the camera coordinate system using the following formula: In the formula, The elevation transformation matrix generated for the DEM is used to map geographic elevations to local terrain surfaces. Its calculation relies on bilinear interpolation or kriging interpolation algorithms. Based on camera tilt angle and azimuth Constructed rotation matrix; offset vector Used to align the optical center of the camera with the origin of the local coordinate system in the satellite image; This represents the three-dimensional coordinates of a point in the WGS84 geodetic coordinate system.

[0009] By combining the camera's internal parameters, three-dimensional points are projected onto the image plane to generate a virtual satellite image that matches the actual camera's oblique angle.

[0010] Furthermore, the terrain enhancement feature extraction method using multi-scale phase fusion described in step 2 includes: Step 2.1. Perform multi-level wavelet decomposition on satellite imagery and tower camera imagery to decompose the input imagery into sub-band images of different scales in order to capture terrain features of different sizes; A series of Log-Gabor filters with different orientation angles are created for sub-band images of various scales. For each scale of sub-band image, each Log-Gabor filter in the corresponding scale filter bank is convolved. The complex convolution kernel captures the amplitude and phase of the same pixel point in both the direction consistent with the filter direction and the direction perpendicular to the filter direction. Anisotropic structure is extracted, and a feature vector is generated for each pixel point. The feature vector is used to describe the local texture orientation distribution of the pixel point, representing the amplitude and phase of each pixel point in each direction and scale. In sub-band images of the same scale, the feature vectors of pixels in each direction need to be multiplied by the direction weights. We obtain the processed feature vectors in each direction, where the weights are... Calculated using the weighting function: in, For the first Angles in all directions For direction-weighted range, Used to traverse all directions It is a weighted function for tilt angle. For the camera angle; Step 2.2. Based on the magnitude and phase of the pixel in each direction at the same scale obtained in Step 2.1, construct a second-order moment tensor for each pixel. Perform feature decomposition on the second-order moment tensor of the pixel at each scale to obtain two eigenvalues ​​(λ1, λ2) and the corresponding eigenvectors (v1, v2), obtain the main direction and shape of the pixel, and identify the main response direction. Step 2.3. Based on the phase information obtained by convolution of the Log-Gabor filter bank in various directions and at various scales, and combined with the local principal direction, construct the phase consistency matrix of each pixel at different scales; Phase consistency matrices obtained at different scales are stacked according to their corresponding pixel positions to form a three-dimensional array. Principal component analysis is then used to reduce the dimensionality of the three-dimensional array and extract the first principal component as a simplified map of the main edge of the image.

[0011] Furthermore, the Log-Gabor filter bank is configured with: 5 wavelet decomposition layers, 22.5° directional spacing, and 30° directional weighting range, indicating that directions within ±30° of the camera direction have higher weights.

[0012] Furthermore, the hierarchical matching framework described in step 3 includes: Initial screening stage: Based on the key points and main direction information in the edge map, the rBRIEF descriptor is dynamically corrected in direction. The distribution of sampling points is adjusted by phase direction rotation. An enhanced rBRIEF descriptor is generated by adopting a non-uniform sampling strategy. Point pairs with weak similarity are excluded by bidirectional nearest neighbor and Hamming distance ratio. The spatial distribution uniformity of the point set is evaluated by using the Poisson distribution model. False responses in texture-mixed areas are removed to form a preliminary matching set. Fine matching stage: Taking each pair of matching points (a,b) from the preliminary matching point set output from the initial screening stage as the center, block images Xa and Xb are extracted respectively, centered on the coordinate positions in the simplified main edge maps of satellite and tower camera images. Block images Xa and Xb are composed of multiple concentric rings, multiple azimuth sectors, or multiple concentric rings and multiple azimuth sectors. Logarithmic polar coordinate transformation is performed to generate block images LXa and LXb respectively. The average gray value of each concentric ring, or each azimuth sector, or each concentric ring and each azimuth sector in block images LXa and LXb is calculated. The average gray values ​​of each concentric ring, or each directional sector, or each concentric ring and each directional sector are arranged in order to form a one-dimensional vector. The average gray values ​​of each concentric ring, or each directional sector, or each concentric ring and each directional sector in block image LXb are arranged in order to form another one-dimensional vector. The Log-Polar feature descriptors of block images LXa and LXb are obtained respectively. The Log-Polar feature descriptors have rotation and scale invariance. Among them, the block images Xa and Xb corresponding to the main edge simplified map of satellite and tower images are divided into the same concentric rings and directional sectors in the same way. Calculate the Pearson correlation coefficient of the Log-Polar feature descriptor corresponding to each pair of matching points (a,b) in the preliminary matching point set. Based on the DEM digital elevation model data or the three-dimensional information of the image itself, evaluate the topographic relief of the area where the matching points are located. For areas with flat terrain and areas with severe topographic relief, different similarity thresholds are used to tolerate feature differences caused by geometric distortion. Output a set of matching point pairs with similarity greater than the similarity threshold. Finally, the set of matching point pairs is used as input, and geometric consistency verification is performed in combination with terrain structure features to output a set of finely matched points with high positioning accuracy.

[0013] Furthermore, the parameter settings for the fine matching stage include: Hamming distance ratio threshold set to 0.7, Log-Polar feature radius set to 32 pixels, number of rings set to 5, number of sectors set to 12, and vector dimension set to 60.

[0014] Furthermore, step 5, constructing the fixed optical center constraint regional network adjustment model, includes: Assuming that the optical center of all images is the same fixed point, the projection relationship between images is simplified to a pure rotation relationship estimation; The first step is to establish an observation equation based on the consistency of direction vectors by taking the interior orientation elements of all images in the benchmark dataset and the selected pairs of corresponding points in adjacent images as input. By transforming the unit direction vector of each pair of corresponding points in the camera coordinate system of each image to a unified coordinate system, the included angle residual is calculated and the initial direction residual vector is output. The second step involves using the initial orientation residual vector as input to construct an adjustment objective function containing three core elements: (1) an orientation consistency residual term, which uses a weighted Huber loss function for robust control; (2) a rotation parameter regularization term, which constrains the deviation between the rotation matrix of each image and the initial pose; and (3) a relative rotation smoothing term, which maintains the continuity of rotation changes between adjacent images. The output is a constrained nonlinear least squares problem statement. The third step involves using the objective function and initial rotation parameters as inputs, and iteratively solving the problem using the Levenberg-Marquardt algorithm. In each iteration, the following are updated synchronously: (1) the three-axis rotation parameters of each image; (2) the adaptive weights of the Huber loss function; (3) the balance coefficients between the regularization term and the residual term; and the optimized precise rotation matrix of each image is output. The fourth step involves using the optimized rotation matrix as input to perform posterior accuracy assessment, calculating the mean square error of the orientation residuals for each pair of corresponding points, eliminating gross observations with residuals greater than 3 times the mean square error, and outputting the final relative orientation parameters after quality control.

[0015] Furthermore, the calculation of the absolute positioning parameters in step 6 includes: The matching point pairs are associated with the corresponding 3D point clouds in the benchmark dataset, and the high-precision positioning information of satellite reference imagery is used to provide a geographic reference control benchmark for pole imagery. Based on the geometric constraints established by the relative orientation parameters calculated in step 5, the entire image network is optimized as a whole. Meter-level control points with geographic coordinates are obtained through spatial forward intersection calculation, and then the absolute position and attitude of each image in a unified coordinate system are solved through back intersection optimization.

[0016] A computer-readable storage medium having a computer program stored thereon, which, when executed, implements the method.

[0017] Compared with the prior art, the present invention has the following advantages: 1. This invention can automatically acquire meter-level control points from high-resolution satellite orthophotos without the need for manual point marking, significantly improving the positioning efficiency of tower cameras.

[0018] 2. This invention improves the matching accuracy of satellite imagery and pole camera imagery by extracting terrain enhancement features through multi-scale phase fusion.

[0019] 3. This method is automated, easy to operate, and has very high overall efficiency. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in this application or related technologies, the accompanying drawings used in the description of the embodiments or related technologies will be briefly introduced below.

[0021] Figure 1 This is a flowchart illustrating the overall design of an embodiment of the present invention; Figure 2 This is a diagram showing the location of the camera; Figure 3 This is a schematic diagram of terrain enhancement features; Figure 4 This is a schematic diagram of cross-view image matching; Figure 5 A schematic diagram of the adjustment of the tower image area network with fixed optical center constraints. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0023] Example With the widespread application of cameras in various monitoring and surveillance scenarios of natural resources, and in response to the need for automatic monitoring and positioning of natural resources, this invention proposes a high-precision positioning method for tower camera images using a space-based collaborative approach. The scheme and process are attached. Figure 1 .

[0024] Step 1: As Figure 2 As shown, a pole-mounted camera in Dengyun Town, Longchuan County, Heyuan City, Guangdong Province was selected for the experiment. High-resolution satellite orthophotos of the area around the pole were obtained as reference images. By combining the ALOS digital elevation model (DEM) of the area with the camera's geometric parameters, the satellite images were virtually obliquely projected onto the oblique coordinate system of the pole-mounted camera. The downward viewing angle of the camera was set to 10 degrees, and the pole-mounted camera images were taken by rotating the downward viewing angle 360° to capture a full circle.

[0025] The specific method for virtual oblique projection of satellite imagery is as follows: First, a 3D point set is constructed in the geographic coordinate system of the satellite imagery, and the elevation value of each pixel is obtained through DEM interpolation. Then, a rotation matrix is ​​constructed using the camera's extrinsic parameters. And through the offset vector Align the optical center with the origin of the local coordinate system to complete the conversion from geographic coordinates to camera coordinates. The conversion formula is as follows: In the formula, The elevation transformation matrix generated for the DEM is used to map geographic elevations to local terrain surfaces. Its calculation relies on bilinear interpolation or kriging interpolation algorithms. Based on camera tilt angle and azimuth Constructed rotation matrix; offset vector Used to align the optical center of the camera with the origin of the local coordinate system in the satellite image; This represents the three-dimensional coordinates of a point in the WGS84 geodetic coordinate system.

[0026] After completing the spatial coordinate transformation, the camera's internal parameters (focal length, principal point coordinates) are further combined to project the three-dimensional points onto the image plane, generating a virtual satellite image that matches the actual camera's oblique angle.

[0027] Step 2: By fusing multi-scale phase consistency analysis and occlusion perception features, simplified edge maps of satellite reference images and tower camera images are generated. Multi-scale phase fusion is then used to extract terrain enhancement features, as shown in the following figure. Figure 3 As shown.

[0028] Multi-scale phase fusion terrain enhancement feature extraction methods include: Step 2.1. First, perform multi-level wavelet decomposition on the satellite imagery and the tower camera imagery to decompose the input imagery into sub-band images of different scales to capture terrain features of varying sizes. For each sub-band image of the same scale, create a series of Log-Gabor filters with different orientation angles. Convolve each Log-Gabor filter in the corresponding filter bank for each scale of the sub-band image to obtain the degree of matching between each pixel and the orientation and scale preferred by each filter. Generate a feature vector for each pixel. The feature vector is used to describe the local texture orientation distribution of the pixel, representing the amplitude and phase of each pixel in each orientation and scale.

[0029] To further enhance edge features consistent with the camera's viewpoint, the feature vectors of pixels in each direction in the sub-band image of the same scale obtained in step 2.1 need to be multiplied by directional weights. We obtain the processed feature vectors in each direction, where the weights are... Calculated using the weighting function: In the formula, For the first Angles in all directions For direction-weighted range, Used to traverse all directions It is a weighted function for tilt angle. This refers to the camera angle.

[0030] By applying weighting coefficients that match the camera's viewing angle to the responses in each direction, significant edges in the tilt direction are additionally enhanced. This weighting method effectively improves the edge response intensity of the target area and avoids feature interference from areas with inconsistent viewing angles.

[0031] In step 2.1, the Log-Gabor filter uses a complex convolution kernel during convolution to extract anisotropic structures. The complex convolution kernel simultaneously captures the amplitude and phase of the same pixel in both the direction aligned with and perpendicular to the filter direction. This corresponds to the center frequency of the Log-Gabor filter in a filter bank of the same scale. The settings are the same. The larger the value, the better the effect of capturing high-frequency information, which is used to capture fine edges and details in images. The smaller the value, the better the effect of capturing low-frequency information. It is used to capture coarse outlines and structures in the image. The filter enhances the edge information of the scale corresponding to its center frequency. In this embodiment, the filter bank is set to 5 groups, corresponding to the sub-band image after wavelet decomposition into 5 levels. Each filter in the same filter bank is set to have the same center frequency but different directional angles. The directional angles are set with the same sampling step size in the range of 0°~180°. In this embodiment, the sampling step size is 22.5°. Therefore, each filter in a filter bank convolves the image in 8 different directions. Thus, the filter parameters in this step are set as follows: wavelet decomposition of 5 levels, directional interval of 22.5°, and directional weighting range of 30°, indicating that directions within ±30° of the camera direction have higher weights.

[0032] Step 2.2. Based on the magnitude and phase of each pixel in each direction at the same scale obtained in Step 2.1, construct a second-order moment tensor for each pixel to characterize the energy distribution trend of each pixel in the spatial frequency domain. For each pixel in the image, at its respective scale, utilize the energy distribution tensor in different directions... The gradient response on the surface is used to construct a second-order moment tensor M. This tensor is a 2×2 symmetric positive semi-definite matrix, and its calculation formula is as follows: in, This represents the direction of the Log-Gabor filter; In direction Above, the gradient component of the pixel in the x-direction (horizontal); In direction Above, the gradient component of the pixel in the y-direction (vertical); It is the direction The weights are typically taken as the response amplitude in that direction.

[0033] For each pixel at each scale, a second-order moment tensor eigenvalue decomposition is performed to obtain two eigenvalues ​​(λ1, λ2) and the corresponding eigenvectors (v1, v2), thus acquiring the principal orientation and shape of the pixel. The eigenvalue decomposition process is as follows: Assumption: in: ; ; .

[0034] Eigenvalues Satisfies the characteristic equation: Right now: =0, expanded, we get: ( ()( )- =0 Solving this quadratic equation yields two eigenvalues: in, The larger eigenvalues ​​are the principal eigenvalues, corresponding to the energy in the principal direction. The smaller eigenvalue, i.e. the secondary eigenvalue, corresponds to the energy in the vertical direction.

[0035] For eigenvalues : Solve the equation Taking the first equation, we obtain the eigenvectors: For eigenvalues Since the eigenvectors of a symmetric matrix are orthogonal, they can be directly taken as follows: The eigenvalue decomposition of the second-order moment tensor reflects the degree of anisotropy of the local image structure, thereby enhancing specific feature characteristics and suppressing irrelevant textures. Specifically, for features with consistent orientation, such as ridges and building boundaries, the gradient directions are highly concentrated in space, resulting in a significantly larger principal feature value (λ1) than the secondary feature value (λ2), and an anisotropy ratio (λ1 / λ2) much greater than 1. These features are enhanced in subsequent processing. Conversely, for regions with scattered orientations and irregular texture changes, such as vegetation and cloud shadows, the gradient directions are statistically evenly distributed, leading to close values ​​for λ1 and λ2, and an anisotropy ratio approaching 1. These regions are naturally suppressed in feature representation, thereby improving the structural stability and clarity of the overall edge information.

[0036] In this step, the eigenvalues ​​of the second-order moment tensor reflect the energy intensity in the principal direction and the anisotropic direction. The principal direction points to the direction of the fastest energy change and corresponds to the eigenvector with a high eigenvalue. The anisotropic direction points to the direction of the slowest energy change and corresponds to the eigenvector with a low eigenvalue. The eigenvectors indicate the principal direction of the local energy distribution. The ratio of the principal direction eigenvalue to the anisotropic direction eigenvalue indicates the degree of anisotropy in the region—the larger the ratio, the more definite the directionality (such as linear edges); if the ratio is close to 1, the directionality is unclear.

[0037] Step 2.3. Obtain the phase consistency matrix of each pixel in the input image at different scales, and stack the phase consistency matrices at different scales according to the corresponding pixel positions to form a three-dimensional array. In order to preserve the main terrain structure while compressing redundant information, this invention uses principal component analysis (PCA) to reduce the dimensionality of the above three-dimensional array. PCA finds the projection direction with the largest data variance, i.e., the first principal component, by performing a linear transformation on the scale dimension. Finally, the first principal component values ​​of all pixels are reconstructed into a two-dimensional image, which is the simplified map of the main edges. Specifically, firstly, a three-dimensional array of size H×W×S (H and W are the image height and width, respectively, and S is the number of scales) is reshaped into a two-dimensional matrix X of size (H×W)×S, so that each row vector contains the phase consistency features of the same pixel at different scales. Next, the matrix X is centered along the column direction, that is, the average value of each scale feature is subtracted to obtain the mean-reduced matrix X'. Then, the S×S covariance matrix C of X' is calculated and eigenvalue decomposition is performed, and the eigenvectors are obtained by sorting them in descending order of eigenvalues. Finally, the eigenvector corresponding to the largest eigenvalue is selected as the projection axis, and X' is projected onto this direction to obtain a one-dimensional vector of length H×W. This vector is then reshaped into a two-dimensional image of size H×W, which is the main edge map.

[0038] This simplified map, while preserving the main information of the terrain structure, significantly reduces redundant texture interference, providing high-quality edge feature support for subsequent image matching or structural analysis.

[0039] Specifically, the phase consistency matrix is ​​constructed based on the phase information obtained by convolving the Log-Gabor filter bank in each direction and at each scale in step 2.1, and combined with the local principal direction (feature vector v1) obtained in step 2.2. This matrix quantifies the consistency of the local phase of a pixel near the dominant direction at different scales, and can effectively reflect the stable structural features across scales in the image.

[0040] The paper "Multi-directional Phase Consistency Local Registration Model for High-Resolution Optical Images" discloses a method for constructing the phase consistency matrix.

[0041] Step 3: Employ a hierarchical matching strategy to perform cross-view matching between satellite reference imagery and tower camera imagery, obtaining matching point pairs. Set the Hamming distance ratio threshold to 0.7, the Log-Polar feature radius to 32 pixels, the number of rings to 5, the number of sectors to 12, and the feature dimension to 60. The matching results are as follows: Figure 4 As shown.

[0042] The process of acquiring matching point pairs between satellite reference imagery and pole imagery includes: Initial screening stage: Feature points were collected from the edge contour maps of satellite and pole images using the FAST algorithm. The principal directions of the feature points in both edge contour maps were calculated using grayscale centroid localization. Standard BRIEF descriptors were generated for each. Based on the principal directions of the feature points in both edge contour maps, the coordinates of the sampling points in the BRIEF descriptors were rotated to align with the principal directions of the feature points in the corresponding edge contour maps, resulting in rBRIEF descriptors for both edge contour maps. The Hamming-Hamm distance between the rBRIEF descriptors of the two edge contour maps was calculated. For each feature point A in the satellite image edge contour map, the feature point B with the closest Hamming distance in the pole image edge contour map was found. Conversely, for feature point B in the pole image edge contour map, the feature point with the closest Hamming distance in the satellite image edge contour map was found. It was verified whether this was A. When A and B were each other's nearest neighbors, this pair of matching points (A, B) was considered a valid candidate matching point pair. The edge map of satellite and tower images is divided into a uniform grid. The number of matching points falling into each grid is counted. Based on the average number of matching points in all grids, a Poisson distribution model is established. If the number of matching points in any grid exceeds the upper limit of the confidence interval of the Poisson distribution (e.g., 95%), the matching point density in that grid is considered abnormal and is collectively removed. False matching points with abnormal distribution in areas with mixed textures are removed. Finally, a preliminary set of matching points that meets the spatial consistency requirements is output.

[0043] The steps for collecting key points and obtaining rBRIEF descriptors in the initial screening stage are similar to those described in CN114170450A, and the method for screening valid candidate matching points is similar to that described in CN114821130A.

[0044] Fine matching stage: Taking the coordinates of each pair of matching points (a,b) in the enhanced feature map (i.e., the simplified main edge map of satellite and tower images) in the preliminary matching point set output from the initial screening stage as the center, block images Xa and Xb are extracted respectively. Block images Xa and Xb are composed of multiple concentric rings, multiple azimuth sectors, or multiple concentric rings and multiple azimuth sectors. Logarithmic polar coordinate transformation is performed to generate block images LXa and LXb respectively. The concentric rings or azimuth sectors in block images LXa and LXb are calculated. The average gray values ​​of each concentric ring or directional sector in block image LXa are arranged sequentially to form a one-dimensional vector. Similarly, the average gray values ​​of each concentric ring or directional sector in block image LXb are arranged sequentially to form another one-dimensional vector. This yields the Log-Polar feature descriptors for block images LXa and LXb, respectively. The Log-Polar feature descriptors are rotation- and scale-invariant. Block images Xa and Xb, corresponding to the simplified main edge maps of satellite and tower images, are divided into the same concentric rings and directional sectors. The Pearson correlation coefficient of the Log-Polar feature descriptor for each pair of matching points (a,b) in the preliminary matching point set is calculated. Generally, a correlation coefficient greater than 0.75–0.85 indicates high similarity, while a correlation coefficient less than 0.5–0.6 indicates low similarity. Based on the 3D information of DEM (Digital Elevation Model) data or imagery, the terrain undulation of the area where the matching points are located is evaluated. For flat areas, a high and strict similarity threshold (e.g., 0.9) is used; for areas with dramatic terrain undulations, the similarity threshold is appropriately reduced (e.g., 0.75) to tolerate feature differences caused by geometric distortion. The set of matching point pairs with a similarity greater than the similarity threshold is output. Finally, this set of matching point pairs is used as input, and geometric consistency is verified by combining terrain structure features to output a final, highly accurate set of precise matching points.

[0045] Step 4: Control the tower camera to automatically adjust the downward viewing angle according to the preset sequence, and perform 360° rotation shooting under each downward viewing angle to collect multi-angle image sequences covering the all-round scene around the tower, which together constitute the original images of the benchmark dataset.

[0046] Step 5: Based on the benchmark dataset, construct a pole image area network adjustment model with fixed optical center constraints, and calculate the relative orientation parameters of adjacent images in the pole camera benchmark dataset. The model optimizes rotation parameters through orientation consistency residuals, uses a weighted Huber loss function to enhance robustness, and introduces a regularization term to constrain attitude deviations and ensure optimization stability. Figure 5 As shown.

[0047] In this invention, the relative orientation parameter of adjacent images (which includes images at the same pitch angle but in different orientations, as well as images at different pitch angles but with a large overlap) specifically refers to the three-axis rotation angle of each image relative to a unified camera coordinate system under a fixed optical center constraint. Since the optical center position of the tower camera remains fixed during shooting, the relative orientation between images is simplified to a pure rotational relationship estimation.

[0048] The construction of a fixed optical center constraint regional network adjustment model shall be carried out according to the following steps: The first step involves establishing an observation equation based on the consistency of orientation vectors, using the interior orientation elements of all images in the benchmark dataset and selected pairs of corresponding points from adjacent images as input. This step transforms the unit orientation vectors of each pair of corresponding points in the camera coordinate systems of each image to a unified coordinate system, calculates their angular residuals, and outputs the initial orientation residual vector. The angular residual refers to the angle difference between the unit orientation vectors of the same pair of corresponding points in the camera coordinate systems of two adjacent images after being transformed to the unified coordinate system using a rotation matrix. The initial orientation residual vector is a column vector composed of the angular residuals of all pairs of corresponding points arranged in order.

[0049] The second step involves using the initial orientation residual vector as input to construct an adjustment objective function containing three core elements: (1) an orientation consistency residual term, which uses a weighted Huber loss function for robust control; (2) a rotation parameter regularization term, which constrains the deviation between the rotation matrix of each image and the initial pose; and (3) a relative rotation smoothing term, which maintains the continuity of rotation changes between adjacent images. The output is a constrained nonlinear least squares problem statement.

[0050] The objective function expression for adjustment constructed in this step is as follows: In the formula, the first term is the orientation consistency residual term, where N is the number of images in the benchmark dataset; This is the set of matching pairs of points with the same name on the k-th image; For the pair of points with the same name, the unit direction vector is given in the camera coordinate system of the k-th image. Let be the rotation matrix of the k-th image to be solved; It is a weighted Huber loss function used to control error, applying smaller weights to larger residuals to improve model robustness.

[0051] The second term is the rotation parameter regularization term, where Let be the initial value of the rotation matrix for the k-th image; The Frobenius norm is used to measure the difference between two matrices. It is a regularization parameter that controls the weight of this term in the overall objective function.

[0052] The third term is the relative rotation smoothing term. It is the set of all image pairs (m,n) that are adjacent in time or space; It is a smoothing parameter that controls the intensity of rotational consistency between adjacent images.

[0053] The third step involves iteratively solving the problem using the aforementioned objective function and initial rotation parameters as input. In each iteration, the following are updated synchronously: (1) the three-axis rotation parameters of each image; (2) the adaptive weights of the Huber loss function; and (3) the balance coefficients between the regularization term and the residual term. The optimized, precise rotation matrices for each image are then output.

[0054] The fourth step involves using the optimized rotation matrix as input for posterior accuracy evaluation. The mean square error of the orientation residuals for each pair of corresponding points is calculated, outlier observations with residuals greater than three times the mean square error are removed, and the final relative orientation parameters, after quality control, are output.

[0055] Step 6: Using the matching point pairs of the tower camera and satellite reference images obtained in Step 3, calculate high-precision meter-level control points and calculate the absolute positioning parameters of each image in the benchmark dataset. The absolute positioning parameters refer to the position and attitude parameters required to transform the tower camera image from its own coordinate system to the globally unified geographic coordinate system. Specifically, they include the three-dimensional coordinates (X, Y, Z) of the image capture center and the rotation angles (pitch, yaw, roll) around the three axes.

[0056] In this step, the benchmark dataset constructed in step 4 and the relative orientation parameters obtained in step 5 together provide the necessary spatial benchmark and geometric constraints for absolute positioning: First, the matching point pairs obtained in step 3 are associated with the corresponding 3D point clouds in the benchmark dataset of step 4. Utilizing the high-precision positioning information from satellite reference imagery, a geographically referenced control benchmark is provided for the tower imagery. The benchmark dataset constructed in step 4 serves as a bridge for spatial coordinate transfer, transmitting the absolute coordinates of the satellite imagery to the tower camera coordinate system.

[0057] Secondly, the relative orientation parameters calculated in step 5 establish the precise rotational relationship between the tower images in the benchmark dataset, forming a regional network with strong geometric constraints. This constraint allows for the optimization of the entire image network as a whole when calculating control points, without having to solve the position parameters independently for each image, thus significantly improving the accuracy and stability of absolute positioning.

[0058] In practice, the aforementioned matching point pairs, benchmark datasets, and relative orientation parameters are used as inputs. Meter-level control points with geographic coordinates are obtained through spatial forward intersection calculations. Then, the absolute position and orientation of each image in a unified coordinate system are solved through back intersection optimization. Finally, tower image data with precise geographic positioning information is output.

[0059] The paper "Geometric Positioning Method for Pan-Tilt-Zone Controlled Camera Assisted by Multi-Source Spatial Data" describes the method of this implementation.

[0060] Step 7: Match the image to be located by the tower camera with the image in the benchmark dataset, and achieve high-precision geolocation based on the known positioning parameters.

[0061] Since the positioning parameters of the images in the benchmark dataset are known, the images to be located can be positioned with high precision based on the matching results.

Claims

1. A high-precision positioning method for tower camera images using a space-based collaborative approach, characterized in that, Includes the following steps: Step 1: Acquire a high-resolution satellite orthophoto of the area surrounding the tower as a reference image. The tower camera rotates 360° from its initial downward viewing angle to capture an initial tower camera image. Combining the digital elevation model of the area with the initial downward viewing angle of the tower camera, a virtual oblique projection is performed on the satellite reference image. The region of interest of the satellite image is mapped to the oblique coordinate system of the tower camera to generate a virtual satellite image that is consistent with the tower camera at the initial downward viewing angle. Step 2: Perform multi-scale phase consistency analysis and occlusion perception feature fusion on the satellite reference image and the tower camera image to generate a simplified edge map of the satellite reference image and the tower camera image, thereby realizing terrain enhancement feature extraction; Step 3: Using a hierarchical matching framework, perform cross-view matching between satellite reference images and tower camera images to obtain matching point pairs; Step 4: Control the tower camera to automatically adjust the downward viewing angle according to the preset sequence, and perform 360° rotation shooting under each downward viewing angle to collect multi-angle image sequences covering the all-round scene around the tower, which together constitute the original images of the benchmark dataset; Step 5: Based on the benchmark dataset, construct a pole image area network adjustment model with fixed optical center constraints, and calculate the relative orientation parameters of adjacent images in the pole camera benchmark dataset; The relative orientation parameters of adjacent images specifically refer to the three-axis rotation angles of each image relative to a unified camera coordinate system under the constraint of a fixed optical center. Step 6: Using the matching point pairs obtained in Step 3, combined with the benchmark dataset in Step 4 and the relative orientation parameters in Step 5, meter-level control points are obtained through spatial forward intersection and backward intersection calculations, and the absolute positioning parameters of the images in the benchmark dataset are calculated. The absolute positioning parameters refer to the position and attitude parameters required to transform the tower camera images from their own coordinate system to the globally unified geographic coordinate system. Step 7: Match the image to be positioned by the tower camera with the images in the reference dataset, and achieve high-precision positioning of the image to be positioned based on the known positioning parameters.

2. The method according to claim 1, characterized in that, The specific steps of the virtual strabismus projection described in step 1 include: A three-dimensional point set is constructed in the geographic coordinate system of satellite imagery, and the elevation value corresponding to each pixel is obtained by interpolation through a digital elevation model; Based on the camera's external parameters, points in the geographic coordinate system are transformed to the camera coordinate system using the following formula: In the formula, The elevation transformation matrix generated for the DEM is used to map geographic elevations to local terrain surfaces. Its calculation relies on bilinear interpolation or kriging interpolation algorithms. Based on camera tilt angle and azimuth Constructed rotation matrix; offset vector Used to align the optical center of the camera with the origin of the local coordinate system in the satellite image; This represents the three-dimensional coordinates of a point in the WGS84 geodetic coordinate system. By combining the camera's internal parameters, three-dimensional points are projected onto the image plane to generate a virtual satellite image that matches the actual camera's oblique angle.

3. The method according to claim 1, characterized in that, The multi-scale phase fusion terrain enhancement feature extraction method described in step 2 includes: Step 2.

1. Perform multi-level wavelet decomposition on satellite imagery and tower camera imagery to decompose the input imagery into sub-band images of different scales in order to capture terrain features of different sizes; A series of Log-Gabor filters with different orientation angles are created for sub-band images of various scales. For each scale of sub-band image, each Log-Gabor filter in the corresponding scale filter bank is convolved. The complex convolution kernel captures the amplitude and phase of the same pixel point in both the direction consistent with the filter direction and the direction perpendicular to the filter direction. Anisotropic structure is extracted, and a feature vector is generated for each pixel point. The feature vector is used to describe the local texture orientation distribution of the pixel point, representing the amplitude and phase of each pixel point in each direction and scale. In sub-band images of the same scale, the feature vectors of pixels in each direction need to be multiplied by the direction weights. We obtain the processed feature vectors in each direction, where the weights are... Calculated using the weighting function: in, For the first Angles in all directions For direction-weighted range, Used to traverse all directions It is a weighted function for tilt angle. For the camera angle; Step 2.

2. Based on the magnitude and phase of the pixel in each direction at the same scale obtained in Step 2.1, construct a second-order moment tensor for each pixel. Perform feature decomposition on the second-order moment tensor of the pixel at each scale to obtain two eigenvalues ​​(λ1, λ2) and the corresponding eigenvectors (v1, v2), obtain the main direction and shape of the pixel, and identify the main response direction. Step 2.

3. Based on the phase information obtained by convolution of the Log-Gabor filter bank in various directions and at various scales, and combined with the local principal direction, construct the phase consistency matrix of each pixel at different scales; Phase consistency matrices obtained at different scales are stacked according to their corresponding pixel positions to form a three-dimensional array. Principal component analysis is then used to reduce the dimensionality of the three-dimensional array and extract the first principal component as a simplified map of the main edge of the image.

4. The method according to claim 3, characterized in that, The Log-Gabor filter bank settings include: 5 wavelet decomposition layers, 22.5° directional spacing, and 30° directional weighting range, indicating that directions within ±30° of the camera direction have higher weights.

5. The method according to claim 1, characterized in that, The hierarchical matching framework described in step 3 includes: Initial screening stage: Based on the key points and main direction information in the edge map, the rBRIEF descriptor is dynamically corrected in direction. The distribution of sampling points is adjusted by phase direction rotation. An enhanced rBRIEF descriptor is generated by adopting a non-uniform sampling strategy. Point pairs with weak similarity are excluded by bidirectional nearest neighbor and Hamming distance ratio. The spatial distribution uniformity of the point set is evaluated by using the Poisson distribution model. False responses in texture-mixed areas are removed to form a preliminary matching set. Fine matching stage: Taking each pair of matching points (a,b) from the preliminary matching point set output from the initial screening stage as the center, block images Xa and Xb are extracted respectively, centered on the coordinate positions in the simplified main edge maps of satellite and tower camera images. Block images Xa and Xb are composed of multiple concentric rings, multiple azimuth sectors, or multiple concentric rings and multiple azimuth sectors. Logarithmic polar coordinate transformation is performed to generate block images LXa and LXb respectively. The average gray value of each concentric ring, or each azimuth sector, or each concentric ring and each azimuth sector in block images LXa and LXb is calculated. The average gray values ​​of each concentric ring, or each directional sector, or each concentric ring and each directional sector are arranged in order to form a one-dimensional vector. The average gray values ​​of each concentric ring, or each directional sector, or each concentric ring and each directional sector in block image LXb are arranged in order to form another one-dimensional vector. The Log-Polar feature descriptors of block images LXa and LXb are obtained respectively. The Log-Polar feature descriptors have rotation and scale invariance. Among them, the block images Xa and Xb corresponding to the main edge simplified map of satellite and tower images are divided into the same concentric rings and directional sectors in the same way. Calculate the Pearson correlation coefficient of the Log-Polar feature descriptor corresponding to each pair of matching points (a,b) in the preliminary matching point set. Based on the DEM digital elevation model data or the three-dimensional information of the image itself, evaluate the topographic relief of the area where the matching points are located. For areas with flat terrain and areas with severe topographic relief, different similarity thresholds are used to tolerate feature differences caused by geometric distortion. Output a set of matching point pairs with similarity greater than the similarity threshold. Finally, the set of matching point pairs is used as input, and geometric consistency verification is performed in combination with terrain structure features to output a set of finely matched points with high positioning accuracy.

6. The method according to claim 5, characterized in that, The parameter settings for the fine matching stage include: Hamming distance ratio threshold set to 0.7, Log-Polar feature radius set to 32 pixels, number of rings set to 5, number of sectors set to 12, and vector dimension set to 60.

7. The method according to claim 1, characterized in that, Step 5, constructing the fixed optical center constraint regional network adjustment model, includes: Assuming that the optical center of all images is the same fixed point, the projection relationship between images is simplified to a pure rotation relationship estimation; The first step is to establish an observation equation based on the consistency of direction vectors by taking the interior orientation elements of all images in the benchmark dataset and the selected pairs of corresponding points in adjacent images as input. By transforming the unit direction vector of each pair of corresponding points in the camera coordinate system of each image to a unified coordinate system, the included angle residual is calculated and the initial direction residual vector is output. The second step involves using the initial orientation residual vector as input to construct an adjustment objective function containing three core elements: (1) an orientation consistency residual term, which uses a weighted Huber loss function for robust control; (2) a rotation parameter regularization term, which constrains the deviation between the rotation matrix of each image and the initial pose; and (3) a relative rotation smoothing term, which maintains the continuity of rotation changes between adjacent images. The output is a constrained nonlinear least squares problem statement. The third step involves using the objective function and initial rotation parameters as inputs, and iteratively solving the problem using the Levenberg-Marquardt algorithm. In each iteration, the following are updated synchronously: (1) the three-axis rotation parameters of each image; (2) the adaptive weights of the Huber loss function; (3) the balance coefficients between the regularization term and the residual term; and the optimized precise rotation matrix of each image is output. The fourth step involves using the optimized rotation matrix as input to perform posterior accuracy assessment, calculating the mean square error of the orientation residuals for each pair of corresponding points, eliminating gross observations with residuals greater than 3 times the mean square error, and outputting the final relative orientation parameters after quality control.

8. The method according to claim 1, characterized in that, The calculation of the absolute positioning parameters in step 6 includes: The matching point pairs are associated with the corresponding 3D point clouds in the benchmark dataset, and the high-precision positioning information of satellite reference imagery is used to provide a geographic reference control benchmark for pole imagery. Based on the geometric constraints established by the relative orientation parameters calculated in step 5, the entire image network is optimized as a whole. Meter-level control points with geographic coordinates are obtained through spatial forward intersection calculation, and then the absolute position and attitude of each image in a unified coordinate system are solved through back intersection optimization.

9. A computer-readable storage medium having a computer program stored thereon, which, when executed, implements the method as described in any one of claims 1 to 8.

Images