A terrestrial ecosystem carbon monitoring satellite laser-assisted regional network adjustment method

Abstract

The application discloses a kind of land ecosystem carbon monitoring satellite laser auxiliary area network adjustment methods, comprising: stereoscopic image connection point generation: the multi-level pyramid of forward-looking, rear view 19° panchromatic image is constructed, and layer-by-layer matching inspection is carried out to obtain in-scene and inter-scene connection points;Free network adjustment: with the connection point as control, gradually tighten threshold to remove gross error, obtain free network adjustment RPC irrelevant to absolute height;Laser control point preparation;Laser-stereoscopic joint adjustment: laser control point and stereoscopic image connection point are jointly used as observation value, and forward intersection is carried out to obtain corrected intersection control point;Final RPC determination: image RPC is corrected again using corrected intersection control point, and final adjustment parameter is obtained after iterative convergence.The application realizes the joint refinement of free network and laser absolute height by "stereoscopic connection point hierarchical optimization, laser height replacement and distance reciprocal weighted coupling correction", which significantly improves the geometric accuracy and reliability of area network adjustment.

Description

Technical Field

[0001] This invention relates to the fields of photogrammetry and remote sensing, and in particular to a satellite laser-assisted regional network adjustment method for monitoring carbon in terrestrial ecosystems. Background Technology

[0002] The Terrestrial Ecosystem Carbon Monitoring Satellite (Jumang) can acquire observation data from multiple perspectives and of different types. This invention relates to the satellite's 19° forward-looking panchromatic imagery and RPC parameters, 19° backward-looking panchromatic imagery and RPC parameters, laser footprint imagery, and auxiliary data. The resolution of the 19° forward-looking and 19° backward-looking panchromatic images is approximately 3.5 meters. The resolution of the laser footprint imagery is approximately 8 meters. Auxiliary data includes the geodetic coordinates of ground laser points and the RPC of the laser footprint imagery. This invention constructs a stereo model from the 19° forward-looking and 19° backward-looking panchromatic images. After regional network adjustment, the footprint imagery is correlated with the stereo imagery to obtain laser control points. Joint adjustment of the laser control points and stereo imagery improves the elevation accuracy of the regional network adjustment. Summary of the Invention

[0003] The purpose of this invention is to provide a satellite laser-assisted regional network adjustment method for monitoring carbon in terrestrial ecosystems, thereby solving the aforementioned problems in the prior art.

[0004] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A method for adjusting a satellite laser-assisted regional network for monitoring carbon in terrestrial ecosystems includes the following steps: S1. Generation of stereoscopic image connection points: Construct a multi-level pyramid for the forward-looking and backward-looking 19° panchromatic images, match layer by layer and perform reverse consistency checks to obtain the connection points between the scene and the scene. S2. Free network adjustment: Using the tie points as control, the image RPC is iteratively corrected, and the threshold is gradually tightened to eliminate gross errors, resulting in a free network adjustment RPC that is independent of the absolute elevation. S3. Laser control point preparation: Analyze satellite laser footprint data, extract local footprint images, and combine free network adjustment RPC to generate laser control points corresponding to the stereo image; S4. Laser-stereo joint adjustment: The laser control points and the connection points of the stereo image are used together as observations to perform forward intersection, and the intersection coordinates are corrected using the laser elevation as the true value to obtain the corrected intersection control points. S5. Final RPC Determination: The image RPC is corrected again using the corrected intersection control points. After iterative convergence, the final adjustment parameters that combine free network geometric consistency and laser absolute elevation accuracy are obtained.

[0005] Furthermore, step S1, generating the stereoscopic image connection points, specifically includes the following steps: S11. Downsample the 19° panchromatic images of both the forward and backward views using the mean DN value of a 3×3 rectangular window, until the image width and height are ≤64 pixels, forming an n-layer pyramid. Then, normalize the corresponding RPC row and column parameters according to... Synchronous scaling yields the top-level pyramid image and its top-level RPC; S12. Project the center pixel of the forward-looking top-level image onto the HEIGHT_OFF elevation surface through the top-level RPC to obtain the geodetic coordinates (L,B,H). Then, back-project the image points of the back-looking top-level RPC to obtain the image points of the back-looking image. Calculate the offsets (roff,coff) between the forward and backward image points: roff=rcB-rcF coff=ccB-ccF S13. Divide the top-level image of the front view into a 7×7 uniform grid. With the center of the grid as the initial point, take an 11×11 pixel window. Calculate the normalized cross-correlation coefficient pixel by pixel in the 7×7 search area of ​​the top-level image of the rear view. Select the point with the largest correlation coefficient as the forward matching point. Then, use this point as the search center to reverse match the front view image. If the distance between the reverse matching point and the initial point is greater than 1 pixel, it is marked as invalid. Otherwise, it is a valid matching point pair. S14. Multiply the coordinates of the effective matching point pairs by 3 and project them to the next layer of the pyramid. Insert one encrypted initial point between two adjacent points. Repeat the forward-backward matching until the original resolution is reached to obtain the in-scene connection points. S15. Calculate the approximate ground distance based on the LONG_OFF and LAT_OFF values ​​of the two image pairs' RPCs. If the distance is less than the threshold, perform projection of the four corner points to determine the overlapping area and correct the LINE_OFF and SAMP_OFF ​​values ​​of the corresponding RPCs, generating images 1-4. Perform layer-by-layer matching on images 1-4 using the same strategy as steps S11-S14, and retain only the four-fold matching point pairs that satisfy r1=r1' and r3=r3', finally obtaining the inter-scene connection points.

[0006] Furthermore, the specific steps for generating the free network adjustment in step S2 include: S21. Perform forward intersection on the in-scene and inter-scene connection points obtained in step S1 to obtain ground coordinates (Lbf, Bbf, Hbf), and combine the ground coordinates with the forward and backward image points in the corresponding connection points to form the intersection control points of each image. S22. Count all intersection control points of each image, and use their image point coordinates as observations to correct the "row and column rational polynomial numerator constants" and "normalized longitude and latitude linear terms" in the RPC of the image to obtain the first-order correction RPC; S23. Using a single-correction RPC, back-project the intersection control points onto the image plane and calculate the distance Δp between the projected image points and the original image points: Δp= ; Δp: Image plane distance, used to measure the geometric deviation between the image point after RPC backprojection and the original image point; rproj: The row coordinates of the "projected image point" obtained from the RPC back projection calculation; robs: Row coordinates of the "observed image point" obtained from the original matching or measurement; cproj: The column coordinates of the "projected image points" obtained from the RPC back projection calculation; cobs: The column coordinates of the "observed image points" obtained from the original matching or measurement; If Δp is greater than a threshold Merr, then mark the control point as invalid; otherwise, retain it. S24. Alternately execute steps S21-S23 a total of 5 times. In the first time, Ndd=5 and Merr=9 pixels. In the subsequent 4 times, Ndd=2 and Merr decreases to 7, 5, 3 and 1 pixels respectively. The final RPC is the free network adjustment RPC.

[0007] Furthermore, step S3, the preparation of the laser control point, specifically includes the following steps: S31. Parse the satellite laser h5 file to obtain the ground coordinates (L,B,H) of the laser point, and obtain the image plane coordinates (rld,cld) by orthographic projection of the footprint image RPC. S32. Extract a rectangular region with height Hld and width Wld centered at (rld,cld) to generate a segmented footprint image. Simultaneously correct the LINE_OFF and SAMP_OFF ​​of the footprint image RPC so that the coordinates of the top left corner image point of the segmented image are (0,0) and the coordinates of the laser point image point within the segmented image are (Hld / 2,Wld / 2). This forms the footprint observation data of the laser control point. S33. Project the same (L,B,H) image onto the RPC of the rearview image, and record the sequence number, row / column coordinates and laser point sequence number of all overlapping rearview images. Generate a laser control point with the same name for each overlap. The total number of control points is greater than or equal to the total number of laser points.

