Semi-automatic ground resolution distance measurement method for on-orbit medium resolution remote sensing image

Abstract

The application discloses a kind of semi-automatic ground resolution distance measurement methods of on-orbit medium resolution remote sensing image.The present application includes the following steps: step S1.natural target selection and target database construction;Step S2.target subpixel extraction and profile coordinate set generation;Step S3.the coordinate set of obtained profile extracts the gray value sequence of normalized image, constructs profile matrix and averages by row, obtains average edge spread function, carries out first-order difference operation to edge spread function, obtains line spread function, further carries out normalization processing to line spread function;Step S4.extract the peak value of normalized line spread function curve and its corresponding half-peak position, determine the distance between left and right half-peak points by linear interpolation, as the measure of width, multiply the measure of width with the ground sampling distance of satellite image, obtain ground resolution distance value;Step S5: result summary output.The present application is applied to the technical field of remote sensing image measurement or surveying and mapping geographic information service.

Description

Technical Field

[0001] This invention relates to the technical field of remote sensing image quality assessment and satellite on-orbit imaging performance evaluation, and particularly to a semi-automatic ground-resolution distance measurement method for medium-resolution remote sensing images on orbit. This method serves remote sensing image quality assessment, calibration, and on-orbit performance evaluation, and can also be applied to other surveying and mapping geographic information services and surveying services other than satellite application services. Background Technology

[0002] Ground Resolved Distance (GRD) represents the smallest ground feature size that a remote sensing satellite imaging system can distinguish under specific imaging conditions. It is a key indicator reflecting the actual spatial resolution performance of the imaging system.

[0003] Existing GRD measurement methods mainly fall into two categories: those based on manually calibrated targets and those based on natural targets. The manual target method, by capturing images of pre-laid high-contrast patterns on the ground and analyzing the imaging results, can accurately determine the system's resolution. However, target production is costly and deployment is complex, limiting its application to high-resolution satellite or aerial imagery and making it unsuitable for medium-resolution satellite on-orbit monitoring applications. In contrast, the natural target method utilizes ground features with clear edges in the imagery (such as building roofs, roads, and land boundaries) as references, estimating the actual resolution using indicators such as the Edge Spread Function (ESF), Line Spread Function (LSF), and Modulation Transfer Function (MTF). This method eliminates the need for manual targets, broadening its applicability. However, in medium-resolution imagery, it is susceptible to the complexity of ground features, noise interference, and edge blurring, making it difficult to consistently obtain high-precision results. Currently, a method for on-orbit resolution measurement that balances efficiency and accuracy is lacking. Existing research has validated the feasibility of natural target assessment on high-resolution imagery, but studies on medium-resolution satellite imagery remain limited, particularly lacking a systematic and repeatable on-orbit GRD measurement method. Furthermore, existing algorithms still suffer from limitations such as reliance on manual edge extraction, resulting in low automation; difficulty in precisely matching profile location with edge direction, affecting computational accuracy; and inconsistent processing procedures for edge spread function and line spread function, leading to significant deviations in Full Width at Half Maxima (FWHM) calculations. Therefore, a suitable method for measuring ground-resolved distances on on-orbit medium-resolution optical remote sensing satellites is urgently needed. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to overcome the shortcomings of the prior art and provide a semi-automatic ground resolution distance measurement method for on-orbit medium resolution remote sensing images. This method serves remote sensing imaging quality assessment, calibration, and on-orbit performance assessment, and can also be applied to other surveying and mapping geographic information services and surveying services other than satellite application services.

[0005] The technical solution adopted in this invention is as follows: This invention includes the following steps:

[0006] Step S1. Natural target selection and target library construction;

[0007] Step S2. Perform high-precision sub-pixel position extraction and line fitting on each natural target sub-image in the target library in step S1, and then generate a series of profiles with a preset number and the same length, and output the coordinate set of each profile.

[0008] Step S3. Based on the coordinate set of the profile obtained in step S2, extract the gray value sequence of the normalized image, construct the profile matrix and calculate the average by row to obtain the average edge spread function. Perform first-order difference operation on the edge spread function to obtain the line spread function. Further normalize the line spread function to normalize its peak value to 1 and eliminate imaging gain differences.

[0009] Step S4. Extract the peak value of the normalized linear diffusion function curve and its corresponding half-peak position. Determine the distance between the left and right half-peak points by linear interpolation as a measure of the system response width. The measure of this width is FWHM, which is used to characterize the spatial response characteristics of the imaging system to edge signals. Then, multiply the obtained FWHM value by the ground sampling distance of the satellite image to obtain the ground resolution distance value.

[0010] Step S5: Summarize the ground resolution distance values ​​calculated for all targets and take the average value as the overall image resolution index.

[0011] Furthermore, in step S3, in order to reduce noise interference and maintain the authenticity of the edge diffusion function curve shape, a linear interpolation method is used to enhance the sampling density between adjacent sampling points, thereby obtaining a smooth and structurally accurate enhanced edge diffusion function. Subsequently, a first-order difference operation is performed on the enhanced edge diffusion function to obtain a more realistic line diffusion function.

