Remote sensing satellite in-orbit relative radiation calibration method and system considering geometric radiation characteristic constraint

Abstract

The invention discloses a remote sensing satellite in-orbit relative radiation calibration method and system considering geometric radiation characteristic constraints. The method comprises the following steps: S1, finishing time sequence image alignment based on geometric information constraints; s2, estimating a time sequence image earth surface truth value based on geometric radiation characteristic constraints among the images; s3, extracting a gray value of each detection element of the area array on each frame sequence; s4, constructing an area array load radiation calibration model; and S5, radiometric calibration model resolving is carried out, in-orbit relative radiometric calibration work of the area array type remote sensing satellite is completed through the steps S1-S4, an in-orbit relative radiometric calibration coefficient is obtained, and input is provided for radiometric correction of an original image. The method is suitable for in-orbit relative radiometric calibration of area array type remote sensing satellites, and has the advantages of high frequency, high efficiency, more convenience and low cost.

Description

Technical Field

[0001] The present invention relates to an on-orbit relative radiometric calibration method for a remote sensing satellite carrying a planar array type payload, and particularly to an on-orbit relative radiometric calibration method and system for a remote sensing satellite that images in a planar array high frame rate mode. Background Art

[0002] In recent years, the trend of diversified development of optical remote sensing satellites has become more obvious, and many remote sensing satellites with different focuses have emerged, such as agile satellites that can highly maneuver to image multiple targets, video satellites that can dynamically monitor the ground, and high-resolution night light remote sensing satellites that can observe human activities at night, etc., providing various types of remote sensing observation images for remote sensing applications. It has become an urgent need to continuously monitor global hot spots and targets using remote sensing satellites to obtain dynamic information. High frame rate remote sensing satellites, such as video satellites, can obtain high-dynamic images of targets within a certain time range, and have the ability to continuously monitor targets, becoming a hot spot in recent years.

[0003] On-orbit relative radiometric calibration is a key technology to ensure the radiation quality of satellites and an indispensable part of the ground processing system after the satellite is on orbit. However, due to the influence of launch vibrations during the satellite launch process, and the drastic changes in physical environments such as temperature, heat, and space environment after the satellite is on orbit, as well as the attenuation problem of the satellite sensor detector elements that changes over time, the on-orbit response state of the satellite sensor has periodic complex changes, resulting in the inability to apply the pre-launch laboratory calibration results of the remote sensing satellite on orbit for a long time, directly reducing the quality of satellite imaging data. In-orbit radiometric calibration is urgently needed to improve the image quality. There have been related studies on the on-orbit relative radiometric calibration of planar array type remote sensing satellites, mostly focusing on using the method of imaging a large area of uniform field on the ground to complete on-orbit calibration. However, such methods have a high dependence on the large area of uniform field on the ground, and the calibration frequency is limited; with the improvement of the resolution of planar array type remote sensing satellites, the large area of uniform field that can meet the calibration requirements will be extremely reduced globally. Therefore, studying the on-orbit high-frequency relative radiometric calibration technology for planar array type remote sensing satellites, getting rid of the dependence on the large area of uniform field on the ground, and improving the satellite radiation quality are very important for ensuring the efficient earth observation application of planar array type remote sensing satellites. Summary of the Invention

[0004] The present invention provides a new on-orbit relative radiometric calibration method for a remote sensing satellite carrying a planar array type payload. First, based on geometric information constraints, the temporal images are aligned to eliminate the geometric distortion of multi-frame sequence images and restore the motion trajectory of the planar array detector elements on the multi-frame sequence images; secondly, based on the temporal image ground truth estimation of the geometric radiation characteristics between images, the relative radiometric calibration reference is obtained, a radiometric calibration model is constructed, and the on-orbit relative radiometric calibration of the planar array type remote sensing satellite is realized.

[0005] The object of the present invention is to provide a method for on-orbit high-frequency relative radiometric calibration of a remote sensing satellite suitable for carrying area array sensors, specifically including the following steps:

[0006] Step 1, complete the alignment of sequential images based on geometric information constraints;

[0007] Step 2, based on the geometric and radiometric characteristics constraints between images and the prior probability estimation method of Bayesian theory, realize the estimation of the true surface value of sequential images

