CN105550639A - Automatic extraction method for Earth observation laser height measurement satellite elevation control points and data processing method - Google Patents

Abstract

The invention relates to an automatic extraction method for Earth observation laser height measurement satellite elevation control points and a data processing method. The laser elevation control point extraction method comprises steps that an effective earth observation laser distance value measurement evaluation method is employed, determined cloudless footprint image blocks are kept, and laser elevation data of determined thin-cloud or thick-cloud footprint image blocks are removed; reflectivity Epsilon smaller than 1 of laser footprint points is taken as a screening parameter, laser footprint points of the kept footprint image blocks are screened, Epsilon=reception pulse energy/emission pulse energy, laser footprint points which have only one wave peak, have the peak value greater than the threshold, have the standard deviation sigma not greater than 3.2ns after waveform fitting are selected from an echo waveform, and parameters used for determining the threshold comprise the emission energy and a reception caliber of a laser device. Through the method, influence of clouds on laser distance measurement can be reduced, laser distance measurement precision can be guaranteed, and accuracy of the laser elevation reference data is effectively improved.

Description

Automatic extraction method and data processing method for elevation control points of earth observation laser altimetry satellite

Technical Field

The invention relates to a satellite data processing method, in particular to a cloud and mist amount detection method of a laser footprint image, an elevation control point automatic extraction method based on echo waveform processing and cloud and mist amount detection, and a satellite data processing method for earth observation laser height measurement.

Background

With the development of economy, three-dimensional geographic information has been widely applied in the aspects of digital earth, urban planning, environmental protection and the like. The rapid development of the remote sensing satellite technology makes the space photography measurement become a means for rapidly acquiring three-dimensional geographic information after the aviation photography measurement, and particularly, the technology for acquiring the three-dimensional remote sensing information makes great progress along with the development of the technologies such as a three-linear array stereo camera and the like in recent years.

For example, resource III is the first high-precision civil three-dimensional mapping satellite in China, the successful launching and the effective application of the three-dimensional mapping satellite break through the dependence of China on foreign high-precision satellite images for a long time, and huge social and economic benefits are generated. However, due to the three-dimensional mapping mode of the optical three-dimensional satellite, the attitude and orbit measurement accuracy of the satellite, the distortion of the camera and the like, the elevation measurement accuracy of the optical three-dimensional satellite is difficult to meet the requirement of high-accuracy mapping under the condition of no ground control point.

The ground control point gcp (group control points) is an important reference data source for geometric correction and geographic positioning of the satellite remote sensing images. In the geometric correction processing process of the remote sensing image, a certain number of ground control points are necessary for achieving a certain correction precision, an image imaging model is constructed through object space coordinates of the control points and corresponding image point coordinates, model parameters are solved or an existing imaging model is optimized, compensated and solved, and compensation parameters are solved, and finally a correct conversion relation between an object space and an image space in the imaging process is established.

The control point in the traditional operation process generally adopts a full field measurement mode, a series of complex procedures such as 'collecting existing control data, surveying area and surveying, selecting points and burying stones, field measurement and field arrangement' need to be carried out, although the workload and complexity of field measurement are greatly reduced along with the development of advanced measurement technologies such as GPS-RTK, the necessary field measurement work of the control point is inevitable.

In addition, since the coverage of the satellite remote sensing image is large (taking the resource three satellite as an example, the coverage is 50 kilometers × 50 kilometers), in order to obtain uniformly distributed ground control points, the measurement is usually required within a range of several hundred kilometers or even thousands of square kilometers, and the huge workload, manpower and material resources consumption are self evident. In addition, in areas with many natural disasters such as earthquake, flood and debris flow, or in regions with rare occurrences such as original forests, swamps and deserts, survey personnel often cannot enter the field to perform measurement.

Even if the control point is obtained by manually selecting the same-name point of the remote sensing image and the topographic map, the problems of low efficiency, difficulty in ensuring precision and the like exist. And the situation of repeated point selection can occur when images of different time phases or different sensors in the same area are corrected.

Ground control points are generally divided into: a level control point, a level control point and an elevation control point. Chinese patent application 201310143369.3 discloses an automatic acquisition method for control points of a multi-source heterogeneous remote sensing image, which can automatically extract remote sensing control image points (control point image slices) from a multi-source heterogeneous image, thereby improving the efficiency and accuracy of control data acquisition, but is still an automatic acquisition method for plane control points in essence.

Under the condition that the data of the ground control points is insufficient, the block adjustment technology of the satellite images can be used as an important means for accurate geometric positioning. For example, chinese patent application 201510191096.9 discloses a satellite image stereo block adjustment method based on satellite-borne laser height measurement data. In the solution disclosed in this patent application, satellite-borne laser altimetry data is used as control data of a generalized elevation control point database.

The development of global mapping is a necessary support for protecting the territorial ownership, related benefits and geographic information safety of China. The satellite remote sensing technology has unique advantages in the aspect of acquiring overseas geographic information. At present, the national geographic information department of surveying and mapping is gradually developing global high-precision three-dimensional mapping around the major and gross war requirements of 'one way taking' and 'high-speed rail going away'. The measurement and the acquisition of global elevation control points are important technical support for carrying out global surveying and mapping. Therefore, in order to smoothly perform the work related to global mapping, the acquisition of global high-precision control points must be performed preferentially, and a corresponding control point database is established. The method for acquiring the global elevation control points by using the laser altimetry satellite is an effective technical means under the current satellite remote sensing technical condition.

Satellite-borne laser height measurement is a ground point elevation measurement technology, a satellite is used as a platform, a laser height measuring instrument is carried to observe the earth from space and time, the distance between the satellite and a measured object is measured in real time at high precision, and information such as landform, vegetation coverage, sea surface form and the like on the earth surface is obtained through data processing and analysis. The satellite-borne laser height measuring technology in foreign countries develops rapidly, and the research of the satellite-borne laser height measuring instrument is carried out in the main developed countries in the world. For example, an ICESat satellite is launched in 2003 in the united states, and a Geoscience Laser Altimeter System (GLAS) carried on the ICESat satellite is the first satellite-borne laser ranging system for continuously observing the earth all over the world, and the main scientific purpose of the system is to measure the elevation and change of an ice cover of a polar region, the distribution characteristics of a cloud layer and aerosol, and the like. A large amount of high-precision elevation data are obtained during the in-orbit operation of an ICESat satellite, the laser footprint plane precision of the ICESat satellite reaches 10m magnitude, and the elevation precision is about 15 cm.

In recent years, satellite-borne laser height measurement research is very important in China. The independently researched and developed moon exploration engineering series satellites 'ChangE I' and 'ChangE II' are provided with high-precision laser altimeters, and three-dimensional images of the surface of the moon are obtained. Except that the 'ChangE' series satellites are provided with laser altimeters, no laser altimeter system for earth observation exists, but the subsequent development planning of earth observation laser altimeter satellites is included, and for example, high-precision civil three-dimensional mapping satellites with a 1:1 ten thousand scale, carbon monitoring satellites of a land ecosystem and the like are provided with laser altimeters for earth observation.