[0012] Furthermore, step S1 includes the following steps:

[0013] S11. Perform multi-level image enhancement preprocessing on the input satellite image to be measured, including anisotropic diffusion filtering, grayscale histogram statistics, conditional bidirectional contrast normalization, fuzzy histogram hyperbolic enhancement, median filtering and Gaussian filtering.

[0014] S12. After completing the multi-level image enhancement preprocessing in S11, perform Canny edge detection and output the binary edge map;

[0015] S13. Apply Hough transform to the edge binary map to detect straight lines and extract line parameters;

[0016] S14. Based on the line strength and length, the lines are selected as candidate lines using a relative threshold method based on line length.

[0017] S15. Use the threshold method to filter out straight lines that do not meet the step edge characteristics on both sides of the light and dark sides, and then manually review the images to remove straight lines with insufficient spectral contrast on both sides of the line.

[0018] S16. The straight lines that have been manually approved are cropped sequentially from the satellite images to be measured as natural target targets, and S11-S15 are repeated to complete the target library construction.

[0019] Furthermore, in step S2, the edge position of each natural target sub-image obtained in step S1 is first estimated with sub-pixel precision. The gray value of each pixel is regarded as the weighted average of the area covered by the edge when it passes through the pixel region and the gray values ​​on both sides, i.e., a partial area model. Assuming that the edge divides the pixel gray value of the natural target image into two parts A and B, the observed gray value F of the pixel is then... ij It can be represented as:

[0020]

[0021] Among them, F ij Let A be the observed value of pixel (i,j), A be the average gray value of the dark side region of the edge, B be the average gray value of the bright side region of the edge, and α be the average gray value of the pixel (i,j). ij ∈[0,1] represents the area ratio of the edge covering the bright side region within a pixel. By analyzing the grayscale distribution of pixels in the local window, the sub-pixel position, direction, and contrast of the edge can be retrieved, thus obtaining the coordinate set of the edge points of the natural target image. That is, target subpixel.

[0022] Then, the least squares method is used to fit a straight line to the coordinate set of edge points of the natural target image, and formula (14) is expressed as follows:

[0023]

[0024] in, , Let m represent a set of edge position coordinates, m be the slope of the fitted line, and b be the intercept; the fitted line is divided into equal parts along the X-axis. Section, denoted as For each segment i, calculate the corresponding point on the fitted line. The slope of the fitted line at point and intercept Formulas for calculating the parameters of each vertical section (15-16):

[0025]

[0026]

[0027] The length of the vertical profile is adaptively set according to the length and width of the input natural target image. After completing the above operations, a series of data and profiles of the same length are generated, and the (x,y) coordinate set of each generated profile is output and stored.

[0028] Furthermore, in step S3, based on the cross-sectional (x,y) coordinate set data obtained in step S2, the grayscale value of each cross-section is extracted and an edge diffusion function is generated to describe the response characteristics of the imaging system to the edge. First, the input natural target image is normalized, mapping the DN value of the image to the [0,1] interval. Then, the normalized grayscale value of the image is extracted using the (x,y) coordinate set of each cross-section. For the grayscale value sequence extracted from each cross-section, an L×N array is constructed. p The ESF matrix, where L is the section length and N is the cross-sectional length. p The total number of profiles is given. Each column in the matrix represents the gray value change sequence of a profile. The mean of the edge diffusion function (ESF) matrix is ​​calculated row by row to obtain the average ESF vector, which is calculated by formula (17):

[0029]

[0030] Among them, DN ik This represents the grayscale value of the k-th profile at the i-th pixel position;

[0031] After generating the edge diffusion function, differentiate it according to formula (18) to obtain the line diffusion function:

[0032] .

[0033] Furthermore, in step S4, the points on the left and right sides of the linear diffusion function curve that intersect with the half-peak value (0.5) are determined by linear interpolation, namely the left midpoint and the right midpoint, and the distance between the two is defined as FWHW, calculated as follows:

[0034]

[0035]

[0036]

[0037] in, , These represent the coordinates of two adjacent sampling points on the left side of the linear diffusion function curve; , This indicates the coordinates of two adjacent sampling points on the right.

[0038] Linear interpolation can accurately determine the left and right intersection points of the line spread function curve at the half-peak, thus obtaining the FWHM value in pixels. Multiplying this value by the ground sampling distance yields the ground resolution distance value, calculated using the following formula (22):

[0039] .

[0040] Furthermore, in step S5, the N effective target GRD values ​​obtained in step S4 are recorded as follows: The average value is calculated as the overall ground resolution distance index of the satellite to be measured, and the formula (23) is as follows:

[0041] (twenty three).

[0042] Furthermore, in step S11, anisotropic diffusion filtering is performed on the main bands of the image to suppress high-frequency noise while preserving significant edge information. The diffusion model of the diffusion filter satisfies the following partial differential equation (1):

[0043]

[0044]