[0008] Step 3, based on Step 1, restore the motion trajectories of area array detector elements in sequential images, and obtain the gray values of each detector element in each frame sequence;

[0009] Step 4, based on the gray values of area array detector elements and the true surface value obtained in Step 2 and Step 3, construct a radiometric calibration model for each detector element of the area array payload;

[0010] Step 5, solve the radiometric calibration model in Step 4 for all detector elements of the area array payload, obtain the relative radiometric calibration coefficients of each detector element of the area array payload, and complete the on-orbit relative radiometric calibration of the area array remote sensing satellite.

[0011] Further, in Step 1, use geometric calibration to complete system error compensation, generate an undistorted RFM model based on the calibrated interior orientation element supplementary model, eliminate the geometric distortion within the scene and the deformation between scenes of the sequential image data, realize the alignment of sequential images, and obtain the one-to-one correspondence between the sequential image coordinates (x, y) and the reference image coordinates (X base , Y base ).

[0012] Further, the specific implementation method of Step 1 is as follows:

[0013] 1.1, System error compensation: The rational polynomial model (RFM) establishes the conversion relationship between the image coordinates (x, y) and the ground coordinates (X, Y, Z) through a geometric polynomial;

[0014] 1.2, Geometric distortion elimination: Use the image of the calibration area to calibrate the interior orientation element error compensation model, and then generate an undistorted RFM based on the following steps:

[0015] 1) Use the solved interior orientation element error compensation model to establish the conversion relationship between the image coordinates and the ground coordinates according to Step 1.1;

[0016] 2) According to the conversion relationship between the image coordinates and the ground coordinates, calculate the control grid (x, y, X, Y, Z) used to generate the undistorted RFM;

[0017] 3) Use the control grid generated in 2) to calculate a new RFM;

[0018] 1.3, Temporal Image Alignment: Based on the distortion-free RFM, a high-precision registration method using the frame-to-frame positioning consistency constraint; the registration steps are as follows:

[0019] ① Use the SIFT algorithm to extract corresponding points (x, y), (X base , Y base ) between the temporal frame images and the reference image;