In the process that a laser beam comes and goes to the atmosphere, the earth observation laser height measurement satellite generates phenomena of scattering, refraction and the like with atmospheric molecules and aerosol, and further causes the problems of attenuation, echo deformation, reduction of distance measurement precision and the like of laser energy.

Disclosure of Invention

According to an aspect of an embodiment of the present invention, there is provided a method for detecting a cloud amount of a laser footprint image, including: carrying out sample training by utilizing a large number of cloud and fog-containing images to obtain an optimal gray threshold value and a texture characteristic value of the cloud and fog; calculating a gray level histogram of the whole footprint image, and preliminarily judging whether cloud exists or not and the content of the cloud based on the preferred gray level threshold; calculating the gray average value of the sub-blocks of the footprint image, judging that the sub-blocks contain cloud when the gray average value of the sub-blocks is higher than a first threshold value, and judging that the sub-blocks do not contain cloud when the gray average value of the sub-blocks is lower than a second threshold value; calculating texture features of the image sub-blocks with the gray mean value between the first threshold value and the second threshold value by utilizing a gray co-occurrence matrix, and comparing the texture features with the cloud and fog texture feature values obtained from the sample training; and counting the total number of the pixels judged as cloud and the number ratio of the pixels in the whole footprint image, and determining the cloud amount.

According to the cloud amount detection method of the laser footprint image, optionally, when calculating the gray histogram of the whole footprint image, histogram equalization enhancement processing is performed on the image, and the equalization formula is as follows:

$s_k=T\left(r_k\right)=\underset{i=0}{\overset{k}{\Σ}}{P}_r\left(r_i\right)=\underset{i=0}{\overset{k}{\Σ}}\frac{n_i}{n}$

0≤r_k≤1,k＝0,1,2,...L-1

wherein s is_kThe gray value of the pixel with the gray value k of the original image is converted to a new gray value; p_r(r_i) Is the pixel frequency with a gray value of i; l is the gray level of the image.

According to the method for detecting the cloud amount of the laser footprint image in the embodiment of the present invention, optionally, 235 is the first threshold, and the second threshold is 80.

According to the cloud amount detection method for the laser footprint image, optionally, the texture features considered when the texture features of the image sub-blocks are calculated by using the gray level co-occurrence matrix include: one or more of angular second moment, homogeneity, contrast and correlation.

According to another aspect of the embodiments of the present invention, there is provided a method for evaluating validity of a ground observation laser ranging value, including: performing system geometric correction and preliminary waveform processing on decoded original data downloaded from a laser altimeter satellite, obtaining a footprint image corresponding to a laser spot according to a geometric corresponding relation, and realizing basic registration of laser and the footprint image based on hardware parameters and geographic coordinates or preliminary calibration of a footprint camera and a laser altimeter; determining the position of the central point of the obtained footprint image corresponding to the laser spot, and determining the area corresponding to the laser spot in the footprint image; cutting the determined image or storing the image into a memory to form a footprint image block; carrying out cloud amount detection on the footprint image block by adopting the cloud amount detection method to obtain a cloud amount value of the footprint image block; determining validity or availability of a laser ranging value for the footprint image block based on the cloud and fog magnitude value.

According to still another aspect of an embodiment of the present invention, there is provided a laser elevation control point extraction method including: by adopting the method for evaluating the effectiveness of the earth observation laser ranging values, the footprint image blocks judged to be cloudless are reserved, and the laser elevation data of the footprint image blocks judged to be thin clouds or thick clouds are removed; reflectivity of laser footprint<1 asA screening parameter to screen laser footprint points of the retained footprint image block, wherein,and selecting a laser footprint point which has only one peak in the echo waveform, the peak value is greater than a threshold value, and the standard deviation sigma after the waveform is fitted is less than or equal to 3.2ns, wherein the parameters for determining the threshold value comprise the transmitting energy and the receiving aperture of the laser.

According to the laser elevation control point extraction method provided by the embodiment of the invention, optionally, a waveform fitting is performed on the echo waveform in a manner of superimposing a plurality of gaussian functions, and a fitting formula is shown as follows:

$w\left(t\right)=\mathrm{\&epsiv;}+\underset{m=1}{\overset{p}{\&Sigma;}}{A}_m{e}^{(-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2})}$

in the above formula, t is time, A_m,t_m,σ_mRespectively, the amplitude, mean and standard deviation of the mth gaussian function,is a waveform noise value.

According to still another aspect of embodiments of the present invention, there is provided a method for processing earth observation laser altimetry satellite data, including: smoothing and denoising the transmitting and receiving waveforms of the laser, extracting waveform characteristic parameters, determining the time corresponding to the gravity centers of the transmitting and receiving waveforms, and calculating the initial distance value of the unidirectional laser transmission according to the laser transmission time interval; calculating rough three-dimensional coordinates of the laser footprint point according to the attitude and orbit parameters of the satellite, the laser emission time, the initial distance value and the laser geometric positioning model; re-sampling to produce a footprint image geometric rough correction product according to the satellite attitude and orbit parameters, the footprint image, the footprint camera parameters and the footprint image geometric positioning model; by adopting the cloud amount detection method, the cloud amount detection is carried out on the foot print image, and the laser foot print point of the foot print image of which the detected cloud amount value does not exceed the preset standard is processed in the next step.

According to the data processing method of the earth observation laser altimetry satellite, optionally, the method further comprises the following steps: obtaining an atmospheric delay correction value according to atmospheric parameters and the rough three-dimensional coordinates of the laser footprint point by adopting an atmospheric delay correction model; calculating an accurate distance value according to the laser ranging systematic error value obtained by geometric calibration, the initial distance value and the atmospheric delay correction value; calculating the three-dimensional coordinates of the laser footprint point according to the attitude and orbit parameters of the satellite, the laser emission time, the accurate distance value and a laser geometric positioning model; and calculating a tide correction value by adopting a tide correction model, and correcting the three-dimensional coordinates of the calculated laser footprint point to obtain the accurate three-dimensional coordinates of the laser footprint point.

According to the data processing method of the earth observation laser altimetry satellite, optionally, the method further comprises the following steps: producing a footprint image orthorectification product by adopting topographic data according to the satellite attitude and orbit parameters, the footprint image, the footprint camera parameters and the footprint image geometric positioning model; judging whether the foot print image with the cloud and fog quantity value not exceeding the preset standard is cloud-free, and judging the characteristic constraint condition of the echo waveform characteristic parameter; and if the characteristic constraint condition is met, combining the laser footprint point accurate three-dimensional coordinate and the footprint image orthorectification product to form a laser elevation control point.

