SAR (Synthetic Aperture Radar) image terrain radiation correction method considering overlay and shadow regions

Abstract

The invention relates to the technical field of SAR (Synthetic Aperture Radar) image processing, and particularly discloses an SAR image terrain radiation correction method considering overlay and shadow regions, which comprises the following steps: acquiring an SAR image, a DEM (Digital Elevation Model) image and a source file in a research region; carrying out multi-view processing on the original SAR image, resampling the DEM image, and carrying out radiometric calibration; geometric positioning is carried out according to the RD model by adopting an indirect positioning method, lookup tables A and B are generated, and images are registered; detecting pixels of the overlay area and the shadow area; calculating an effective scattering area corresponding to each SAR pixel; and performing radiation correction according to the effective scattering area, repairing a back scattering coefficient value of the shadow area, and finally performing geometric correction. According to the method, the topographic effect caused by an SAR imaging mechanism and topographic relief is eliminated to a certain extent, the radiation correction effect of an area with large topographic relief is improved, and the corrected product can be widely applied to the fields of ground feature inversion and classification, deformation detection, target detection and the like.

Description

Technical Field

[0001] This invention relates to the field of SAR image processing technology, specifically to a SAR image topographic radiometric correction method that takes into account overlapping and shadowed areas. Background Technology

[0002] Synthetic Aperture Radar (SAR) imagery, due to its all-weather, all-time, high-resolution, and high-precision characteristics, has been widely used in remote sensing fields, including geological mapping, resource exploration, disaster monitoring, agricultural forecasting, and forestry management. With continuous technological advancements and the expanding application areas, the performance indicators of SAR technology have significantly improved, leading to its increasingly widespread and in-depth application in various fields. However, the complex imaging geometry of SAR imagery, especially the problem of topographic radiation distortion, has become a significant factor limiting its application. The imaging principle of SAR imagery dictates that its representation in complex terrain exhibits significant geometric distortions, such as perspective contraction, shadows, and overlays. These distortions negatively impact the accuracy of quantitative analysis and classification. In mountainous or hilly areas, the changes in echo signals caused by terrain undulations result in a deviation between the radiometric values ​​of SAR imagery and the true surface reflectivity, thus affecting image interpretation and application. Furthermore, the polarization characteristics of SAR imagery are also affected by terrain, further exacerbating the radiation distortion problem. To improve the application effect of SAR imagery in areas with complex terrain, geometric and radiometric corrections must be performed on the imagery. Geometric correction is typically achieved using a digital elevation model (DEM) to eliminate geometric distortions caused by terrain; while radiometric correction adjusts the radiometric intensity of an image to make it closer to the actual surface reflection.

[0003] Topographic radiometric correction is a crucial method for addressing geometric and radiometric distortion issues in SAR images. SAR image topographic radiometric correction consists of two parts: geometric correction and radiometric correction. Geometric correction refers to the process of adjusting and correcting geometric distortions in synthetic aperture radar (SAR) images to match them with the actual geographic coordinate system. The purpose of geometric correction is to improve the spatial accuracy of the image, making the positions of ground targets in the image more precise. Radiometric correction aims to eliminate radiometric distortions or aberrations caused by external factors (such as atmospheric conditions, sensor response, and topographic effects), ensuring that images acquired at different times and angles have the same radiometric standard. The goal of radiometric correction is to improve image quality and interpretability, enabling reliable quantitative analysis of SAR images acquired at different times or by different sensors.

[0004] Currently, commonly used SAR image geometric correction models include the Range-Doppler Model (RD model), the Polynomial Correction Model, the Rational Function Model (RFM), and the Affine Transformation Model. Among these, the Polynomial Correction Model, the Rational Function Model, and the Affine Transformation Model are empirical or semi-empirical models and are not specifically designed for SAR image geometric correction. The RD model, however, utilizes the imaging geometric parameters of the SAR system (such as orbit, attitude, and slant range) to establish a mathematical relationship between image coordinates and geographic coordinates; it is a physical model used only for SAR image geocoding or geometric correction. Commonly used SAR image topographic radiometric correction models include the Local Incident Angle Model, the Projection Angle Model, and the Area Integral Model. Ultimately, all three models calculate the effective scattering area of ​​the image and then perform radiometric correction.

[0005] To address the above issues, some scholars have proposed local incident angle correction models and projection angle correction models based on the RD model. Both models utilize angles to perform radiometric correction on images based on geometric correction. While both methods offer some correction effectiveness, this is only applicable to plains or relatively flat areas. In areas with significant terrain undulations, the correction results of both methods are less than satisfactory. The reason for this is that these angle-based radiometric correction methods do not consider the overlapping and shadowing issues in SAR images, only addressing the radiometric distortion caused by pixel scaling.