[0008] Furthermore, step S4, laser-stereo joint adjustment, includes: S41. Using the same laser control point obtained in step S33 as the center, crop out 1000×1000 pixel local images of the front and back views, and use free network adjustment RPC to calculate the local DSM of the area. S42. Project the laser ground coordinates (L,B,H) onto the ground using the backsight adjustment RPC to obtain the backsight image point (Rb,Cb); then use the backsight adjustment RPC and the local DSM to project back onto the ground to obtain the ground coordinates (L1,B1,H1). S43. Project (L1,B1,H1) onto the forward-looking image point (Rf,Cf) using the forward-looking adjustment RPC orthographic projection. Then replace H1 with the original laser elevation H to form the final ground coordinates (L1,B1,H). S44. A laser control point is formed by (L1,B1,H), the back view image point (Rb,Cb), and the front view image point (Rf,Cf). Repeat steps S41-S44 to complete the preparation of all laser control points.

[0009] Furthermore, step S5, the final RPC determination, specifically includes the following steps: S51. Using the free network adjustment RPC as the initial value, perform forward intersection on the connection points of the stereo image to obtain the ground coordinates (Lbf, Bbf, Hbf) of the intersection control point, and combine them with the corresponding forward and backward image points to form the intersection control point. S52. For each laser control point, intersect its forward image point (Rf, Cf) and backward image point (Rb, Cb) to obtain the intersection coordinates (Ljs, Bjs, Hjs), and calculate the difference between these coordinates and the laser coordinates (L1, B1, H). Loffi=L1-Ljs Boffi=B1-Bjs Hoffi=H-Hjs After multiplying Loffi and Boffi by 100000 to convert them to meters, the error threshold Merrld∈{10,8,6,4,2}m was successively tightened using RANSAC to eliminate gross errors, and the remaining points were the effective laser control points. S53. Calculate the correction amount for each intersection control point using a weighted average of the inverse distance: Lall=Add(Loffi / Di) / Add(1 / Di); Ball=Add(Boffi / Di) / Add(1 / Di); Hall=Add(Hoffi / Di) / Add(1 / Di); Where Add represents cumulative calculation, Di (i=0,…, Nldyx-1) is the ground distance between the intersection control point and the effective laser control point, and (Lbf, Bbf, Hbf) is added to (Lall, Ball, Hall) to obtain the corrected intersection control point; S54. Using the corrected intersection control points, the row and column numerator constants and the normalized latitude and longitude linear terms of the image RPC are corrected again. Steps S51-S54 are executed iteratively a total of 5 times. The RPC obtained in the last iteration is the final adjustment parameter.

[0010] Furthermore, step S1, generating the stereoscopic image connection points, also includes the following steps: S16. After downsampling each layer of the pyramid, wavelet detail injection is used to enhance the texture of the 3.5m panchromatic image, and the high-frequency components of the RPC normalization parameters are corrected simultaneously, specifically including: S161. Perform Daubechies-4 wavelet decomposition on the downsampled image to obtain the low-frequency subband LL and the high-frequency subbands LH, HL, and HH. S162. Multiply the high-frequency subband coefficients by the enhancement coefficient α = 0.8 to obtain the enhanced high-frequency subbands LH′, HL′, and HH′. S163. Perform inverse wavelet transform on LL, LH′, HL′, and HH′ to obtain texture-enhanced images. The texture information enhancement amount ΔT = |I_enhanced - I_original| / I_original, ΔT ≥ 12%; Where I_enhanced represents the image grayscale value after wavelet detail injection; I_original represents the original grayscale value after downsampling and without enhancement. S164. Synchronously correct the high-frequency components of the RPC normalization parameters: multiply the coefficients of the higher-order terms of the RPC row and column rational polynomials by the compensation factor β = 1 / (1+α) to maintain the geometric consistency between the enhanced image and the RPC and improve the cross-scale matching success rate; where β: is the compensation factor for the higher-order terms of the RPC; α: is the wavelet high-frequency enhancement coefficient, which is 0.8. S17. Introduce laser waveform gradient weights in the reverse consistency check, specifically including: S171. Calculate the gradient dI / dx along the direction of the footprint image to form the gradient map G; S172. Normalize the gradient graph G to obtain the waveform gradient weight w = |dI / dx| / Σ|dI / dx|, w∈[0,1]; Where dI / dx is the grayscale gradient of the footprint image along the row direction; Σ|dI / dx|: The sum of the absolute values ​​of the gradients of all pixels within the current search window, used for normalization; S173. Use w as the matching quality weight. If w < 0.3, mark the matching point as invalid to reduce false matching of vegetation edges. S174: Weights are applied only to the top of the pyramid, reducing computation by 30% while maintaining matching reliability.

[0011] Furthermore, in step S2, during the RPC correction stage, ridge estimation regularization is adopted, with a regularization coefficient λ = 0.005 × (base-to-height ratio)², to suppress the ill-conditioned normal equation caused by the ultra-large base-to-height ratio of 19°+19°. In step S3, the 8m footprint image is reconstructed with a super-resolution of 0.5m. The laser waveform half-width at half-maximum is used as the reconstruction constraint to improve the effective resolution of the footprint image to 2.5m, and then it is matched with the 3.5m stereo image.

[0012] Furthermore, during local DSM matching, a finer grid of 0.3 pixels × 0.3 pixels is used for mountainous areas with a slope > 15°, while a grid of 0.5 pixels × 0.5 pixels is maintained for plains areas to reduce the failure rate of matching in mountainous areas. In the joint adjustment, the laser control points are weighted according to NDVI: the weight is multiplied by 0.7 when NDVI>0.6 and by 1.0 when NDVI≤0.6, in order to suppress forest canopy bias. In the distance reciprocal weighted correction formula, an elevation accuracy factor σH=|Hcanopy-Hground| / 2 is introduced, with a weight of 1 / (Di×σH_i), to reduce the vegetation-ground double elevation unmixing error; Where Hcanopy is the elevation value of the first peak of the laser echo; Hground represents the elevation value of the final peak of the laser echo; |Hcanopy−Hground| represents the elevation difference between the two peaks, often referred to as the "canopy-ground elevation difference" or "laser beamwidth"; σH_i is the σH value corresponding to the i-th effective laser control point; Di is the horizontal distance on the ground from the intersection control point to the laser point.

[0013] Furthermore, this method is only enabled when jointly processing 19° forward / backward stereo images from the terrestrial ecosystem carbon monitoring satellite with 8m laser footprint data to form a dedicated algorithmic closed loop.

[0014] The beneficial effects of this invention are: This invention provides a method for adjusting a satellite laser-assisted regional network for terrestrial ecosystem carbon monitoring. It divides stereo image tie points into intra-scene and inter-scene tie points, and obtains reliable stereo tie points through pyramid hierarchical matching, progressive refinement, and reverse matching. Through multiple iterative calculations and progressively reducing the error threshold, gross errors in tie points are eliminated, achieving free network adjustment. Corresponding stereo image points to laser points are obtained through local DSM extraction and point projection calculation in the laser point's region, thereby acquiring laser control points and avoiding direct matching failures. By modifying the elevation in the intersection coordinates of the corresponding stereo image points to the laser point's elevation, the reliable elevation accuracy of laser points is utilized, avoiding the instability of laser point planar accuracy. Based on the distribution of object-space residual errors of laser points, multiple iterations and progressive reduction of the gross error threshold are performed. Gross errors in laser points are eliminated using the object-space RANSAC algorithm, ensuring adjustment accuracy. By modifying the intersection coordinates of tie points using a distance-inverse weighted calculation of the laser point object-space coordinate correction, absolute accuracy can be improved with the help of laser control points without compromising the accuracy of the free network. Attached Figure Description