According to the cloud amount detection method for the earth observation laser footprint image, the extraction method for the earth observation laser elevation control point based on cloud amount detection and echo waveform processing, and the laser height measurement satellite data processing method, the influence of a cloud layer on laser distance measurement can be reduced, the precision of the laser distance measurement is ensured, and the accuracy of laser elevation reference data is effectively improved.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings of the embodiments will be briefly described below, and it is apparent that the drawings in the following description only relate to some embodiments of the present invention and are not limiting on the present invention.

FIG. 1 is a flow chart of a method for detecting cloud amount of a laser footprint image according to an embodiment of the present invention;

FIG. 2 schematically illustrates a laser elevation control point fine extraction flow according to an embodiment of the invention;

FIG. 3A is a schematic diagram of a rigorous geometric positioning model of a laser altimetry satellite, and FIG. 3B is a schematic diagram of an included angle between a laser emission direction and a body coordinate system;

FIGS. 4A and 4B schematically illustrate the transmit and receive waveforms, respectively, of a laser altimeter;

FIG. 5 schematically illustrates a process flow of satellite borne laser altimetry data according to one embodiment of the invention;

FIG. 6 schematically illustrates a hierarchy of laser altimetry satellite data products;

FIG. 7 shows a schematic production flow of a base product;

FIG. 8 shows a schematic production flow of a standard product;

fig. 9 shows a schematic production flow of a special/high-grade product.

Detailed Description

In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the drawings of the embodiments of the present invention. It is to be understood that the embodiments described are only a few embodiments of the present invention, and not all embodiments. All other embodiments, which can be derived by a person skilled in the art from the described embodiments of the invention without any inventive step, are within the scope of protection of the invention.

Unless defined otherwise, technical or scientific terms used herein shall have the ordinary meaning as understood by one of ordinary skill in the art to which this invention belongs. The use of "first," "second," and similar terms in the description and claims of the present application do not denote any order, quantity, or importance, but rather the terms are used to distinguish one element from another. Also, the use of the terms "a" or "an" and the like do not denote a limitation of quantity, but rather denote the presence of at least one.

A1: 1 ten thousand high-precision civil three-dimensional mapping satellite planned to be launched in 2018 in China is to carry a set of laser height measuring system. The laser height measurement system is used as an auxiliary load and is used for adjusting the height in combination with an optical image shot by a CCD camera, and the purpose of improving the high-resolution satellite uncontrolled positioning precision is achieved. The set of laser height measuring system is also provided with a footprint camera, and a footprint image formed by the pulse reaching the ground can be acquired, so that the real position of the footprint can be determined, and the problem of large plane positioning error of the laser height measuring system is solved. For example, the laser height measurement system adopts laser pulses with a wavelength of 1064nm, the wavelength is in an atmospheric window band and has good penetrability to atmosphere, however, when the system is in a cloud shielding or fog situation, the laser cannot penetrate the cloud fog to reach the ground, so that the obtained footprint position may be the position of the cloud fog surface, and such data cannot meet the observation requirement and needs to be rejected.

Alternatively, the measured data with the difference exceeding a certain threshold value is judged to be the measured footprint point on the cloud fog by comparing the obtained height measurement data with the existing Digital Elevation Model (DEM), so that the measured footprint point is eliminated. However, the method has certain limitation, on one hand, the method has certain requirement on the DEM precision, and on the other hand, when the laser can partially penetrate through the cloud mist to reach the ground and can receive ground echo signals, the algorithm is easy to lose effectiveness and cannot accurately reflect the reliability of the laser ranging value.

According to the laser ranging data extraction method provided by the embodiment of the invention, the laser ranging data is screened based on the detection result of the cloud amount. The cloud and fog detection and judgment are carried out by utilizing the footprint image and adopting the idea of remote sensing image classification, the laser echo waveform characteristics can be further combined, and if the cloud and fog exist on the footprint path, and the echo waveform peak value is smaller and the broadening is larger, the measured data is considered to have larger error or be unusable.

The method for detecting the cloud amount using the footprint image will be described in detail below. The cloud amount detection of the remote sensing image is essentially the classification of the image, and the cloud is regarded as one of a plurality of ground objects. Generally speaking, cloud amount detection relies on two characteristics of the cloud: spectral features and textural features. The spectral characteristics reflect radiation reflection conditions of the cloud cluster at different wave bands, including radiation characteristics of the cloud and reflection spectral characteristics of the cloud and ground objects. Texture features can be thought of as the spatial distribution and spatial interaction between image gray levels over a local range.

Since the footprint image of the laser altimetry system does not contain multispectral information, according to the invention, a cloud and fog detection method using a footprint panchromatic image is considered, optionally a threshold method based on a grey-level histogram and a detection method based on a grey-level co-occurrence matrix are used to detect cloud and fog. Furthermore, in order to improve the accuracy of cloud and mist detection, the two methods can be integrated to judge the cloud and mist.

The gray level histogram reflects the frequency of occurrence of each gray level pixel in an image, takes the gray level as an abscissa and the frequency of the gray level as an ordinate, is an important feature of the image, and reflects the gray level distribution of the image. When the gray level histogram of the whole footprint image has obvious bimodality, the method indicates that a brighter area and a darker area of the image can be well separated, and the two peak values are used as references to set a threshold value, so that a good binary image processing effect can be obtained.

And training a certain amount of remote sensing image samples to obtain an optimal threshold value of the gray value characteristics of cloud and fog, and judging that cloud exists when the threshold value is at a certain frequency in a gray level histogram of the whole remote sensing image. The cloud content can be judged overall by using the gray histogram. The image is divided into a plurality of sub-blocks, the gray average value of each sub-block is counted, and cloud and fog judgment is preliminarily realized according to the size of the average value.

The cloud has uniform gray distribution and small jumping degree on the image, and the texture is thick and fuzzy and is similar to the characteristics of partial ground objects, for example, the texture characteristics of desert and stratosphere are very close. In order to increase the accuracy of cloud and fog detection, histogram equalization enhancement processing can be performed on the image, so that the implicit texture details are highlighted. The equation for equalization is as follows:

$s_k=T\left(r_k\right)=\underset{i=0}{\overset{k}{\&Sigma;}}{P}_r\left(r_i\right)=\underset{r=0}{\overset{k}{\&Sigma;}}\frac{n_i}{n}$

0≤r_k≤1,k＝0,1,2,...L-1

in the above formula, s_kThe gray value of the pixel with the gray value k of the original image is converted to a new gray value; p_r(r_i) Is the pixel frequency with a gray value of i; l is the gray level of the image.