[0045] Where I represents the original image and I' represents the image after diffusion filtering. For image gradient, Let represent the Laplacian operator of the image, i.e., the sum of the second-order partial derivatives, and c be the adaptive diffusion coefficient dependent on the local gradient of the image, used to control the diffusion intensity. The gradient of the diffusion coefficient c is represented by the gradient operator applied to c. The resulting vector reflects the rate of change of the diffusion coefficient c at various locations in the image space. K is the diffusion sensitivity factor, used to adjust the sensitivity of the diffusion model to gradient changes. Its value can be adaptively set according to the noise level of the image.

[0046] Perform grayscale histogram statistics on the anisotropically diffused filtered image I'(x,y) and calculate the brightness deviation coefficient. , formula (3):

[0047]

[0048] in, The average value of the image. This is the median gray value;

[0049] according to The value is used to perform bidirectional contrast normalization on the image, as expressed in formula (4) as follows:

[0050]

[0051] in, This is the brightness adjustment scaling factor;

[0052] After adjustment, result I was then... n Perform linear normalization to limit the output grayscale value to the range [0,1].

[0053] The image after conditional bidirectional contrast normalization is first normalized pixel by pixel using formula (5), then scaled using formula (6), and finally multiplied by the scaling formula and the normalization result using formula (7) to obtain the final enhanced pixel value. The entire image is denoted as . After this step, the local contrast of the image is enhanced, and the structural edges of low-contrast areas are significantly highlighted.

[0054]

[0055]

[0056]

[0057] in, This represents the grayscale value located in the i-th column and j-th row. Let L and λ be the minimum and maximum gray values ​​of the image, respectively, where L is the number of gray levels and λ is the scaling function. This refers to the final enhanced pixel value;

[0058] For the enhanced The image is filtered and smoothed. First, m×n median filtering is used for smoothing to remove isolated noise and unstructured details. The calculation formula (8) is as follows:

[0059]

[0060] in, This represents a local neighborhood window centered at pixel (x, y) with a size of m×n, where median represents the median grayscale value within that neighborhood window.

[0061] Subsequently, the median filtering results were analyzed. Gaussian smoothing is performed to suppress residual high-frequency noise, using an m×n Gaussian kernel. 1. Set the standard deviation σ, and its calculation formula (9) is as follows:

[0062]

[0063] In equation (9) This represents the convolution operation;

[0064] In step S12, the Canny edge detection uses the Canny operator to perform edge detection on the image. This algorithm combines gradient thresholding and non-maximum suppression strategies to extract edge results with strong structure and low noise interference; the output result is a single-band edge binary image. Its grayscale value takes only 0 and 1, representing non-edge and edge regions respectively;

[0065] In step S13, Hough transform is used to detect edges in a binary image. Edges with significant linear characteristics are used to obtain a set of parameterized lines, and their calculation formula (10) is as follows:

[0066]

[0067] Where ρ represents the distance from the line to the origin, and θ is the angle between the normal and the horizontal direction;

[0068] In step S14, all detected straight lines are sorted from longest to shortest. Given a total number of lines, the longest line is selected as a candidate. The calculation formula (11) is as follows:

[0069]

[0070] Where β is the relative threshold coefficient, L max The longest straight line length N top The settings are adaptively adjusted based on image features and detection sensitivity.

[0071] In step S15, a threshold method is used to filter candidate lines, eliminating targets that do not possess step edge features. First, the average gray value and gradient magnitude of the regions on both sides of the candidate line are calculated, and a threshold is set. and As a constraint on spectral contrast and gradient, the following formula (12) is satisfied:

[0072]

[0073] in, and These represent the average gray values ​​of the regions on either side of the line. The maximum value of the gradient magnitude is used. If the above conditions are not met, the target is determined to be an atypical step edge feature and is removed. After the automatic screening is completed, the remaining candidate edges are manually reviewed. The manual review mainly combines the image to visually verify the candidate edge targets to prevent the existence of atypical step edge feature targets. After automatic threshold filtering and manual review, the straight lines are retained as valid candidate edge targets.

[0074] In step S16, the effective candidate edge targets processed in step S15 are cropped sequentially on the satellite image to be measured. With the center of the straight line and the normal direction as references, a local sub-image region containing complete edge information is extracted. This sub-image region is a natural target sample with a size of H×W pixels. After cropping, the natural target images are uniformly numbered and stored in the target database. Then, steps S11 to S15 are repeated to realize the target extraction and screening of multiple satellite images to be measured, and finally construct a natural target library covering different land cover types and different brightness conditions.

[0075] The beneficial effects of this invention are: 1. By using algorithms for natural target edge detection, sub-pixel fitting, profile extraction, and linear enhancement, high-precision evaluation of the spatial resolution of satellite on-orbit images is achieved; 2. This invention can complete ground resolution distance measurement in a semi-automatic manner in natural ground feature images without standard targets, which can effectively improve the reliability and automation level of ground resolution distance measurement; 3. Stable and repeatable resolution measurement results can be obtained through a parameter-controllable edge spread function-line spread function analysis process. Attached Figure Description