[0006] To address the shortcomings of the aforementioned radiometric correction algorithms, another radiometric correction method emerged—the area integration method. This method fully considers the topological relationship between DEM and SAR images. For both one-to-many and many-to-one relationships between DEM and SAR images, the area integration method is used to calculate the effective scattering area, significantly improving the radiometric correction effect. Furthermore, other studies have categorized the factors causing radiometric distortion in SAR images into two types: area effect and local incident angle effect. After removing the area effect, a step to eliminate the local incident angle effect is added. This method greatly improves the radiometric correction effect of SAR images in mountainous and plateau regions with significant topographic relief, and it remains widely used for topographic radiometric correction of SAR images today. However, this method still has some limitations:

[0007] 1) Although the one-to-many and many-to-one relationships between DEM and SAR images were considered, the impact of shadowed pixels on the overlapping area was not taken into account, so the radiometric correction of pixels in the overlapping area still has problems.

[0008] 2) Although radiometric correction was performed on the image, no special processing was done on the shadow pixels in the image, i.e., shadow pixel restoration.

[0009] Therefore, a method for topographic radiometric correction of SAR images that takes into account overlapping and shadowed areas has become an urgent problem to be solved. Summary of the Invention

[0010] The purpose of this invention is to solve the problems of unsatisfactory correction effect of existing SAR image topographic radiometric correction methods in areas with large topographic relief and lack of repair of shadowed pixels. This invention provides a SAR image topographic radiometric correction method that takes into account overlapping and shadowed areas. This method takes into account the topological relationship between DEM image and SAR image, performs radiometric correction on SAR image by calculating the real and effective surface scattering area corresponding to each SAR pixel, and repairs shadowed areas.

[0011] To achieve the above objectives, the technical solution provided by this invention is: a SAR image topographic radiometric correction method that takes into account overlapping and shadowed areas, comprising the following steps:

[0012] S1. Data Acquisition: Acquire SAR and DEM images within the study area, and obtain the source files of the SAR images;

[0013] S2. Image preprocessing: Perform multi-view processing on the original SAR image according to the resolution of the SAR image and DEM image, and resample the DEM image to obtain SAR image and DEM image with relatively consistent resolution. At the same time, perform radiometric calibration on the SAR image according to the radiometric calibration file or calibration constant.

[0014] S3. Geometric positioning: After preprocessing, the DEM image and SAR image are geometrically positioned according to the RD model using the indirect positioning method. The correspondence between the DEM image and the calculated simulated SAR image, as well as the correspondence between the simulated SAR image and the real SAR image, are obtained. Lookup table A and lookup table B are generated. At the same time, the real SAR image is registered onto the simulated SAR image according to lookup table B.

[0015] S4. Overlapping and Shadow Area Pixel Detection: Obtain overlapping area pixels on the DEM image based on the generated lookup table A, and then detect shadow area pixels on the DEM image based on the shadow area imaging principle of SAR image.

[0016] S5. Calculate the effective scattering area: Calculate the surface area of ​​each pixel in the DEM image, project it onto the isophase surface using the projection angle, and calculate its area; then calculate the effective scattering area corresponding to each SAR pixel based on the overlapping area pixels, the shadow area pixels, and lookup table A.

[0017] S6. Geometric and radiometric correction: Perform radiometric correction on the SAR image based on the calculated effective scattering area, and repair the backscattering coefficient values ​​of the SAR image in the shadowed area based on the detected overlapping pixels and shadowed pixels; after the radiometric correction is completed, perform geometric correction on the SAR image according to lookup table A.

[0018] Furthermore, in step S2, the specific method for image preprocessing is as follows:

[0019] S21. The DEM resampling rate is obtained by resampling the DEM image:

[0020] Where θ refers to the angle of incidence, and δ dem This refers to the resolution of the DEM image, δ out This refers to the slant range resolution of the output analog image, and f refers to the minimum resampling rate of the DEM.

[0021] S22. To balance the azimuth resolution of the DEM image, multi-look processing is performed on the SAR image to make the pixels of the SAR image as square as possible. The number of looks is calculated as follows:

[0022] Where, δ r For slant range resolution, δ a This refers to the azimuth resolution.

[0023] S23. Based on the radiometric calibration factor given in the SAR image source file, perform absolute radiometric calibration on the SAR image, converting the original DN value into radar brightness β. 0 .

[0024] Furthermore, in step S3, the geometric positioning method is as follows:

[0025] S31. Based on the positioning parameters in the SAR image source file, use the RD model to generate a simulated SAR image using the indirect positioning method, and store the correspondence between DEM image pixels and simulated SAR image pixels as a lookup table A.

[0026] During the generation of SAR simulated images, the projection angle of each pixel in the DEM image is calculated, and a projection angle file is generated. The formula for calculating the projection angle is:

[0027] Where η is the projection angle, Let be the normal vector of a pixel in the DEM image. This is the normal vector of the imaging surface corresponding to the pixel;

[0028] S32. Based on the generated simulated SAR image, use the Harris corner detection algorithm to establish the pixel correspondence between the simulated SAR image and the real SAR image, and generate lookup table B. At the same time, register the real SAR image onto the simulated SAR image.

[0029] Furthermore, in step S4, the specific methods for detecting overlapping and shadow region pixels are as follows:

[0030] S41. Overlapping area detection: When multiple DEM image pixels correspond to one SAR image pixel, it is called overlapping. Therefore, according to the mapping relationship in lookup table A, we search for the case in lookup table A where one SAR image pixel corresponds to one DEM image pixel, directly detect the overlapping area pixels of the DEM image, and generate an overlapping area mask file.

[0031] S42. Shadow area detection: Perform inverse geometric correction on the DEM image based on the coarse track data to convert it into an image with coordinates approximately the same as the uncorrected SAR image, and then detect shadow areas based on pixel occlusion.

[0032] Furthermore, the method for detecting the shadow area is as follows:

[0033] S421. Perform inverse geometric correction. Based on the satellite's orbit equation, calculate the angle between the satellite and the due north direction. Use this angle to convert the DEM image into an image with the same direction as the SAR image.

[0034] S422. Detect the image line by line from left to right, set a 1×n detection window with a step size of 1 and n as the base number; take the current pixel as the detection center, if the elevation value of the pixel is greater than the elevation values ​​of the pixels on the left and right sides, then the pixel is the mountain peak point, otherwise it is not the mountain peak point; detect the pixels of the entire row in sequence to obtain the position of the mountain peak point in each row of the image.

[0035] S423. Calculate the incident angle θ corresponding to each pixel and save it as an incident angle file. The formula for calculating the incident angle is:

[0036]

[0037] in, This is the satellite's position vector in the WGS84 or CGCS2000 coordinate system at this moment, with the direction pointing from the satellite towards the coordinate center. The distance vector between the satellite and the point at this moment, with the direction from the satellite to the ground point;

[0038] S424. Detect the image line by line. When the first mountain peak is detected, check the incident angle θ corresponding to that point in the incident angle file. Detect the next point to the right and check the incident angle θ1. If θ1 is greater than θ, it is not a shadow point. Continue to the right to find the next mountain peak. If θ1 is less than or equal to θ, the point is a shadow point. Mark it as a shadow point and continue to detect the next point until a non-shadow point is detected. Then the shadow area occluded by that mountain peak is detected, and the detection of the next shadow area begins.

[0039] S425. After all images have been inspected, a shadow area mask file is generated.

[0040] Furthermore, the method for inverse geometric correction is as follows:

[0041] S4211. Calculate the velocity vector of the satellite at the reference point during its transit time, based on the satellite orbit formula:

[0042]

[0043] S4212. Convert the satellite's spatial rectangular coordinates at the reference point to geodetic coordinates;

[0044] S4213. Construct the basis vectors of the northeast celestial coordinate system:

[0045] (East direction);

[0046] (North direction);

[0047] (Heavenly direction);

[0048] Where φ is the latitude of the satellite at the reference point, and λ is the longitude of the satellite at the reference point;

[0049] S4214. Calculate the satellite's velocity vector in the east direction ΔE and the north direction ΔN in the northeast celestial coordinate system:

[0050]

[0051] S4215. Calculate the angle between the satellite's velocity direction and true north, i.e., the azimuth angle:

[0052]

[0053] S4216. Using the center point of the DEM image as the rotation center, rotate it according to the azimuth angle to make it an image with coordinates approximately the same as the uncorrected SAR image.

[0054] S4217. Adjust the rotated DEM image again according to the satellite's flight direction; when the satellite flies from south to north, flip the DEM image vertically; when the satellite flies from north to south, flip the DEM image horizontally; this completes the inverse geometric correction of the DEM image.

[0055] Furthermore, in step S5, the specific method for calculating the effective scattering area is as follows:

[0056] S51. Use the triangle method to calculate the land area of ​​DEM image pixels, that is, divide each DEM image pixel into two triangles, and calculate the sum of the areas of the two triangles, which is the land area of ​​the DEM image pixel.

[0057] S52. Based on the projection angle file, calculate the projected area of ​​each DEM image pixel on the equiphase surface, and generate a projected area file for the equiphase surface. The calculation formula is: γ dem (i,j)=σ dem (i,j)·sinη(i,j);

[0058] Where (i,j) is the row and column number corresponding to the DEM image pixel, σ dem η is the surface area of ​​the DEM image cell, and η is the projection angle of the DEM image cell; γ dem This represents the projected area of ​​the DEM image pixel on the equiphase surface;

[0059] S53. The obtained overlay mask file and shadow mask file are superimposed to obtain shadow pixels in the overlay area. Although these pixels are located in the overlay area, they do not theoretically contribute to the backscattering coefficient. The shadow mask file is subtracted from these shadow pixels in the overlay area to generate a new shadow area file.

[0060] S54. Based on the shadow pixels in the overlay area obtained in step S3, change the area of ​​these shadow pixels in the projection area file of the equiphase surface to 0.

[0061] Furthermore, in step S6, the specific methods for geometric and radiometric corrections are as follows:

[0062] S61. Combining lookup table A and the projected area file of the equiphase surface, calculate the effective scattering area corresponding to each SAR image pixel. The calculation formula is as follows:

[0063] Where (r,a) represents the row and column number in the SAR image. γ represents the effective scattering area corresponding to the pixel at that location in the SAR image, (i,j) represents the row and column number of the pixel in the DEM, and γ dem It is the effective projected area of ​​the pixel point of the DEM image on the equiphase surface;

[0064] S62. Based on the backscattering coefficient expression of the SAR image imaging surface, perform topographic radiometric correction on the polarimetric SAR image, as shown in the following formula:

[0065] Where, β 0 The original radar brightness values ​​for the SAR image. γ is the scattering area corresponding to the imaging surface, i.e., the pixel size of the SAR image. 0 (r,a) represents the corrected SAR image pixel brightness value, i.e., the backscattering coefficient value;

[0066] S63. Based on the new shadow area file detected by DEM and lookup table A, obtain the corresponding shadow area pixels in the SAR image; calculate the average backscattering coefficient value of the row where the shadow area pixel is located, and fill and repair the shadow area pixel.

[0067] S64. According to lookup table A, perform geometric correction on the topographically radiometrically corrected SAR image, project the SAR image onto geospatial coordinates, and obtain the geocoded SAR product after topographic radiometric correction.

[0068] The advantages of this invention compared to the prior art are:

[0069] This invention fully considers the impact of overlapping and shadowed regions on the radiometric characteristics of SAR images. By accurately detecting pixels in overlapping and shadowed regions and performing special processing on these regions when calculating the effective scattering area and performing radiometric correction, the accuracy of radiometric correction is improved.

[0070] In the image preprocessing process, this invention maintains consistent resolution between SAR and DEM images through reasonable resampling and multi-view processing, laying a solid data foundation for subsequent geometric positioning and radiometric correction.

[0071] This invention employs the RD model and indirect positioning method for geometric positioning, combined with the Harris corner detection algorithm for image registration, which improves the accuracy of geometric positioning and ensures the reliable correspondence between DEM images and SAR images.

[0072] The calculation of the effective scattering area in this invention takes into account the projection angle of the terrain and the effects of overlapping and shadowed areas, making the calculation results more consistent with the actual situation and providing accurate parameters for radiometric correction.

[0073] This invention fills in and repairs the backscattering coefficient values ​​in shadowed areas, improving the integrity and usability of SAR images and providing more reliable data support for subsequent image interpretation and applications.

[0074] This invention eliminates, to a certain extent, the terrain effect caused by the SAR imaging mechanism and topographic relief, significantly improving the radiometric correction effect of SAR images in areas with large terrain undulations. Products corrected for terrain radiometry and geometry can be widely used in fields such as ground feature inversion and classification, deformation detection, and target detection. Attached Figure Description

[0075] Figure 1 This is a flowchart of a SAR image topographic radiometric correction method that takes into account overlapping and shadowed areas according to the present invention.

[0076] Figure 2 This is a flowchart of inverse geometric correction.

[0077] Figure 3 This is a flowchart of the shadow area detection method. Detailed Implementation

[0078] Various exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be noted that, unless otherwise specifically stated, the relative arrangement, numerical expressions, and values ​​of the components and steps set forth in these embodiments do not limit the scope of the present invention.

[0079] The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the invention or its application or use.

[0080] Techniques, methods, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, such techniques, methods, and equipment should be considered part of the specification.

[0081] In all examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0082] The following detailed description, in conjunction with the accompanying drawings, provides a method for SAR image topographic radiometric correction that takes into account overlapping and shadowed areas.

[0083] Combined with appendix Figure 1-3 The specific implementation process of the SAR image topographic radiometric correction method that takes into account overlapping and shadowed areas of the present invention is as follows:

[0084] A method for topographic radiometric correction of SAR images that takes into account overlapping and shadowed areas includes the following steps:

[0085] S1. Obtain Data

[0086] Acquire SAR and DEM images of the study area, and obtain the source files of the SAR images. The SAR image source files contain important information such as the SAR image positioning parameters and radiometric calibration factors, providing a data foundation for subsequent preprocessing, geometric positioning, and radiometric correction steps; the DEM images are used to reflect the topographic relief of the study area.

[0087] S2, Image Preprocessing

[0088] The original SAR image undergoes multi-look processing based on its resolution, and the DEM image is resampled to obtain SAR and DEM images with relatively consistent resolution. Simultaneously, the SAR image is radiometrically calibrated according to a radiometric calibration file or calibration constant. The specific methods are as follows:

[0089] S21. The DEM resampling rate is obtained by resampling the DEM image:

[0090] Where θ refers to the angle of incidence, and δ dem This refers to the resolution of the DEM image, δ out This refers to the slant range resolution of the output analog image, and f refers to the minimum resampling rate of the DEM. By calculating this resampling rate, it can be ensured that the resampled DEM image matches the resolution of subsequent processing, thus providing a guarantee for accurate geometric positioning and radiometric correction.

[0091] S22. To balance the azimuth resolution of the DEM image, multi-look processing is performed on the SAR image to make the pixels of the SAR image as square as possible. The number of looks is calculated as follows:

[0092] Where, δ r For slant range resolution, δ a For azimuth resolution; multi-view processing can reduce speckle noise in SAR images and adjust pixel shape, which is beneficial for subsequent image registration and processing.