On the other hand, texture can be thought of as the spatial distribution between image gray levels and their spatial interrelation in a local context. The gray level co-occurrence matrix is obtained by counting the specific gray levels of two pixels which are kept at a certain distance from an image, can reflect the distribution characteristic of the gray levels, can reflect the position distribution characteristic among the pixels with the same gray level or the pixels close to the same gray level, and is a method based on the image gray level joint probability matrix, namely a second-order statistical characteristic about the change of the image gray level. Because cloud detection is based on the detection of the gray value of the footprint image, the cloud detection scheme can use the gray co-occurrence matrix as the basis for feature extraction.

The gray level co-occurrence matrix is to count the frequency P (i, j, d, theta) of the simultaneous occurrence of pixels (x + a, y + b) with the distance d and the gray level j from the pixel (x, y) with the gray level i of the image, and the mathematical expression is as follows:

P(i,j,d,θ)＝{[(x,y),(x+a,y+b)|f(x,y)＝i；f(x+a,y+b)＝j]}

where θ is the generation direction of the gray level co-occurrence matrix. In short, d and θ determine the positions and directions of two pixels, and the gray level co-occurrence matrix P can be obtained under certain conditions:

$P=\left(\begin{array}{ccc}P(0,0)& ...& P(0,L-1)\\ ...& ...& ...\\ P(L-1,0)& ...& P(L-1,L-1)\end{array}\right)$

where each element P (i, j) represents the number of times gray i and gray j occur at a given position and orientation, and L represents the number of gray levels. Only 4 directional values need to be found to represent all texture information.

The following 4 feature quantities f may be used₁，f₂，f₃，f₄As texture features:

(1) angular second moment f₁

$f_1=ASM=\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}{P}^2(i,j)$

The uniformity of image distribution and the thickness degree of texture are reflected: when most elements are distributed close to the main diagonal, the gray level distribution of the image is uniform, the texture is thicker, and therefore f₁The value of (a) will be large; otherwise f₁Smaller, indicating finer texture.

(2) Homogeneity f₂

$f_2=IDM=\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}\frac{P(i,j)}{1+{(i-j)}^2}$

The inverse difference moment reflects the homogeneity of the image texture, and if the local texture of the image lacks variation and is very uniform, the value is large.

(3) Contrast f₃

$f_3=\underset{n=1}{\overset{L}{\&Sigma;}}{n}^2\left(\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}P\right(i,j\left)\right)$

Where i-j is n, the image contrast reflects the sharpness of the texture. In an image, if the gray difference of a pixel pair is larger, the contrast of textures is stronger, and the visual effect is more obvious. Otherwise, it indicates that the texture is not apparent. The contrast of rough texture is often small, so the contrast may also reflect the thickness of the texture.

(4) Correlation f₄

$f_4=COR=\frac{\underset{i=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{L}{\&Sigma;}}(i-j)P(i,j)-{\mathrm{\&mu;}}_x{\mathrm{\&mu;}}_y}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}$

Wherein,

${\mathrm{\&mu;}}_x=\underset{i=1}{\overset{L}{\&Sigma;}}i\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j),{\mathrm{\&mu;}}_y=\underset{i=1}{\overset{L}{\&Sigma;}}j\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j)$

${\mathrm{\&sigma;}}_x=\underset{i=1}{\overset{L}{\&Sigma;}}{(i-{\mathrm{\&mu;}}_x)}^2\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j),{\mathrm{\&sigma;}}_y=\underset{i=1}{\overset{L}{\&Sigma;}}{(j-{\mathrm{\&mu;}}_y)}^2\underset{j=1}{\overset{L}{\&Sigma;}}P(i,j)$

the correlation reflects the direction of the texture, which measures the similarity of the gray level co-occurrence matrix elements in the row or column direction. The longer a certain gradation extends in a certain direction, the larger the correlation value in that direction. When similar texture regions in an image have a certain directionality, the correlation value is large.

And (3) acquiring the texture features of the cloud and fog by using a large number of images containing the cloud and fog as training samples, and preferably adopting an angular second moment as the texture features as a threshold value for judging the amount of the cloud and fog.

Then, combining the two methods of the gray level histogram and the gray level co-occurrence matrix to detect the cloud and fog of the footprint image, the steps are as follows:

carrying out sample training by utilizing a large number of cloud and fog-containing images to obtain an optimal gray threshold value and a texture characteristic value of the cloud and fog;

calculating a gray level histogram of the whole footprint image, preliminarily judging whether cloud and fog exist and the content of the cloud and fog, and optionally performing gray level equalization enhancement on the footprint image;

calculating the average value of the gray levels of the sub-blocks of the footprint image, judging whether cloud fog is contained or not according to the average value, for example, the gray level is higher than 235 and represents that the cloud fog is contained, and the gray level is lower than 80 and does not contain the cloud fog, and further detecting the sub-blocks of the image between the two threshold values;

calculating the texture characteristics of the image subblocks between the two threshold values by utilizing a gray level co-occurrence matrix, and comparing the texture characteristics with cloud and fog texture characteristic values obtained by sample training;

and counting the total number of the pixels judged to be cloud and the number ratio (percentage) of the pixels in the footprint image to determine the cloud amount.

Fig. 1 schematically shows the steps of the cloud amount detection method for laser footprint images.

Based on the laser footprint image cloud amount detection method, according to the embodiment of the invention, a laser ranging value effectiveness evaluation method based on footprint image cloud detection is provided.

The basic flow of the laser ranging value validity evaluation method is as follows:

and performing system geometric correction and primary waveform processing on decoded original data downloaded from the laser height measurement satellite, obtaining a footprint image corresponding to a laser spot according to the geometric correspondence, and realizing basic registration of the laser and the footprint image based on hardware parameters and geographic coordinates. If the preliminary calibration of the footprint camera and the laser altimeter is performed, the calibration result can also be directly used for the registration of the laser and the footprint image.

Since the image obtained by the footprint camera is much larger than the actual ground footprint size of the laser, some clipping is required. Determining the central point position of the obtained footprint image corresponding to the laser spot, and determining the area S ═ Pi R corresponding to the laser spot in the footprint image²Wherein R is the radius of the laser footprint on the ground, and R can be doubled in consideration of the laser pointing precision so as to ensure the full coverage of the footprint image on the laser footprint point. And cutting the determined image or storing the determined image into a memory to form a footprint image block.

The cloud amount detection algorithm is adopted to perform cloud amount detection on the footprint image block to obtain a cloud amount value of the footprint image block, and the cloud condition is judged based on the cloud amount value, for example, if the cloud amount is less than 10%, the footprint image block is judged to be cloudless, if the cloud amount is greater than 10% and less than 30%, the footprint image block is judged to be thin cloud, and if the cloud amount is greater than 30%, the footprint image block is judged to be thick cloud.