[0015] Figure 1 This is a flowchart of a satellite laser-assisted regional network adjustment method for terrestrial ecosystem carbon monitoring according to the present invention. Detailed Implementation

[0016] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0017] Reference Figure 1 The method for adjusting a satellite laser-assisted regional network for monitoring carbon in terrestrial ecosystems, as shown, includes the following steps: S1. Generation of stereoscopic image connection points: Construct a multi-level pyramid for the forward-looking and backward-looking 19° panchromatic images, match layer by layer and perform reverse consistency checks to obtain the connection points between the scene and the scene. S2. Free network adjustment: Using the tie points as control, the image RPC is iteratively corrected, and the threshold is gradually tightened to eliminate gross errors, resulting in a free network adjustment RPC that is independent of the absolute elevation. S3. Laser control point preparation: Analyze satellite laser footprint data, extract local footprint images, and combine free network adjustment RPC to generate laser control points corresponding to the stereo image; S4. Laser-stereo joint adjustment: The laser control points and the connection points of the stereo image are used together as observations to perform forward intersection, and the intersection coordinates are corrected using the laser elevation as the true value to obtain the corrected intersection control points. S5. Final RPC Determination: The image RPC is corrected again using the corrected intersection control points. After iterative convergence, the final adjustment parameters that combine free network geometric consistency and laser absolute elevation accuracy are obtained.

[0018] Furthermore, step S1, generating the stereoscopic image connection points, specifically includes the following steps: S11. Downsample the 19° panchromatic images of both the forward and backward views using the mean DN value of a 3×3 rectangular window, until the image width and height are ≤64 pixels, forming an n-layer pyramid. Then, normalize the corresponding RPC row and column parameters according to... Synchronous scaling yields the top-level pyramid image and its top-level RPC.

[0019] This step specifically involves: reducing the height and width of the 19-degree panchromatic image to one-third of its original size by averaging the DN values ​​within a 3×3 rectangular window, thus obtaining the first-layer pyramid image. This process is repeated iteratively until both the height and width of the pyramid image are less than or equal to 64. Assume this process has been repeated n times, resulting in n layers of pyramid images. At this point, the normalized row and column coordinate parameters of the RPC corresponding to the smallest pyramid image (called the top-level pyramid image) are 1.0 / 3^n times the parameters corresponding to the RPC of the 19-degree panchromatic image.

[0020] The same method was used to obtain the pyramid images of each layer and the top RPC of the 19-degree panchromatic rearview image.

[0021] Hereinafter, the forward-looking 19-degree panchromatic image is referred to as the forward-looking image, and the rear-looking 19-degree panchromatic image is referred to as the rear-looking image.

[0022] S12. Project the center pixel of the forward-looking top-level image onto the HEIGHT_OFF elevation surface through the top-level RPC to obtain the geodetic coordinates (L,B,H). Then, back-project the image points of the back-looking top-level RPC to obtain the image points of the back-looking image. Calculate the offsets (roff,coff) between the forward and backward image points: roff=rcB-rcF coff=ccB-ccF; The specific steps are as follows: The center coordinates (rcF, ccf) of the foreseeable top pyramid image are projected onto an elevation surface with elevation H equal to the HEIGHT_OFF parameter of the top pyramid image via the foreseeable top RPC, yielding geodetic coordinates (L, B, H). These geodetic coordinates (L, B, H) are then projected onto the backseeable top pyramid image via the backseeable top RPC, yielding coordinates (rcB, ccB). The offset of the image points in the backseeable top pyramid image relative to the image points in the foreseeable top pyramid image is (roff, coff): roff=rcB-rcF coff = ccB - ccF.

[0023] S13. Divide the top-level image of the front view into a 7×7 uniform grid. With the center of the grid as the initial point, take an 11×11 pixel window. Calculate the normalized cross-correlation coefficient pixel by pixel in the 7×7 search area of ​​the top-level image of the rear view. Select the point with the largest correlation coefficient as the forward matching point. Then, use this point as the search center to reverse match the front view image. If the distance between the reverse matching point and the initial point is greater than 1 pixel, it is marked as invalid. Otherwise, it is a valid matching point pair. The specific steps are as follows: Divide the top-level pyramid image in the front view into 7×7 uniform grids, and take the center point of each grid as the initial point for coordinate matching (rsF_i_j, csF_i_j) (i=0,…,6; j=0,…,6). The initial point corresponding to the pyramid image in the back view is (rsB_i_j, csB_i_j), then: rsB_i_j = rsB_i_j + roff csB_i_j = csB_i_j + coff For each initial point in the forward-looking top-level pyramid image, an 11×11 window of DN values ​​is taken. The corresponding initial point in the backward-looking top-level pyramid image is used as the search center. Within a 7×7 search window at this center, the corresponding 11×11 window DN values ​​are taken pixel by pixel. The correlation coefficient between the forward-looking and backward-looking 11×11 window DN values ​​is calculated, resulting in 7×7 correlation coefficients. The point in the backward-looking top-level pyramid image with the largest correlation coefficient is taken as the optimal matching point (hereinafter referred to as the backward-looking matching point). This process is called forward matching. Using the forward-looking search center point as the search center, the corresponding 11×11 window DN values ​​are taken pixel by pixel within a 7×7 search window at this center. The correlation coefficient between the backward-looking and forward-looking 11×11 window DN values ​​is calculated, and the optimal matching point (hereinafter referred to as the forward-looking reverse matching point) is calculated. This process is called reverse matching. If the row and column coordinates of the forward-looking reverse matching point are more than one pixel away from the row and column coordinates of the forward-looking search center point, the point is marked as an invalid matching point. Otherwise, it is marked as a valid matching point. This will create a maximum of 7×7 matching point pairs between the forward-looking top-level pyramid image and the rear-looking top-level pyramid image.

[0024] S14. Multiply the coordinates of the effective matching point pairs by 3 and project them to the next layer of the pyramid. Insert one encrypted initial point between two adjacent points. Repeat the forward-backward matching until the original resolution is reached to obtain the in-scene connection points. The specific steps are as follows: Valid matching point pairs from the front-view and back-view top-level pyramid images are projected onto the next pyramid level by multiplying the coordinates of the matching point pairs by 3. An encrypted initial point is interpolated between two adjacent projected points to obtain the encrypted initial point for searching the front-view and back-view pyramid images. Matching point pairs for this pyramid level are obtained through forward and reverse matching of the pyramid images.

[0025] Repeat the above point encryption and forward and reverse matching process until matching point pairs of the original resolution layer are obtained.

[0026] S15. Calculate the approximate ground distance based on the LONG_OFF and LAT_OFF values ​​of the two image pairs' RPCs. If the distance is less than the threshold, perform projection of the four corner points to determine the overlapping area and correct the LINE_OFF and SAMP_OFF ​​values ​​of the corresponding RPCs, generating images 1-4. Perform layer-by-layer matching on images 1-4 using the same strategy as steps S11-S14, and retain only the four-fold matching point pairs that satisfy r1=r1' and r3=r3', finally obtaining the inter-scene connection points.

[0027] The specific operation of this step is as follows: First, calculate the approximate ground distance between the two images based on the LONG_OFF and LAT_OFF parameters of the RPC of the two forward-looking panchromatic images. If the distance in the longitude and latitude directions is less than a certain threshold, the overlapping area is calculated; otherwise, it is considered that the two images do not overlap.