[0093] S23. Based on the radiometric calibration factor given in the SAR image source file, perform absolute radiometric calibration on the SAR image, converting the original DN value into radar brightness β. 0 Radiometric calibration can eliminate the influence of the sensor itself, giving the gray values ​​of SAR images physical meaning, and facilitating comparison and analysis between images acquired at different times and by different sensors.

[0094] S3, Geometric Positioning

[0095] After preprocessing, DEM and SAR images are used to perform geometric positioning of the SAR image based on the RD model using an indirect positioning method. This obtains the correspondence between the DEM image and the calculated simulated SAR image, as well as the correspondence between the simulated SAR image and the real SAR image, generating lookup tables A and B. Simultaneously, the real SAR image is registered onto the simulated SAR image according to lookup table B. The specific method is as follows:

[0096] S31. Based on the positioning parameters in the SAR image source file, use the RD model to generate a simulated SAR image using the indirect positioning method, and store the correspondence between DEM image pixels and simulated SAR image pixels as a lookup table A.

[0097] During the generation of SAR simulated images, the projection angle of each pixel in the DEM image is calculated, and a projection angle file is generated. The formula for calculating the projection angle is:

[0098] Where η is the projection angle, Let be the normal vector of a pixel in the DEM image. The normal vector of the imaging surface corresponding to the pixel; the projection angle reflects the angular relationship between the terrain surface and the imaging surface, and is an important parameter for calculating the effective scattering area.

[0099] S32. Based on the generated simulated SAR image, the Harris corner detection algorithm is used to establish the pixel correspondence between the simulated SAR image and the real SAR image, and a lookup table B is generated. Simultaneously, the real SAR image is registered onto the simulated SAR image. The Harris corner detection algorithm has good stability and accuracy, and can effectively extract corner features from the image, providing reliable control points for image registration.

[0100] S4, Overlapping and Shadow Region Pixel Detection

[0101] The overlapping pixels in the DEM image are obtained from the generated lookup table A. Then, based on the imaging principle of shadow areas in SAR images, the shadow area pixels in the DEM image are detected. The specific method is as follows:

[0102] S41. Overlapping Area Detection: When multiple DEM image pixels correspond to one SAR image pixel, it is called overlapping. Therefore, based on the mapping relationship in lookup table A, we search for the case in lookup table A where one SAR image pixel corresponds to one DEM image pixel, directly detect the overlapping area pixels of the DEM image, and generate an overlapping area mask file. The existence of overlapping areas will cause information confusion on the SAR image. Accurate detection of overlapping areas is a prerequisite for effective correction.

[0103] S42. Shadow Area Detection: The DEM image is subjected to inverse geometric correction based on the coarse track data, transforming it into an image with coordinates approximately the same as the uncorrected SAR image. Then, shadow areas are detected based on pixel occlusion. The method for shadow area detection is as follows:

[0104] S421. Perform inverse geometric correction. Based on the satellite's orbit equations, calculate the angle between the satellite and true north. Use this angle to convert the DEM image into an image aligned with the SAR image. The method for inverse geometric correction is as follows:

[0105] S4211. Calculate the velocity vector of the satellite at the reference point during its transit time, based on the satellite orbit formula:

[0106]

[0107] S4212. Convert the satellite's spatial rectangular coordinates at the reference point to geodetic coordinates;

[0108] S4213. Construct the basis vectors of the northeast celestial coordinate system:

[0109] (East direction);

[0110] (North direction);

[0111] (Heavenly direction);

[0112] Where φ is the latitude of the satellite at the reference point, and λ is the longitude of the satellite at the reference point;

[0113] S4214. Calculate the satellite's velocity vector in the east direction ΔE and the north direction ΔN in the northeast celestial coordinate system:

[0114]

[0115] S4215. Calculate the angle between the satellite's velocity direction and true north, i.e., the azimuth angle:

[0116]

[0117] S4216. Using the center point of the DEM image as the rotation center, rotate it according to the azimuth angle to make it an image with coordinates approximately the same as the uncorrected SAR image.

[0118] S4217. Based on the satellite's flight direction, the rotated DEM image is adjusted again; when the satellite flies from south to north, the DEM image is flipped vertically; when the satellite flies from north to south, the DEM image is flipped horizontally; this completes the inverse geometric correction of the DEM image. Through inverse geometric correction, the DEM image and the SAR image are kept consistent in coordinate direction, creating conditions for subsequent shadow area detection.

[0119] S422. Detect the image row by row from left to right, setting a 1×n detection window with a step size of 1 and n as the base number; take the current pixel as the detection center, if the elevation value of the pixel is greater than the elevation values ​​of the pixels on the left and right sides, then the pixel is the mountain peak point, otherwise it is not the mountain peak point; detect all pixels in the entire row in sequence to obtain the position of the mountain peak point in each row of the image; the mountain peak point is a key position that may produce shadows, and accurately identifying the mountain peak point helps to improve the accuracy of shadow area detection.