The validity or availability of the laser ranging value for the footprint image block may then be determined from the cloud magnitude value. For example, the laser range measurement value of the footprint image block determined to be cloudless is marked as good, the laser range measurement value of the footprint image block determined to be thin cloud is marked as good, and the laser range measurement value of the footprint image block determined to be thick cloud is marked as bad. Further, the laser ranging value marked as excellent is fully usable, the laser ranging value marked as good is considered to be standby, the laser ranging value marked as poor is not usable at all, and the laser ranging value marked as poor is not involved in the production of subsequent products.

An important application of a laser height measurement satellite is that high-precision elevation control point data in a global range can be obtained, however, the final elevation precision of a ground footprint point of the laser height measurement satellite is influenced by the distance measurement precision of a laser, and is also related to factors such as the satellite position and attitude measurement precision, the atmospheric influence, the ground fluctuation size and the like, so that laser elevation reference points subjected to precise geometric calibration, atmospheric correction, tide correction and the like are required to be screened in combination with related auxiliary parameters, and reliable laser elevation reference data are guaranteed to be finally reserved as special products. By combining the laser ranging value effectiveness evaluation method based on footprint image cloud detection, according to the embodiment of the invention, a multi-parameter constrained laser elevation control point extraction method is further provided, so that the effective production of the special product of the laser elevation control point is realized.

The basic steps of the multi-parameter constrained laser elevation control point extraction method are as follows:

reserving the laser ranging value of the footprint image block judged to be cloud-free, and eliminating the laser elevation data of the footprint image block judged to be thin cloud or thick cloud;

since the received pulse energy is susceptible to atmospheric backscattering, solar veiling glare, overexposure, etc., the reflectivity of the laser footprint is determined<1 (wherein) As a screening parameter, screening laser footprint points in the laser elevation data;

selecting a laser footprint point which has only one peak in the echo waveform, namely nPeaks is 1, the peak E is larger than a certain threshold V (the specific threshold is determined according to parameters such as the transmitting energy and the receiving caliber of a laser), and the standard deviation sigma after waveform fitting is less than or equal to 3.2 ns;

and taking the laser footprint points obtained after screening as laser elevation control points, and combining footprint images to construct a laser elevation control point database.

FIG. 2 schematically illustrates a laser elevation control point fine extraction flow according to an embodiment of the invention. In the process shown in fig. 2, the topographic reference data and the footprint image orthorectification product are further combined to improve the accuracy of the laser elevation control point extraction.

The method of defining and fitting the earth observation laser altimetry satellite echo waveform is described as follows.

In order to reflect the vertical distribution of the ground features more accurately and reveal the geometric and physical properties of the ground features, a laser altimetry satellite for ground observation generally carries out sampling recording on backscattering energy formed by the action of transmitted laser pulses and a target at small equal sampling intervals (such as 1GHz) in time sequence to obtain an echo signal which changes along with time, and the echo signal is called echo waveform data. The echo waveform includes characteristic information such as amplitude, pulse width, number of peaks, backscatter cross section, etc., and waveform characteristic analysis is performed based on the echo waveform data, and may also constitute an important step in laser height measurement data processing.

The basic principle of the laser height measurement satellite is as follows: laser beams emitted by a satellite are received by the satellite after being reflected by the ground, the time interval t between the emission and the reception of the laser is calculated, the propagation speed of light is c, the laser one-way transmission distance p is c t/2, and then the three-dimensional ground coordinates of the laser footprint point can be obtained by combining the satellite position and attitude information obtained by a GPS (global positioning system) and a star sensor carried on the satellite. Its rigorous geometric model is shown in FIG. 2, where P_laserReference point for laser emission, P_GPSIs the GPS antenna phase center, O_BodyAs the center of mass of the satelliteOr origin of the body coordinate system, P_GroundIs a laser ground footprint point.

The original attitude parameters of the satellite are all defined under an international celestial sphere reference system ICRS (International celestial Reference System), and an international celestial sphere reference frame ICRF (International celestial Reference frame) is a specific implementation of the ICRS, and currently, a J2000 coordinate system is mainly adopted. The generally described ground three-dimensional coordinates are defined under the international earth reference system ITRS (international terrestrialreferencesystem), such as the WGS84 coordinate system that best fits the ITRF2000, which is a specific implementation of ITRS. The coordinate transformation matrix from the international celestial sphere reference frame ICRF to the international earth reference frame ITRF is related to the earth rotation, the time difference, the nutation and the like,wherein P (t) is a time matrix, N (t) is a nutation matrix, R (t) is a rotation matrix of the earth, and W (t) is a polar shift matrix.

The satellite body coordinate system is defined as: the center of mass of the satellite is the origin, the X-axis points to the direction of the satellite flight, the Z-axis points to the zenith direction, the Y-axis is perpendicular to the satellite orbit plane, and forms a right-hand coordinate system with the X, Z-axis, as shown in fig. 3A. When the laser is emitted, a certain included angle exists between the pointing direction and the satellite body coordinate system, and if the negative included angle between the laser pointing direction and the Z axis of the body coordinate system is θ, the included angle between the projection on the XOY plane and the positive direction of the X axis is α, as shown in fig. 3B.

If the distance measurement value of the laser is rho, the coordinate of the laser reference point in the body coordinate system is as follows:the coordinates of the laser footprint point in the satellite body coordinate system are:

$\stackrel{\&OverBar;}{\mathrm{\&rho;}}=\left(\begin{array}{c}Lx\\ Ly\\ Lz\end{array}\right)+{\left(\begin{array}{c}\mathrm{\&rho;}sin\mathrm{\&theta;}cos\mathrm{\&alpha;}\\ \mathrm{\&rho;}sin\mathrm{\&theta;}sin\mathrm{\&alpha;}\\ -\mathrm{\&rho;}cos\mathrm{\&theta;}\end{array}\right)}_{Body}$

the phase center of the GPS generally does not completely coincide with the center of mass of the satellite, and a certain deviation exists between the two, and the coordinate of the phase center of the GPS in the satellite body coordinate system is assumed as follows:the rotation matrix between the satellite star sensor body coordinate system and the satellite body coordinate isThe star sensor measures the rotation matrix from the star sensitive body system to the J2000 coordinate system

Therefore, the strict geometric positioning formula of the laser altimetry satellite is as follows:

${\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_{ITRF}=R_{J2000}^{ITRF}\&lsqb;{\left(\begin{array}{c}{X}_{GPS}\\ Y_{GPD}\\ Z_{GPS}\end{array}\right)}_{J2000}+R_{Star}^{J2000}{R}_{Body}^{Star}\&lsqb;\left(\begin{array}{c}Lx\\ Ly\\ Lz\end{array}\right)+\left(\begin{array}{c}(\mathrm{\&rho;}+\mathrm{\&Delta;}\mathrm{\&rho;})\sin \mathrm{\&theta;}\cos \mathrm{\&alpha;}\\ (\mathrm{\&rho;}+\mathrm{\&Delta;}\mathrm{\&rho;})\sin \mathrm{\&theta;}\sin \mathrm{\&alpha;}\\ -(\mathrm{\&rho;}+\mathrm{\&Delta;}\mathrm{\&rho;})\cos \mathrm{\&theta;}\end{array}\right)-\left(\begin{array}{c}Dx\\ Dy\\ Dz\end{array}\right)\&rsqb;\&rsqb;$

where Δ ρ is the range correction due to atmospheric refraction and hardware errors.

The correspondence between the spatial rectangular coordinates (X, Y, Z) and the latitude and longitude (B, L) and the ground height H in the earth frame of reference coordinate system ITRF is shown in the following formula.

$\left\{\begin{array}{c}X=(N+H)\cos B\cos L\\ Y=(N+H)\cos B\sin L\\ Z=\&lsqb;N(1-e^2)+H\&rsqb;\sin B\end{array}\right.$

In the above formula, N is the radius of the unitary point-and-mortise circle, and e is the first eccentricity of the reference ellipsoid. The formula for calculating N is:

$N=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}$

wherein a is the major semi-axis of the reference ellipsoid.

The emitted waveform of the laser altimeter can be approximate to Gaussian pulse, and the echo waveform of the laser reflected by the ground surface can be approximately regarded as superposition of one or more Gaussian pulses. Fig. 4A and 4B schematically show a transmission waveform and a reception waveform of the laser altimeter, respectively.

And for the multi-echo waveform data, fitting by adopting a mode of superposing a plurality of Gaussian functions, wherein a fitting formula is shown as the following formula.

$w\left(t\right)=\mathrm{\&epsiv;}+\underset{m=1}{\overset{p}{\&Sigma;}}{A}_m{e}^{(-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2})}---\left(2.14\right)$

In the above formula, t is time, A_m,t_m,σ_mThe amplitude, the mean and the standard deviation of the mth Gaussian function are respectively waveform noise values.

The 3 parameters of the Gaussian function model of the echo waveform can be represented in a vector form as

c_m＝[A_m,t_m,σ_m](2.15)

The unknowns can be represented as X by a vector and the dimension of X is 3p + 1.

X＝[,c₁,c₂,…,c_p]^T

Setting the number of samples of the echo waveform as N, taking GLAS as an example, where N is 544, the actual echo value of each sample point is expressed as: r ═ R₁,r₂,...,r_N]^T(2.16)

The value of the waveform model corresponding to each sampling point calculated by equation (2.14) is

W＝[w₁,w₂,...,w_N]^T(2.17)

Defining the error l of the sampled value and the calculated value of the waveform model as

F(C)＝W-R(2.18)

To determine the unknown parameter vector C at the optimization of the model, i.e. f (C) ═ 0;

for the waveform model of (2.14), the unknown parameter vector C partial derivative is obtained:

$A=\frac{\&part;W}{\&part;C}=\left[\begin{array}{cccc}\frac{\&part;w_1}{\&part;\mathrm{\&epsiv;}}& \frac{\&part;w_1}{\&part;c_1}& \mathrm{...}& \frac{\&part;w_1}{\&part;c_p}\\ \frac{\&part;w_2}{\&part;\mathrm{\&epsiv;}}& \frac{\&part;w_2}{\&part;c_1}& \mathrm{...}& \frac{\&part;w_2}{\&part;c_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \frac{\&part;w_N}{\&part;\mathrm{\&epsiv;}}& \frac{\&part;w_N}{\&part;c_1}& \mathrm{...}& \frac{\&part;w_N}{\&part;c_p}\end{array}\right]---\left(2.19\right)$

the waveform model is developed using the first order taylor equation:

$w_n^{j+1}=w_n^j+\frac{\&part;w_n^j}{\&part;\mathrm{\&epsiv;}}{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}+\underset{m=1}{\overset{p}{\&Sigma;}}\frac{\&part;w_n^j}{\&part;c_m}{\mathrm{\&delta;}}_{c_m}---\left(2.20\right)$

wherein:

$\frac{\&part;w_n}{\&part;\mathrm{\&epsiv;}}=1$

$\frac{\&part;w_n}{\&part;c_m}=\&lsqb;\frac{\&part;w_n}{\&part;L},\frac{\&part;w_n}{\&part;t_m},\frac{\&part;w_n}{\&part;{\mathrm{\&sigma;}}_m}\&rsqb;$

$\frac{\&part;w_n}{\&part;A_m}=\exp \&lsqb;-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2}\&rsqb;$

$\frac{\&part;w_n}{\&part;t_m}=\frac{A_m(t-t_m)}{{{\mathrm{\&sigma;}}_m}^2}\exp \&lsqb;-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2}\&rsqb;=\frac{A_m(t-t_m)}{{{\mathrm{\&sigma;}}_m}^2}\frac{\&part;w_n}{\&part;A_m}$

$\frac{\&part;w_n}{\&part;{\mathrm{\&sigma;}}_m}=\frac{A_i{(t-t_i)}^2}{{{\mathrm{\&sigma;}}_m}^3}\exp \&lsqb;-\frac{{(t-t_i)}^2}{2{{\mathrm{\&sigma;}}_m}^2}\&rsqb;=\frac{(t-t_m)}{{\mathrm{\&sigma;}}_m}\frac{\&part;w_n}{\&part;t_m}---\left(2.21\right)$

with respect to the formula (2.18),

f (c) ═ w (c) + Ax-R,

wherein W (C) is an echo approximate value calculated according to the parameter C to be estimated and the model (2.14).

The constraint conditions are as follows: f (C) is 0, then

Ax＝L

L＝R-W(C)(2.22)

Then x is ═ a^TA]^-1A^TL, whereinConsidering the weights and the prior values, the formula can be written as:

x＝[A^TPA+V₀]^-1A^TPL(2.23)

wherein: p is a weight matrix, V₀Is a matrix of prior values, P_ij＝p_{i ij}，[V₀]_jk＝p_{ck jk}，p_iIs the weight of the ith sample value, p_ckIs the k-th vector c_kA priori value weight of.

Wherein,_ijas a function of kronecker:

${\mathrm{\&delta;}}_{ij}=\left\{\begin{array}{c}0,i\&NotEqual;j\\ 1,i=j\end{array}\right.$

after x is calculated, new parameters can be obtained:

Cⁱ⁺¹＝Cⁱ+x(2.24)

and judging | x | < d, wherein d is a set threshold value. If so, stopping iteration and outputting a parameter C; if not, the new parameter C is adopted, and the next iteration operation is carried out by returning to (2.18).

By combining the footprint image cloud detection and the echo waveform data processing, the data processing flow of the earth observation laser height measurement satellite can be described as follows:

1) smoothing and denoising the transmitting and receiving waveforms of the laser, extracting waveform characteristic parameters, determining the time corresponding to the gravity centers of the transmitting and receiving waveforms, and calculating the initial distance value of the unidirectional laser transmission according to the laser transmission time interval;