[0076] Figure 1 This is the logic flowchart of the present invention;

[0077] Figure 2 This is the binary edge image after Canny edge detection is completed;

[0078] Figure 3 This is the result image of precise sub-pixel extraction and profile generation of the target;

[0079] Figure 4 This is a normalized grayscale value sequence curve for each profile line;

[0080] Figure 5 This is a graph showing the results of the enhanced edge spread function after a three-fold enhancement;

[0081] Figure 6 This is a graph showing the results of the average edge spread function;

[0082] Figure 7 This is a graph showing the results of the line diffusion function;

[0083] Figure 8This is a graph showing the correspondence between the line spread function and the FWHM value. Detailed Implementation

[0084] like Figures 1 to 8 As shown, in this embodiment, the present invention mainly comprises five parts: ① natural target selection and target library construction; ② target sub-pixel extraction and profile generation; ③ edge spread function construction and line spread function generation; ④ calculation of full width at half maximum (FWHM) and ground resolution distance; ⑤ result statistics and accuracy evaluation.

[0085] This invention uses the CMOS2 sensor (400-1000nm, 32 bands, GSD 10m) of the Chinese OHS-3A satellite as an example of the satellite to be measured, in order to implement and test the scheme of this invention. The specific steps are as follows:

[0086] Step S1. Natural target selection and target library construction;

[0087] Step S1 includes the following steps:

[0088] S11. The input satellite image to be measured is subjected to multi-level image enhancement preprocessing in sequence, including anisotropic diffusion filtering, grayscale histogram statistics, conditional bidirectional contrast normalization, fuzzy histogram hyperbolic enhancement, median filtering, and Gaussian filtering: Specifically, the input satellite image to be measured is 50560 pixels × 50560 pixels. The main bands of the image (such as red, green, and near-infrared) are selected to perform anisotropic diffusion filtering in order to suppress high-frequency noise while retaining significant edge information. The diffusion model of the diffusion filter satisfies the following partial differential equation (1):

[0089]

[0090]

[0091] Where I represents the original image and I' represents the image after diffusion filtering. For image gradient, Let represent the Laplacian operator of the image, i.e., the sum of the second-order partial derivatives, and c be the adaptive diffusion coefficient dependent on the local gradient of the image, used to control the diffusion intensity. The gradient of the diffusion coefficient c is represented by the gradient operator applied to c. The resulting vector reflects the rate of change (including the direction and magnitude of change) of the diffusion coefficient c at various locations in the image space. K is the diffusion sensitivity factor, used to adjust the sensitivity of the diffusion model to gradient changes. Its value can be adaptively set according to the noise level of the image.

[0092] Perform grayscale histogram statistics on the anisotropically diffused filtered image I'(x,y) and calculate the brightness deviation coefficient. , formula (3):

[0093]

[0094] in, The average value of the image. This is the median gray value.

[0095] according to The value is used to perform bidirectional contrast normalization on the image, as expressed in formula (4) as follows:

[0096]

[0097] in, The brightness adjustment scaling factor is taken as [value] in this implementation example. ;

[0098] After adjustment, result I was then... n Linear normalization is performed to limit the output grayscale value to the range of [0,1]. This operator can achieve bidirectional high and low contrast balance and avoid edge loss caused by strong reflection or shadow.

[0099] The image after conditional bidirectional contrast normalization is subjected to fuzzy histogram hyperbolic enhancement to improve local brightness contrast and enhance weak edge response. First, the pixel-by-pixel normalization formula (5) is applied, then the image is scaled using formula (6), and finally the scaling formula is multiplied by the normalization result using formula (7) to obtain the final enhanced pixel value. The entire image is denoted as . After this step, the local contrast of the image is enhanced, and the structural edges of low-contrast areas (such as light-colored buildings and dark roads) are significantly highlighted.

[0100]

[0101]

[0102]

[0103] in, This represents the grayscale value located in the i-th column and j-th row. Let L and λ be the minimum and maximum gray values ​​of the image, respectively, where L is the number of gray levels and λ is the scaling function. This refers to the final enhanced pixel value;

[0104] For the enhanced The image is filtered and smoothed. First, m×n median filtering is used for smoothing to remove isolated noise and unstructured details. The calculation formula (8) is as follows:

[0105]

[0106] in, This represents a local neighborhood window centered at pixel (x, y) with a size of m×n, where median represents the median grayscale value within that neighborhood window.

[0107] Subsequently, the median filtering results were analyzed. Gaussian smoothing is performed to suppress residual high-frequency noise, using an m×n Gaussian kernel. 1. Set the standard deviation σ, and its calculation formula (9) is as follows:

[0108]

[0109] In equation (9) This represents the convolution operation.

