Time normalization method for thermal infrared remote sensing data of unmanned aerial vehicle

Abstract

The invention discloses a time normalization method for thermal infrared remote sensing data of an unmanned aerial vehicle, and belongs to the technical field of unmanned aerial vehicle remote sensing. The method comprises the following steps: performing temperature drift removal processing on a thermal infrared remote sensing image sequence, and then splicing the thermal infrared remote sensing image sequence into a brightness temperature orthographic mosaic graph; estimating equivalent acquisition time information of each pixel in the brightness temperature orthographic mosaic graph; performing earth surface classification on the measurement area, and obtaining earth surface classification sub-images in one-to-one correspondence with the images in the original thermal infrared image sequence; using the surface classification sub-images as masks, extracting brightness temperature differences between dominant ground features and non-dominant ground features in the thermal infrared images, and constructing a ground feature brightness temperature difference change curve covering the duration of the whole flight mission through an interpolation method; and correcting the preliminarily obtained brightness temperature orthographic mosaic diagram by adopting a temperature difference compensation strategy so as to obtain brightness temperature data with consistent time. The temperature difference information in the original thermal infrared observation data is fully utilized, the requirement of introducing a large amount of additional auxiliary data is effectively avoided, and errors are reduced.

Description

Technical Field

[0001] This application relates to the field of UAV remote sensing, and particularly to a method for temporal normalization of UAV thermal infrared remote sensing data. Background Art

[0002] UAV thermal infrared remote sensing technology is a technical means of capturing the thermal radiation information emitted by the ground and objects through a UAV platform equipped with a thermal infrared sensor. With the high mobility, wide operation range and low cost - effectiveness of UAVs, this technology has shown great application potential in many fields such as environmental monitoring, urban heat island effect analysis, forest fire warning and agricultural crop health assessment, and has realized the rapid and accurate measurement of the surface temperature distribution and changes.

[0003] In field operations, UAVs usually need to execute multiple flight missions to obtain a complete thermal infrared image sequence covering the measurement area. However, this data acquisition process often takes a long time, ranging from dozens of minutes to several hours. During this period, the surface temperature is likely to change significantly, especially during the day or when the weather changes violently. On the other hand, users usually need to use a mosaic tool to stitch the collected thermal infrared image sequence into an ortho - mosaic image and further convert it into a temperature map for subsequent applications. However, due to the difference in the acquisition time of each image in the image sequence, the temperatures of each part (i.e., pixels) of the stitched ortho - mosaic image are not consistent in time, making it difficult to directly compare the temperature values of each part. Users usually expect to obtain ortho - temperature image data with consistent time to support relevant decisions. However, there is currently a lack of a method for temporal normalization of UAV thermal infrared remote sensing data during a single complete flight mission in the current research field, which may lead to large uncertainties in the data and thus restrict the further development of the technology. Summary of the Invention

[0004] The invention objective of this application is: aiming at the problem that the existing UAV thermal infrared remote sensing data is not unified in the time dimension, to provide a method for temporal normalization of UAV thermal infrared remote sensing data, so as to obtain a UAV thermal infrared ortho - mosaic image with consistent time, improve the comparability of the temperature values of each part of the image and reduce errors.

[0005] The technical solution adopted in this application is as follows:

[0006] A method for temporal normalization of UAV thermal infrared remote sensing data, the method includes the following steps:

[0007] Step 1: Perform preliminary normalization on the input original thermal infrared image sequence (denoted as image sequence 1) to remove temperature drift;

[0008] Step 2: Stitch the thermal infrared image sequence after removing temperature drift to generate a brightness temperature ortho-mosaic image covering the entire target area (denoted as ortho- Figure 1 );

[0009] Step 3: Estimate the equivalent acquisition time of each pixel in the ortho- Figure 1 to provide time parameters for subsequent time normalization;

[0010] Step 4: Obtain the complete land surface classification image of the target area (denoted as ortho- Figure 2 ). Taking the spatial reference of ortho- Figure 1 as the standard, georeference this complete land surface classification image to ensure their spatial consistency, and resample the georeferenced ortho- Figure 2 to match the resolution of ortho- Figure 1 ;

[0011] Step 5: Based on the position and attitude parameters of the original thermal infrared image sequence and the ortho- Figure 2 processed in Step 4, obtain the land surface classification sub-images corresponding to each frame image in the image sequence 1 input in Step 1 (denoted as image sequence 2);

[0012] Step 6: Based on the ratio of the area occupied by each land cover type in the target area to the total area of the target area, classify the land cover type with the highest area ratio as the dominant land cover type, and the remaining land cover types as secondary land cover types (i.e., non-dominant land cover types);