2) calculating rough three-dimensional coordinates of the laser footprint point according to the attitude and orbit parameters of the satellite, the laser emission time, the initial distance value and the laser geometric positioning model;

3) re-sampling to produce a footprint image geometric rough correction product according to the satellite attitude and orbit parameters, the footprint image, the footprint camera parameters and the footprint image geometric positioning model;

wherein, the waveform characteristic parameters, the rough three-dimensional coordinates of the laser footprint points of the laser initial distance values and the rough correction products of the geometry of the footprint images obtained in the steps 1 to 3 are basic products, all laser data points need to be processed,

4) processing the footprint image by adopting a cloud and fog detection model to obtain a footprint image cloud and fog detection product;

5) if the cloud amount in the footprint image cloud detection product exceeds the standard (if the ratio exceeds 20%), returning to the step 4 to process the next laser point; if the cloud amount does not exceed the standard, entering the step 6;

6) obtaining an atmospheric delay correction value according to atmospheric parameters and the rough three-dimensional coordinates of the laser footprint point in the step 2 by adopting an atmospheric delay correction model;

7) calculating an accurate distance value according to the laser ranging systematic error value obtained by geometric calibration, the initial distance value in the step 1 and the atmospheric delay correction value in the step 6;

8) calculating the three-dimensional coordinates of the laser footprint point according to the attitude and orbit parameters of the satellite, the laser emission time, the accurate distance value and a laser geometric positioning model;

9) calculating a tide correction value by adopting a tide correction model, and correcting the result in the step 8 to obtain the accurate three-dimensional coordinates of the laser footprint point;

10) producing a footprint image orthorectification product by adopting topographic data according to the satellite attitude and orbit parameters, the footprint image, the footprint camera parameters and the footprint image geometric positioning model;

wherein the footprint image cloud and fog detection product, the atmosphere delay correction value, the accurate distance value, the tide correction value, the accurate three-dimensional coordinate of the laser footprint point and the footprint image orthorectification product obtained in the 4-10 steps are all standard products,

11) judging the footprint image cloud and fog detection product with the cloud and fog amount not exceeding the standard in the step 5, and entering the step 12 if no cloud exists; otherwise, entering the step 17;

12) judging the characteristic constraint conditions of the echo waveform characteristic parameters (for example, adopting the echo waveform fitting), and if the characteristic constraint conditions are met, entering the step 13; otherwise, entering the step 16;

13) combining the laser footprint point accurate three-dimensional coordinates and the footprint image orthorectification products in the 9 th step and the 10 th step to form laser elevation control points;

14) performing elevation change analysis by using laser elevation reference points obtained by multi-period measurement in the same area to obtain a polar ice cover change monitoring product;

15) constructing a database by using laser elevation reference points obtained by multi-track measurement to form a laser elevation control point database product;

16) extracting and analyzing forestry parameters from the laser echo waveform characteristic parameters to obtain forestry special-subject products;

17) and (5) analyzing the cloud cover, and combining the results of the 9 th step and the 10 th step to obtain other special products.

Wherein, the polar ice cover change monitoring products, laser elevation control point library products, forestry special products and other special products generated in the steps 14 to 17 belong to special high-grade products generated by further processing.

Fig. 5 schematically shows a process flow of the satellite-borne laser altimetry data according to one embodiment of the invention. For example, a footprint image after system geometric correction is established according to the satellite attitude and orbit parameters and the basic parameters of the camera, and the approximate position of the laser spot in the footprint image is determined according to the satellite attitude and orbit parameters, the laser pointing parameters and the laser initial ranging value.

Figure 6 schematically shows the hierarchy of laser altimetry satellite data products. The concrete description is as follows:

level _0 (Level 0) raw data product: the decoded original data downloaded from the laser height measurement satellite comprise laser footprint images, satellite attitude measurement data, satellite orbit measurement data, laser emission waveform data, laser echo waveform data, relevant hardware parameters and the like. The raw data products are not generally provided externally and are only used by the ground data processing system.

Level _1 (Level 1) base product: the method is a product which is subjected to system geometric correction and primary waveform processing aiming at an original data product, and comprises characteristic parameters of a waveform, rough laser distance parameters, a rough footprint image geometric correction product, rough laser footprint three-dimensional coordinates and the like, wherein the attitude and orbit parameters after post-precision processing are classified as basic products. The basic product is not provided for the outside generally, and is only used by a ground data processing system.

A schematic production flow of the basic product is shown in fig. 7, for example. The production flow shown in fig. 7, the obtained product is: waveform characteristic parameters, laser initial distance values, laser footprint point rough three-dimensional coordinates and footprint image geometric rough correction products. In the production stage, according to the production requirements of the product, as long as the format is complete and the data is not lost, the product is produced, only the quality inspection is added later, and the mark with poor data quality is made and is not issued.

Level _2 (Level 2) standard product: the product is processed by accurate geometric calibration, atmospheric correction, tide correction and the like, and comprises atmospheric correction parameters, tide correction parameters, footprint image cloud detection products, footprint image orthorectification products and laser footprint accurate three-dimensional coordinate products. Standard products may be released to the user.

A schematic production flow of a standard product is shown in fig. 8, for example. As shown in the production flow of fig. 8, the obtained product is: atmospheric delay correction parameters, tide correction parameters, accurate distance values, laser footprint three-dimensional coordinates, footprint image cloud detection products and footprint image orthorectification products.

Level _3 (Level 3) thematic product/high-Level product: the high-grade laser height measurement satellite products are developed aiming at specific requirements of users on the basis of standard products, and comprise laser height control point library products, forestry special-purpose products, polar ice cover monitoring products, other special-purpose products and the like. Topical products may be released to users.

A schematic production flow of the special/high-grade product is shown in fig. 9, for example. As shown in the production flow of fig. 9, the obtained product is: polar region ice cover monitoring products, laser elevation control point library products, forestry special products and other special products.

The laser ranging value validity evaluation method can be used for subsequent processing of level 1 and level 2 data (products); the foregoing multi-parameter constrained laser elevation control point extraction method can be used for subsequent processing of level 2 data (product).

The research of utilizing the laser altimetry satellite to acquire the global elevation control points is developed, and a technical means for acquiring the overseas land elevation control points can be provided for China at the present stage; providing a data guarantee of elevation control points for global surveying and mapping in China; in addition, the method can be used for acquiring and improving the elevation control precision of ground control points in China in a large range and providing calibration data for other remote sensing equipment; and finally, the positioning precision of the domestic remote sensing satellite data can be further improved, and the application field of the domestic satellite data is widened.