[0028] The image plane coordinates of the four corner points of the foreseeable panchromatic image (hereinafter referred to as the master image) in one image pair are projected onto the foreseeable panchromatic image (hereinafter referred to as the projected image) in another image pair using the elevation plane of the master image RPC with the HEIGHT_OFF parameter, resulting in four projected image plane coordinates. If the row or column coordinate of the four projected image plane coordinates is less than 0, it is set to 0; if the row coordinate of the four projected image plane coordinates is greater than the number of rows of the projected image Ht, it is set to Ht; if the column coordinate of the four projected image plane coordinates is greater than the number of columns of the projected image Wt, it is set to Wt. Then, the range of the four defined ranges of the projected image plane coordinates (hereinafter referred to as the projected image range) is calculated.

[0029] Swap the order of the master image and the projected image, i.e., treat the projected image as the master image and the master image as the projected image. Calculate the image plane coordinate range of the projected image.

[0030] The same method is used to determine the overlap range of the main image and the rear view image corresponding to the projected image.

[0031] Assuming the top-left row and column coordinates of the overlapping area of ​​an image are (ro, co), subtracting ro from the LINE_OFF parameter and co from the SAMP_OFF ​​parameter of this image's RPC yields the RPC of the overlapping area. This results in the overlapping area images and corresponding RPC parameters of the main image, the corresponding rear-view image, the projected image, and the corresponding rear-view image. These four overlapping area images are hereinafter referred to as Image 1, Image 2, Image 3, and Image 4.

[0032] Based on method 1.1, pyramid images corresponding to the back-view images of images 1, 2, 3, and 4 are calculated. Based on method 1.2, matching point pairs between images 1 and 3 (denoted as (r1, c1), (r3, c3)), between images 1 and 2 (denoted as (r1', c1'), (r2, c2)), and between images 3 and 4 (denoted as (r3', c3'), (r4, c4)) are calculated. Only valid matching point pairs that meet the condition (r1 equals r1' and r3 equals r3') are retained at each pyramid level and passed to the next level. This yields fourfold connectivity points between two image pairs.

[0033] Furthermore, the specific steps for generating the free network adjustment in step S2 include: S21. Perform forward intersection on the in-scene and inter-scene connection points obtained in step S1 to obtain ground coordinates (Lbf, Bbf, Hbf), and combine the ground coordinates with the forward and backward image points in the corresponding connection points to form the intersection control points of each image. The specific operation of this step is as follows: The ground coordinates (Lbf, Bbf, Hbf) are obtained by intersecting the front and rear view stereoscopic image connection points (including intra-scene connection points and inter-scene connection points). If it is an intra-scene connection point, the ground coordinates and the coordinates of the rear view image points in the corresponding connection point constitute a rear view image control point. At the same time, the ground coordinates and the coordinates of the front view image points in the corresponding connection point constitute a front view image control point. If it is an inter-scene connection point, the ground coordinates and all the image points corresponding to the connection point constitute a corresponding control point, which is called the intersection control point.

[0034] S22. Count all intersection control points of each image, and use their image point coordinates as observations to correct the "row and column rational polynomial numerator constants" and "normalized longitude and latitude linear terms" in the RPC of the image to obtain the first-order correction RPC; The specific steps in this process are as follows: Calculate the control points corresponding to all connection points, and statistically analyze all control points corresponding to each image (forward-looking or backward-looking image). Based on the control points of an image, correct the row and column constant terms of the rational polynomial numerators and the normalized longitude and latitude linear terms in the RPC parameters of that image. This yields the corrected RPC parameters.

[0035] S23. Using a single-correction RPC, back-project the intersection control points onto the image plane and calculate the distance Δp between the projected image points and the original image points: Δp= ; Δp: Image plane distance, used to measure the geometric deviation between the image point after RPC backprojection and the original image point; rproj: The row coordinates of the "projected image point" obtained from the RPC back projection calculation; robs: Row coordinates of the "observed image point" obtained from the original matching or measurement; cproj: The column coordinates of the "projected image points" obtained from the RPC back projection calculation; cobs: The column coordinates of the "observed image points" obtained from the original matching or measurement; If Δp is greater than a threshold Merr, then mark the control point as invalid; otherwise, retain it. The specific operation of this step is as follows: Steps S21 and S22 are executed alternately for a total of Ndd times. The corrected RPC of the intersection control point is projected onto the forward and backward images to obtain the coordinates of the projected image point. The image plane distance between the coordinates of the projected image point and the coordinates of the image point recorded by the control point is calculated. If the distance is greater than a threshold Merr, the point is set as an invalid point; otherwise, it is a valid point.

[0036] S24. Alternately execute steps S21-S23 a total of 5 times. In the first time, Ndd=5 and Merr=9 pixels. In the subsequent 4 times, Ndd=2 and Merr decreases to 7, 5, 3 and 1 pixels respectively. The final RPC is the free network adjustment RPC.

[0037] The specific steps are as follows: First, execute step S23, setting Ndd to 5 and Merr to 9 pixels. Then, repeat step S23 four times, setting Ndd to 2 and Merr to 7, 5, 3, and 1 respectively each time. The final corrected RPC is the RPC after free network adjustment.

[0038] Furthermore, step S3, the preparation of the laser control point, specifically includes the following steps: S31. Parse the satellite laser h5 file to obtain the ground coordinates (L,B,H) of the laser point, and obtain the image plane coordinates (rld,cld) by orthographic projection of the footprint image RPC. The specific steps are as follows: The laser data from the terrestrial ecosystem carbon monitoring satellite includes a footprint image in TIFF format, a corresponding RPC file for the footprint image, and an h5 file containing the ground coordinates of the laser points. First, the ground coordinates (L, B, H) of the laser points are parsed from the h5 file according to the file description. The (L, B, H) coordinates are then projected onto the image plane using the corresponding footprint image RPC to obtain the row and column coordinates (rld, cld). The local area of ​​the footprint image corresponding to this laser point is set, namely, the height Hld and the width Wld. A rectangular area image with a starting position of (rld-Hld / 2, cld-Wld / 2) from the top left corner, and a height and width of (Hld, Wld) is read and saved as a separate TIFF image, called a segmented footprint image. The RPC obtained by subtracting (rld-Hld / 2) from the LINE_OFF and (cld-Wld / 2) from the SAMP_OFF ​​of the entire footprint image RPC is the RPC of this segmented footprint image. The complete data for a laser control point consists of a segmented footprint image, a segmented footprint image RPC, ground coordinates (L, B, H), and the image point coordinates (Hld / 2, Wld / 2) of the laser point on the segmented footprint image. The complete data for each laser point is obtained sequentially using the same method.

[0039] S32. Extract a rectangular region with height Hld and width Wld centered at (rld,cld) to generate a segmented footprint image. Simultaneously correct the LINE_OFF and SAMP_OFF ​​of the footprint image RPC so that the coordinates of the top left corner image point of the segmented image are (0,0) and the coordinates of the laser point image point within the segmented image are (Hld / 2,Wld / 2). This forms the footprint observation data of the laser control point. S33. Project the same (L,B,H) image onto the RPC of the rearview image, and record the sequence number, row / column coordinates and laser point sequence number of all overlapping rearview images. Generate a laser control point with the same name for each overlap. The total number of control points is greater than or equal to the total number of laser points.

[0040] The specific operation of this step is as follows: Project the ground coordinates of the laser point onto the rear-view image, and record the corresponding image sequence number, row / column coordinates, and the sequence number of the laser point. Overlapping rear-view images are recorded separately. That is, a laser point and its corresponding point on a pair of forward and rear-view images constitute a laser control point. Several overlapping image pairs in the area corresponding to a laser point constitute several laser control points with the same ground coordinates. The number of laser control points is greater than or equal to the number of laser points.

[0041] Furthermore, step S4, laser-stereo joint adjustment, includes: S41. Using the same laser control point obtained in step S33 as the center, crop out 1000×1000 pixel local images of the front and back views, and use free network adjustment RPC to calculate the local DSM of the area. The local DSM is calculated based on the front and back view local images (both 1000 pixels in height and width) corresponding to the laser points and the RPC after free mesh adjustment (hereinafter referred to as the adjusted RPC).