[0013] Using image sequence 2 as a mask, extract the brightness temperature corresponding to each land cover type in image sequence 1 and calculate the brightness temperature difference between the dominant land cover type and each secondary land cover type in each thermal infrared image (i.e., the initially calculated brightness temperature difference);

[0014] Step 7: Perform interpolation on the initially calculated brightness temperature difference to obtain the brightness temperature difference change curves of different land cover types covering the UAV flight mission period;

[0015] Step 8: Use the temperature difference compensation strategy to compensate the temperature values of the areas belonging to non-dominant land cover type pixels in the ortho- Figure 1 to obtain a brightness temperature ortho-mosaic image that completely covers the target area and has unified time (denoted as ortho- Figure 3 ).

[0016] Furthermore, in Step 1, during the preliminary normalization, estimate the temperature drift amount of each frame image in the original thermal infrared image sequence based on the probability density function of the thermal infrared image gray value.

[0017] Furthermore, in Step 1, the preliminary normalization is specifically:

[0018] According to the formula RDN i= LM(DN o,i ) calculates the gray value corresponding to the location with the maximum probability density of the upper dominant feature type in any frame image of the image sequence 1, where LM() represents a function for finding the local probability maximum of the dominant feature type, i represents the image serial number, and DN o,i represents the original gray value matrix corresponding to the i-th frame image;

[0019] Randomly select a reference image from the original thermal infrared image sequence, and calculate the gray value RDN of each non-reference image in the original thermal infrared image sequence j and the gray value RDN of the reference image * of the gray difference ΔRDN j , where j is the image serial number of each non-reference image in the original thermal infrared image sequence;

[0020] For each non-reference image, based on its original gray value matrix DN o,j and the gray difference ΔRDN j to obtain its gray value matrix DN after temperature drift normalization r,j , and directly use the original gray value matrix of the reference image as the gray value matrix DN after temperature drift normalization r * ;

[0021] Based on the ground observation data, determine the temperature drift amount DN of the reference image base ;

[0022] For each frame image of the original thermal infrared image sequence, according to the formula DN c,i = DN r,i - DN base to obtain the gray value matrix DN of each frame image after temperature drift correction c,i , to obtain the thermal infrared image sequence after removing the temperature drift, where DN r,i represents the gray value matrix after temperature drift normalization of the i-th frame image in the original thermal infrared image sequence.

[0023] Further, in step 1, use the first frame image in the original thermal infrared image sequence as the reference image.

[0024] Further, in step 1, determine the temperature drift amount DN of the reference image based on the ground observation data base Specifically:

[0025] For the original gray value matrix of the reference image, according to the conversion formula between gray and bright temperature, convert it to the corresponding bright temperature value to obtain the bright temperature value of the reference image;

[0026] The measured brightness temperature of the reference image is obtained based on the ground observation data, and the temperature drift value is obtained based on the difference between the measured brightness temperature and the brightness temperature value of the reference image;

[0027] According to the conversion formula of grayscale and brightness temperature, the temperature drift value is converted into the corresponding grayscale drift value to obtain the temperature drift value DN of the reference image. base .

[0028] Furthermore, the conversion formula between grayscale and brightness temperature is: b =k1×DN+k0, where k0 and k1 represent the offset value and gain value in the temperature conversion coefficient respectively, T b represents brightness temperature, DN represents the brightness temperature T b The gray value of .

[0029] Furthermore, in step 3, the orthophoto is estimated Figure 1 The equivalent acquisition time of each pixel is as follows:

[0030] Ortho Figure 1 For any pixel of the image, each thermal infrared image in the thermal infrared image sequence corresponding to the pixel is used as the contribution image of the current pixel;

[0031] According to the formula Calculate the equivalent acquisition time of the current pixel, where t j Indicates the acquisition time of the j-th contribution image and the weight of the j-th contribution image x and y represent the distance between the current pixel and the image center of the current contribution image in the horizontal and vertical coordinate directions, W and H represent the width and height of the current contribution image, and N represents the total number of contribution images corresponding to the current pixel.

[0032] Furthermore, in step 4, a complete surface classification image of the target area can be obtained based on a deep learning algorithm;

[0033] Furthermore, in step 4, the types of objects in the target area can be set to five types: wetland, water body, road, grassland, and building.

[0034] Furthermore, in step 6, the brightness temperature difference between the dominant ground object type and each secondary ground object type in each thermal infrared image is calculated as follows:

[0035] Calculate the representative values ​​of brightness temperature of the dominant land object type and each secondary land object type according to the formula;

[0036]

[0037] Among them, T b-dom-i represents the brightness temperature of the ith pixel belonging to the dominant ground object type; Md represents the total number of pixels of the dominant ground object type in the current thermal infrared image, Tb (LC_dom) represents the representative value of the brightness temperature of the dominant land cover type in the current thermal infrared image; T b-non-i represents the brightness temperature of the i-th pixel belonging to a certain secondary land cover type; Mn represents the total number of pixels of this secondary land cover type in the current thermal infrared image; T b (LC_non) represents the representative value of the brightness temperature of this secondary land cover type in the current thermal infrared image;

