Elevation inversion method for complex terrain area based on small unmanned aerial vehicle borne interferometric SAR
By setting up multi-track flight paths to acquire images in a small UAV-borne interferometric SAR system, performing image offset calculation and registration, and using the Goldstein algorithm to filter the phase map, elevation inversion in complex terrain areas was achieved, solving the problems of long baseline phase ambiguity and insufficient unwinding accuracy, and generating a high-precision digital elevation model.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- CHONGQING INNOVATION CENTER OF BEIJING INSTITUTE OF TECHNOLOGY
- Filing Date
- 2023-09-04
- Publication Date
- 2026-06-23
AI Technical Summary
Existing technologies for elevation inversion in complex terrain areas using small UAV-borne interferometric SAR suffer from long baselines leading to coherence failure and phase ambiguity, insufficient phase unwrapping accuracy, and difficulty in handling strong discontinuities and steep slopes in complex terrain.
By setting N orbits for small UAVs to acquire N SAR images, calculating image offsets and performing image registration, using the Goldstein algorithm to filter the interferometric phase map, and combining multi-baseline phase unwrapping, the failure of the phase continuity assumption is overcome, and elevation inversion is achieved.
It ensures the generation of high-precision digital elevation models in complex terrain areas, solves the phase ambiguity problem caused by surface discontinuities and steep slopes, and improves the accuracy of phase unwinding.
Smart Images
Figure CN117372482B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of image processing technology, and specifically to a method for elevation inversion in complex terrain areas based on small unmanned aerial vehicle-borne interferometric SAR. Background Technology
[0002] When using small UAV-borne interferometric SAR for terrain elevation inversion, according to the interferometric SAR terrain elevation inversion formula, the longer the baseline, the higher the inversion accuracy. However, long baselines can cause problems such as decoherence between small UAV-borne SAR image pairs, leading to interferometric measurement failure. Furthermore, existing interferometric techniques are mainly based on a single baseline; during phase unwrapping, the interferometric phase wraps to -π to π, requiring phase unwrapping processing.
[0003] Existing single-baseline interferometric SAR phase unwrapping methods mostly use Itoh as a constraint, assuming that the absolute phase difference between adjacent pixels is less than π, meaning that the terrain surface of the measurement area has no steep slopes or strong discontinuities. However, for complex mountainous areas, the terrain surface has many strong discontinuous steep slopes. When using small UAV-borne interferometric SAR technology to invert its elevation, the phase difference between adjacent pixels in the interferometric phase will be higher than π, resulting in interferometric phase ambiguity and leading to errors in the phase unwrapping results.
[0004] Therefore, the current challenges are how to overcome the failure of the interferometric phase continuity assumption caused by strong surface discontinuities or steep slopes, and how to solve the insufficient unwrapping accuracy of multi-baseline phase in small UAV-borne interferometric SAR. Summary of the Invention
[0005] To address the shortcomings of existing technologies, this invention provides a method for elevation inversion in complex terrain areas based on small UAV-borne interferometric SAR, which solves the technical problem of how to overcome the failure of the interferometric phase continuity assumption caused by strong surface discontinuities or steep slopes, and solves the technical problem of insufficient multi-baseline phase unwinding accuracy of small UAV-borne interferometric SAR.
[0006] This invention provides a method for elevation inversion in complex terrain areas based on small unmanned aerial vehicle-borne interferometric SAR, including:
[0007] S1. Set the flight path of N-track small UAVs, where the X-axis and Y-axis coordinates of the start and end points of each track of small UAVs remain unchanged, and the value of the Z-axis increases step by step along the vertical direction. The airborne SAR of the small UAVs acquires N-track echo data according to the predetermined flight path, and uses the BP algorithm to process the data into images to obtain N SAR images with N baselines.
[0008] S2. Obtain the offset of the SAR image of a small UAV with an adjacent short baseline;
[0009] S3. Calculate the offset of the SAR image of the small UAV based on the offset of the small UAV, and perform image registration on the long baseline SAR image.
[0010] S4. Process the SAR image after image registration to obtain the corresponding interferometric phase map, and use the Goldstein algorithm to filter the interferometric phase map;
[0011] S5. Based on the interferometric phase of the interferometric phase map, perform multi-baseline phase unwrapping to complete the elevation inversion of complex terrain.
[0012] Optionally, obtaining N SAR images from N baselines includes:
[0013] The SAR image is represented as S i , i = 1, 2, 3…N.
[0014] Optionally, acquiring the offset of the SAR image of a small UAV with an adjacent short baseline includes:
[0015] SAR images acquired by small UAVs with adjacent short baselines i As the primary graphic element, S n-i As the first image, the first main image S i High signal-to-noise ratio pixels are used as control points. Integer offset and fractional offset estimations are performed with the coordinates of these control points as the center. The sum of the integer and fractional offsets is then used as the offset of the control point, which is expressed as:
[0016] Wherein, the integer offset is based on the first main image S i Centered on the coordinates of the control points, with a matching window size of 5*5, the first image S n-i Centered on the corresponding coordinates, with a search window size of 12*12, the integer offset is estimated using the multiple correlation method.
[0017] The decimal offset is based on the first main image S. i Centered on the control point coordinates, with a matching window size of 7*7, the first image S n-i Centered on the corresponding coordinates, the search window size is 15*15, and 16 times upsampling is performed. The fractional offset is estimated using the multiple correlation method.
[0018] Optionally, the calculation of the offset of the SAR image based on the small UAV's SAR image with a long baseline includes:
[0019] SAR images acquired by small drones with long baselines iS1 serves as the second primary image, S n For n≥3 images, the second image is resampled according to the control point coordinates of the second main image, and the offset of adjacent short baseline control points is also considered. Iterative accumulation corresponds to obtaining the second main image S1 and the second secondary image S. n The offset is expressed as:
[0020]
[0021] Optionally, the image registration of the long-baseline SAR image includes:
[0022] Using the control point coordinates and corresponding point coordinates of the second main image S1, in the second sub-image S n The offset is calculated using the least squares method to obtain the second image S. n The offset was fitted, and the second image S was completed using bilinear interpolation. n Resampling processing.
[0023] Optionally, the SAR image processing after image registration to obtain a corresponding interferometric phase map, and the filtering processing of the interferometric phase map using the Goldstein algorithm, includes:
[0024] The second main image S1 and the second sub-image S n Conjugate multiplication yields the corresponding interferometric phase map, which is then filtered using the Goldstein algorithm. Furthermore, the trajectory information from a small UAV is used to remove the ground-level phase from the interferometric phase. Therefore, the interferometric phase expressions for the target point in the monitoring area under different baselines are:
[0025]
[0026] Where, φ n Indicates S1 and S n Interference phase, λ represents the corresponding vertical baseline length, h represents the ground elevation of the target point P, λ represents the wavelength, r represents the shortest slant distance from the scene center to the antenna phase center of the track, and θ represents the antenna incident angle from the scene center to the track at the shortest slant distance.
[0027] Optionally, the multi-baseline phase unwrapping based on the interferometric phase map to complete the elevation inversion of complex terrain includes:
[0028] Transform Equation 2, making the baseline as The target elevation at that time can be expressed as:
[0029]
[0030] Based on the premise that the elevation h of the target is not affected by the baseline, Equation 3 can be expressed as follows in the case of multiple baselines:
[0031]
[0032] Formula 4 can be transformed into:
[0033]
[0034] According to the winding phase φ n The absolute phase can be expressed as Formula 5 can then be reformulated as:
[0035]
[0036] in, k represents the winding phase at different baseline lengths. n This represents the blur number of the same pixel in each wrapped phase image, if the length ratio of any two baselines in N-1 interferograms satisfies:
[0037]
[0038] Among them, Γ i (i = 1, 2, ..., N-1) are pairwise coprime positive integers, and according to the Chinese Remainder Theorem, the baseline solution is expressed as:
[0039]
[0040] Where C0 is any rational number that is not zero;
[0041] The combined formula 6-8 can be expressed as:
[0042]
[0043] Formula 9 can be transformed as follows:
[0044]
[0045] make
[0046]
[0047] Where · represents the integer division operation;
[0048] Substituting Formula 11 into Formula 10, Formula 10 can be expressed as:
[0049] 2π(ξ1+ζ1+Γ1k1)=2π(ξ2+ζ2+Γ2k2)=…=2π(ξ N-1 +ζ N-1 +ΓN-1 k N-1 (12)
[0050] Taking the remainder of each term in Equation 12 with respect to 2π, the system of congruence equations in Equation 12 can be expressed as:
[0051]
[0052] The value of x is obtained using the Chinese Remainder Theorem. If we assume the solution value is p, then the blur number of the same pixel in each wrapped phase image is represented as:
[0053]
[0054] Therefore, the absolute phase of the same pixel in a multi-baseline interferometric phase map can be expressed as:
[0055]
[0056] Substituting Formula 15 into Formula 3, the elevation h of any target point P within the region can be obtained. P , can be represented as:
[0057]
[0058] Compared with the prior art, the present invention has the following beneficial effects:
[0059] By using long-baseline SAR image registration by cyclically accumulating the offsets of airborne SAR images from adjacent short-baseline small UAVs, and by performing multi-baseline phase unwrapping through the interferometric phase of the interferometric phase map, the problem of the failure of the interferometric phase continuity assumption caused by strong surface discontinuities or steep slopes is overcome, ensuring the accuracy of the digital elevation model generated by InSAR. Attached Figure Description
[0060] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention.
[0061] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0062] Figure 1 This is a schematic diagram of the process of the present invention;
[0063] Figure 2 This is a schematic diagram of multi-baseline phase unwrapping of small unmanned aerial vehicle-borne interferometric SAR in this invention;
[0064] Figure 3 This is a schematic diagram of interpolation and resampling of control points for different baseline SAR images in this invention. Detailed Implementation
[0065] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application. Functional units with the same reference numerals in the examples of this invention have the same and similar structures and functions.
[0066] See Figure 1 This invention provides a method for elevation inversion in complex terrain areas based on small unmanned aerial vehicle-borne interferometric SAR, comprising:
[0067] S1. Set the flight path of N-track small UAVs, where the X-axis and Y-axis coordinates of the start and end points of each track of small UAVs remain unchanged, and the value of the Z-axis increases step by step along the vertical direction. The airborne SAR of the small UAVs acquires N-track echo data according to the predetermined flight path, and uses the BP algorithm to process the data into images to obtain N SAR images with N baselines.
[0068] S2. Obtain the offset of the SAR image of a small UAV with an adjacent short baseline;
[0069] S3. Calculate the offset of the SAR image of the small UAV based on the offset of the small UAV, and perform image registration on the long baseline SAR image.
[0070] S4. Process the SAR image after image registration to obtain the corresponding interferometric phase map, and use the Goldstein algorithm to filter the interferometric phase map;
[0071] S5. Based on the interferometric phase of the interferometric phase map, perform multi-baseline phase unwrapping to complete the elevation inversion of complex terrain.
[0072] In this embodiment, S1, N-track small UAV flight paths are set, wherein the X-axis and Y-axis coordinates of the start and end points of each track of small UAV remain unchanged, and the Z-axis value increases step by step along the vertical direction. The airborne SAR of the small UAV acquires N-track echo data according to the predetermined flight path, and the BP algorithm is used to process the data to obtain N SAR images.
[0073] S2. Obtain the offset of the SAR image of a small UAV with an adjacent short baseline.
[0074] See Figure 1 and Figure 2 SAR images of adjacent short-baseline small UAVs i As the first primary image, S n-i As the first image. In the first main image S i High signal-to-noise ratio pixels are selected as control points to complete the offset of the main and secondary SAR image pairs.
[0075] Offset estimation is divided into integer offset estimation and fractional offset estimation. For integer offset estimation, the control point coordinates of the first main image are used as the center, with a matching window size of 5*5. The corresponding coordinates in the first sub-image are used as the center, with a search window size of 12*12. The offset estimation is completed using the multiple correlation method.
[0076] Similarly, for decimal offset estimation, the first image is first offset by an integer offset. Then, using the control point coordinates as the center, the matching window size in the first main image is 7*7, and the search window size is 15*15. A 16x upsampling process is performed, and the offset is estimated using the multiple correlation method. The integer and decimal offsets of the control point are added together to obtain the final offset for that control point.
[0077] Similarly, for adjacent short-baseline small UAV-borne SAR images S i+1 and S i+2 The same method is used to complete the offset estimation. Repeat the above steps until the small UAV-borne SAR image S is complete. n-1 and S n Offset estimation. Assume the offset of the control point being filtered in each iteration of the loop is...
[0078] S3. Calculate the offset of the SAR image of the small UAV based on the offset of the small UAV, and perform image registration on the long baseline SAR image.
[0079] The primary task is to obtain their relative offsets. Due to the variation in baselines, the control points differ between SAR images of different baseline lengths, meaning that the offset of a long-baseline SAR image cannot be simply obtained by summing the offsets of short-baseline SAR images.
[0080] Therefore, see Figure 3 SAR images acquired by small drones with long baselines i S1 serves as the second primary image, S n For n≥3 images, the second image is resampled according to the control point coordinates of the second main image, and the offset of adjacent short baseline control points is also considered. Iterative accumulation corresponds to obtaining the second main image S1 and the second secondary image S.n The offset is expressed as:
[0081]
[0082] Using the control point coordinates and corresponding point coordinates of the second main image S1, the second sub-image S n The offset is calculated using the least squares method to obtain the second image S. n The offset was fitted, and the second image S was completed using bilinear interpolation. n Resampling processing.
[0083] S4. Process the SAR image after image registration to obtain the corresponding interferometric phase map, and use the Goldstein algorithm to filter the interferometric phase map.
[0084] The second main image S1 and the second sub-image S n Conjugate multiplication yields the corresponding interferometric phase map, which is then filtered using the Goldstein algorithm. Furthermore, the trajectory information from a small UAV is used to remove the ground-level phase from the interferometric phase. Therefore, the interferometric phase expressions for the target point in the monitoring area under different baselines are:
[0085]
[0086] Where, φ n Indicates S1 and S n Interference phase, λ represents the corresponding vertical baseline length, h represents the ground elevation of target point P, λ represents the wavelength, r represents the shortest slant distance from the scene center to the antenna phase center of track 1, and θ represents the antenna incident angle from the scene center to track 1 at the shortest slant distance.
[0087] S5. Based on the interferometric phase of the interferometric phase map, perform multi-baseline phase unwrapping to complete the elevation inversion of complex terrain.
[0088] Transform Equation 2, making the baseline as The target elevation at that time can be expressed as:
[0089]
[0090] Based on the premise that the elevation h of the target is not affected by the baseline, Equation 3 can be expressed as follows in the case of multiple baselines:
[0091]
[0092] Formula 4 can be transformed into:
[0093]
[0094] According to the winding phase φ n The absolute phase can be expressed as Formula 5 can then be reformulated as:
[0095]
[0096] in, k represents the winding phase at different baseline lengths. n This represents the blur number of the same pixel in each wrapped phase image, if the length ratio of any two baselines in N-1 interferograms satisfies:
[0097]
[0098] Among them, Γ i (i = 1, 2, ..., N-1) are pairwise coprime positive integers, and according to the Chinese Remainder Theorem, the baseline solution is expressed as:
[0099]
[0100] Where C0 is any rational number that is not zero;
[0101] The combined formula 6-8 can be expressed as:
[0102]
[0103] Formula 9 can be transformed as follows:
[0104]
[0105] make
[0106]
[0107] Where · represents the integer division operation;
[0108] Substituting Formula 11 into Formula 10, Formula 10 can be expressed as:
[0109] 2π(ξ1+ζ1+Γ1k1)=2π(ξ2+ζ2+Γ2k2)=…=2π(ξ N-1 +ζ N-1 +Γ N-1 k N-1 (12)
[0110] Taking the remainder of each term in Equation 12 with respect to 2π, the system of congruence equations in Equation 12 can be expressed as:
[0111]
[0112] The value of x is obtained using the Chinese Remainder Theorem. If we assume the solution value is p, then the blur number of the same pixel in each wrapped phase image is represented as:
[0113]
[0114] Therefore, the absolute phase of the same pixel in a multi-baseline interferometric phase map can be expressed as:
[0115]
[0116] Substituting Formula 15 into Formula 3, the elevation h of any target point P within the region can be obtained. P , can be represented as:
[0117]
[0118] This invention overcomes the problem of the failure of the interferometric phase continuity assumption caused by strong surface discontinuities or steep slopes by performing long baseline SAR image registration by cyclically accumulating the offsets of airborne SAR images from adjacent short baseline small UAVs, and by performing multi-baseline phase unwrapping through the interferometric phase of the interferometric phase map, thus ensuring the accuracy of the digital elevation model generated by InSAR.
[0119] It should be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.
[0120] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features claimed herein.
Claims
1. A method for elevation inversion in complex terrain areas based on small unmanned aerial vehicle-borne interferometric SAR, characterized in that, include: S1. Set N-track flight paths for small UAVs, where the X and Y coordinates of the start and end points of each UAV track remain unchanged, and the Z-axis value increases progressively along the vertical direction. The onboard SAR of the small UAVs acquires N-track echo data according to the predetermined path, and performs image processing using the BP algorithm to obtain N SAR images with N baselines. The SAR images are represented as follows: ; S2. Obtain the offset of the SAR image of a small UAV with an adjacent short baseline; S3. Calculate the offset of the SAR image of the long-baseline small UAV based on the offset of the SAR image of the small UAV, where the SAR image is acquired by the long-baseline small UAV. middle As the second main image, As a second image, the second image is resampled according to the control point coordinates of the second main image, and the offsets of adjacent short baseline control points are iteratively accumulated to obtain the corresponding second main image. Second image The offset is determined, and image registration is performed on long-baseline SAR images; S4. Process the SAR image after image registration to obtain the corresponding interferometric phase map, and use the Goldstein algorithm to filter the interferometric phase map. The interferometric phase expression for the target point in the monitoring area under different baselines is as follows: （2） in, express and Interference phase, This indicates the corresponding vertical baseline length. Indicates the target point Surface elevation, Indicates wavelength. This represents the shortest slant distance from the center of the scene to the center of the antenna phase of the flight path. This represents the antenna incident angle from the scene center to the track at the shortest slant range; S5. Based on the interferometric phase of the aforementioned interferometric phase map, perform multi-baseline phase unwrapping to complete the elevation inversion of complex terrain. Specifically, this includes: transforming Formula 2 to obtain the baseline... The target elevation at that time can be expressed as: （3） Based on the elevation of the target h Under the premise of not being affected by baselines, Equation 3 can be expressed as follows in the case of multiple baselines: （4） Formula 4 can be transformed into: （5） According to the winding phase The absolute phase can be expressed as Then formula 5 can be rewritten as: （6） in, This indicates the winding phase at different baseline lengths. This represents the blur number of the same pixel in each wrapped phase map. The length ratio of any two baselines in an interferogram satisfies: （7） in, Given pairwise coprime positive integers, and according to the Chinese Remainder Theorem, the baseline solution is expressed as: （8） in, Let be any rational number that is not zero; The combined formula 6-8 can be expressed as: （9） Formula 9 can be transformed as follows: （10） make （11） in, For integer operations; Substituting Formula 11 into Formula 10, Formula 10 can be expressed as: （12） Adjust the terms in formula 12 Taking the remainder, the system of congruence equations in formula 12 can be expressed as: （13） Calculated using the Chinese Remainder Theorem The value of , if we assume the solution value is Then, the blur number of the same pixel in each wrapped phase image is represented as: （14） Therefore, the absolute phase of the same pixel in a multi-baseline interferometric phase map can be expressed as: （15） Substituting Formula 15 into Formula 3 will yield any target point within the region. elevation , can be represented as: （16）。 2. The elevation inversion method for complex terrain areas based on small unmanned aerial vehicle-borne interferometric SAR as described in claim 1, characterized in that, The method of acquiring the offset of SAR images of small UAVs with adjacent short baselines includes: SAR images acquired by small UAVs with adjacent short baselines As the first main image As the first image, the first main image High signal-to-noise ratio pixels are used as control points. Integer offset and fractional offset estimations are performed with the coordinates of these control points as the center. The sum of the integer and fractional offsets is then used as the offset of the control point, which is expressed as: ; Wherein, the integer offset is based on the first main image. Centered on the control point coordinates, with a matching window size of 5*5, the first image... Centered on the corresponding coordinates, with a search window size of 12*12, the integer offset is estimated using the multiple correlation method. The decimal offset is based on the first main image. Centered on the control point coordinates, with a matching window size of 7*7, the first image... Centered on the corresponding coordinates, the search window size is 15*15, and 16 times upsampling is performed. The fractional offset is estimated using the multiple correlation method.
3. The elevation inversion method for complex terrain areas based on small unmanned aerial vehicle-borne interferometric SAR as described in claim 2, characterized in that, Second main image Second image The offset is expressed as: （1） in, This indicates the offset between adjacent short baseline control points.
4. The elevation inversion method for complex terrain areas based on small UAV-borne interferometric SAR as described in claim 3, characterized in that, The process of image registration for long-baseline SAR images includes: With the second main image The coordinates of the control points and the corresponding point coordinates in the second image The offset was used to complete the second image using the least squares method. The offset was fitted, and the second image was completed using bilinear interpolation. Resampling processing.
5. The elevation inversion method for complex terrain areas based on small unmanned aerial vehicle-borne interferometric SAR as described in claim 4, characterized in that, The SAR image processing after image registration obtains the corresponding interferometric phase map, and the Goldstein algorithm is used to filter the interferometric phase map, including: The second main image Second image Conjugate multiplication yields the corresponding interferometric phase map, which is then filtered using the Goldstein algorithm. Furthermore, the trajectory information from a small UAV is used to remove the ground-level phase from the interferometric phase. Therefore, the interferometric phase expressions for the target point in the monitoring area under different baselines are: （2） in, express and Interference phase, This indicates the corresponding vertical baseline length. Indicates the target point Surface elevation, Indicates wavelength. This represents the shortest slant distance from the center of the scene to the center of the antenna phase of the flight path. This represents the antenna incident angle from the scene center to the track at the shortest slant range.