[0120] S423. Calculate the incident angle θ corresponding to each pixel and save it as an incident angle file. The formula for calculating the incident angle is:

[0121]

[0122] in, This is the satellite's position vector in the WGS84 or CGCS2000 coordinate system at this moment. This is the distance vector between the satellite and the point at this moment; the incident angle is an important factor affecting the backscattering characteristics of SAR images and a key parameter for shadow area detection.

[0123] S424. Detect the image line by line, up to the first mountain peak. Check the incident angle θ corresponding to this point in the incident angle file. Detect the next point to the right and check the incident angle θ1. If θ1 is greater than θ, it is not a shadow point, and continue to the right to find the next mountain peak. If θ1 is less than or equal to θ, the point is a shadow point. Mark it as a shadow point and continue to detect the next point until a non-shadow point is detected. Then the shadow area occluded by that mountain peak is detected, and the detection of the next shadow area begins. This method based on incident angle comparison can effectively detect the range of the shadow area.

[0124] S425. After all images have been inspected, a shadow area mask file is generated.

[0125] S5. Calculate the effective scattering area

[0126] The surface area of ​​each pixel in the DEM image is calculated by projecting it onto an equiphase surface using the projection angle, and then calculating its area. Next, the effective scattering area corresponding to each SAR pixel is calculated based on the overlapping area pixels, the shadow area pixels, and lookup table A. The specific method is as follows:

[0127] S51. The surface area of ​​DEM image pixels is calculated using the triangle method, which involves dividing each DEM image pixel into two triangles and calculating the sum of the areas of the two triangles as the surface area of ​​the DEM image pixel. The triangle method can accurately calculate the area of ​​irregular terrain surfaces, providing a basis for subsequent calculation of effective scattering area.

[0128] S52. Based on the projection angle file, calculate the projected area of ​​each DEM image pixel on the equiphase surface, and generate a projected area file for the equiphase surface. The calculation formula is: γ dem (i,j)=σ dem (i,j)·sinη(i,j);

[0129] Where (i,j) is the row and column number corresponding to the DEM image pixel, σ dem η is the surface area of ​​the DEM image cell, and η is the projection angle of the DEM image cell; γ dem This represents the projected area of ​​the DEM image pixels on the equiphase surface; the projected area on the equiphase surface reflects the effective contribution area of ​​the terrain to radar wave scattering.

[0130] S53. The obtained overlay mask file and shadow mask file are superimposed to obtain shadow pixels in the overlay area. Although these pixels are located in the overlay area, they do not theoretically contribute to the backscattering coefficient. The shadow mask file is subtracted from these shadow pixels in the overlay area to generate a new shadow area file. Through this process, the region that contributes to the backscattering coefficient can be determined more accurately.

[0131] S54. Based on the shadow pixels in the overlay region obtained in step S3, change the area of ​​these shadow pixels in the projected area file of the equiphase surface to 0. This is because the shadow pixels in the overlay region do not participate in backscattering, and setting their area to 0 can improve the accuracy of the effective scattering area calculation.

[0132] S6, Geometric and Radiometric Corrections

[0133] Radiometric correction is performed on the SAR image based on the calculated effective scattering area, and the backscattering coefficient values ​​of the SAR image in the shadow area are repaired based on the detected overlapping pixels and shadow pixels. After radiometric correction is completed, geometric correction is performed on the SAR image according to lookup table A.

[0134] The specific method is as follows:

[0135] S61. Combining lookup table A and the projected area file of the equiphase surface, calculate the effective scattering area corresponding to each SAR image pixel. The calculation formula is as follows:

[0136] Where (r,a) represents the row and column number in the SAR image. γ represents the effective scattering area corresponding to the pixel at that location in the SAR image, (i,j) represents the row and column number of the pixel in the DEM, and γ dem It is the effective projected area of ​​the DEM image pixel on the equiphase surface; the effective scattering area is a key parameter for radiometric correction, and its accuracy directly affects the correction effect.

[0137] S62. Based on the backscattering coefficient expression of the SAR image imaging surface, perform topographic radiometric correction on the polarimetric SAR image, as shown in the following formula:

[0138] Where, β 0 The original radar brightness values ​​for the SAR image. γ is the scattering area corresponding to the imaging surface, i.e., the pixel size of the SAR image. 0 (r,a) represents the corrected SAR image pixel brightness value, i.e., the backscattering coefficient value; this formula takes into account the influence of terrain factors on the backscattering coefficient and can effectively eliminate radiation distortion caused by terrain undulation.

[0139] S63. Based on the new shadow area file detected by DEM and lookup table A, obtain the corresponding shadow area pixels in the SAR image; calculate the average backscattering coefficient value of the row where the shadow area pixel is located, and fill and repair the shadow area pixel; the backscattering coefficient value of the shadow area is usually low and inaccurate, and filling and repair can improve the integrity and usability of the image.

[0140] S64. According to lookup table A, perform geometric correction on the topographically radiometrically corrected SAR image, project the SAR image onto geospatial coordinates, and obtain the geocoded SAR product after topographic radiometric correction. Geometric correction can eliminate the geometric distortion of SAR images, giving them accurate geographic coordinate information, which facilitates fusion and analysis with other geographic data.