[0038] Based on the difference between T b (LC_dom) and each T b (LC_non) in each thermal infrared image, the brightness temperature difference between the dominant land cover type and each secondary land cover type in each thermal infrared image is obtained.

[0039] Furthermore, in step 7, the linearly interpolated brightness temperature difference is interpolated by using the piecewise linear interpolation method.

[0040] Furthermore, step 8 specifically includes:

[0041] Adopt the temperature difference compensation strategy to compensate the temperature value of the area belonging to the non-dominant land cover pixels in the ortho Figure 1 to obtain an ortho-rectified mosaic image of brightness temperature that completely covers the target area and has a unified time (denoted as ortho Figure 3 ).

[0042] Step 8: Obtain the ortho-rectified mosaic image of brightness temperature after time normalization through temperature difference compensation;

[0043] Calculate the brightness temperature of a certain secondary land cover type pixel after removing temperature drift and temperature difference compensation according to the formula:

[0044] T b-norm (t eq →t tg ,LC_non)

[0045] =T b-r (t eq →t tg ,LC_non)+[δT b (t tg ,LC_non)-δT b (t eq ,LC_non)]

[0046] Wherein, t tg represents the target time after time normalization, t eq represents the equivalent acquisition time of the current pixel, T b-r (t eq →t tg, (LC_non) represents the brightness temperature of a certain secondary feature type pixel after temperature drift correction (i.e., the brightness temperature obtained after conversion from DN r,i ), and δT b (t tg , LC_non), δT b (t eq , LC_non) respectively represent the brightness temperature difference between a certain secondary feature type pixel and the dominant feature type at the target time t tg and the equivalent acquisition time t eq , that is, the representative value of the brightness temperature calculated in step 6.

[0047] The technical solution provided by this application at least brings the following beneficial effects:

[0048] A time normalization method for UAV thermal infrared remote sensing data provided by this application can perform time normalization on the thermal infrared ortho-mosaic map, so that the temperature values of each pixel in the image reflect the surface temperature state at the same moment. After time normalization correction, the comparability of the temperature data is significantly improved, and the error is significantly reduced. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] The above and / or additional aspects and advantages of this application will become obvious and easy to understand from the following description of the embodiments in conjunction with the drawings, where:

[0050] Figure 1 is a schematic flow chart of a time normalization method for UAV thermal infrared remote sensing data provided by an embodiment of this application;

[0051] Figure 2 is the brightness temperature ortho-mosaic map provided by an embodiment of this application. Among them, (a) is the image before removing temperature drift, and (b) is the image after removing temperature drift. The rectangular parts (A, B, C, D) in the figure are the selected sub-regions for showing the processing effect;

[0052] Figure 3 is the estimated result map of the equivalent acquisition time of pixels provided by an embodiment of this application;

[0053] Figure 4 is the surface classification image corresponding to the brightness temperature ortho-mosaic map provided by an embodiment of this application;

[0054] Figure 5 is the schematic diagram of the surface classification image sequence provided by an embodiment of this application;

[0055] Figure 6 is the brightness temperature difference curve of different surface coverage types provided by an embodiment of this application;

[0056] Figure 7The orthorectified mosaic map of brightness temperature after time normalization provided by the embodiment of the present application; where (a) is the temperature compensation map, (b) is the time normalization result, and the rectangular parts (A, B, C, D) in the figure are the selected sub-regions for demonstrating the processing effect;

[0057] Figure 8 The box plot of temperature compensation for different ground object types provided by the embodiment of the present application;

[0058] Figure 9 The verification result map based on ground observation data provided by the embodiment of the present application; where (a) is the change curve of the brightness temperature difference between roads and wetland types, and (b) is the change curve of the brightness temperature difference between grasslands and wetlands. Detailed implementation manners

[0059] To enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions of the embodiments of the present application will be described in detail and completely below in conjunction with the accompanying drawings in the embodiments of the present application. Obviously, the embodiments described by referring to the accompanying drawings are exemplary and are intended to explain the present application, rather than being construed as a limitation to the present application.

[0060] At present, the orthorectified mosaic map of surface temperature obtained by unmanned aerial vehicle (UAV) thermal infrared remote sensing usually does not consider the problem of non-uniformity in the time dimension of each part of the image, making the temperature values of each part of the image lack comparability and have a large uncertainty. Starting from the original thermal infrared image sequence, the present application considers the influence of time variation on the temperature of different ground objects through a temperature difference compensation strategy and realizes time normalization.