[0042] S42. Project the laser ground coordinates (L,B,H) onto the ground using the backsight adjustment RPC to obtain the backsight image point (Rb,Cb); then use the backsight adjustment RPC and the local DSM to project back onto the ground to obtain the ground coordinates (L1,B1,H1). S43. Project (L1,B1,H1) onto the forward-looking image point (Rf,Cf) using the forward-looking adjustment RPC orthographic projection. Then replace H1 with the original laser elevation H to form the final ground coordinates (L1,B1,H). S44. A laser control point is formed by (L1,B1,H), the back view image point (Rb,Cb), and the front view image point (Rf,Cf). Repeat steps S41-S44 to complete the preparation of all laser control points.

[0043] The ground coordinates (L, B, H) of the laser point are used to calculate the back-view image point coordinates (Rb, Cb) using back-view image adjustment RPC. The back-view image point coordinates (Rb, Cb) are then used to calculate the ground coordinates (L1, B1, H1) using back-view image adjustment RPC and local DSM. The ground coordinates (L1, B1, H1) are then used to calculate the forward-view image point coordinates (Rf, Cf) using forward-view image adjustment RPC. (L1, B1, H1) is then modified to (L1, B1, H). (L1, B1, H), along with (Rb, Cb) and (Rf, Cf), constitute a laser control point. All laser control points are calculated based on the above method.

[0044] Furthermore, step S5, the final RPC determination, specifically includes the following steps: S51. Using the free network adjustment RPC as the initial value, perform forward intersection on the connection points of the stereo image to obtain the ground coordinates (Lbf, Bbf, Hbf) of the intersection control point, and combine them with the corresponding forward and backward image points to form the intersection control point. The adjusted RPC parameters calculated in step S2 are used as the initial values ​​of the RPC for the front and rear view stereo images. The ground coordinates (Lbf, Bbf, Hbf) are obtained by intersecting the front and rear view stereo image connection points (including intra-scene connection points and inter-scene connection points). If it is an intra-scene connection point, the ground coordinates and the coordinates of the rear view image points in the corresponding connection point constitute a rear view image control point. At the same time, the ground coordinates and the coordinates of the front view image points in the corresponding connection point constitute a front view image control point. If it is an inter-scene connection point, the ground coordinates and all the image points corresponding to the connection point constitute a corresponding intersection control point.

[0045] S52. For each laser control point, intersect its forward image point (Rf, Cf) and backward image point (Rb, Cb) to obtain the intersection coordinates (Ljs, Bjs, Hjs), and calculate the difference between these coordinates and the laser coordinates (L1, B1, H). Loffi=L1-Ljs Boffi=B1-Bjs Hoffi=H-Hjs After multiplying Loffi and Boffi by 100000 to convert them to meters, the error threshold Merrld∈{10,8,6,4,2}m was successively tightened using RANSAC to eliminate gross errors, and the remaining points were the effective laser control points. The specific operation of this step is as follows: Assuming there are Nld laser control points, the ground coordinates (Ljs, Bjs, Hjs) are obtained by intersecting the front and rear view stereo image points corresponding to the laser control points. The difference (Loffi, Boffi, Hoffi) (i=0,…,Nld-1) between these ground coordinates and the ground coordinates (L1, B1, H) of the laser control points calculated in step S4.2 is then calculated: Loffi=L1-Ljs Boffi=B1-Bjs Hoffi=H-Hjs Set an error threshold Merrld in meters. Since Loffi and Boffi are in degrees, multiply Loffi and Boffi by 100,000 to approximate them in meters. This is used for gross error removal; that is, the RANSAC algorithm is used to remove laser control points that do not meet the error threshold. After removing gross errors, Loffi and Boffi are restored to their original values. The remaining laser control points are called effective laser control points, and their number is assumed to be Nldyx.

[0046] S53. Calculate the correction amount for each intersection control point using a weighted average of the inverse distance: Lall=Add(Loffi / Di) / Add(1 / Di); Ball=Add(Boffi / Di) / Add(1 / Di); Hall=Add(Hoffi / Di) / Add(1 / Di); Where Add represents cumulative calculation, Di (i=0,…, Nldyx-1) is the ground distance between the intersection control point and the effective laser control point, and (Lbf, Bbf, Hbf) is added to (Lall, Ball, Hall) to obtain the corrected intersection control point; The specific steps in this process are as follows: Calculate the distance Di (i=0,…, Nldyx-1) between the ground coordinates of the intersection control points (Lbf, Bbf, Hbf) corresponding to the connection point and the ground coordinates of all valid laser control points. Calculate the object-space coordinate corrections (Lall, Ball, Hall) for the ground coordinates of the intersection control points (Lbf, Bbf, Hbf): Lall=Add(Loffi / Di) / Add(1 / Di) Ball=Add(Boffi / Di) / Add(1 / Di) Hall=Add(Hoffi / Di) / Add(1 / Di) Where Add represents cumulative calculation, i=0,…, Nldyx-1, which means that the correction of the ground coordinates of the intersection control point is calculated by weighting the distances between the ground coordinates of the intersection control point and the ground coordinates of all valid laser control points by the reciprocal.

[0047] The corrected ground coordinates of the intersection control points are obtained by adding their object coordinate corrections to the ground coordinates of all intersection control points. The corresponding intersection control points are called corrected intersection control points.

[0048] S54. Using the corrected intersection control points, the row and column numerator constants and the normalized latitude and longitude linear terms of the image RPC are corrected again. Steps S51-S54 are executed iteratively a total of 5 times. The RPC obtained in the last iteration is the final adjustment parameter.

[0049] For each image, the corrected intersection control points are statistically analyzed. Based on these corrected intersection control points, the row and column constant terms of the rational polynomials in the RPC parameters of that image are corrected, along with the normalized longitude and latitude linear terms. The corrected RPC parameters are then obtained.

[0050] The iteration steps S51-S54 are executed 5 times, with Merrld set to 10, 8, 6, 4, and 2 respectively. The RPC after the last modification is the final adjustment RPC.

[0051] Furthermore, step S1, generating the stereoscopic image connection points, also includes the following steps: S16. After downsampling each layer of the pyramid, wavelet detail injection is used to enhance the texture of the 3.5m panchromatic image, and the high-frequency components of the RPC normalization parameters are corrected simultaneously, specifically including: S161. Perform Daubechies-4 wavelet decomposition on the downsampled image to obtain the low-frequency subband LL and the high-frequency subbands LH, HL, and HH. S162. Multiply the high-frequency subband coefficients by the enhancement coefficient α = 0.8 to obtain the enhanced high-frequency subbands LH′, HL′, and HH′. S163. Perform inverse wavelet transform on LL, LH′, HL′, and HH′ to obtain texture-enhanced images. The texture information enhancement amount ΔT = |I_enhanced−I_original| / I_original, ΔT ≥ 12%; Where I_enhanced represents the image grayscale value after wavelet detail injection; I_original represents the original grayscale value after downsampling and without enhancement. S164. Synchronously correct the high-frequency components of the RPC normalization parameters: multiply the coefficients of the higher-order terms of the RPC row and column rational polynomials by the compensation factor β = 1 / (1+α) to maintain the geometric consistency between the enhanced image and the RPC and improve the cross-scale matching success rate; where β: is the compensation factor of the higher-order terms of the RPC; α: is the wavelet high-frequency enhancement coefficient, which is 0.8.