[0020] ② Use the distortion-free RFM of the reference image to calculate the corresponding ground coordinates (X base , Y base ) of the matching corresponding points (X base , Y base ) on the reference image, and form control points (X

[0021] ③ Use the control points in ② to solve for the distortion-free RFM orientation parameters of the temporal frame image;

[0022] ④ Generate a registered result image of the same size as the reference image;

[0023] ⑤ For any point (x res , y res ) in the result image, use the distortion-free RFM of the reference image to calculate its ground coordinates (X res , Y res , Z res );

[0024] ⑥ Based on the distortion-free RFM orientation results of other sequence data of the temporal image in ③, calculate the pixel coordinates (x, y) of each image sequence corresponding to (X res , Y res , Z res ), calculate the gray value of (x, y), and assign it to (x res , y res );

[0025] ⑦ Repeat ⑤ - ⑥ to generate the registered result image of the temporal image, and obtain the one-to-one correspondence relationship (x, y), (X base , Y base ) between each pixel of the temporal frame image and the reference image, that is, the precise relative position relationship of the temporal image.

[0026] Furthermore, the conversion relationship between the image coordinates (x, y) and the ground coordinates (X, Y, Z) in step 1.1 is given by the following formula:

[0027]

[0028] Num x (X, Y, Z) = a0 + a1Z + a2Y + a3X + a4ZY + a5ZX + a6YX + a7Z2 +a8Y 2 +a9X 2

[0029] +a 10 ZYX+a 11 Z 2 Y+a 12 Z 2 X+a 13 Y 2 Z+a 14 Y 2 X+a 15 ZX 2 +a 16 YX 2

[0030] +a 17 Z 3 +a 18 Y 3 +a 19 X 3

[0031] Den x (X,Y,Z) = b0 + b1Z + b2Y + b3X + b4ZY + b5ZX + b6YX + b7Z 2 +b8Y 2 +b9X 2

[0032] +b 10 ZYX+b 11 Z 2 Y+b 12 Z 2 X+b 13 Y 2 Z+b 14 Y 2 X+b 15 ZX 2 +b 16 YX 2

[0033] +b 17 Z 3 +b 18 Y 3 +b 19 X 3

[0034] Num y (X,Y,Z) = c0 + c1Z + c2Y + c3X + c4ZY + c5ZX + c6YX + c7Z 2 +c8Y 2 +c9X 2

[0035] +c 10 ZYX + c 11 Z 2 Y + c 12 Z 2 X + c 13 Y 2 Z + c 14 Y 2 X + c 15 ZX 2 +c 16 YX 2

[0036] +c 17 Z 3 +c 18 Y 3 +c 19 X 3

[0037] Den y (X, Y, Z) = d0 + d1Z + d2Y + d3X + d4ZY + d5ZX + d6YX + d7Z 2 +d8Y 2 +d9X 2

[0038] +d 10 ZYX + d 11 Z 2 Y + d 12 Z 2 X + d 13 Y 2 Z + d 14 Y 2 X + d 15 ZX 2 +d 16 YX 2

[0039] +d 17 Z 3 +d 18 Y 3 +d 19 X 3

[0040] where a i , b i , c i , d i are polynomial coefficients; the RFM geometric positioning error comes from the satellite imaging geometric parameters, that is, the positioning error (x e , y e ) caused by the exterior orientation elements and the positioning error (x i , y i ) caused by the interior orientation elements, and there is:

[0041]

[0042] Among them, Num x , Den x , Num y , Den y are intermediate quantities.

[0043] Furthermore, the specific implementation of step 2 includes the following sub-steps:

[0044] 2.1, Establish the following model relationship between the true surface image X and the image Y obtained by the remote sensing payload imaging:

[0045] Y i = D i B i F i X + n i

[0046] In the formula, i is the sequence image number, i = 1, 2,..., N, D is the motion parameter matrix, B is the remote sensing payload diffusion function matrix, F is the resampling matrix, and n is the noise matrix;

[0047] 2.2, Obtain X in the above formula. According to the probability estimation theory, under the constraint of prior knowledge, the objective function for estimating the true ground object image can be written in the following form:

[0048]

[0049] In the formula, i is the sequence image frame number, Y i is the i-th frame image. Therefore, the first part of the function ‖Y i - D i B i F i X‖ p reflects the brightness consistency between the sequence images and the estimated true image, and p is the order of the norm; the second part ‖Y i -(D i - D i+1 )Y i+1 ‖ q reflects the brightness consistency between adjacent sequence images, and q is the order of the norm; the third part is the regularization term R(X) and the regularization parameter λ;

[0050] 2.3, Select the M-estimation in robust estimation as the first term of the objective function for the motion parameter matrix D between sequence images, then:

[0051]

[0052] In the formula, N is the number of sequence images, M is the number of detector elements of the area array payload, e i,nis the residual between the digital quantization value DN of the nth pixel in the ith image and the DN value of the nth pixel in the estimated image. ρ(e) is a loss function, which is symmetric, greater than or equal to 0, and monotonic when e is greater than or less than 0; j i (n) is the detector element number that images the nth pixel in the ith image, DN i (j i (n)) is the DN value of the j i (n)th detector element in the ith image of the area array payload, DN True (n) is the DN value of the nth pixel in the estimated image;

[0053] 2.4, the ρ(e) loss function is defined as follows:

[0054]

[0055] a is a parameter;

[0056] 2.5, select bilateral filtering as the regularization term, and the final objective function can become:

[0057]

[0058] where w is the weight function, P is the one-dimensional bilateral filter kernel parameter, and the filter kernel size is 2 / P + 1, respectively representing a translation of l pixels in the horizontal direction and a translation of h pixels in the vertical direction;

[0059] 2.6, finally realize the estimation of the ground truth and obtain the ground truth DN true (x, y).

[0060] Furthermore, the calculation method of the parameter a is as follows:

[0061]

[0062] where E i represents the average residual on the ith image, and Max and Min are the maximum and minimum value functions respectively.

[0063] Furthermore, in step 3, the method for extracting the gray value of the area array detector element on each frame sequence is as follows:

[0064] 3.1, it is known that the coordinates of the area array detector element D on the reference image are (X base_D , Y base_D );

[0065] 3.2, the coordinates (x D , y D ) of the detector element D on a certain frame image are obtained from step 1;

[0066] 3.3. According to the coordinates (x D , y D ) of the detection element D on a certain frame of image, the gray value DN D can be extracted.

[0067] Further, in step 4, the radiation calibration model of each detection element of the area array payload is as follows:

[0068] DN true (j) = a i × DN i,j + b i

[0069] where a i , b i are the relative radiation calibration coefficients of the i-th detection element of the area array payload, j is the frame sequence number, y ≤ N, and N is the number of frame sequences.

[0070] Further, in step 5, the least squares method is used to solve the radiation calibration model in step 4.

[0071] The present invention also provides a relative on-orbit radiation calibration system for a remote sensing satellite considering geometric radiation characteristics constraints, including:

[0072] A processor and a memory, the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the relative on-orbit radiation calibration method for a remote sensing satellite considering geometric radiation characteristics constraints as described in the above technical solution.

[0073] Compared with the prior art, the present invention has the following characteristics and beneficial effects:

[0074] (1) A relative on-orbit radiation calibration method for area array remote sensing satellites based on time-series frame images is proposed.

[0075] (2) On-orbit calibration does not require imaging of a uniformly distributed ground object.

[0076] (3) The satellite does not need to have on-board calibration processing capabilities.

[0077] (4) Any frame sequence can be used to complete calibration, greatly improving the on-orbit calibration frequency of area array remote sensing satellites. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] Figure 1 is a schematic diagram (partial) of the radiation calibration reference extracted in the embodiment of the present invention;