[0061] The embodiment of the present application discloses a time normalization method for UAV thermal infrared remote sensing data. The method includes: First, perform temperature drift removal processing on the thermal infrared remote sensing image sequence and splice these images into an orthorectified mosaic map of brightness temperature; Then, estimate the equivalent acquisition time information of each pixel in the orthorectified mosaic map of brightness temperature as the time reference for normalization; Then, perform surface classification on the measurement area and obtain the surface classification sub-images corresponding to each image in the original thermal infrared image sequence one by one; Subsequently, use these surface classification sub-images as masks to extract the brightness temperature differences between the dominant ground objects and non-dominant ground objects in each thermal infrared image, and construct the change curve of the brightness temperature difference of the ground objects covering the entire flight mission duration through interpolation; Finally, adopt a temperature difference compensation strategy to correct the initially obtained orthorectified mosaic map of brightness temperature to obtain time-consistent brightness temperature data. The method proposed in the embodiment of the present application makes full use of the temperature difference information in the original thermal infrared observation data and avoids the need to introduce a large amount of additional auxiliary data. This makes the UAV thermal infrared remote sensing temperature data more comparable and reduces errors, providing a more reliable basic data source for related applications.

[0062] In one embodiment, asFigure 1 As shown in Figure 1 , a method for temporal normalization of thermal infrared remote sensing data for unmanned aerial vehicles provided by an embodiment of the present application includes eight steps: removing temperature drift and stitching the original thermal infrared image sequence; obtaining a brightness temperature ortho-mosaic map after removing temperature drift; estimating the equivalent acquisition time of image pixels; obtaining a complete surface classification image of the measurement area; obtaining a surface classification sub-image corresponding to each frame of thermal infrared image; calculating the brightness temperature difference between the dominant land cover type and each secondary land cover type in each thermal infrared image; interpolating to obtain a brightness temperature difference change curve covering the UAV flight mission period; and obtaining a brightness temperature ortho-mosaic map after temporal normalization through temperature difference compensation. The specific implementation of each step is as follows:

[0063] Step 1: Removing temperature drift and stitching the original thermal infrared image sequence;

[0064] UAVs usually carry uncooled thermal imagers for temperature observation. However, such thermal imagers are prone to interference from the external environment during operation, resulting in the observed temperature values deviating significantly from the true temperature of the target. This phenomenon is called temperature drift.

[0065] Previous studies have shown that the temperature drift between thermal infrared images can be regarded as a set of constants. To remove this temperature drift, in the embodiment of the present application, the temperature drift amount is first estimated through the probability density function of the gray value of the thermal infrared image. Then, a reference image is selected, and the temperature drift levels of other images in the image sequence are normalized to the same level as the reference image. Finally, correction is performed using ground observation data to complete the removal of temperature drift from the original thermal infrared image sequence. This process can be expressed by formula (1):

[0066]

[0067] In the formula, LM() represents a function for finding the local probability maximum of the dominant land cover; i represents the image serial number; DN o,i represents the original gray value matrix corresponding to the i-th thermal infrared image; RDN i represents the gray value corresponding to the location where the probability density of the dominant land cover type on the i-th thermal infrared image is the largest; RDN1 represents the gray value corresponding to the location where the probability density of the dominant land cover type on the selected reference thermal infrared image (i.e., the reference image for temperature drift normalization) is the largest; ΔRDN i represents the difference between the RDN value of the i-th thermal infrared image and the RDN value of the reference thermal infrared image; DN r,i represents the gray value matrix of the i-th thermal infrared image after temperature drift normalization; DN base represents the temperature drift amount of the reference thermal infrared image determined by ground observation data (such as an infrared radiometer) (quantified by the gray value, which can be determined by subsequent formula (2)); DN c,iIt represents the grayscale value matrix of the i-th thermal infrared image after removing the temperature drift.

[0068] Step 2: Obtain the bright temperature ortho-mosaic map after removing the temperature drift;

[0069] The field of view of a single thermal infrared image is limited. To serve subsequent applications, the image sequence needs to be stitched into an ortho-mosaic map that completely covers the measurement area. Common automated mosaic software includes Pix4D mapper and Metashape, etc. In the embodiment of this application, Metashape software is taken as an example to automate the stitching of thermal infrared images. This software can automatically complete operations such as image alignment, point cloud extraction, digital elevation model construction, ortho-image generation, and index image conversion, and finally output a bright temperature ortho-mosaic map with geographic information. The conversion from the original grayscale value to the bright temperature value can be carried out through formula (2):

[0070] T b = k1×DN + k0 (2)

[0071] In the formula, k0 and k1 respectively represent two temperature conversion coefficients provided by the manufacturer. Among them, k0 represents the offset value, and k1 represents the gain value; T b represents the bright temperature (K).