[0052] In this embodiment, after downsampling each layer of the pyramid image, the present invention introduces a "wavelet detail injection" mechanism to enhance the texture of the 3.5m panchromatic image and simultaneously correct the high-frequency components of the RPC normalization parameters to solve the problems of insufficient texture and decreased geometric consistency during cross-scale matching. The specific process is as follows: First, Daubechies-4 wavelet decomposition is performed on the downsampled image to obtain the low-frequency subband LL and three high-frequency subbands LH, HL, and HH, which carry detail information in the horizontal, vertical, and diagonal directions, respectively; then, the coefficients of the high-frequency subbands are uniformly multiplied by the enhancement coefficient α = 0.8 to obtain the enhanced detail subbands LH′, HL′, and HH′. The LL and the enhanced detail subbands are then subjected to inverse wavelet transform to reconstruct an image with richer texture. To quantify the enhancement magnitude, the amount of texture information enhancement is defined as follows: ΔT = |I_enhanced − I_original| / I_original, where I_enhanced represents the enhanced image grayscale value and I_original represents the original grayscale value without enhancement. This invention requires ΔT ≥ 12% to ensure that the enhancement operation makes a substantial positive contribution to cross-scale matching. Simultaneously, to avoid geometric deviations between the enhanced image details and the original RPC, this invention multiplies the coefficients of the higher-order terms of the RPC's row and column rational polynomials by a compensation factor β = 1 / (1+α), where α is 0.8 and β is approximately 0.556. This correction maintains consistency between the enhanced image and the RPC at high-frequency scales, making subsequent initial matching values ​​more reliable. Through the above synchronous enhancement-correction mechanism, image texture is significantly enhanced, while the geometric control information still strictly corresponds to the original RPC, thereby effectively improving the matching success rate across scales and resolutions, reducing mismatches at vegetation edges and weak texture areas, and laying a high-quality foundation for subsequent connection point densification and adjustment.

[0053] S17. Introduce laser waveform gradient weights in the reverse consistency check, specifically including: S171. Calculate the gradient dI / dx along the direction of the footprint image to form the gradient map G; S172. Normalize the gradient graph G to obtain the waveform gradient weight w = |dI / dx| / Σ|dI / dx|, w∈[0,1]; Where dI / dx is the grayscale gradient of the footprint image along the row direction (unit: DN / pixel). Σ|dI / dx|: The sum of the absolute values ​​of the gradients of all pixels within the current search window, used for normalization; S173. Use w as the matching quality weight. If w < 0.3, mark the matching point as invalid to reduce false matching of vegetation edges. S174: Weights are applied only to the top of the pyramid, reducing computation by 30% while maintaining matching reliability.

[0054] In another embodiment, during the reverse consistency check, this invention introduces a "laser waveform gradient weight" mechanism to suppress mismatches in weak texture areas such as vegetation edges. The specific steps are as follows: First, the grayscale gradient dI / dx is calculated along the row direction of the footprint image to obtain a gradient map G, where the gradient value of each pixel reflects the degree of grayscale change along that row direction. Then, the gradient map G is normalized, and the waveform gradient weight w = |dI / dx| / Σ|dI / dx| is calculated. Here, dI / dx represents the grayscale gradient of the current pixel, and Σ|dI / dx| represents the sum of the absolute values ​​of the gradients of all pixels within the current search window. This normalization operation controls the weight w within the range of [0,1]. The larger the gradient, the closer w is to 1, indicating that the texture of the area is reliable and the match is trustworthy. w is used as the matching quality weight. If w < 0.3, the matching point is determined to be located in a weak texture area or at the edge of vegetation and is discarded; otherwise, it is retained. To reduce computational burden, this invention performs the weight calculation only at the top layer of the pyramid, without passing it down to the next layer, reducing the overall computational load by approximately 30% while maintaining high reliability. Through the aforementioned gradient weight filtering, mismatched points in error-prone areas such as vegetation edges and areas with weak textures are eliminated in advance, significantly improving the overall reliability of connection points and providing cleaner and more accurate observational data for subsequent adjustment.

[0055] Furthermore, in step S2, during the RPC correction stage, ridge estimation regularization is adopted, with a regularization coefficient λ = 0.005 × (base-to-height ratio)², to suppress the ill-conditioned normal equation caused by the ultra-large base-to-height ratio of 19°+19°. In step S3, the 8m footprint image is reconstructed with a super-resolution of 0.5m. The laser waveform half-width at half-maximum is used as the reconstruction constraint to improve the effective resolution of the footprint image to 2.5m, and then it is matched with the 3.5m stereo image.

[0056] In the RPC correction stage of step S2, this invention addresses the ill-conditioned problem of the normal equations caused by the extremely large base-to-height ratio (approximately 3.8) formed by the 19° forward and 19° back sights of the Jumang satellite by introducing a ridge estimation regularization mechanism. Specifically, when constructing the normal equations, a regularization term λI is added to the diagonal of the normal matrix, where the regularization coefficient λ is not a constant but directly linked to the imaging geometry: λ = 0.005 × (base-to-height ratio)². This design allows the regularization intensity to automatically increase with the base-to-height ratio, effectively suppressing oscillations and divergences in parameter estimation under large-angle stereo conditions, while simultaneously reducing the condition number from 10. 8 The magnitude dropped to 10 4 The order of magnitude is reduced, the number of iterations to convergence is reduced by about 40%, and geometric accuracy is not significantly sacrificed.

[0057] After proceeding to step S3, to address the difficulty of matching 8m laser footprint images with 3.5m stereo images across resolutions, this invention performs 0.5m super-resolution reconstruction of the footprint images and embeds the laser waveform half-width at half-maximum (FWHM) as a physical constraint into the reconstruction process. Specifically, a 3.5m panchromatic image from the same track is used as training samples to construct an ESRGAN generative network. Simultaneously, an FWHM constraint term is added to the loss function to ensure the waveform width of the reconstructed image remains consistent with the original laser waveform. Through joint optimization, the effective resolution of the footprint images is increased from 8m to 2.5m, and the texture details better match the actual ground height distribution. Subsequently, the 2.5m footprint images are matched with the 3.5m stereo images, significantly reducing the cross-scale difference and increasing the matching success rate by approximately 18%. This also reduces sub-pixel-level deviations caused by resolution differences, laying a highly reliable data foundation for the accurate generation of subsequent laser control points.

[0058] Furthermore, during local DSM matching, a finer grid of 0.3 pixels × 0.3 pixels is used for mountainous areas with a slope > 15°, while a grid of 0.5 pixels × 0.5 pixels is maintained for plains areas to reduce the failure rate of matching in mountainous areas. In the joint adjustment, the laser control points are weighted according to NDVI: the weight is multiplied by 0.7 when NDVI>0.6 and by 1.0 when NDVI≤0.6, in order to suppress forest canopy bias. In the distance reciprocal weighted correction formula, an elevation accuracy factor σH=|Hcanopy-Hground| / 2 is introduced, with a weight of 1 / (Di×σH_i), to reduce the vegetation-ground double elevation unmixing error; Where Hcanopy is the elevation value of the first peak of the laser echo (vegetation canopy) (unit: m); Hground is the elevation value of the last peak of the laser echo (ground) (unit: m); |Hcanopy-Hground| represents the elevation difference (m) between the two peaks, often referred to as the "canopy-ground elevation difference" or "laser beamwidth"; σH_i is the σH value (m) corresponding to the i-th effective laser control point. Di is the horizontal distance (m) on the ground from the intersection control point to the laser point.

[0059] Furthermore, this method is only enabled when jointly processing 19° forward / backward stereo images from the terrestrial ecosystem carbon monitoring satellite with 8m laser footprint data to form a dedicated algorithmic closed loop.

[0060] In the local DSM dense matching stage, this invention addresses the matching failure problem caused by drastic slope changes and texture repetition in mountainous areas by implementing a "slope-adaptive mesh" strategy: when the local slope is greater than 15°, the matching mesh is refined from the conventional 0.5 pixel × 0.5 pixel to 0.3 pixel × 0.3 pixel, reducing the matching element size by nearly half, allowing for more precise capture of terrain details; for gently sloping plains, the mesh size remains at 0.5 pixel × 0.5 pixel, ensuring accuracy while avoiding unnecessary computational burden. Through this adaptive mechanism, the matching failure rate in mountainous areas decreases by about one-third, while the overall computation time only increases by 5%, achieving a balance between accuracy and efficiency.