[0110] S12. After completing the multi-level image enhancement preprocessing in S11, perform Canny edge detection and output a binary edge image. Specifically, after multi-level preprocessing, the Canny operator is used to perform edge detection on the image. This algorithm combines gradient thresholding and non-maximum suppression strategies to extract edge results with strong structure and low noise interference; the output result is a single-band binary edge image. Its grayscale value takes only 0 and 1, representing non-edge and edge regions respectively, such as... Figure 2 The image shown is the result after Canny edge detection processing.

[0111] S13. Apply Hough Transform to the binary edge image to detect lines and extract line parameters. Specifically, use Hough Transform (Duda and Hart) to detect edge images. Edges with significant linear characteristics are used to obtain a set of parameterized lines, and their calculation formula (10) is as follows:

[0112]

[0113] Where ρ represents the distance from the line to the origin, and θ is the angle between the normal and the horizontal direction;

[0114] By adjusting the Hough transform accumulator threshold and angular resolution parameters, the system can adapt to different ground feature characteristics. This step can identify natural structures with strong step edges, such as road boundaries, roof outlines, and runways.

[0115] S14. Sort the lines according to their intensity and length, and use a length-based relative threshold filtering method to retain only the most significant lines. The specific steps are as follows:

[0116] Sort all detected straight lines in descending order of length. Assume a total number of lines, and select the longest line as a candidate. The calculation formula (11) is as follows:

[0117]

[0118] Where β is the relative threshold coefficient, L max The longest straight line length N top The method is adaptively set based on image features and detection sensitivity; it can effectively preserve the main straight lines representing the edges of image structures and eliminate weak edges and noise interference.

[0119] S15. A thresholding method is used to filter straight lines that do not meet the step edge characteristics on both sides of the light and dark sides. Then, the images are manually reviewed to remove straight lines with insufficient spectral contrast on both sides. Specifically, a thresholding method is used to filter candidate straight lines and remove targets that do not have step edge characteristics. First, the average gray value and gradient magnitude of the regions on both sides of the candidate straight line are calculated, and a threshold is set. and As a constraint on spectral contrast and gradient, the following formula (12) is satisfied:

[0120]

[0121] in, and These represent the average gray values ​​of the regions on either side of the line. The maximum value of the gradient magnitude is used. If the above conditions are not met, the target is determined to be an atypical step edge feature and is removed. After the automatic screening is completed, the remaining candidate edges are manually reviewed. The manual review mainly combines the image to visually verify the candidate edge targets to prevent the existence of atypical step edge feature targets. After automatic threshold filtering and manual review, the straight lines are retained as valid candidate edge targets.

[0122] S16. The lines that have been manually approved are cropped sequentially from the satellite images to be measured as natural target objects. S11-S15 are repeated to complete the construction of the target library. Specifically, the effective candidate edge targets processed in step S15 are cropped sequentially on the satellite images to be measured. With the center of the line and the normal direction as references, local sub-image regions containing complete edge information are extracted. This sub-image region is a natural target sample with a size of H×W pixels (e.g., 256×256). After cropping, the natural target images are uniformly numbered and stored in the target database. Then, steps S11 to S15 are repeated to realize the target extraction and screening of multiple satellite images to be measured, and finally a natural target library covering different land cover types and different brightness conditions is constructed.

[0123] Step S2. Target sub-pixel extraction and profile generation: Perform high-precision sub-pixel position extraction and line fitting on each natural target sub-image in the target library of Step S1, and then generate a series of profiles with a preset number and the same length, and output the coordinate set of each profile.

[0124] Step S2 specifically involves: performing high-precision edge location extraction and sub-pixel-level fitting on each natural target sub-image in the target library to obtain high-precision edge data suitable for resolution measurement. First, sub-pixel-precision edge location estimation is performed on each natural target image obtained in step S1. The key to this step is: treating the gray value of each pixel as the weighted average of the area covered by the edge when it crosses the pixel region and the gray values ​​on both sides, i.e., a partial area model. Assuming the edge divides the pixel gray value of the natural target image into two parts, A and B, then the observed gray value F of that pixel... ij It can be represented as:

[0125]

[0126] Among them, F ij Let A be the observed value of pixel (i,j), A be the average gray value of the dark side region of the edge, B be the average gray value of the bright side region of the edge, and α be the average gray value of the pixel (i,j). ij ∈[0,1] represents the area ratio of the edge covering the bright side region within a pixel. By analyzing the grayscale distribution of pixels in the local window, the sub-pixel position, direction, and contrast of the edge can be retrieved, thus obtaining the coordinate set of the edge points of the natural target image. That is, target subpixel.

[0127] Then, the least squares method is used to fit a straight line to the coordinate set of edge points of the natural target image, and formula (14) is expressed as follows:

[0128]

[0129] in, , Let m represent a set of edge position coordinates, m be the slope of the fitted line, and b be the intercept; the fitted line is divided into equal parts along the X-axis. Section, denoted as For each segment i, calculate the corresponding point on the fitted line. The slope of the fitted line at point and intercept Formulas for calculating the parameters of each vertical section (15-16):

[0130]

[0131]