[0072] As Figure 2 shown in (a) of [], the bright temperature ortho-mosaic map generated by using the original thermal infrared image sequence has obvious patch effects and poor image quality; while Figure 2 the bright temperature ortho-mosaic map after removing the temperature drift shown in (b) of [] has significantly improved quality, the patch effect is significantly weakened, and the overall temperature distribution is also more uniform.

[0073] Step 3: Estimate the equivalent acquisition time of the image pixels;

[0074] Since the time normalization process involves the quantization of the time information of each pixel, and the stitching process of the ortho-mosaic map usually does not consider the processing of time information. Therefore, in the embodiment of this application, the index of "equivalent acquisition time" is used to quantify the time information of the bright temperature mosaic image pixels. The calculation process of this index mainly considers the weighted average of the acquisition times of all "contributing images" of a certain pixel on the mosaic map, and the determination of the weight considers the actual fusion method of the images. The process can be described by formula (3):

[0075]

[0076] In the formula, j represents the serial number of the contributing image; N represents the total number of contributing images; w j () represents the weight of the j-th contributing image; t jIt represents the acquisition time (h) of the j-th contribution image; x and y respectively represent the distances (unit: m) between the current pixel on the brightness temperature mosaic map and the center of the current contribution image in the X direction and the Y direction; EAT represents the equivalent acquisition time (unit: h); W represents the width (unit: m) of the current contribution image; H represents the height (unit: m) of the current contribution image.

[0077] Figure 3 shows the estimation results of the equivalent acquisition time for the Figure 2 brightness temperature ortho-mosaic map in [reference], and it can be seen that the time generally shows the characteristics of "slow longitudinal change and obvious lateral change", which is consistent with the actual flight direction of the UAV, and the numerical range of the image also corresponds to the start and end times of the flight mission.

[0078] Step 4: Obtain the complete surface classification image of the measurement area;

[0079] Since the time normalization process involves the quantification of the temperature changes of various ground object types, surface classification is essential. The surface classification can be carried out in a supervised classification manner, using the maximum likelihood criterion or some deep learning algorithms. In the embodiment of the present application, the measurement area is classified by using a deep learning algorithm through Python tools. In order to utilize richer ground object spectral information, the algorithm first uses principal component transformation to fuse the multi-spectral data and thermal infrared data of the measurement area to generate a fused ortho-mosaic map; then manually annotates various ground object types on the mosaic map, including wetland, water body, road, grassland, and building, a total of five ground object types; finally, the Deeplv3 deep learning framework is used for high-precision ground object classification. If there is no corresponding multi-spectral image of the measurement area, the visible light image synchronously acquired by the thermal imager can also be used for classification, and the classification results can be manually post-processed subsequently to correct the obviously misclassified areas.

[0080] Figure 4 shows the surface classification results of the measurement area, and it can be seen that the wetland occupies the largest area proportion (being the dominant ground object type in this area), the water body type ranks second in terms of area proportion, while the road, grassland, and building have relatively small proportions. The overall classification result is good and conforms to the actual situation of the measurement area.

[0081] Step 5: Obtain the surface classification sub-images corresponding to each frame of thermal infrared images;

[0082] Since it is necessary to describe the difference in brightness temperature between non-dominant ground object types and dominant ground object types at multiple moments, it is necessary to extract the brightness temperature information of each ground object type on each frame of the thermal infrared image; the premise for this process is to construct a ground object type mask for each image in the thermal infrared image sequence, that is, a classified sub-image. In the embodiments of the present application, by combining the azimuth and attitude data of each thermal infrared image recorded by the thermal imager, a sub-classified image corresponding to each image in the thermal infrared image sequence is "intercepted" from the complete classified image of the measurement area through a rotation transformation. The image rotation transformation can be represented by formula (4):

[0083]

[0084] In the formula, θ represents the yaw angle of the unmanned aerial vehicle; x0 and y0 respectively represent the abscissa and ordinate before transformation; x' and y' respectively represent the abscissa and ordinate after transformation.

[0085] Figure 5 Part of the original thermal infrared images and their corresponding surface classified sub-images are shown. It can be seen that the classified images accurately describe the ground object coverage of each infrared image and maintain the same orientation.

[0086] Step 6: Calculate the difference in brightness temperature between the dominant ground object type and each secondary ground object type in each thermal infrared image;

[0087] By performing a dot product operation on the gray value matrix corresponding to the original thermal infrared image and the sub-classified image matrix obtained in step 5, a set of gray values corresponding to the dominant ground object type and non-dominant ground object types on each thermal infrared image can be obtained, and combined with formula (2), a set of brightness temperature values can be obtained. In order to suppress the errors introduced by the classification process and different observation azimuths of the thermal imager, the arithmetic mean of the brightness temperature values of each ground object type is taken as the representative brightness temperature value of this type of ground object type on this thermal infrared image. This process can be represented by formula (5):

[0088]