[0141] The present invention and its embodiments have been described above. This description is not restrictive, and the accompanying drawings are only one embodiment of the present invention; the actual structure is not limited thereto. In conclusion, if those skilled in the art are inspired by this description and design similar structures and embodiments without departing from the spirit of the invention, such designs should fall within the protection scope of the present invention.

Claims

1. A method for topographic radiometric correction of SAR images that takes into account overlapping and shadowed areas, characterized in that, Includes the following steps: S1. Data Acquisition: Acquire SAR and DEM images within the study area, and obtain the source files of the SAR images; S2. Image preprocessing: Perform multi-view processing on the original SAR image according to the resolution of the SAR image and DEM image, and resample the DEM image to obtain SAR image and DEM image with relatively consistent resolution. At the same time, perform radiometric calibration on the SAR image according to the radiometric calibration file or calibration constant. S3. Geometric positioning: After preprocessing, the DEM image and SAR image are geometrically positioned according to the RD model using the indirect positioning method. The correspondence between the DEM image and the calculated simulated SAR image, as well as the correspondence between the simulated SAR image and the real SAR image, are obtained. Lookup table A and lookup table B are generated. At the same time, the real SAR image is registered onto the simulated SAR image according to lookup table B. S4. Overlapping and Shadow Area Pixel Detection: Obtain overlapping area pixels on the DEM image based on the generated lookup table A, and then detect shadow area pixels on the DEM image based on the shadow area imaging principle of SAR image. S5. Calculate the effective scattering area: Calculate the surface area of ​​each pixel in the DEM image, project it onto the equiphase surface using the projection angle, and calculate its area. Then, the effective scattering area corresponding to each SAR pixel is calculated based on the overlapping area pixels, the shadow area pixels, and lookup table A; S6. Geometric and radiometric correction: Perform radiometric correction on the SAR image based on the calculated effective scattering area, and repair the backscattering coefficient values ​​of the SAR image in the shadowed area based on the detected overlapping pixels and shadowed pixels; after the radiometric correction is completed, perform geometric correction on the SAR image according to lookup table A.

2. The SAR image topographic radiometric correction method considering overlapping and shadowed areas according to claim 1, characterized in that: In step S2, the specific method for image preprocessing is as follows: S21. The DEM resampling rate is obtained by resampling the DEM image: Where θ refers to the angle of incidence, and δ dem This refers to the resolution of the DEM image, δ out This refers to the slant range resolution of the output analog image, and f refers to the minimum resampling rate of the DEM. S22. To balance the azimuth resolution of the DEM image, multi-look processing is performed on the SAR image to make the pixels of the SAR image as square as possible. The number of looks is calculated as follows: Where, δ r For slant range resolution, δ a This refers to the azimuth resolution. S23. Based on the radiometric calibration factor given in the SAR image source file, perform absolute radiometric calibration on the SAR image, converting the original DN value into radar brightness β. 0 .

3. The SAR image topographic radiometric correction method considering overlapping and shadowed areas according to claim 2, characterized in that: In step S3, the geometric positioning method is as follows: S31. Based on the positioning parameters in the SAR image source file, use the RD model to generate a simulated SAR image using the indirect positioning method, and store the correspondence between DEM image pixels and simulated SAR image pixels as a lookup table A. During the generation of SAR simulated images, the projection angle of each pixel in the DEM image is calculated, and a projection angle file is generated. The formula for calculating the projection angle is: Where η is the projection angle, Let be the normal vector of a pixel in the DEM image. This is the normal vector of the imaging surface corresponding to the pixel; S32. Based on the generated simulated SAR image, use the Harris corner detection algorithm to establish the pixel correspondence between the simulated SAR image and the real SAR image, and generate lookup table B. At the same time, register the real SAR image onto the simulated SAR image.

4. The SAR image topographic radiometric correction method considering overlapping and shadowed areas according to claim 3, characterized in that: In step S4, the specific methods for detecting overlapping and shadow region pixels are as follows: S41. Overlapping area detection: When multiple DEM image pixels correspond to one SAR image pixel, it is called overlapping. Therefore, according to the mapping relationship in lookup table A, we search for the case in lookup table A where one SAR image pixel corresponds to one DEM image pixel, directly detect the overlapping area pixels of the DEM image, and generate an overlapping area mask file. S42. Shadow area detection: Perform inverse geometric correction on the DEM image based on the coarse track data to convert it into an image with coordinates approximately the same as the uncorrected SAR image, and then detect shadow areas based on pixel occlusion.