According to the cloud amount detection method for the earth observation laser footprint image, the earth observation laser elevation control point extraction method based on cloud amount detection and echo waveform processing and the laser height measurement satellite data processing method, the influence of a cloud layer on laser distance measurement can be reduced, the precision of the laser distance measurement is ensured, and the accuracy of laser elevation reference data is effectively improved.

The above description is intended to be illustrative of the present invention and not to limit the scope of the invention, which is defined by the claims appended hereto.

Claims (10)

1. A cloud amount detection method for a laser footprint image is characterized by comprising the following steps:

carrying out sample training by utilizing a large number of cloud and fog-containing images to obtain an optimal gray threshold value and a texture characteristic value of the cloud and fog;

calculating a gray level histogram of the whole footprint image, and preliminarily judging whether cloud exists or not and the content of the cloud based on the preferred gray level threshold;

calculating the gray average value of the sub-blocks of the footprint image, judging that the sub-blocks contain cloud when the gray average value of the sub-blocks is higher than a first threshold value, and judging that the sub-blocks do not contain cloud when the gray average value of the sub-blocks is lower than a second threshold value;

calculating texture features of the image sub-blocks with the gray mean value between the first threshold value and the second threshold value by utilizing a gray co-occurrence matrix, and comparing the texture features with the cloud and fog texture feature values obtained from the sample training;

and counting the total number of the pixels judged as cloud and the number ratio of the pixels in the whole footprint image, and determining the cloud amount.

2. The cloud amount detection method of claim 1, wherein when calculating the gray histogram of the entire footprint image, histogram equalization enhancement processing is performed on the image, and the equalization formula is as follows:

$s_k=T\left(r_k\right)=\underset{i=0}{\overset{k}{\&Sigma;}}{P}_r\left(r_i\right)=\underset{i=0}{\overset{k}{\&Sigma;}}\frac{n_i}{n}$

0≤r_k≤1,k＝0,1,2,...L-1

wherein s is_kThe gray value of the pixel with the gray value k of the original image is converted to a new gray value; p_r(r_i) Is the pixel frequency with a gray value of i; l is the gray level of the image.

3. The cloud amount detection method of claim 1, wherein said first threshold is 235 and said second threshold is 80.

4. The cloud amount detection method of claim 1, wherein the texture features considered when computing the texture features of the image sub-blocks using a gray level co-occurrence matrix comprises: one or more of angular second moment, homogeneity, contrast and correlation.

5. A method for evaluating the effectiveness of a ground observation laser ranging value is characterized by comprising the following steps:

performing system geometric correction and preliminary waveform processing on decoded original data downloaded from a laser altimeter satellite, obtaining a footprint image corresponding to a laser spot according to a geometric corresponding relation, and realizing basic registration of laser and the footprint image based on hardware parameters and geographic coordinates or preliminary calibration of a footprint camera and a laser altimeter;

determining the position of the central point of the obtained footprint image corresponding to the laser spot, and determining the area corresponding to the laser spot in the footprint image;

cutting the determined image or storing the image into a memory to form a footprint image block;

carrying out cloud amount detection on the footprint image block by adopting the cloud amount detection method according to any one of claims 1-4 to obtain a cloud amount value of the footprint image block;

determining validity or availability of a laser ranging value for the footprint image block based on the cloud and fog magnitude value.

6. A laser elevation control point extraction method is characterized by comprising the following steps:

the method for evaluating the effectiveness of the earth observation laser ranging values according to claim 5 is adopted, the footprint image blocks judged to be cloudless are reserved, and the laser elevation data of the footprint image blocks judged to be thin clouds or thick clouds are removed;

by reflection of laser foot-printsRate of change<1 as a screening parameter, screening the laser footprint points of said retained footprint image block, wherein,

and selecting a laser footprint point which has only one peak in the echo waveform, the peak value is greater than a threshold value, and the standard deviation sigma after the waveform is fitted is less than or equal to 3.2ns, wherein the parameters for determining the threshold value comprise the transmitting energy and the receiving aperture of the laser.

7. The laser elevation control point extraction method according to claim 6,

and performing waveform fitting on the echo waveform by adopting a mode of superposing a plurality of Gaussian functions, wherein a fitting formula is shown as the following formula:

$w\left(t\right)=\mathrm{\&epsiv;}+\underset{m=1}{\overset{p}{\&Sigma;}}{A}_m{e}^{(-\frac{{(t-t_m)}^2}{2{{\mathrm{\&sigma;}}_m}^2})}$

in the above formula, t is time, A_m,t_m,σ_mThe amplitude, the mean and the standard deviation of the mth Gaussian function are respectively waveform noise values.

8. A method for processing data of a ground observation laser altimetry satellite is characterized by comprising the following steps:

smoothing and denoising the transmitting and receiving waveforms of the laser, extracting waveform characteristic parameters, determining the time corresponding to the gravity centers of the transmitting and receiving waveforms, and calculating the initial distance value of the unidirectional laser transmission according to the laser transmission time interval;

calculating rough three-dimensional coordinates of the laser footprint point according to the attitude and orbit parameters of the satellite, the laser emission time, the initial distance value and the laser geometric positioning model;

re-sampling to produce a footprint image geometric rough correction product according to the satellite attitude and orbit parameters, the footprint image, the footprint camera parameters and the footprint image geometric positioning model;

the cloud amount detection method according to any one of claims 1 to 4 is adopted to carry out cloud amount detection on the footprint image, and the laser footprint point of the footprint image of which the detected cloud amount value does not exceed a preset standard is subjected to further processing.

9. The earth observation laser altimetry satellite data processing method of claim 8, further comprising:

obtaining an atmospheric delay correction value according to atmospheric parameters and the rough three-dimensional coordinates of the laser footprint point by adopting an atmospheric delay correction model;

calculating an accurate distance value according to the laser ranging systematic error value obtained by geometric calibration, the initial distance value and the atmospheric delay correction value;

calculating the three-dimensional coordinates of the laser footprint point according to the attitude and orbit parameters of the satellite, the laser emission time, the accurate distance value and a laser geometric positioning model;

and calculating a tide correction value by adopting a tide correction model, and correcting the three-dimensional coordinates of the calculated laser footprint point to obtain the accurate three-dimensional coordinates of the laser footprint point.

10. The earth observation laser altimetry satellite data processing method of claim 9, further comprising:

producing a footprint image orthorectification product by adopting topographic data according to the satellite attitude and orbit parameters, the footprint image, the footprint camera parameters and the footprint image geometric positioning model;

if the foot print image with the cloud and fog quantity value not exceeding the preset standard is judged to be cloud-free, characteristic constraint condition judgment is carried out on the echo waveform characteristic parameters;

and if the characteristic constraint condition is met, combining the laser footprint point accurate three-dimensional coordinate and the footprint image orthorectification product to form a laser elevation control point.