[0132] The length of the vertical profile is adaptively set based on the length and width of the input natural target image. After completing the above operations, a series of data and profiles of the same length are generated, such as... Figure 3 The image shown is the target sub-pixel extraction and profile result after calculation by formula (13)-formula (16) in sequence. The (x,y) coordinate set of each generated profile is output and stored.

[0133] Step S3. Edge diffusion function construction and line diffusion function generation: Based on the coordinate set of the profile obtained in step S2, extract the gray value sequence of the normalized image, such as... Figure 4 The image shows the normalized grayscale value sequence curves for each profile line. A profile matrix is ​​constructed, and the average edge spread function is obtained by averaging the rows. To reduce noise interference and maintain the authenticity of the edge spread function curve shape, linear interpolation is used to increase the sampling density between adjacent sampling points, thereby obtaining a smooth and structurally accurate enhanced edge spread function, as shown below. Figure 5 The image shows the result of the enhanced edge spread function after three-fold enhancement. Then, a first-order difference operation is performed on the enhanced edge spread function to obtain a more realistic line spread function. The line spread function is then normalized to make its peak value 1, thus eliminating the difference in imaging gain.

[0134] Step S3 specifically involves: based on the profile (x,y) coordinate set data obtained in step S2, extracting the grayscale value of each profile and generating an edge diffusion function to describe the response characteristics of the imaging system to edges. First, the input natural target image is normalized, mapping the DN value of the image to the [0,1] interval. Then, the normalized grayscale value of the image is extracted using the (x,y) coordinate set of each profile. For the grayscale value sequence extracted from each profile, an L×N array is constructed. p The ESF matrix, where L is the section length and N is the cross-sectional length. p The total number of profiles is represented by a column in the matrix, where each column represents a sequence of grayscale value changes for one profile. The mean of the edge spread function (ESF) matrix is ​​calculated row-wise, as shown below. Figure 6 As shown, the average ESF vector is obtained, and its calculation formula is (17):

[0135]

[0136] Among them, DN ik This represents the grayscale value of the k-th profile at the i-th pixel position;

[0137] After generating the ESF, differentiate it according to formula (18), as follows: Figure 7 As shown, the line spread function is obtained:

[0138]

[0139] To preserve the true shape of the line spread function (LSF), this invention uses linear interpolation results for derivative calculations.

[0140] Step S4. Calculation of full width at half maximum (FWHM) and ground resolution distance: Extract the peak value of the normalized linear diffusion function curve from step S3 and its corresponding half-peak position. Determine the distance between the left and right half-peak points through linear interpolation, which is used as a measure of the system response width. This width is called FWHM and is used to characterize the spatial response characteristics of the imaging system to edge signals. Then, multiply the obtained FWHM value by the ground sampling distance GSD of the satellite image to obtain the ground resolution distance value GRD.

[0141] Step S4 specifically involves the following: The normalized line spread function (FWHW) from step S3 typically does not possess a continuous Gaussian distribution. Due to the resolution limitations of medium-resolution optical satellites, its noise and edge structure characteristics lead to local discontinuities or asymmetries in the FWHW curve. Therefore, this invention uses linear interpolation to determine the points on the left and right sides of the FWHW curve that intersect with the half-peak value (0.5), namely the left midpoint (LMP) and the right midpoint (RMP), and defines the distance between them as FWHW. The calculation formula is as follows:

[0142]

[0143]

[0144]

[0145] in, , These represent the coordinates of two adjacent sampling points on the left side of the LSF curve; , This indicates the coordinates of two adjacent sampling points on the right.

[0146] Linear interpolation can accurately determine the position of the left and right intersection points of the LSF curve at the half-peak, thus obtaining the FWHM value in pixels. Figure 8 As shown, the FWHM value calculated by formula (19)-formula (21) is 1.1954 pixels. Multiplying it by the ground sampling distance GSD, the ground resolution distance value GRD can be obtained. The calculation formula (22) is as follows:

[0147]

[0148] Step S5. Results Summary and Output: Summarize the ground resolution distance values ​​(GRD) calculated for all targets and take the average value as the overall image resolution index;

[0149] Step S5 specifically involves: Recording the N effective target ground resolution distance values ​​GRD obtained in step S4 as... The average value is calculated as the overall ground resolution distance index of the satellite to be measured, and the formula (23) is as follows:

[0150] (twenty three).

[0151] Although the embodiments of the present invention are described with reference to actual solutions, they do not constitute a limitation on the meaning of the present invention. Modifications to the embodiments and combinations with other solutions based on this specification will be obvious to those skilled in the art.

Claims