[0089] In the formula, T b-dom-i represents the brightness temperature (unit: K) of the i-th pixel belonging to the dominant ground object type; Md represents the total number of pixels of the dominant ground object type in this thermal infrared image; T b (LC_dom) represents the representative brightness temperature value (unit: K) of the dominant ground object type (LC_dom) in this thermal infrared image; T b-non-i represents the brightness temperature (unit: K) of the i-th pixel belonging to a certain non-dominant ground object type; Mn represents the total number of pixels of this non-dominant ground object type in this thermal infrared image; T b (LC_non) represents the representative brightness temperature value (unit: K) of this non-dominant ground object type (LC_non) in this thermal infrared image.

[0090] By performing the above operations on all thermal infrared images in the thermal infrared image sequence, the sequence of representative brightness temperature values of each ground object type can be solved, and then the sequence of differences in representative brightness temperature values between each non-dominant ground object type and the dominant ground object type can be obtained.

[0091] Step 7: Interpolate to obtain the change curve of the difference in brightness temperature between different ground objects covering the UAV flight mission period;

[0092] Since the coverage status of each ground object type in the measurement area is different, its frequency of appearance in the thermal infrared image is also different; considering that subsequent time normalization calculations may involve the difference in brightness temperature at any point during the flight process, it is necessary to interpolate and expand the difference sequence initially obtained in Step 7. Considering that the shooting interval of the UAV thermal imager is usually small, the piecewise linear interpolation method is used, and this process can be expressed by formula (6):

[0093]

[0094] In the formula, t start and t end respectively represent the start and end times of the interpolation interval; δT b (t, LC_non) represents the difference in brightness temperature between a certain non-dominant ground object type and the dominant ground object type at time t (t start , t end ).

[0095] In order to suppress extreme values, the interpolation sequence is processed by spline curve smoothing. Figure 6 Shows the difference curves of brightness temperature between different non-dominant ground object types (water body, road, grassland, building) and the dominant ground object type (wetland). The missing parts in the figure represent the intervals when the UAV stops flying. It can be seen that the interpolated curve better reflects the change trend of the temperature difference between different ground object types, and the temperature difference fluctuations between the water body, road, grassland and wetland are relatively smaller, while the temperature difference fluctuation between the building and the wetland is larger, which is jointly caused by the geometric characteristics, emissivity characteristics of different ground object types and their distribution in the measurement area.

[0096] Step 8: Obtain the ortho-mosaic map of brightness temperature after time normalization through temperature difference compensation;

[0097] In the process of removing temperature drift, a certain thermal infrared image in the thermal infrared image sequence is used as a reference, and the deviation correction method is adopted. This method has little impact on the accuracy of the dominant land cover type (wetland type in the embodiments of the present application), but for the remaining non-dominant land cover types, the difference in the surface temperature change rate will further introduce uncertainty in the correction process. Theoretically, the impact brought by the temperature change rate can be corrected by compensating for the difference in temperature change between the non-dominant land cover type from the current moment (i.e., the equivalent acquisition moment) to the target moment (i.e., the target moment of time normalization / temperature drift correction). This process can be represented by formula (7):

[0098]

[0099] In the formula, t tg represents the target moment of time normalization; t eq represents the equivalent acquisition moment of the current pixel; T b-r (t eq →t tg , LC_non) represents the brightness temperature of a pixel of a certain non-dominant land cover type after temperature drift correction; T b-norm (t eq →t tg , LC_non) represents the brightness temperature of a pixel of a certain non-dominant land cover type after time normalization (i.e., removing temperature drift and performing temperature difference compensation).

[0100] Figure 7 (a) in shows the temperature compensation status of each part on the brightness temperature ortho-mosaic map after removing temperature drift in the embodiments of the present application. It can be found that the temperature compensation value is closely related to the land cover type. Figure 7 (b) in shows the brightness temperature ortho-mosaic map after time normalization in the embodiments of the present application. It can be found that while maintaining good image quality, the temperatures of the same land cover type are closer, and the temperature difference is significantly reduced. Further observing Figure 7 sub-regions A, B, C, and D in (b) in, it can be found that the temperature change of the water body type is more continuous and natural, and the features of lakes and rivers that are discontinuous in the image after only removing temperature drift become clearly visible; in addition, the temperatures of the building type are also more consistent.

[0101] Figure 8 Further shows the box plot of the temperature compensation values of various land cover types. It can be found that the distribution status and dispersion degree of the temperature compensation values of different land cover types are also different. The water body type is mainly negative compensation (the average compensation is about -1.6K), while other land cover types are positive compensation (the average compensation is about 0.8 - 1.7K).