[0079] Figure 2 is a schematic diagram of the radiation calibration model of each detection element of the area array payload in the embodiment of the present invention;

[0080] Figure 3Effect diagrams before and after relative radiometric calibration in the embodiments of the present invention, where (a) and (c) are before relative radiometric calibration, and (b) and (d) are after relative radiometric calibration. Detailed implementation manners

[0081] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

[0082] A method for on-orbit relative radiometric calibration of a remote sensing satellite considering geometric and radiometric characteristics constraints provided by the present invention specifically includes the following steps:

[0083] S1. Complete the alignment of sequential images based on geometric information constraints, and the specific implementation is as follows:

[0084] Use geometric calibration to complete system error compensation, generate an undistorted RFM model based on the calibrated interior orientation element supplementary model, eliminate the geometric distortion within the scene and the deformation between scenes of the sequential image data, and achieve the alignment of sequential images, obtaining the one-to-one correspondence between the sequential image coordinates (x, y) and the reference image coordinates (X base , Y base ), where x, y ∈ [n, m], and n and m are the width and height of the image. The further specific steps are as follows:

[0085] 1.1 System error compensation: The Rational Function Model (RFM) establishes the conversion relationship between the image coordinates (x, y) and the ground coordinates (X, Y, Z) through a geometric polynomial, which is given by the following formula.

[0086]

[0087]

[0088] In the formula, a i , b i , c i , d i , i ≤ 19 are polynomial coefficients, Num x , Den x , Num y , Den y are intermediate quantities. The geometric positioning error of RFM mainly comes from the satellite imaging geometric parameters, that is, the positioning error (x e , y e ) caused by the exterior orientation elements and the positioning error (x i , y i ) caused by the interior orientation elements, and there is:

[0089]

[0090] Among them, the calibration method of internal and external orientation elements can be implemented by using the geometric calibration method in relevant public literature.

[0091] 1.2, Geometric distortion elimination: Use the calibration image of the calibration area to calibrate the internal orientation element error compensation model, and then generate an undistorted RFM based on the following steps:

[0092] 1) Use the solved internal orientation element error compensation model to establish the conversion relationship between image coordinates and ground coordinates according to step 1.1:

[0093] 2) Calculate the control grid (x, y, X, Y, Z) for generating the undistorted RFM according to the conversion relationship between image coordinates and ground coordinates;

[0094] 3) Use the control grid generated in 2) to calculate a new RFM.

[0095] The RFM generated through 1)-3) eliminates the influence of high-order distortion of internal orientation elements.

[0096] 1.3, Temporal image alignment: Based on the undistorted RFM, adopt a high-precision registration method with frame-to-frame positioning consistency constraints. The registration steps are as follows:

[0097] ① Use the SIFT algorithm to extract corresponding points (x, y), (X base , Y base ) between the temporal frame images and the reference image;