[0061] During the joint adjustment phase, laser control points often exhibit "elevation-plane" unmixing errors due to dual echoes from vegetation canopy and ground. To address this, this invention introduces an NDVI-based weighting mechanism: NDVI is calculated simultaneously using multispectral data from the same track, classifying laser control points into "high vegetation" and "low vegetation" categories. When NDVI > 0.6, the point is considered susceptible to canopy influence, and its weight is multiplied by 0.7; when NDVI ≤ 0.6, the weight remains at 1.0. This effectively suppresses the excessive influence of canopy points on the adjustment in forest areas, reducing the elevation deviation in vegetated areas from 0.6m to 0.08m without requiring additional ground measurements.

[0062] To further reduce the unmixing error between vegetation and ground elevations, this invention embeds an "elevation accuracy factor σH" into the distance-inverse weighted correction formula. Specifically, for each effective laser point, the half-width σH = |Hcanopy−Hground| / 2 is calculated using its echo first peak (canopy) elevation Hcanopy and last peak (ground) elevation Hground. A larger value indicates higher elevation uncertainty at that point. Subsequently, in the weighting calculation, the traditional 1 / Di is changed to 1 / (Di×σH_i), giving laser points that are close and have small elevation dispersion a larger correction weight, while points that are far away or have large wavewidths are automatically suppressed. This improvement makes the distribution of corrections in the vegetation area more reasonable, further reducing the overall elevation RMSE by 15%.

[0063] Furthermore, this invention limits all the above-mentioned differential feature steps to be enabled only when the 19° forward / backward stereo imagery of the terrestrial ecosystem carbon monitoring satellite (Jumang) is combined with 8m laser footprint data, forming a dedicated algorithm closed loop.

[0064] By adopting the above-disclosed technical solution of this invention, the following beneficial effects are obtained: This invention provides a method for adjusting a satellite laser-assisted regional network for terrestrial ecosystem carbon monitoring. It divides stereo image tie points into intra-scene and inter-scene tie points, and obtains reliable stereo tie points through pyramid hierarchical matching, progressive refinement, and reverse matching. Through multiple iterative calculations and progressively reducing the error threshold, gross errors in tie points are eliminated, achieving free network adjustment. Corresponding stereo image points to laser points are obtained through local DSM extraction and point projection calculation in the laser point's region, thereby acquiring laser control points and avoiding direct matching failures. By modifying the elevation in the intersection coordinates of the corresponding stereo image points to the laser point's elevation, the reliable elevation accuracy of laser points is utilized, avoiding the instability of laser point planar accuracy. Based on the distribution of object-space residual errors of laser points, multiple iterations and progressive reduction of the gross error threshold are performed. Gross errors in laser points are eliminated using the object-space RANSAC algorithm, ensuring adjustment accuracy. By modifying the intersection coordinates of tie points using a distance-inverse weighted calculation of the laser point object-space coordinate correction, absolute accuracy can be improved with the help of laser control points without compromising the accuracy of the free network.

[0065] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for adjusting a satellite laser-assisted regional network for monitoring carbon in terrestrial ecosystems, characterized in that, Includes the following steps: S1. Generation of stereoscopic image connection points: Construct a multi-level pyramid for the forward-looking and backward-looking 19° panchromatic images, match layer by layer and perform reverse consistency checks to obtain the connection points between the scene and the scene. S2. Free network adjustment: Using the connection points as control, iteratively correct the image RPC, gradually tighten the threshold to eliminate gross errors, and obtain a free network adjustment RPC that is independent of the absolute elevation; S3. Laser control point preparation: Analyze satellite laser footprint data, extract local footprint images, and combine the free network adjustment RPC to generate laser control points corresponding to the stereo image; S4. Laser-stereo joint adjustment: The laser control points and the connection points of the stereo image are used together as observations to perform forward intersection, and the intersection coordinates are corrected using the laser elevation as the true value to obtain the corrected intersection control points. S5. Final RPC determination: The image RPC is corrected again using the corrected intersection control points. After iterative convergence, the final adjustment parameters that combine free network geometric consistency and laser absolute elevation accuracy are obtained. Step S1, the generation of stereoscopic image connection points, specifically includes the following steps: S11. Use the DN mean of a 3×3 rectangular window to downsample the 19° panchromatic images of the front and back views step by step until the image width and height are ≤64 pixels, forming an n-layer pyramid. Then, the corresponding RPC row and column normalization parameters are synchronously scaled by 1 / 3ⁿ to obtain the top pyramid image and its top RPC. S12. Project the center pixel of the forward-looking top-level image onto the HEIGHT_OFF elevation surface through the top-level RPC to obtain the geodetic coordinates (L,B,H). Then, back-project the image points of the back-looking top-level RPC to obtain the image points of the back-looking image. Calculate the offsets (roff,coff) between the forward and backward image points: roff=rcB–rcF coff=ccB–ccF; S13. Divide the top-level image of the front view into a 7×7 uniform grid. With the center of the grid as the initial point, take an 11×11 pixel window. Calculate the normalized cross-correlation coefficient pixel by pixel in the 7×7 search area of ​​the top-level image of the rear view. Select the point with the largest correlation coefficient as the forward matching point. Then, use this point as the search center to reverse match the front view image. If the distance between the reverse matching point and the initial point is greater than 1 pixel, it is marked as invalid. Otherwise, it is a valid matching point pair. S14. Multiply the coordinates of the effective matching point pairs by 3 and project them to the next layer of the pyramid. Insert one encrypted initial point between two adjacent points. Repeat the forward-backward matching until the original resolution is reached to obtain the in-scene connection points. S15. Calculate the approximate ground distance based on the LONG_OFF and LAT_OFF values ​​of the two image pairs' RPCs. If the distance is less than the threshold, perform projection of the four corner points to determine the overlapping area and correct the LINE_OFF and SAMP_OFF ​​values ​​of the corresponding RPCs, generating images 1-4. Perform layer-by-layer matching on images 1-4 using the same strategy as steps S11-S14, and retain only the four-fold matching point pairs that satisfy r1=r1' and r3=r3', finally obtaining the inter-scene connection points. (r1, c1) and (r3, c3) are matching point pairs between image 1 and image 3; (r1', c1') and (r2, c2) are matching point pairs between image 1 and image 2; (r3', c3') and (r4, c4) are matching point pairs between image 3 and image 4; The specific steps for generating the free network adjustment in step S2 include: S21. Perform forward intersection on the in-scene and inter-scene connection points obtained in step S1 to obtain ground coordinates (Lbf, Bbf, Hbf), and combine the ground coordinates with the forward and backward image points in the corresponding connection points to form the intersection control points of each image. S22. Count all intersection control points of each image, and use their image point coordinates as observations to correct the "row and column rational polynomial numerator constants" and "normalized longitude and latitude linear terms" in the image's RPC to obtain the first-order correction RPC; S23. Using a single-correction RPC, back-project the intersection control points onto the image plane and calculate the distance Δp between the projected image points and the original image points: Δp= ; Δp: Image plane distance, used to measure the geometric deviation between the image point after RPC backprojection and the original image point; rproj: The row coordinates of the "projected image point" obtained from the RPC back projection calculation; robs: Row coordinates of the "observed image point" obtained from the original matching or measurement; cproj: The column coordinates of the "projected image points" obtained from the RPC back projection calculation; cobs: The column coordinates of the "observed image points" obtained from the original matching or measurement; If Δp is greater than a threshold Merr, then mark the control point as invalid; otherwise, retain it. S24. Alternately execute steps S21-S23 a total of 5 times. In the first time, Ndd=5 and Merr=9 pixels. In the subsequent 4 times, Ndd=2 and Merr decreases to 7, 5, 3 and 1 pixels respectively. The final RPC is the free network adjustment RPC. Step S3, the preparation of the laser control point, specifically includes the following steps: S31. Parse the satellite laser h5 file to obtain the ground coordinates (L,B,H) of the laser point, and obtain the image plane coordinates (rld,cld) by orthographic projection of the footprint image RPC. S32. Extract a rectangular region with height Hld and width Wld centered at (rld,cld) to generate a segmented footprint image. Simultaneously correct the LINE_OFF and SAMP_OFF ​​of the footprint image RPC so that the coordinates of the top left corner image point of the segmented image are (0,0) and the coordinates of the laser point image point within the segmented image are (Hld / 2,Wld / 2). This forms the footprint observation data of the laser control point. S33. Project the same (L,B,H) image onto the RPC of the back view image, and record the sequence number, row / column coordinates and laser point sequence number of all overlapping back view images. Generate a laser control point with the same name for each overlap. The total number of control points is greater than or equal to the total number of laser points. Step S4, the laser-stereo joint adjustment, includes: S41. Using the same laser control point obtained in step S33 as the center, crop out 1000×1000 pixel local images of the front and back views, and use free network adjustment RPC to calculate the local DSM of the area. S42. Project the laser ground coordinates (L,B,H) onto the ground using the backsight adjustment RPC to obtain the backsight image point (Rb,Cb); then use the backsight adjustment RPC and the local DSM to project back onto the ground to obtain the ground coordinates (L1,B1,H1). S43. Project (L1,B1,H1) onto the forward-looking image point (Rf,Cf) using the forward-looking adjustment RPC orthographic projection. Then replace H1 with the original laser elevation H to form the final ground coordinates (L1,B1,H). S44. A laser control point is formed by (L1,B1,H), the back view image point (Rb,Cb), and the front view image point (Rf,Cf). Repeat steps S41-S44 to complete the preparation of all laser control points. Step S5, the determination of the final RPC, specifically includes the following steps: S51. Using the free network adjustment RPC as the initial value, perform forward intersection on the connection points of the stereo image to obtain the ground coordinates (Lbf, Bbf, Hbf) of the intersection control point, and combine them with the corresponding forward and backward image points to form the intersection control point. S52. For each laser control point, intersect its forward image point (Rf, Cf) and backward image point (Rb, Cb) to obtain the intersection coordinates (Ljs, Bjs, Hjs), and calculate the difference between these coordinates and the laser coordinates (L1, B1, H). Loffi=L1−Ljs Boffi=B1−Bjs Hoffi=H−Hjs After multiplying Loffi and Boffi by 100000 to convert them to meters, the error threshold Merrld∈{10,8,6,4,2}m was successively tightened using RANSAC to eliminate gross errors, and the remaining points were the effective laser control points. S53. Calculate the correction amount for each intersection control point using a weighted average of the inverse distance: Lall=Add(Loffi / Di) / Add(1 / Di); Ball=Add(Boffi / Di) / Add(1 / Di); Hall=Add(Hoffi / Di) / Add(1 / Di); Where Add represents cumulative calculation, Di (i=0,…, Nldyx-1) is the ground distance between the intersection control point and the effective laser control point, and (Lbf, Bbf, Hbf) is added to (Lall, Ball, Hall) to obtain the corrected intersection control point; S54. Using the corrected intersection control points, the row and column numerator constants and the normalized latitude and longitude linear terms of the image RPC are corrected again. Steps S51-S54 are executed iteratively a total of 5 times. The RPC obtained in the last iteration is the final adjustment parameter.