[0102] Figure 9(a) in it shows the verification results of the difference in brightness temperature between the road and the wetland based on the observation data of the ground infrared radiometer. It can be found that the temperature difference change of the uncorrected image and the image corrected only for temperature drift is a horizontal straight line, which completely fails to reflect the change of temperature difference of different surface types over time; while the temperature difference curve after time normalization is in good agreement with the measured data (the correlation coefficient is about 0.73), and can reflect the change of surface temperature difference to a large extent. Similarly, Figure 9 (b) in it shows the verification results of the difference in brightness temperature between the grassland and the wetland. The temperature differences of the uncorrected and the one corrected only for temperature drift are also constant values, while the curve after time normalization is in better agreement with the actual observed value (the correlation coefficient is about 0.79).

[0103] A time normalization method for unmanned aerial vehicle (UAV) thermal infrared remote sensing data provided by an embodiment of the present application is applicable to the time normalization of UAV thermal infrared data in multiple flight areas, and has good scalability and practicability.

[0104] In the description of this specification, the description referring to terms such as "one embodiment", "some embodiments", "example", "specific example", or "some examples" means that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present application. In this specification, the schematic expressions of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in a suitable manner in any one or more embodiments or examples. In addition, without contradiction, those skilled in the art can combine and combine the different embodiments or examples described in this specification and the features of different embodiments or examples.

[0105] In addition, the descriptions such as "first", "second", etc. are only for descriptive purposes, and cannot be understood as indicating or implying relative importance or implicitly indicating the quantity of the indicated technical features. Thus, the features defined with "first", "second", etc. can explicitly or implicitly include at least one of the features.

[0106] Any process or method description shown in the flowchart or described in other ways in this specification can be understood as representing a module, segment, or part of code including one or more executable instructions for implementing a customized logic function or process. The scope of the preferred embodiments of the present application includes additional implementations, where the functions can be executed in a manner that is not shown or discussed in the order, including in a substantially simultaneous manner according to the involved functions or in the reverse order, which should be understood by those skilled in the art of the embodiments of the present application.

[0107] It should be understood that each part of the present application can be implemented by hardware, software, firmware, or a combination thereof. In the above embodiments, multiple steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, any one or a combination of the following techniques well known in the art can be used: discrete logic circuits with logic gate circuits for implementing logic functions on data signals, application specific integrated circuits with appropriate combinational logic gate circuits, programmable gate arrays (PGAs), field programmable gate arrays (FPGAs), etc.

[0108] Those of ordinary skill in the art can understand that all or part of the steps carried by the method of implementing the above embodiments can be completed by instructing relevant hardware through a program, and the program can be stored in a computer-readable storage medium. When the program is executed, it includes one or a combination of the steps of the method embodiments.

[0109] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements on some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A time normalization method for unmanned aerial vehicle thermal infrared remote sensing data, characterized in that: The following steps are involved: Step 1: Perform preliminary normalization on the input raw thermal infrared image sequence to remove temperature drift; Step 2: Stitch the thermal infrared image sequence after removing the temperature drift to generate a brightness temperature ortho-mosaic image covering the entire target area; Step 3: Estimate the equivalent acquisition time of each pixel of the brightness temperature orthomosaic image obtained in step 2; Step 4: Obtain a complete land surface classification image of the target area, georeference the complete land surface classification image based on the spatial reference of the brightness temperature orthomosaic image, and resample the georeferenced complete land surface classification image to match the resolution of the brightness temperature orthomosaic image; Step 5: Based on the position and attitude parameters of the original thermal infrared image sequence and the complete surface classification image processed in step 4, obtain the surface classification sub-image corresponding to each frame image in the image sequence 1 input in step 1; Step 6: Based on the ratio of the area occupied by each feature type in the target area to the total area of ​​the target area, the feature type with the largest area share is classified as the dominant feature type, and the remaining feature types are classified as secondary feature types; Using the surface classification sub-image as a mask, the brightness temperature corresponding to each surface cover type in the original thermal infrared image sequence is extracted, and the brightness temperature difference between the dominant surface object type and each secondary surface object type in each thermal infrared image is calculated to obtain the preliminary brightness temperature difference; Step 7: Interpolate the preliminary brightness temperature difference to obtain the brightness temperature difference curve of different objects during the UAV flight mission; Step 8: Use the temperature difference compensation strategy to compensate the regional temperature values ​​of the pixels belonging to the secondary object type in the brightness temperature ortho-mosaic image to obtain a brightness temperature ortho-mosaic image that completely covers the target area and has a unified time.

2. The method according to claim 1, characterized in that In step 1, during the initial normalization, the temperature drift of each frame in the original thermal infrared image sequence is estimated based on the probability density function of the thermal infrared image grayscale value.