[0098] ② Use the undistorted RFM of the reference image to calculate the corresponding ground coordinates (X, Y, Z) of the matching corresponding points (X base , Y base ) on the reference image, and form control points (X base , Y base , X, Y, Z);

[0099] ③ Use the control points in ② to solve the undistorted RFM orientation parameters of the temporal frame image;

[0100]

[0101] ④ Generate a registered result image of the same size as the reference image;

[0102] ⑤ For any point (x res , y res ) in the result image, use the undistorted RFM of the reference image to calculate its ground coordinates (X res , Y res , Z res );

[0103] ⑥ Based on the undistorted RFM orientation results of other sequence data of the temporal image in ③, calculate (Xres , Y res , Z res ) For each pixel coordinate (x, y) of the corresponding image sequences, calculate the gray value of (x, y) and assign it to (x res , y res );

[0104] ⑦ Repeat ⑤ - ⑥ to generate the time - series image registration result image, and obtain the one - to - one correspondence (x, y), (X base , y base ) between each pixel of the time - series frame image and the reference image, that is, the precise relative position relationship of the time - series images.

[0105] S2. Estimation of the ground truth of time - series images based on the geometric - radiation characteristics constraint between images: Considering the geometric - radiation characteristics constraint between time - series images, based on the Bayesian theory prior - probability estimation method, realize the estimation of the ground truth DN t1ue (x, y). The specific steps are as follows:

[0106] 2.1. Establish the model relationship between the ground - truth image X and the image Y obtained by the remote - sensing payload imaging as follows:

[0107] Y i = D i B i F i X + n i (11)

[0108] In the formula, i is the time - series image serial number (i = 1, 2,..., N), D is the motion - parameter matrix, B is the remote - sensing payload diffusion - function matrix, F is the resampling matrix, and n is the noise matrix.

[0109] 2.2. Obtain X in the above formula. According to the probability - estimation theory, under the above prior - knowledge constraints, the objective function for estimating the ground - truth image can be written in the following form:

[0110]

[0111] In the formula, i is the time - series image - frame serial number (i = 1, 2,..., N), Y i is the i - th frame image, X is the ground - truth object image. Therefore, the first part of the function ‖Y i - D i B i F i X‖ p reflects the brightness consistency between the time - series image and the estimated ground - truth image, and p is the norm order; the second part ‖Y i - (D i - D i+1 )Y i+1 ‖ qIt reflects the brightness consistency between adjacent temporal images, where q is the order of the norm; the third part is the regularization term R(X) and the regularization parameter λ.

[0112] 2.3. If the motion parameter matrix D between temporal images selects the M-estimation in robust estimation as the first term of the objective function, then:

[0113]

[0114] In the formula, N is the number of temporal images, M is the number of detector elements of the area array payload, and e i,n is the residual between the digital quantization value (Digital Number, DN) of the nth pixel in the i-th image and the DN value of the nth pixel in the estimated image. ρ(e) is the loss function. e is a fuzzy value, which is symmetric, greater than or equal to 0, and has monotonicity when e is greater than or less than 0; j i (n) is the serial number of the detector element that images the nth pixel in the i-th image, and DN i (j i (n)) is the DN value of the i-th imaging of the j i (n)th detector element of the area array payload, and DN True (n) is the DN value of the nth pixel in the estimated image.

[0115] 2.4. To take into account small errors while reducing the impact of large errors, the loss function ρ(e) is defined as follows:

[0116]

[0117] In the formula, the parameter a is automatically determined by the magnitude of the residual error. Affected by factors such as observation conditions, alignment errors, geometric distortions, and noise level differences, the residuals on different sequences of images are different. Based on the contribution of each error, the size parameter a of each residual of the temporal images is determined based on the following formula:

[0118]

[0119]

[0120] In the formula, E i represents the average residual on the i-th image, and Max and Min are the maximum and minimum value functions respectively.

[0121] 2.5. Selecting bilateral filtering as the regularization term, the final objective function can become:

[0122]

[0123] In the formula, w is the weight function, P is the kernel parameter of the one-dimensional bilateral filter, and the filter kernel size is 2 / P + 1. They respectively represent a translation of l pixels in the horizontal direction and a translation of h pixels in the vertical direction.