5. A SAR image topographic radiometric correction method considering overlapping and shadowed areas according to claim 4, characterized in that: The method for detecting the shadow area is as follows: S421. Perform inverse geometric correction. Based on the satellite's orbit equation, calculate the angle between the satellite and the due north direction. Use this angle to convert the DEM image into an image with the same direction as the SAR image. S422. Detect the image line by line from left to right, set a 1×n detection window with a step size of 1 and n as the base number; take the current pixel as the detection center, if the elevation value of the pixel is greater than the elevation values ​​of the pixels on the left and right sides, then the pixel is the mountain peak point, otherwise it is not the mountain peak point; detect the pixels of the entire row in sequence to obtain the position of the mountain peak point in each row of the image. S423. Calculate the incident angle θ corresponding to each pixel and save it as an incident angle file. The formula for calculating the incident angle is: in, This is the satellite's position vector in the WGS84 or CGCS2000 coordinate system at this moment, with the direction pointing from the satellite towards the coordinate center. The distance vector between the satellite and the point at this moment, with the direction from the satellite to the ground point; S424. Detect the image line by line. When the first mountain peak is detected, check the incident angle θ corresponding to that point in the incident angle file. Detect the next point to the right and check the incident angle θ1. If θ1 is greater than θ, it is not a shadow point. Continue to the right to find the next mountain peak. If θ1 is less than or equal to θ, the point is a shadow point. Mark it as a shadow point and continue to detect the next point until a non-shadow point is detected. Then the shadow area occluded by that mountain peak is detected, and the detection of the next shadow area begins. S425. After all images have been inspected, a shadow area mask file is generated.

6. A SAR image topographic radiometric correction method considering overlapping and shadowed areas according to claim 5, characterized in that: The method for inverse geometric correction is as follows: S4211. Calculate the velocity vector of the satellite at the reference point during its transit time, based on the satellite orbit formula: S4212. Convert the satellite's spatial rectangular coordinates at the reference point to geodetic coordinates; S4213. Construct the basis vectors of the northeast celestial coordinate system: Where φ is the latitude of the satellite at the reference point, and λ is the longitude of the satellite at the reference point; S4214. Calculate the satellite's velocity vector in the east direction ΔE and the north direction ΔN in the northeast celestial coordinate system: S4215. Calculate the angle between the satellite's velocity direction and true north, i.e., the azimuth angle: S4216. Using the center point of the DEM image as the rotation center, rotate it according to the azimuth angle to make it an image with coordinates approximately the same as the uncorrected SAR image. S4217. Adjust the rotated DEM image again according to the satellite's flight direction; when the satellite flies from south to north, flip the DEM image vertically; when the satellite flies from north to south, flip the DEM image horizontally; this completes the inverse geometric correction of the DEM image.

7. A SAR image topographic radiometric correction method considering overlapping and shadowed areas according to claim 6, characterized in that: In step S5, the specific method for calculating the effective scattering area is as follows: S51. Use the triangle method to calculate the land area of ​​DEM image pixels, that is, divide each DEM image pixel into two triangles, and calculate the sum of the areas of the two triangles, which is the land area of ​​the DEM image pixel. S52. Based on the projection angle file, calculate the projected area of ​​each DEM image pixel on the equiphase surface, and generate a projected area file for the equiphase surface. The calculation formula is: γ dem (i,j)=σ dem (i,j)·sinη(i,j); in (i,j) is the row and column number corresponding to the pixel in the DEM image, σ dem η is the surface area of ​​the DEM image cell, and η is the projection angle of the DEM image cell; γ dem This represents the projected area of ​​the DEM image pixel on the equiphase surface; S53. The obtained overlay mask file and shadow mask file are superimposed to obtain shadow pixels in the overlay area. Although these pixels are located in the overlay area, they do not theoretically contribute to the backscattering coefficient. The shadow mask file is subtracted from these shadow pixels in the overlay area to generate a new shadow area file. S54. Based on the shadow pixels in the overlay area obtained in step S3, change the area of ​​these shadow pixels in the projection area file of the equiphase surface to 0.

8. The SAR image topographic radiometric correction method considering overlapping and shadowed areas according to claim 1, characterized in that: In step S6, the specific methods for geometric and radiometric corrections are as follows: S61. Combining lookup table A and the projected area file of the equiphase surface, calculate the effective scattering area corresponding to each SAR image pixel. The calculation formula is as follows: Where (r,a) represents the row and column number in the SAR image. γ represents the effective scattering area corresponding to the pixel at that location in the SAR image, (i,j) represents the row and column number of the pixel in the DEM, and γ dem It is the effective projected area of ​​the pixel point of the DEM image on the equiphase surface; S62. Based on the backscattering coefficient expression of the SAR image imaging surface, perform topographic radiometric correction on the polarimetric SAR image, as shown in the following formula: Where, β 0 The original radar brightness values ​​for the SAR image. γ is the scattering area corresponding to the imaging surface, i.e., the pixel size of the SAR image. 0 (r,a) represents the corrected SAR image pixel brightness value, i.e., the backscattering coefficient value; S63. Based on the new shadow area file detected by DEM and lookup table A, obtain the corresponding shadow area pixels in the SAR image; calculate the average backscattering coefficient value of the row where the shadow area pixel is located, and fill and repair the shadow area pixel. S64. According to lookup table A, perform geometric correction on the topographically radiometrically corrected SAR image, project the SAR image onto geospatial coordinates, and obtain the geocoded SAR product after topographic radiometric correction.

Images