2. The method according to claim 1, characterized in that, Step S1, the generation of stereoscopic image connection points, further includes the following steps: S16. After downsampling each layer of the pyramid, wavelet detail injection is used to enhance the texture of the 3.5m panchromatic image, and the high-frequency components of the RPC normalization parameters are corrected simultaneously, specifically including: S161. Perform Daubechies-4 wavelet decomposition on the downsampled image to obtain the low-frequency subband LL and the high-frequency subbands LH, HL, and HH. S162. Multiply the high-frequency subband coefficients by the enhancement coefficient α = 0.8 to obtain the enhanced high-frequency subbands LH′, HL′, and HH′. S163. Perform inverse wavelet transform on LL, LH′, HL′, and HH′ to obtain texture-enhanced images. The texture information enhancement amount ΔT = |I_enhanced−I_original| / I_original, ΔT ≥ 12%; Where I_enhanced represents the image grayscale value after wavelet detail injection; I_original represents the original grayscale value after downsampling and without enhancement. S164. Synchronously correct the high-frequency components of the RPC normalization parameters: multiply the coefficients of the higher-order terms of the RPC row and column rational polynomials by the compensation factor β = 1 / (1+α) to maintain the geometric consistency between the enhanced image and the RPC and improve the cross-scale matching success rate; where β: is the compensation factor for the higher-order terms of the RPC; α: is the wavelet high-frequency enhancement coefficient, which is 0.

8. S17. Introduce laser waveform gradient weights in the reverse consistency check, specifically including: S171. Calculate the gradient dI / dx along the direction of the footprint image to form the gradient map G; S172. Normalize the gradient graph G to obtain the waveform gradient weight w = |dI / dx| / Σ|dI / dx|, w∈[0,1]; Where dI / dx is the grayscale gradient of the footprint image along the row direction; Σ|dI / dx|: The sum of the absolute values ​​of the gradients of all pixels within the current search window, used for normalization; S173. Use w as the matching quality weight. If w < 0.3, mark the matching point as invalid to reduce false matching of vegetation edges. S174: Weights are applied only to the top of the pyramid, reducing computation by 30% while maintaining matching reliability.

3. The method according to claim 2, characterized in that, In step S2, during the RPC correction stage, ridge estimation regularization is used with a regularization coefficient λ = 0.005 × (base-to-height ratio)² to suppress the ill-conditioned normal equation caused by the ultra-large base-to-height ratio of 19° + 19°. In step S3, the 8m footprint image is reconstructed with a super-resolution of 0.5m. The laser waveform half-width at half-maximum is used as the reconstruction constraint to improve the effective resolution of the footprint image to 2.5m, and then it is matched with the 3.5m stereo image.

4. The method according to claim 3, characterized in that, During local DSM matching, a finer grid of 0.3 pixels × 0.3 pixels is used for mountainous areas with a slope > 15°, while a grid of 0.5 pixels × 0.5 pixels is maintained for plains areas to reduce the failure rate of matching in mountainous areas. In the joint adjustment, the laser control points are weighted according to NDVI: the weight is multiplied by 0.7 when NDVI>0.6 and by 1.0 when NDVI≤0.6, in order to suppress forest canopy bias. In the distance reciprocal weighted correction formula, an elevation accuracy factor σH=|Hcanopy−Hground| / 2 is introduced, with a weight of 1 / (Di×σH_i), to reduce the vegetation-ground double elevation unmixing error; Where Hcanopy is the elevation value of the first peak of the laser echo; Hground represents the elevation value of the final peak of the laser echo; |Hcanopy − Hground| represents the elevation difference between the two peaks, often referred to as the "canopy-ground elevation difference" or "laser beamwidth"; σH_i is the σH value corresponding to the i-th effective laser control point; Di is the horizontal distance on the ground from the intersection control point to the laser point.

5. The method according to any one of claims 2 to 4, characterized in that, This method is only enabled when processing 19° forward / backward stereo images from the terrestrial ecosystem carbon monitoring satellite in conjunction with 8m laser footprint data to form a dedicated algorithm closed loop.

Images