[0124] 2.6, and finally the ground truth is estimated to obtain the ground truth DN true (x, y).

[0125] S3. Based on the motion trajectories of the area array detector elements in the time-series images restored in Step 1, the gray values DN of each detector element in each frame sequence are obtained i , i ∈ [1, M], where M is the number of area array detector elements; taking a certain detector element D of the area array as an example, the method for extracting the gray value of the area array detector element in each frame sequence is described as follows:

[0126] 3.1. It is known that the coordinates of the area array detector element D on the reference image are (X base_D , Y base_D );

[0127] 3.2. The coordinates (x D , y D ) of the detector element D on a certain frame image can be known from Step 1.

[0128] (x D , y D ) = f(X base_D , Y base_D ) (18)

[0129] f refers to the functional relationship, that is, the function of the precise relative position relationship of the time-series images;

[0130] 3.3. According to the coordinates (x D , y D ) of the detector element D on a certain frame image, the gray value DN can be extracted D .

[0131] S4. Based on the gray values of the area array detector elements and the ground truth obtained in Step 2 and Step 3, a radiation calibration model for each detector element of the area array payload is constructed as follows:

[0132] DN true (j) = a i × DN i,j + b i (19)

[0133] In the formula, a i , b i are the relative radiation calibration coefficients of the i-th detector element of the area array payload, j is the frame sequence number, y ≤ N, and N is the number of frame sequences.

[0134] S5. Using the least squares method to solve the formula in Step 4 for all detector elements of the area array payload, that is, the relative radiation calibration coefficients of each detector element of the area array payload are obtained, and the on-orbit relative radiation calibration of the area array type remote sensing satellite is completed.

[0135] In specific implementation, the above process can be automatically run by using computer software technology. The hardware for running and implementing this method should also be within the protection scope of the present invention.

[0136] On the other hand, the embodiment of the present invention also provides a relative on-orbit radiometric calibration system for remote sensing satellites considering geometric radiation characteristic constraints, including:

[0137] A processor and a memory, where the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the relative on-orbit radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints as described in the above technical solution.

[0138] In a specific embodiment, 51 or more frame sequences of data are used for relative radiometric calibration experiments. The alignment accuracy of the temporal images is completed based on geometric information constraints and measured by the matching accuracy of homologous points between frames. In the verification example, the maximum matching error of homologous points in the frame sequence does not exceed 0.4 pixels, and the inter-frame matching accuracy is better than 0.1 pixel (1σ), as shown in Table 1.

[0139] Table 1 Matching accuracy of homologous points in multi-frame sequences

[0140]

[0141] After completing the estimation of the ground truth of the temporal images, the ground truth DN true image is obtained, as Figure 1 shown. A radiometric calibration model for each detector element of the area array payload is constructed, as Figure 2 shown. The calibration coefficients are applied to the area array image, and the radiometrically corrected image is evaluated with the image signal-to-noise ratio (SNR) as the objective evaluation index. The objective evaluation results are shown in Table 2, and the visual correction effect is as Figure 3 shown.

[0142] Table 2 Statistics of SNR improvement of the image after radiometric correction

[0143]

[0144] The specific embodiments described herein are merely illustrative of the spirit of the present invention. Those skilled in the technical field to which the present invention pertains can make various modifications or supplements to the described specific embodiments or use similar methods for substitution, but will not deviate from the spirit of the present invention or exceed the scope defined by the appended claims.

Claims

1. A method for on-orbit relative radiometric calibration of remote sensing satellites considering geometric radiation characteristic constraints, characterized in that It includes the following steps: Step 1: Complete the alignment of sequential images based on geometric information constraints; Step 2: Realize the estimation of the true surface value of sequential images based on the geometric and radiometric characteristic constraints between images and the prior probability estimation method of Bayesian theory Step 3: Based on Step 1, restore the motion trajectories of the area array detector elements in sequential images, and obtain the gray values of each detector element in each frame sequence; Step 4: Based on the gray values of the area array detector elements and the true surface values obtained in Step 2 and Step 3, construct the radiometric calibration model for each detector element of the area array payload; Step 5: Solve the radiometric calibration model in Step 4 for all detector elements of the area array payload, obtain the relative radiometric calibration coefficients for each detector element of the area array payload, and complete the on-orbit relative radiometric calibration of the area array remote sensing satellite.

2. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints according to claim 1, characterized in that: In Step 1, geometric calibration is used to complete system error compensation. Based on the calibrated interior orientation element supplementary model, a distortion-free RFM model is generated to eliminate intra-scene geometric distortion and inter-scene deformation in the time-series image data, achieve time-series image alignment, and obtain the one-to-one correspondence between the time-series image coordinates (x, y) and the reference image coordinates (X base , Y base ).

3. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints as claimed in claim 1 or 2, wherein: The specific implementation method of Step 1 is as follows: 1.1 System error compensation: The rational polynomial model (RFM) establishes the conversion relationship between the image coordinates (x, y) and the ground coordinates (X, Y, Z) through geometric polynomials; 1.2 Geometric distortion elimination: Use the calibration internal orientation element error compensation model of the calibration area image, and then generate the distortion-free RFM based on the following steps: 1) Use the solved internal orientation element error compensation model to establish the conversion relationship between the image coordinates and the ground coordinates according to Step 1.1; 2) Calculate the control grid (x, y, X, Y, Z) for generating the distortion-free RFM according to the conversion relationship between the image coordinates and the ground coordinates; 3) Calculate the new RFM using the control grid generated in 2); 1.3 Alignment of sequential images: Based on the distortion-free RFM, adopt a high-precision registration method with frame-to-frame positioning consistency constraints; The registration steps are as follows: ① Use the SIFT algorithm to extract corresponding points (x, y), (X base , Y base ) between the sequential frame images and the reference image; ② Using the undistorted RFM of the reference image, calculate the corresponding ground coordinates (X, Y, Z) of the matched homologous points (X base , Y base ) on the reference image, and form control points (X base , Y base , X, Y, Z); ③ Use the control points in ② to solve the distortion-free RFM orientation parameters of the sequential frame images; ④ Generate a registered result image of the same size as the reference image; ⑤For any point (x res , y res ) in the resulting image, calculate its ground coordinates (X res , Y res , Z res ) using the undistorted RFM of the reference image; ⑥Based on the undistorted RFM orientation results of other sequence data in the time-series images in ③, calculate the pixel coordinates (x, y) of each image sequence corresponding to (X res , Y res , Z res ), calculate the gray value of (x, y), and assign it to (x res , y res ); ⑦ Repeat steps ⑤ - ⑥ to generate the time-series image registration result image, and obtain the one-to-one correspondence (x, y), (X base , Y base ) of each pixel between the time-series frame image and the reference image, which is the precise relative position relationship of the time-series images.

4. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints according to claim 3, characterized in that: The conversion relationship between the image coordinates (x, y) and the ground coordinates (X, Y, Z) in Step 1.1 is given by the following formula: Num x (X, Y, Z) = a0 + a1Z + a2Y + a3X + a4ZY + a5ZX + a6YX + a7Z 2 + a8Y 2 + a9X 2 + a 10 ZYX + a 11 Z 2 Y + a 12 Z 2 X + a 13 Y 2 Z + a 14 Y 2 X + a 15 ZX 2 + a 16 YX 2 + a 17 Z 3 + a 18 Y 3 + a 19 X 3 Den x (X,Y,Z) = b0 + b1Z + b2Y + b3X + b4ZY + b5ZX + b6YX + b7Z 2 + b8Y 2 + b9X 2 + b 10 ZYX + b 11 Z 2 Y + b 12 Z 2 X + b 13 Y 2 Z + b 14 Y 2 X + b 15 ZX 2 + b 16 YX 2 + b 17 Z 3 + b 1a Y 3 + b 19 X 3 Num y (X,Y,Z) = c0 + c1Z + c2Y + c3X + c4ZY + c5ZY + c6YX + c7Z 2 + c8Y 2 + c9X 2 + c 10 ZYX + c 11 Z 2 Y + c 12 Z 2 X + c 13 Y 2 Z + c 14 Y 2 X + c 15 ZX 2 + c 16 YX 2 + c 17 Z 3 + c 18 Y 3 + c 19 X 3 Den y (X,Y,Z) = d0 + d1Z + d2Y + d3X + d4ZY + d5ZX + d6YX + d7Z 2 + d8Y 2 + d9X 2 + d 10 ZYX + d 11 Z 2 Y + d 12 Z 2 X + d 13 Y 2 Z + d 14 Y 2 X + d 15 ZX 2 + d 16 YX 2 + d 17 Z 3 + d 18 Y 3 + d 19 X 3 where a i , b i , c i , d i are polynomial coefficients; the RFM geometric positioning error comes from the satellite imaging geometric parameters, i.e., the positioning error (x e , y e ) caused by the exterior orientation elements and the positioning error (x i , y i ) caused by the interior orientation elements, and there is: Among them, Num x , Den x , Num y , Den y are intermediate quantities.

5. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints according to claim 1, wherein: The specific implementation of Step 2 includes the following sub-steps: 2.1 Establish the model relationship between the true surface image X and the image Y obtained by the remote sensing payload imaging as follows: Y i = D i B i F i X + n i where i is the sequential image number, i = 1, 2,..., N, D is the motion parameter matrix, B is the remote sensing payload diffusion function matrix, F is the resampling matrix, and n is the noise matrix; 2.2 Solve for X in the above formula. According to the probability estimation theory, under the constraint of prior knowledge, the objective function for estimating the true image of the ground object can be written in the following form: where i is the serial number of the sequential image frame, and Y i is the i-th frame image. Therefore, the first part of the function ‖Y i - D i B i F i X‖ p reflects the luminance consistency between the sequential image and the estimated true image, and p is the order of the norm; the second part ‖Y i -(D i - D i+1 )Y i+1 ‖ q reflects the luminance consistency between adjacent sequential images, and q is the order of the norm; the third part is the regularization term R(X) and the regularization parameter λ; 2.3 Select the M-estimation in robust estimation as the first term of the objective function for the motion parameter matrix D between sequential images, then: where N is the number of sequential images, M is the number of detector elements of the area array payload, and e i,n is the residual between the digital quantization value DN of the nth pixel in the ith image and the DN value of the nth pixel in the estimated image. ρ(e) is the loss function, which is symmetric, greater than or equal to 0, and monotonic when e is greater than or less than 0; j i (n) is the detector element number that images the nth pixel in the ith image, DN i (j i (n)) is the DN value of the j i (n)th detector element for the ith image of the area array payload, DN True (n) is the DN value of the nth pixel in the estimated image; 2.4 ρ(e) is the loss function defined as follows: a is a parameter; 2.5 Select bilateral filtering as the regularization term, and the final objective function can be changed to: where \(w\) is the weight function, \(P\) is the kernel parameter of the one-dimensional bilateral filter, and the kernel size of the filter is \(2 / P + 1\). They represent translating \(l\) pixels in the horizontal direction and \(h\) pixels in the vertical direction, respectively. 2.

6. Finally, estimate the true surface value to obtain the true surface DN value true (x, y).

6. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints as claimed in claim 5, wherein: The calculation method of the parameter a is as follows: where E i represents the average residual on the i-th image, and Max and Min are the maximum and minimum value functions respectively.

7. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints as described in claim 1, characterized in that: In Step 3, the method for extracting the gray values of the area array detector elements in each frame sequence is as follows: 3.

1. Given that the coordinates of the planar array detection element D on the reference image are (X base_D , Y base_D ); 3.2, The coordinates (x D , y D ) of the exploration element D on a certain frame of image are obtained by step 1; 3.

3. According to the coordinates (x D , y D ) of the exploration element D on a certain frame of image, the gray value DN D can be extracted.

8. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints according to claim 1, wherein: In Step 4, the radiometric calibration model for each detector element of the area array payload is as follows: DN true (j) = a i × DN i,j + b i where a i , b i is the relative radiometric calibration coefficient of the i-th detector element of the area array load, j is the frame sequence number, y ≤ N, and N is the number of frame sequences.

9. The on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints as claimed in claim 1, characterized in that: In Step 5, the least squares method is used to solve the radiometric calibration model in Step 4.

10. A on-orbit relative radiometric calibration system for remote sensing satellites considering geometric radiation characteristic constraints, characterized in that, It includes: A processor and a memory, the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the on-orbit relative radiometric calibration method for remote sensing satellites considering geometric radiation characteristic constraints according to any one of claims 1-9.

Images