A small unmanned aerial vehicle InSAR image partition registration method
Patent Information
- Application Number
- CN202410235823.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2024-03-01
- Publication Date
- 2026-09-25
- Estimated Expiration
- 2044-03-01
AI Technical Summary
对于复杂地形下的InSAR图像,其偏移量大且空变,采取全局多项式拟合的传统InSAR图像配准算法难以很好地解决
[0031]相比于现有技术,本发明的优点及有益效果在于:本发明能够通过构建沿高度向等间隔的多个成像平面,并采用BP算法对无人机载SAR获取的两轨回波数据在多个成像平面中进行成像处理,得到多组包括主图像和辅图像的主辅图像;计算多组主辅图像在相同成像平面下主辅图像间的相干系数,生成相干系数图;将成像区域均匀划分为等大的子块,获取每个成像平面对应的主辅图像之间子块的最大平均相干系数,将对应的成像平面高度作为子块区域类DEM的平均高度,并按照子块对应位置进行组合和插值处理,得到整个成像区域的类DEM;计算目标点在主辅图像中的偏移量,基于偏移量设定偏移量阈值,并结合各个像素点的位置坐标获取高程门限,基于高程门限对类DEM进行区域分割;在子块中选取多个参考点,采用滑窗计算干涉图质量评价指标的方式,估计多个参考点的二维偏移量,并基于多项式参数模型,估计得到子块中所有像素点的偏移量;对所有子块像素点的偏移量进行全局融合,得到全局偏移量,基于全局偏移量对辅图像进行图像重采样,完成主辅图像的配准;基于多成像平面成像结果,进行相干系数计算并构建类DEM,对成像区域进行地形分割,通过对子块偏移量的估计和全局偏移量的融合,最终实现了复杂地形下的高精度图像匹配,确保了后续干涉测量的性能。
Smart Images