1. A semi-automatic ground resolution distance measuring method for medium resolution remote sensing images in orbit, characterized in that, It comprises the following steps: Step S1. Natural target selection and target library construction; Step S2. Perform high-precision sub-pixel position extraction and straight line fitting on each natural target sub-image in the target library of step S1, and then generate a series of profiles with a preset number and the same length, and output the coordinate set of each profile; Step S3. Extract the gray value sequence of the normalized image according to the coordinate set of the profile obtained in step S2, construct a profile matrix and average it by row to obtain an average edge spread function, perform first-order difference operation on the edge spread function to obtain a line spread function, and further normalize the line spread function to normalize its peak value to 1, eliminating the difference in imaging gain; Step S4. Extract the peak value of the normalized line spread function curve and its corresponding half-peak position, determine the distance between the left and right half-peak points by linear interpolation as the measurement value of the system response width, which is FWHM, used to represent the spatial response characteristics of the imaging system to the edge signal, then multiply the obtained FWHM value with the ground sampling distance (GSD) of the satellite image to obtain the ground resolution distance value (GRD); Step S5. Collect all the ground resolution distance values (GRD) calculated from the targets to obtain the average value as the overall resolution index of the image; The step S2 firstly carries out the edge position estimation of sub-pixel accuracy to each natural target sub-image obtained in the step S1, and the gray value of each pixel is regarded as the weighted average value of the area covered by the edge passing through the pixel area and the gray value on both sides, that is, the partial area model. Assuming that the edge divides the pixel gray value of the natural target image into two parts A and B, the observed gray value F of the pixel is ij which can be expressed as: ， where F ij is the observed value of pixel (i,j), A is the average gray value of the dark side region of the edge, B is the average gray value of the bright side region of the edge, and a ij ∈[0,1] represents the proportion of the area of the bright side region covered by the edge in the pixel, by analyzing the gray distribution of the pixels in the local window, the sub-pixel position, direction and contrast of the edge can be inverted, and the edge point coordinate set of the natural target image is obtained , that is, the target sub-pixel Then, the least squares method is used to perform straight line fitting on the edge point coordinate set of the natural target image, and the formula (14) is as follows: ， wherein, , represents a set of edge position coordinates, m is the slope of the fitted line, and b is the intercept; the fitted line is equally divided into segments along the X-axis, denoted as For each segment i, the slope and intercept of the fitted line at the corresponding point on the fitted line are calculated, and the parametric equations (15-16) of each vertical section are calculated. ， ， The length of the vertical profile is adaptively set according to the length and width of the input natural target image. After the above operation, a series of data and profiles with the same length are generated, and the (x, y) coordinate set of each generated profile is output and stored.

2. The method of claim 1, wherein the method further comprises: In step S3, in order to reduce noise interference and maintain the authenticity of the edge spread function curve shape, the linear interpolation method is used to enhance the sampling density between adjacent sampling points, so as to obtain a smooth and structure-accurate enhanced edge spread function, and then a first-order difference operation is performed on the enhanced edge spread function to obtain a more real line spread function.

3. The method of claim 1, wherein the method further comprises: The step S1 comprises the following steps: S11. Perform anisotropic diffusion filtering, gray histogram statistics, conditional bidirectional contrast normalization, fuzzy histogram hyperbolic enhancement, median filtering and multi-level image enhancement preprocessing on the input satellite image to be measured; S12. After completing the multi-level image enhancement preprocessing of S11, perform Canny edge detection to output an edge binary image; S13. Apply Hough transform to detect straight lines and extract straight line parameters on the edge binary image; S14. Sort the straight lines according to the line intensity and length, and use the relative threshold method based on the line length to filter the straight lines as candidate straight lines; S15. Filter the straight lines that do not meet the step edge feature on both light and dark sides by threshold method, and then manually review the image to remove the straight lines with insufficient spectral contrast on both sides of the straight lines; S16. Cut the straight lines that pass the manual review from the satellite image to be measured as natural target objects, and repeat S11-S15 to complete the target library construction.