3. The method according to claim 1, characterized in that In step 1, the initial normalization is as follows: According to the formula RDN i =LM(DN o,i ) calculates the gray value corresponding to the maximum probability density of the dominant ground object type in any frame of the original thermal infrared image sequence, where LM() represents the function of finding the local probability maximum value of the dominant ground object type, i represents the image number, and DN o,i Represents the original gray value matrix corresponding to the i-th frame image; A reference image is randomly selected from the original thermal infrared image sequence, and the grayscale value RDN of each non-reference image in the original thermal infrared image sequence is calculated. j The gray value RDN of the reference image * Grayscale difference ΔRDN j , where j is the image sequence number of each non-reference image in the original thermal infrared image sequence; For each non-reference image, based on its original gray value matrix DN o,j Grayscale difference ΔRDN j The gray value matrix DN after the temperature drift is normalized is obtained by the difference r,j , and directly use the original gray value matrix of the reference image as the gray value matrix after temperature drift normalization Determine the temperature drift DN of the reference image based on ground observation data base ; For each frame of the original thermal infrared image sequence, according to the formula DN c,i =DN r,i -DN base Get the gray value matrix DN of each frame image after temperature drift correction c,i , we get the thermal infrared image sequence after removing the temperature drift, where DN r,i Represents the gray value matrix after normalization of the temperature drift of the i-th frame image in the original thermal infrared image sequence.

4. The method according to claim 3, characterized in that In step 1, the first frame image in the original thermal infrared image sequence is used as a reference image.

5. The method according to claim 3, characterized in that In step 1, the temperature drift DN of the reference image is determined from the ground observation data. base Specifically: The original gray value matrix of the reference image is converted into the corresponding brightness temperature value according to the conversion formula between gray value and brightness temperature to obtain the brightness temperature value of the reference image; The measured brightness temperature of the reference image is obtained based on the ground observation data, and the temperature drift value is obtained based on the difference between the measured brightness temperature and the brightness temperature value of the reference image; According to the conversion formula of grayscale and brightness temperature, the temperature drift value is converted into the corresponding grayscale drift value to obtain the temperature drift value DN of the reference image. base .

6. The method according to claim 1, characterized in that In step 3, the equivalent acquisition time of each pixel of the brightness temperature orthomosaic image is estimated as follows: For any pixel of the brightness temperature orthomosaic image, each thermal infrared image in the thermal infrared image sequence corresponding to the pixel is used as the contribution image of the current pixel; According to the formula Calculate the equivalent acquisition time of the current pixel, where t j Indicates the acquisition time of the j-th contribution image and the weight of the j-th contribution image x and y represent the distance between the current pixel and the image center of the current contribution image in the horizontal and vertical coordinate directions, W and H represent the width and height of the current contribution image, and N represents the total number of contribution images corresponding to the current pixel.

7. The method according to claim 1, characterized in that In step 4, a complete surface classification image of the target area can be obtained based on the deep learning algorithm.

8. The method according to claim 1, characterized in that In step 6, the brightness temperature difference between the dominant ground object type and each secondary ground object type in each thermal infrared image is calculated as follows: Calculate the representative values ​​of brightness temperature of the dominant land object type and each secondary land object type according to the formula; Among them, T b-dom-i represents the brightness temperature of the ith pixel belonging to the dominant ground object type; Md represents the total number of pixels of the dominant ground object type in the current thermal infrared image, T b (LC_dom) represents the brightness temperature representative value of the dominant ground object type in the current thermal infrared image; T b-non-i represents the brightness temperature of the i-th pixel belonging to a certain secondary object type; Mn represents the total number of pixels of this secondary object type in the current thermal infrared image; T b (LC_non) represents the brightness temperature representative value of this type of secondary features in the current thermal infrared image; Based on T in each thermal infrared image b (LC_dom) and each T b The difference between the brightness temperatures of the dominant land object type and each secondary land object type in each thermal infrared image is obtained by calculating the difference between the brightness temperatures of the dominant land object type and each secondary land object type in each thermal infrared image.

9. The method according to claim 1, characterized in that In step 7, the initially calculated brightness temperature difference is interpolated using piecewise linear interpolation.

10. The method according to claim 1, characterized in that Step 8 specifically includes: The brightness temperature of a certain type of pixel of a minor land object is calculated by removing temperature drift and temperature difference compensation according to the formula: T b-norm (t eq →t tg ,LC_non) =T b-r (t eq →t tg ,LC_non)+[δT b (t tg ,LC_non)-δT b (t eq ,LC_non)] Among them, t tg represents the time-normalized target moment, t eq Indicates the equivalent acquisition time of the current pixel, T b-r (t eq →t tg ,LC_non) represents the brightness temperature of a certain primary feature type pixel after temperature drift correction, δT b (t tg ,LC_non)、δT b (t eq , LC_non) represent the target time t tg , equivalent acquisition time t eq The brightness temperature difference between a pixel of a minor land object type and the dominant land object type.

Images