Figure CN118314174B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of synthetic aperture radar technology, and more particularly to a method for partitioning and registering small unmanned aerial vehicles (UAVs) carrying InSAR images. Background Technology
[0002] Interferometric Synthetic Aperture Radar (InSAR), as a novel radar system, is widely used in remote sensing fields such as landslide monitoring and topographic mapping. With increasing demands for monitoring flexibility, miniaturized and lightweight unmanned aerial vehicles (UAVs) provide a more flexible and efficient observation platform for InSAR, and have gradually become a research hotspot.
[0003] InSAR image registration involves geometrically correcting two images acquired from different viewpoints / baselines to align them and obtain an accurate interferometric phase map. Small UAV-borne InSAR platforms are limited by flight altitude, and the terrain elevation within the measurement area cannot be ignored. For images with spatial baselines, the viewpoints acquired vary greatly, and the same target may shift significantly between two passing radar images, causing image mismatch and resulting in changes in image coherence. In complex terrain, the pixel offset between image pairs can vary considerably.
[0004] Traditional InSAR image registration methods estimate pixel offsets using a sliding window search with a measure function. Based on the extracted offsets, a polynomial fitting model is constructed using the pixel coordinates, and the model parameters are estimated using the least squares method to establish a mapping relationship from the main image to the auxiliary image. However, for InSAR images with complex terrain, the offsets are large and spatially variable, making it difficult for traditional InSAR image registration algorithms using global polynomial fitting to handle them effectively.
[0005] Therefore, there is an urgent need for a high-precision image registration method suitable for small UAV-borne InSAR to ensure the performance of subsequent interferometric measurements. Summary of the Invention
[0006] Therefore, it is necessary to provide a small UAV-borne InSAR image partitioning and registration method to address the above-mentioned technical problems, so as to achieve high-precision matching of InSAR images under complex terrain conditions.
[0007] A method for partitioning and registering InSAR images on a small unmanned aerial vehicle (UAV) includes the following steps: constructing multiple imaging planes at equal intervals along the height direction; based on the multiple imaging planes, using the BP algorithm to perform imaging processing on the two-track echo data acquired by the UAV-borne SAR to obtain multiple sets of primary and secondary images, wherein the primary and secondary images include a primary image and a secondary image; calculating the coherence coefficient between the primary and secondary images under the same imaging plane for the multiple sets of primary and secondary images, generating a coherence coefficient map; uniformly dividing the imaging region into equally sized sub-blocks, obtaining the maximum average coherence coefficient between the sub-blocks corresponding to the primary and secondary images of each imaging plane, using the height of the imaging plane corresponding to the maximum average coherence coefficient as the average height of the DEM-like structure of the sub-block region, and dividing the imaging region into equal-sized sub-blocks; and finally, dividing the imaging region into equal-sized sub-blocks. The average height is combined and interpolated according to the corresponding positions of the sub-blocks to obtain a DEM-like image of the entire imaging area. The offset of the target point in the main and auxiliary images is calculated, an offset threshold is set based on the offset, and an elevation threshold is obtained by combining the position coordinates of each pixel. The DEM-like image is then segmented based on the elevation threshold. Multiple reference points in the sub-blocks are selected, and the two-dimensional offsets of the multiple reference points are estimated by using a sliding window method to calculate the interferogram quality evaluation index. Based on a polynomial parameter model, the offsets of all pixels in the sub-blocks are estimated. The offsets of all sub-block pixels are globally fused to obtain the global offset. The auxiliary image is resampled based on the global offset to complete the registration of the main and auxiliary images.
[0008] In one embodiment, the construction of multiple equally spaced imaging planes along the altitude direction, and the use of the BP algorithm to image the two-track echo data acquired by the UAV-borne SAR based on these multiple imaging planes, yields multiple sets of primary and secondary images. This includes: setting the UAV's flight altitude to H, constructing equally spaced imaging planes along the altitude direction with an altitude interval of dh, and different altitude imaging planes being L. n n = 1, 2, ..., N, where N = H / dh represents the number of imaging planes constructed, and the imaging plane height h n Represented as:
[0009] h n =dh·(n-1) (1)
[0010] The constructed plane L n The two-track echo data acquired by the UAV-borne SAR are sequentially used as imaging planes, and the BP algorithm is used to sequentially image the two-track echo data on the imaging planes to obtain N sets of main and auxiliary images S. 1n S 2n , among which, S 1n S 2n This represents the BP imaging result of the constructed nth imaging plane.
[0011] In one embodiment, the step of calculating the coherence coefficient between the primary and secondary images under the same imaging plane and generating a coherence coefficient map for the multiple sets of primary and secondary images includes: for the primary and secondary images S 1n S 2n The coherence coefficient between the primary and secondary images under the same imaging plane is calculated using the following formula:
[0012]
[0013] In the formula, γ n The coherence coefficient between the primary and secondary images corresponding to the nth image plane is given. The conjugate representation of the auxiliary image is used, and E is the energy intensity of the calculated image. A coherence coefficient map is generated based on all the calculated coherence coefficients.
[0014] In one embodiment, the step of uniformly dividing the imaging region into equally sized sub-blocks, obtaining the maximum average coherence coefficient between the main and auxiliary images corresponding to each imaging plane, using the imaging plane height corresponding to the maximum average coherence coefficient as the average height of the sub-block region's DEM-like structure, and combining and interpolating the average heights of each sub-block according to their corresponding positions to obtain the DEM-like structure of the entire imaging region includes: calculating the mismatch and decoherence of target points within the region between the main and auxiliary images, using the following formula:
[0015]
[0016] In the formula, γ coreg For mismatch and incoherence, μ r For the offset pixel, ρ r For distance resolution, (x p ,y p ,z p Let ) represent the position coordinates of the target point P, and Δz p =|z p -h| represents the relative elevation between the target elevation and the imaging plane elevation, where h is the imaging plane elevation and B is the baseline length. The size of the imaging region is represented as X*Y and uniformly divided into k*l equal-sized sub-blocks. The average coherence coefficient between the main and auxiliary images corresponding to each imaging plane is obtained. The maximum average coherence coefficient is compared to obtain the maximum average coherence coefficient, and the imaging plane elevation corresponding to the maximum average coherence coefficient is taken as the average elevation of the sub-block DEM, expressed as:
[0017]
[0018] In the formula, The height of the imaging plane corresponding to the maximum average coherence coefficient of the i-th sub-block. For the i-th sub-block in the imaging plane h nThe corresponding average coherence coefficient; the average height of each sub-block is combined according to the corresponding position of the sub-block to obtain the DEM-like of the sub-block, and the DEM-like of the sub-block is interpolated to obtain the DEM-like of the entire imaging area, with the interpolation factor being (X / k, Y / l).
[0019] In one embodiment, the calculation of the target point's offset in the primary and secondary images, setting an offset threshold based on the offset, and obtaining an elevation threshold by combining the position coordinates of each pixel, and performing region segmentation on the DEM-like image based on the elevation threshold, includes: using the imaging plane height 0 as a reference, calculating the offset of the target point P in the primary and secondary images within the imaging region during two flyby interferometric measurements of the UAV-borne SAR, using the following formula:
[0020]
[0021] Based on the offset of the target point between the main and auxiliary images, an offset threshold is set for the DEM-like image; the corresponding elevation threshold is obtained by combining the position coordinates of each pixel point, and the different terrain elevation areas in the DEM-like image are segmented based on the elevation threshold, into elevation areas below the elevation threshold and elevation areas above the elevation threshold.
[0022] In one embodiment, the step of selecting multiple reference points in a sub-block, using a sliding window method to calculate the interferogram quality evaluation index, estimating the two-dimensional offset of the multiple reference points, and estimating the offset of all pixels in the sub-block based on a polynomial parameter model includes: selecting a high signal-to-noise ratio point in the sub-block as a reference point; setting a matching window centered on the reference point (i,j) in the main image and setting a search window at the same position in the auxiliary image; moving pixel by pixel in the search window in both row and column directions; obtaining the two-dimensional offset of the reference point when the interferogram quality evaluation index between the search window and the matching window is maximized; using the two-dimensional offset of the reference point as a benchmark, estimating the offset parameters of the sub-block using a polynomial parameter model and the least squares method; and combining the offset parameters with the position information of all pixels in the sub-block to estimate the offset of all pixels in the sub-block.
[0023] In one embodiment, the polynomial parameter model is:
[0024]
[0025] In the formula, △x i,j , △y i,j Let be the two-dimensional offset of pixel (i,j) in the sub-block, and let a0, a1, a2 and b0, b1, b2 be the parameters to be estimated.
[0026] In one embodiment, the evaluation metric is the coherence coefficient, the average fluctuation function, or the spectral function.
[0027] In one embodiment, the global fusion of the offsets of all sub-block pixels to obtain the global offset includes: globally fusing the estimated offsets of all sub-block pixels, wherein the offsets of non-edge region pixels in a sub-block remain unchanged, and for pixel A in the edge region of an adjacent sub-block, the offset is represented by a weighted linear combination, as follows:
[0028]
[0029] In the formula, f1(A) and f2(A) are the offsets of pixel A between two adjacent sub-blocks, l1 is the width of the edge region, and r1 and r2 are the distances from pixel A to the boundaries of the two sub-blocks.
[0030] In one embodiment, the resampling method is a linear interpolation or a bilinear interpolation method.
[0031] Compared with existing technologies, the advantages and beneficial effects of this invention are as follows: This invention can construct multiple imaging planes at equal intervals along the height direction, and use the BP algorithm to process the two-track echo data acquired by UAV-borne SAR in multiple imaging planes to obtain multiple sets of main and auxiliary images including a main image and an auxiliary image; calculate the coherence coefficient between the main and auxiliary images in the same imaging plane to generate a coherence coefficient map; uniformly divide the imaging area into equally sized sub-blocks, obtain the maximum average coherence coefficient between the sub-blocks corresponding to the main and auxiliary images of each imaging plane, use the height of the corresponding imaging plane as the average height of the DEM-like image of the sub-block area, and perform combination and interpolation processing according to the corresponding positions of the sub-blocks to obtain the DEM-like image of the entire imaging area; calculate the offset of the target point in the main and auxiliary images, and set the target point based on the offset. Offset thresholds are used, and elevation thresholds are obtained by combining the position coordinates of each pixel. Based on the elevation thresholds, the DEM-like region is segmented. Multiple reference points are selected in the sub-blocks, and the two-dimensional offsets of multiple reference points are estimated by using a sliding window method to calculate the interferogram quality evaluation index. Based on a polynomial parameter model, the offsets of all pixels in the sub-blocks are estimated. The offsets of all sub-block pixels are globally fused to obtain the global offset. Based on the global offset, the auxiliary image is resampled to complete the registration of the main and auxiliary images. Based on the imaging results of multiple imaging planes, coherence coefficients are calculated and a DEM-like region is constructed. The imaging area is segmented into terrain. By estimating the sub-block offsets and fusing the global offsets, high-precision image matching under complex terrain is finally achieved, ensuring the performance of subsequent interferometry. Attached Figure Description
[0032] Figure 1 This is a flowchart illustrating a method for partitioning and registering InSAR images on a small unmanned aerial vehicle (UAV) in one embodiment.
[0033] Figure 2 The images shown in one embodiment are the original image and the interferometric phase map and coherence coefficient map processed by MCC cross-correlation registration and the partition registration method of this application, wherein (a) is the original interferometric phase map, (b) is the original coherence coefficient map, (c) is the cross-correlation interferometric phase map, (d) is the cross-correlation correlation coefficient map, (e) is the partition registration interferometric phase map, and (f) is the partition registration coherence coefficient map.
[0034] Figure 3 This is a histogram of the original image and the coherence coefficients processed using MCC cross-correlation registration and the partition registration method of this application, as shown in one embodiment. Detailed Implementation
[0035] Before describing the specific embodiments of the present invention, the overall concept of the present invention will be explained as follows:
[0036] This invention is mainly based on the InSAR image registration process of small UAVs. Traditional InSAR image registration methods have large and variable offsets when registering InSAR images in complex terrain, and cannot achieve high-precision image registration in complex terrain.
[0037] Therefore, this invention proposes a small UAV-borne InSAR image partitioning and registration method. This method constructs multiple equally spaced imaging planes along the altitude direction and uses the BP algorithm to process the two-track echo data acquired by the UAV-borne SAR in these multiple imaging planes, obtaining multiple sets of master and slave images. The method calculates the coherence coefficients between the master and slave images on the same imaging plane, generating a coherence coefficient map. The imaging region is uniformly divided into equally sized sub-blocks. The maximum average coherence coefficient between the sub-blocks corresponding to the master and slave images on each imaging plane is obtained. The height of the corresponding imaging plane is used as the average height of the DEM-like structure of the sub-block region. Combination and interpolation processing are performed according to the corresponding positions of the sub-blocks to obtain the DEM-like structure of the entire imaging region. The offset of the target point in the master and slave images is calculated, and based on the offset... An offset threshold is set, and an elevation threshold is obtained by combining the position coordinates of each pixel. Based on the elevation threshold, the DEM-like region is segmented. Multiple reference points are selected in the sub-blocks, and the two-dimensional offsets of multiple reference points are estimated by using a sliding window method to calculate the interferogram quality evaluation index. Based on a polynomial parameter model, the offsets of all pixels in the sub-blocks are estimated. The offsets of all sub-block pixels are globally fused to obtain the global offset. Based on the global offset, the auxiliary image is resampled to complete the registration of the main and auxiliary images. Based on the imaging results of multiple imaging planes, the coherence coefficient is calculated and a DEM-like region is constructed. The imaging area is segmented into terrain. By estimating the sub-block offsets and fusing the global offsets, high-precision image matching under complex terrain is finally achieved, ensuring the performance of subsequent interferometry.
[0038] Having introduced the overall concept of the present invention, to make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below through specific embodiments in conjunction with the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative of the present invention and are not intended to limit the present invention.
[0039] In one embodiment, such as Figure 1 As shown, a method for partitioning and registering InSAR images on a small unmanned aerial vehicle (UAV) is provided, including the following steps:
[0040] Step S101: Construct multiple imaging planes at equal intervals along the height direction. Based on the multiple imaging planes, use the BP algorithm to process the two-track echo data acquired by the UAV-borne SAR to obtain multiple sets of main and auxiliary images, which include a main image and an auxiliary image.
[0041] Specifically, two-track echo data are acquired by UAV-borne SAR, and multiple imaging planes with equal intervals along the altitude direction are constructed. The obtained two-track echo data are processed by the BP (Error Back Propagation) algorithm through multiple imaging planes to obtain multiple sets of main and auxiliary images. The main and auxiliary images include a one-to-one corresponding main image and auxiliary image.
[0042] Step S101 includes: setting the UAV's flight altitude to H, constructing equally spaced imaging planes along the altitude direction with an altitude interval of dh, and different altitude imaging planes being L. n n = 1, 2, ..., N, where N = H / dh represents the number of imaging planes constructed, and the imaging plane height h n Represented as:
[0043] h n =dh·(n-1) (1)
[0044] The constructed plane L n The two-track echo data acquired by the UAV-borne SAR are sequentially used as imaging planes, and the BP algorithm is used to sequentially image the two-track echo data on the imaging planes to obtain N sets of main and auxiliary images S. 1n S 2n , among which, S 1n S 2n This represents the BP imaging result of the constructed nth imaging plane.
[0045] Specifically, assuming the UAV flies at an altitude of H, multiple imaging planes with equal intervals along the altitude direction are constructed, with an altitude interval of dh, resulting in N imaging planes at different altitudes. Based on the obtained multiple imaging planes, the BP algorithm is used sequentially to process the two-track echo data acquired by the UAV-borne SAR on the imaging planes, resulting in N sets of primary and secondary images. This facilitates high-precision image matching under loaded terrain based on the imaging results of multiple imaging planes.
[0046] Step S102: For multiple sets of primary and secondary images, calculate the coherence coefficient between the primary and secondary images under the same imaging plane, and generate a coherence coefficient map.
[0047] Specifically, based on the obtained multiple sets of primary and secondary images, the coherence coefficient between the primary and secondary images under each imaging plane is calculated, and a coherence coefficient map is generated according to all the obtained coherence coefficients, so as to make an intuitive judgment on the coherence of the primary and secondary images.
[0048] Step S102 includes: for the main and auxiliary images S 1n S 2n The coherence coefficient between the primary and secondary images under the same imaging plane is calculated using the following formula:
[0049]
[0050] In the formula, γ n The coherence coefficient between the primary and secondary images corresponding to the nth image plane is given. The conjugate representation of the auxiliary image is used, and E is the energy intensity of the calculated image. A coherence coefficient map is generated based on all the calculated coherence coefficients.
[0051] Specifically, based on the obtained primary and secondary images, the coherence coefficient between the primary and secondary images under each imaging plane is calculated using formula (2), and a coherence coefficient map can be generated based on all the obtained coherence coefficients, so as to intuitively judge the coherence in the image.
[0052] Step S103: Divide the imaging area into equally sized sub-blocks, obtain the maximum average coherence coefficient between the sub-blocks of the main and auxiliary images corresponding to the imaging plane, take the height of the imaging plane corresponding to the maximum average coherence coefficient as the average height of the sub-block region DEM, combine and interpolate the average height of each sub-block according to the corresponding position of the sub-block to obtain the DEM of the entire imaging area.
[0053] Specifically, in the imaging plane, the imaging area is uniformly divided into multiple sub-blocks of the same size. Combined with the coherence coefficient calculated in step S102, the average coherence coefficient between the sub-blocks of the main and auxiliary images corresponding to each imaging plane is obtained. The maximum average coherence coefficient is obtained by comparison. The height of the imaging plane corresponding to the maximum average coherence coefficient is taken as the average height of the sub-block region's DEM (Digital Elevation Model). The average heights of all the sub-blocks are combined according to the corresponding positions of the sub-blocks to obtain the sub-block's DEM. The sub-block's DEM is interpolated to obtain the DEM of the entire imaging area. By constructing the DEM, the digital simulation of the ground terrain based on the terrain elevation data is realized, which facilitates subsequent high-precision image registration.
[0054] Step S103 includes: calculating the mismatch and decoherence of target points within the region between the primary and secondary images, using the following formula:
[0055]
[0056] In the formula, γ coreg For mismatch and incoherence, μ r For the offset pixel, ρ r For distance resolution, (x p y p , z p Let ) represent the position coordinates of the target point P, and Δz p =|z p -h| indicates the relative elevation between the target elevation and the imaging surface height, where h is the imaging surface height and B is the baseline length;
[0057] The imaging region is represented as X*Y and uniformly divided into k*l equal-sized sub-blocks. The average coherence coefficient between the main and auxiliary images corresponding to each imaging plane is obtained. The maximum average coherence coefficient is obtained by comparison, and the height of the imaging plane corresponding to the maximum average coherence coefficient is taken as the average height of the sub-block region DEM, expressed as:
[0058]
[0059] In the formula, The height of the imaging plane corresponding to the maximum average coherence coefficient of the i-th sub-block. For the i-th sub-block in the imaging plane h n The corresponding average coherence coefficient; the average height of each sub-block is combined according to the corresponding position of the sub-block to obtain the DEM-like of the sub-block, and the DEM-like of the sub-block is interpolated to obtain the DEM-like of the entire imaging area, with the interpolation factor being (X / k, Y / l).
[0060] Specifically, for the coherence coefficient of the target point in the region generated on different imaging planes, the factors affecting the change of the coherence coefficient are mainly the decoherence caused by the mismatch between the main and auxiliary images. The mismatch decoherence is caused by the positional shift of the target point between the main and auxiliary images.
[0061] For a target point P within the region, as the imaging height h gradually approaches the target elevation z... p The relative elevation Δz between the target and the imaging plane p Decrease, γ coreg As the value increases, the target coherence coefficient gradually increases. When the coherence coefficient generated by a certain imaging plane is the largest, the imaging plane height is optimal.
[0062] The imaging region is uniformly divided into sub-blocks of the same size, and the average coherence coefficient between the primary and secondary images corresponding to each imaging plane is obtained. The height of the imaging plane corresponding to the maximum average coherence coefficient of the sub-block is taken as the optimal imaging plane height. The optimal imaging surface height is used as the average height of the sub-block region's DEM. The average heights of each sub-block are combined according to their corresponding positions to obtain the sub-block's DEM. The sub-block's DEM is then interpolated using an interpolation factor to obtain the DEM of the entire imaging region. The interpolation factor is calculated based on the imaging region size and the number of sub-blocks.
[0063] Step S104: Calculate the offset of the target point in the main and auxiliary images, set the offset threshold based on the offset, and obtain the elevation threshold by combining the position coordinates of each pixel point. Perform region segmentation on the DEM-like image based on the elevation threshold.
[0064] Specifically, the offset of the target point in the main and auxiliary images is calculated, an offset threshold is set based on the calculated offset, and the corresponding elevation threshold is obtained by combining the position coordinates of each pixel. Based on the elevation threshold, different terrain elevation areas in the DEM-like image are segmented to facilitate the regional processing of complex terrain and improve the accuracy of image matching.
[0065] Step S104 includes: using the imaging plane height 0 as a reference, calculating the offset of target point P in the main and auxiliary images within the imaging area during two flyby interferometric measurements of the UAV-borne SAR, using the formula:
[0066]
[0067] Based on the offset of the target point between the main and auxiliary images, an offset threshold is set for the DEM-like image; the corresponding elevation threshold is obtained by combining the position coordinates of each pixel point, and the different terrain elevation areas in the DEM-like image are segmented based on the elevation threshold, into elevation areas below the elevation threshold and elevation areas above the elevation threshold.
[0068] Specifically, taking the imaging plane height of 0 as a reference, during two interferometric measurements by the UAV, the target point P in the imaging area has a large offset in the main and auxiliary images. The corresponding offset is calculated, and the offset threshold of the DEM-like image is set based on this offset. The corresponding elevation threshold is obtained by combining the position coordinates of each pixel. The DEM-like image is segmented according to the elevation threshold into different elevation regions below and above the elevation threshold, so as to facilitate regional processing. The offset threshold can be set to 0.5 pixel units.
[0069] Step S105: Select multiple reference points in the sub-block, use the sliding window method to calculate the interferogram quality evaluation index, estimate the two-dimensional offset of multiple reference points, and estimate the offset of all pixels in the sub-block based on the polynomial parameter model.
[0070] Specifically, after segmenting different elevation regions, multiple reference points are selected in the sub-blocks. The two-dimensional offset of multiple reference points is estimated by using a sliding window method to calculate the interferogram quality evaluation index. At the same time, a polynomial parameter model is used to estimate the offset of all pixels in the sub-block based on the two-dimensional offset of the reference points, thereby realizing the offset estimation of the sub-block region to facilitate subsequent partitioning and registration operations.
[0071] Step S105 includes: selecting a high signal-to-noise ratio point in the sub-block as a reference point; setting a matching window centered on the reference point (i,j) in the main image and setting a search window at the same position in the auxiliary image; moving pixel by pixel in the search window along the row and column directions respectively; obtaining the two-dimensional offset of the reference point when the interferogram quality evaluation index between the search window and the matching window is maximized; using the two-dimensional offset of the reference point as a benchmark, estimating the offset parameters of the sub-block using a polynomial parameter model and the least squares method; and combining the offset parameters and the position information of all pixels in the sub-block to estimate the offset of all pixels, thereby obtaining the offset of all pixels in the sub-block.
[0072] The polynomial parametric model is as follows:
[0073]
[0074] In the formula, △x i,j , △y i,j Let be the two-dimensional offset of pixel (i,j) in the sub-block, and let a0, a1, a2 and b0, b1, b2 be the parameters to be estimated.
[0075] Specifically, after segmenting different elevation regions, points with relatively high signal-to-noise ratios in the sub-regions are selected as reference points. The corresponding signal-to-noise ratio threshold can be set according to the required number of reference points to filter out high signal-to-noise ratio points as reference points. The two-dimensional offset of the reference points is estimated by using a sliding window to calculate the interferogram quality evaluation index. The evaluation index can be the coherence coefficient, the average fluctuation function, or the spectral function, etc.
[0076] A matching window is set in the main image with the reference point as the center, and a search window is set in the same position in the auxiliary image. In the search window, the pixels are moved one by one in the row and column directions respectively. When the interferogram quality evaluation index of the two windows is maximized, the two-dimensional offset of the reference point is obtained.
[0077] Using the two-dimensional offset of the reference point as a benchmark, the offset parameters of the sub-block are estimated by using a polynomial parameter model and the least squares method. The offset parameters are then combined with the position information of all pixels in the sub-block to complete the offset estimation of all pixels and obtain the offset of all pixels in the sub-block.
[0078] Step S106: Globally fuse the offsets of all sub-block pixels to obtain the global offset.
[0079] Specifically, the offsets of all sub-block pixels in the image are globally fused to obtain global offsets. These global offsets are then fused to facilitate high-precision image matching in complex terrains.
[0080] Step S106 includes: globally fusing the estimated offsets of all sub-block pixels, wherein the offsets of non-edge region pixels in a sub-block remain unchanged, and the offsets of pixel A in the edge region of an adjacent sub-block are represented by a weighted linear combination.
[0081]
[0082] In the formula, f1(A) and f2(A) are the offsets of pixel A between two adjacent sub-blocks, l1 is the width of the edge region, and r1 and r2 are the distances from pixel A to the boundaries of the two sub-blocks.
[0083] Specifically, the estimated offsets of all sub-block pixels in the image are globally fused. In order to make the two-dimensional registration offsets of pixels in the image continuous, the offsets of non-edge region pixels in the sub-block remain unchanged, and the pixel A in the edge region of adjacent sub-blocks is processed by linear weighted combination, thereby realizing the global fusion of offsets and enabling high-precision partition registration under complex terrain based on the global offsets.
[0084] Step S107: Resample the auxiliary image based on the global offset to complete the registration of the main and auxiliary images.
[0085] Specifically, based on the global offset obtained after fusion, the auxiliary image is resampled to complete the registration of the main and auxiliary images, thus achieving high-precision image matching under complex terrain.
[0086] The resampling method can be either linear interpolation or bilinear interpolation.
[0087] In this embodiment, multiple imaging planes are constructed at equal intervals along the height direction, and the backpropagation (BP) algorithm is used to process the two-track echo data acquired by the UAV-borne SAR in multiple imaging planes to obtain multiple sets of main and auxiliary images, including a main image and an auxiliary image. The coherence coefficient between the main and auxiliary images in the same imaging plane is calculated to generate a coherence coefficient map. The imaging area is uniformly divided into equally sized sub-blocks, and the maximum average coherence coefficient between the sub-blocks corresponding to the main and auxiliary images of each imaging plane is obtained. The height of the corresponding imaging plane is used as the average height of the DEM-like structure of the sub-block area, and the sub-blocks are combined and interpolated according to their corresponding positions to obtain the DEM-like structure of the entire imaging area. The offset of the target point in the main and auxiliary images is calculated, an offset threshold is set based on the offset, and combined with... The elevation threshold is obtained by determining the position coordinates of each pixel, and the DEM-like region is segmented based on the elevation threshold. Multiple reference points are selected in the sub-block, and the two-dimensional offset of multiple reference points is estimated by calculating the interferogram quality evaluation index using a sliding window method. Based on a polynomial parameter model, the offset of all pixels in the sub-block is estimated. The offsets of all sub-block pixels are globally fused to obtain the global offset. Based on the global offset, the auxiliary image is resampled to complete the registration of the main and auxiliary images. Based on the imaging results of multiple imaging planes, the coherence coefficient is calculated and a DEM-like region is constructed. The imaging area is segmented into terrain. By estimating the sub-block offsets and fusing the global offsets, high-precision image matching under complex terrain is finally achieved, ensuring the performance of subsequent interferometry.
[0088] In one embodiment, such as Figure 2 As shown, these are the interferometric phase diagrams and coherence coefficient diagrams after processing using different methods. The unprocessed original image is shown below. Figure 2 As shown in (a) and (b), the images processed using the existing MCC cross-correlation registration method are as follows: Figure 2 As shown in (c) and (d) of this paper, and the image processed by the partition registration method proposed in this invention, as shown in the figure. Figure 2 As shown in (e) and (f).
[0089] like Figure 3As shown, the original image has the largest fluctuation in coherence coefficient, so the calculated average coherence coefficient is the largest; the MCC cross-correlation registration method has smaller fluctuations compared to the original image, so the average coherence coefficient is smaller compared to the original image; while the partition registration method proposed in this invention has the smallest fluctuation in coherence coefficient, so the obtained average coherence coefficient will be the smallest.
[0090] The evaluation results of the original image and the image processed by the MCC cross-correlation method and the partition registration method of the present invention are shown in Table 1:
[0091] Table 1. Evaluation Results of Interferograms
[0092] Original image 0.28 459045 MCC cross-correlation 0.35 437530 Partition registration 0.39 415784
[0093] As shown in Table 1, the original image has the highest average coherence coefficient and the lowest residual points when unregistered. Conversely, the partitioned registration method proposed in this application yields the highest average coherence coefficient and the lowest residual points. Since a higher coherence coefficient results in clearer interference fringes, and a smaller residual point indicates a higher degree of model fit optimization, it can be concluded that the partitioned registration method proposed in this application has better registration results than the unprocessed original image and the existing MCC cross-correlation registration method, thus ensuring the performance of subsequent interferometric measurements.
[0094] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. The storage medium can be a magnetic disk, optical disk, read-only memory (ROM), or random access memory (RAM), etc.
[0095] It will be apparent to those skilled in the art that the steps of the present invention described above can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using device-executable program code, thereby storing them in a computer storage medium (ROM / RAM, magnetic disk, optical disk) for execution by the computing device. Furthermore, in some cases, the steps shown or described can be performed in a different order than those presented herein, or they can be fabricated as separate integrated circuit modules, or implemented as a single integrated circuit module. Therefore, the present invention is not limited to any particular hardware and software combination.
[0096] The above description, in conjunction with specific embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such deductions or substitutions should be considered within the scope of protection of the present invention.
Claims
1. A method for partitioning and registering InSAR images on a small unmanned aerial vehicle (UAV), characterized in that, Includes the following steps: Multiple imaging planes are constructed at equal intervals along the height direction. Based on the multiple imaging planes, the BP algorithm is used to perform imaging processing on the two-track echo data acquired by the UAV-borne SAR to obtain multiple sets of main and auxiliary images, which include a main image and an auxiliary image. For the multiple sets of primary and secondary images, calculate the coherence coefficient between the primary and secondary images under the same imaging plane, and generate a coherence coefficient map; The imaging area is uniformly divided into equal-sized sub-blocks. The maximum average coherence coefficient between the main and auxiliary images corresponding to each imaging plane is obtained. The height of the imaging plane corresponding to the maximum average coherence coefficient is used as the average height of the sub-block region-like DEM. The average height of each sub-block is combined and interpolated according to the corresponding position of the sub-block to obtain the DEM of the entire imaging area. Calculate the offset of the target point in the main and auxiliary images, set an offset threshold based on the offset, and obtain an elevation threshold by combining the position coordinates of each pixel. Perform region segmentation on the DEM-like image based on the elevation threshold. Multiple reference points are selected within a sub-block. A sliding window method is used to calculate the interferogram quality evaluation index to estimate the two-dimensional offset of these reference points. Based on a polynomial parameter model, the offset of all pixels within the sub-block is estimated, including: selecting high signal-to-noise ratio points within the sub-block as reference points, and using these reference points in the main image (…). i,j A matching window is set as the center, and a search window is set at the same position in the auxiliary image. Pixels are moved one by one in the search window in both row and column directions. When the interferogram quality evaluation index between the search window and the matching window is maximized, the two-dimensional offset of the reference point is obtained. Based on the two-dimensional offset of the reference point, the offset parameters of the sub-block are estimated using a polynomial parameter model and the least squares method. Combined with the offset parameters and the position information of all pixels in the sub-block, the offset of all pixels is estimated to obtain the offset of all pixels in the sub-block. The offsets of all sub-block pixels are globally fused to obtain the global offset. Based on the global offset, the auxiliary image is resampled to complete the registration of the main and auxiliary images.
2. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle (UAV) according to claim 1, characterized in that, The process involves constructing multiple imaging planes at equal intervals along the altitude direction. Based on these imaging planes, a backpropagation (BP) algorithm is used to process the two-track echo data acquired by the UAV-borne SAR, resulting in multiple sets of primary and secondary images, including: Set the drone's flight altitude to H Construct imaging planes that are equally spaced along the height direction, with a height interval of [missing information]. dh Imaging planes at different heights are L n , n=1,2,L,N, where, N=H / dn Indicates the number of imaging planes constructed, and the height of the imaging planes. h n Represented as: ;(1) The constructed plane L n The two-track echo data acquired by the UAV-borne SAR are sequentially used as imaging planes, and the BP algorithm is applied to these imaging planes to sequentially perform imaging processing on the two-track echo data. N Group of primary and secondary images S 1n , S 2n ,in, S 1n , S 2n Indicates the first construction n BP imaging results for each imaging plane.
3. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle (UAV) according to claim 2, characterized in that, The step of calculating the coherence coefficient between the primary and secondary images under the same imaging plane and generating a coherence coefficient map for the multiple sets of primary and secondary images includes: For primary and secondary images S 1n , S 2n The coherence coefficient between the primary and secondary images under the same imaging plane is calculated using the following formula: ;(2) In the formula, For the first n The coherence coefficient between the primary and secondary images corresponding to the image plane. The conjugate representation of the auxiliary image. E To calculate the energy intensity of the image; A coherence coefficient map is generated based on all the calculated coherence coefficients.
4. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle according to claim 2, characterized in that, The process involves uniformly dividing the imaging region into equally sized sub-blocks, obtaining the maximum average coherence coefficient between the primary and secondary images corresponding to each imaging plane, using the imaging plane height corresponding to the maximum average coherence coefficient as the average height of the sub-block region's DEM-like structure, and combining and interpolating the average heights of each sub-block according to their corresponding positions to obtain the DEM-like structure of the entire imaging region. This includes: The formula for calculating the mismatch and decoherence of target points within the region between the primary and secondary images is as follows: ;(3) In the formula, This is due to mismatch and incoherence. For offset pixels, For distance resolution, For target point P Location coordinates, This indicates the relative elevation of the target to the image plane. h The height of the imaging plane, B Baseline length; The size of the imaging region is represented as and evenly divided into For each equally sized sub-block, the average coherence coefficient between the primary and secondary images corresponding to each imaging plane is obtained. The maximum average coherence coefficient is compared to obtain the average height of the imaging plane corresponding to the maximum average coherence coefficient, and the height of the imaging plane corresponding to the maximum average coherence coefficient is taken as the average height of the sub-block DEM, expressed as: ;(4) In the formula, For the first i The imaging plane height corresponding to the maximum average coherence coefficient of each sub-block For the first i Each sub-block on the imaging plane h n The corresponding average coherence coefficient; The average heights of the obtained sub-blocks are combined according to their corresponding positions to obtain the DEM-like structure of each sub-block. The DEM-like structure of the sub-blocks is then interpolated to obtain the DEM-like structure of the entire imaging area, with an interpolation factor of [value missing]. .
5. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle according to claim 4, characterized in that, The calculation involves determining the offset of the target point in the primary and secondary images, setting an offset threshold based on the offset, and obtaining an elevation threshold by combining the position coordinates of each pixel. Based on the elevation threshold, region segmentation is performed on the DEM-like image, including: Using the imaging plane height 0 as a reference, the offset of target point P in the main and auxiliary images within the imaging area is calculated during two flyby interferometric measurements of the UAV-borne SAR. The formula is as follows: ;(5) Based on the offset of the target point between the main and auxiliary images, an offset threshold is set for the DEM-like image. The elevation threshold is obtained by combining the position coordinates of each pixel. Based on the elevation threshold, different terrain elevation regions in the DEM are segmented into elevation regions below the elevation threshold and elevation regions above the elevation threshold.
6. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle according to claim 1, characterized in that, The polynomial parameter model is: ;(6) In the formula, For pixels in a sub-block The two-dimensional offset, a 0 , a 1 , a 2 and b 0 , b 1 , b 2 These are the parameters to be estimated.
7. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle according to claim 1, characterized in that, The evaluation index is the coherence coefficient, average fluctuation function, or spectral function.
8. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle according to claim 2, characterized in that, The global offset of all sub-block pixels is obtained by global fusion, including: The estimated offsets of all sub-block pixels are globally fused. The offsets of pixels in non-edge regions within a sub-block remain unchanged. For pixels A in the edge regions of adjacent sub-blocks, the offsets are represented by a weighted linear combination: ;(7) In the formula, f 1 (A) , f 2 (A) These represent the offsets of pixel A between its two adjacent sub-blocks. The width of the edge region. r 1 , r 2 This is the distance from pixel A to the boundary between the two sub-blocks.
9. The method for partitioning and registering InSAR images on a small unmanned aerial vehicle according to claim 1, characterized in that, The resampling method is either linear interpolation or bilinear interpolation.
Citation Information
Patent Citations
Fast correlation coefficient method for interferometric synthetic aperture radar image precise registration
CN102955157A
A DEM-aided SAR image registration method with high accuracy
CN109035312A