4. The method of claim 1, wherein the method further comprises: determining a ground resolution distance of the image; and determining a ground resolution distance of the image based on the ground resolution distance of the image and the scale factor. The step S3 is based on the profile (x, y) coordinate set data obtained in the step S2, extracts the gray value of each profile and generates an edge spread function to describe the response characteristics of the imaging system to the edge. First, the input natural target image is normalized to map the DN value of the image to the interval [0, 1]. Then, the gray value of the normalized image is extracted by using each profile (x, y) coordinate set. For the gray value sequence extracted for each profile, an LxN p ESF matrix is constructed, L is the length of the profile, N p is the total number of profiles, and each column in the matrix represents a gray value change sequence of a profile. The average ESF vector is obtained by calculating the average value of the edge spread function ESF matrix by row, and the calculation formula (17) is as follows. ， where DN ik represents the gray value of the kth profile at the i th pixel position; After generating the edge spread function, derive it according to formula (18) to obtain the line spread function: 。 5. The method of claim 4, wherein the method further comprises: The step S4 determines the left and right intersection points of the line spread function curve with the half peak (0.5) by linear interpolation, i.e. the left midpoint (LMP) and the right midpoint (RMP), and defines the FWHM as the distance between the two, and the calculation formula is as follows: ， ， ， wherein, , respectively represent the coordinates of the left two adjacent sampling points on the curve of the line spread function; , represent the coordinates of the right two adjacent sampling points. The left and right intersection points of the line spread function curve with the half peak can be accurately determined by linear interpolation, so as to obtain the FWHM value in pixel units, and the ground resolution distance (GRD) can be obtained by multiplying the ground sampling distance (GSD), and the calculation formula (22) is as follows: 。 6. The method of claim 5, wherein: The step S5 records the N effective target GRD values obtained in the step S4 as , and calculates the average value as the overall ground resolution distance value index of the satellite to be measured, and the formula (23) is as follows: （23）。 7. The method of claim 3, wherein the method further comprises: determining a ground resolution distance of the image; and determining a ground resolution distance of the image based on the ground resolution distance of the image and the scale factor. In the step S11, anisotropic diffusion filtering is performed on the selected image main band, so as to retain significant edge information while suppressing high-frequency noise, and the diffusion model of the diffusion filtering satisfies the following partial differential equation (1): ， ， where I is the original image, I' is the image after diffusion filtering, is the image gradient, is the Laplacian of the image, i.e., the sum of the second-order partial derivatives, c is an adaptive diffusion coefficient depending on the local gradient of the image, used to control the diffusion strength, is the gradient of the diffusion coefficient c, K is a diffusion sensitivity factor, used to adjust the sensitivity of the diffusion model to the gradient variation; The gray scale histogram of the anisotropic diffusion filtered image I'(x, y) is counted to calculate the brightness deviation coefficient , formula (3): ， wherein is the mean value of the image, is the median gray value; According to The conditional bidirectional contrast normalization is performed on the image, and formula (4) is as follows: ， wherein, is a luminance adjustment scale factor; After adjustment, the results I are linearly normalized so that the output gray scale values are limited to the range [0, 1]. n After adjustment, the results I are linearly normalized so that the output gray scale values are limited to the range [0, 1]. The image processed by the conditional bidirectional contrast normalization is firstly normalized by formula (5) pixel by pixel, then scaled by formula (6), and finally the scaled formula is multiplied by the normalized result to obtain the final enhanced pixel value, and the whole image is denoted as After the processing, the local contrast of the image is enhanced, and the structure edges in the low-contrast region are significantly highlighted. ， ， ， wherein, represents the gray value at the i-th column and j-th row, respectively the minimum and maximum gray value of the image, L is the number of gray levels, and λ is a scaling function, is the final enhanced pixel value. For the enhanced The image is filtered and smoothed. First, m×n median filtering is used for smoothing to remove isolated noise and unstructured details. The calculation formula (8) is as follows: ， wherein, represents a local neighborhood window centered at pixel (x, y) with size m x n, median denotes taking the median value of the gray levels in the neighborhood window, Subsequently, the median filtering results were analyzed. Gaussian smoothing is performed to suppress residual high-frequency noise, using an m×n Gaussian kernel.

1. Set the standard deviation σ, and its calculation formula (9) is as follows: ， (9) in the formula denotes a convolution operation; In step S12, the Canny edge detection uses the Canny operator to perform edge detection on the image. This algorithm combines gradient thresholding and non-maximum suppression strategies to extract edge results with strong structure and low noise interference; the output result is a single-band edge binary image. Its grayscale value takes only 0 and 1, representing non-edge and edge regions respectively; The step S13 employs the Hough transform to detect the edge binary image The edges with significant linear features in the image are obtained, and a set of parameterized straight lines is obtained, and the calculation formula (10) is as follows: ， Wherein, ρ represents the distance of a straight line to the origin, and θ is the included angle between the normal and the horizontal direction; In the step S14, all detected straight lines are sorted in descending order of length, and the total number is set, and the longest straight line is taken as a candidate, and the calculation formula (11) is as follows: ， wherein β is a relative threshold coefficient, L max is the longest straight line length N top According to the image features and the detection sensitivity, the threshold is adaptively set. The threshold method is used in the step S15 to filter the candidate straight lines, and the targets without the step edge feature are removed. First, the average gray value and the gradient amplitude of the regions on both sides of the candidate straight line are calculated, and the threshold is set With As the spectral contrast and gradient constraint condition, when the following formula (12) is satisfied: ， wherein, and respectively represent the average gray value of the two regions on both sides of the straight line, is the maximum value of the gradient amplitude, if the above condition is not met, it is determined as a non-typical step edge feature target and is rejected, after the automatic screening is completed, manual visual review is performed on the remaining candidate edges, manual review is a visual review of the candidate edge targets combined with the image to prevent the existence of non-typical step edge feature targets, after the automatic threshold filtering and manual double review, the remaining straight line is used as an effective candidate edge target; In the step S16, the effective candidate edge targets processed in the step S15 are sequentially cropped on the satellite image to be measured, the local sub-image region containing complete edge information is extracted with the straight line center and the normal direction as the reference, the sub-image region is a natural target sample, the size is set as HxW pixels, the cropped natural target image is uniformly numbered and stored into the target database, and then the steps S11 to S15 are repeatedly executed, so that the target extraction and screening of multiple scenes of satellite images to be measured are realized, and finally the natural target database covering different ground object types and different brightness conditions is constructed.

Images