A universal method for acquiring and modeling BRDF information from multi-angle remote sensing images from UAVs

Abstract

The present invention proposes a method for extracting and modeling the bidirectional reflectance distribution function (BRDF) of multi-angle images from unmanned aerial vehicles. The method is as follows: first, radiometric correction is performed on the unmanned aerial vehicle images acquired in any multi-angle observation mode; second, the overlapping images are converted into dense point clouds, digital surface models, and orthophotos; then, a grid of a certain size is set for the DSM, and the three-dimensional coordinates of each corner point in the grid are obtained in the earth projection coordinate system. Based on the projection transformation, the three-dimensional coordinates of the grid corner points are converted into the two-dimensional pixel coordinates of the original image; then, the observation zenith angle, azimuth angle, and average value BRF are calculated grid by grid based on the waypoint position of the original image and the three-dimensional coordinates of the grid corner points to obtain the multi-angle BRF of each grid; finally, the multi-angle BRF of each grid is fitted band by band based on the kernel-driven model. This method is applicable to frame-type remote sensing images acquired in any multi-angle observation mode of unmanned aerial vehicles, and provides important technical support for the flexible and convenient acquisition of BRDF information of ground objects.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing image preprocessing, and in particular to a general method for acquiring and modeling multi-angle remote sensing BRDF information of unmanned aerial vehicles. Background Art

[0002] Multi-angle remote sensing technology is a fundamental means of acquiring multi-angle reflectance data of ground objects. Fine-scale, multi-angle data plays a crucial role in optimizing canopy BRDF modeling and improving vegetation parameter inversion capabilities. BRDF is defined as the ratio of the reflected radiance of a ground object in a given direction to the incident irradiance, describing the anisotropic characteristics of the ground object's reflectance. This parameter is difficult to measure directly, so the bidirectional reflectance factor (BRF) is often measured in the field. This is defined as the ratio of the reflected radiance of the ground object to the reflected radiance of a standard reference surface under the same conditions. Satellite multi-angle sensors (such as MODIS, MISR, and POLDER) typically have resolutions ranging from hundreds of meters to kilometers, making it impossible to observe canopy BRFs at a resolution of ten meters, let alone directly measure the underlying surface. Aerial multi-angle remote sensing platforms typically offer high resolution and a wide range of angular sampling, making them commonly used for validating canopy BRDF models. However, their high operating costs limit their widespread application. Ground-based multi-angle platforms can provide multi-angle observations of low-lying vegetation, but only over very small areas. The development of drones and miniaturized imaging sensor technology has spurred the emergence of multi-angle remote sensing platforms. The flexible and controllable flight characteristics of multi-rotor drones enable them to provide multi-angle data from a wider range of observation angles and at higher spatial resolution. Ultra-high-resolution drone remote sensing imagery also provides direct observation of the underlying surface. Therefore, multi-angle and multispectral drone remote sensing imagery has the potential to enable flexible and efficient acquisition of ground feature BRF information.

[0003] The technical means of UAV multi-angle remote sensing is to obtain information from multiple observation directions of the same surface target from different positions of the UAV in three-dimensional space. The main features of this new type of multi-angle remote sensing are: (1) The flexible and controllable flight mode of multi-rotor UAVs enables them to provide multi-angle images with high angle sampling rate and centimeter-level spatial resolution; (2) The relatively fast flight speed of UAVs enables them to complete multi-angle observations in a short time to ensure that the illumination geometry does not change too much; (3) The easy deployment of UAVs enables them to achieve long-term series observations in different areas; (4) UAVs can achieve multi-angle observations of tall vegetation such as forests at a sample scale, which is not available in ground-based multi-angle remote sensing. In summary, UAV multi-angle remote sensing can be achieved through a rotating gimbal, overlapping images, or multiple tilt cameras. Currently, there have been many studies on extracting BRDF information of ground objects based on UAV multi-angle observations. However, the relevant technical methods are not very universal and are often limited to special types of observation data. There are also many assumptions. Therefore, this patent proposes a universal UAV multi-angle remote sensing BRDF information acquisition and modeling method. Summary of the Invention

[0004] In order to achieve the above-mentioned object of the present invention, the present invention provides a general method for obtaining and modeling BRDF information of multi-angle remote sensing of unmanned aerial vehicles, comprising the following steps:

[0005] (1) Obtain multi- / hyperspectral frame images from multi-angle remote sensing observations by UAVs;

[0006] (2) Perform radiometric correction on the image to obtain the reflectivity factor of the bottom of the atmosphere;

[0007] (3) Import the image into computer vision software (such as Metashape or Pix4Dmapper) and use the structure-from-motion (SfM) algorithm to perform 3D reconstruction to obtain dense point clouds and DSM;

[0008] (4) Back-projecting the three-dimensional coordinates of the grid corner points set for the DSM into the two-dimensional pixel coordinates of the original image based on the projection transformation (or collinearity condition equation);

[0009] (5) Calculate the observation geometry of the four corner points of each grid and the average BRF of all pixels in the grid, and finally obtain the multi-angle BRF of each grid;

[0010] (6) A semi-empirical linear kernel driven model is used to construct a kernel function equation for each aggregated pixel, and the constructed kernel function equation is used to estimate the BRDF of the unknown angle. The combination of the kernel function is the Ross-Thick-Maignan kernel that characterizes volume scattering and the Li-Sparse-Reciprocal kernel that characterizes geometric optical scattering, that is, the RTMLSR kernel function model.

[0011] In a preferred embodiment of this method, when performing radiometric correction on an image, since both direct light and diffuse light from the sky exist in a field environment, the radiation received by the ground object includes both direct solar radiation and diffuse radiation from the sky. Therefore, strictly speaking, the drone camera sensor obtains the hemispherical directional reflectance factor (HDRF). If the contribution of diffuse light is ignored at this time, it can be approximately considered as BRF. When performing radiometric correction, the sensor noise, lens vignetting effect, camera exposure time and gain must first be corrected to obtain a standardized DN value. (Formula 1); then the standardized DN value is converted into BRF (BRF) using the single-point correction method (using only one standard reflectance Lambertian plate, Formula 2) or the empirical line correction method (using two or more standard reflectance Lambertian plates, Formula 3) target ), in this process, implicit atmospheric correction is also achieved.

[0012]

[0013] Among them, DN raw is the DN value of the original image obtained, DN noise is the global noise of the camera sensor (or pixel-by-pixel noise value), t exposure is the exposure time of the camera when shooting, g iso is the sensitivity, r is the distance between the pixel point and the center pixel of the image (in pixels), k i BRF is the correction coefficient value of the polynomial model for the vignetting effect. ref is the known reflectivity value of the standard reflectivity Lambertian plate, Normalized DN values ​​of the Lambertian plate obtained for UAV imagery.

[0014] In a preferred embodiment of this method, since the multi-angle images have a high degree of overlap, they can be imported into computer vision software (such as Metashape or Pix4Dmapper) and used to perform 3D reconstruction using the structure-from-motion (SfM) algorithm to obtain dense point clouds and DSMs, and then orthorectify the image to obtain orthophotos. Note that since the quality of 3D reconstruction is highly dependent on image overlap, the quality of the generated dense point clouds and DSMs is not high for areas with low image overlap. Therefore, slight cropping is required to retain only the point clouds and DSMs in high-quality 3D reconstruction areas. To facilitate subsequent calculations, the geographic coordinate system of the DSM (i.e., longitude, latitude, and elevation) needs to be converted to a projected coordinate system (with the units of the x, y, and z axes all in meters). Note that the reference ellipsoid (usually WGS1984), projection mode (usually UTM projection), and projection zone partitioning need to be specified. This conversion can be directly implemented in the above software.

[0015] In a preferred embodiment of the method, the three-dimensional coordinates of the grid corner points set for the DSM are back-projected into the two-dimensional pixel coordinates of the original image based on the projection transformation (or the collinearity condition equation) (4).

[0016]

[0017] In a preferred embodiment of this method, the observation geometry needs to be solved pixel by pixel or grid by grid. In order to calculate the observation geometry at all observation positions covered by each grid, the actual shooting point coordinates of all images containing the grid (not the waypoint coordinates, but the camera position coordinates in the external parameter matrix) are extracted, and then the observation zenith angle θ (Equation 5) and azimuth angle of each of the four corner points of each grid are calculated based on the spatial triangulation relationship. (6), then the observation zenith angle of the four corner points is averaged to obtain the observation zenith angle of the grid, the observation azimuth angle of the four corner points is median to obtain the observation azimuth angle of the grid, and the average BRF of all pixels in the grid is calculated to finally obtain the multi-angle BRF of each grid.

[0018]

[0019]

[0020] Where x, y, and z are the easting and northing components and elevation (all in meters) in the UTM projected coordinate system. The subscripts p and c represent the ground point and camera point corresponding to a pixel in the image, respectively.

[0021] In a preferred embodiment of this method, a semi-empirical linear kernel driven BRDF algorithm based on a kernel driven model has strong fitting ability, can ensure inversion accuracy and speed, and has been widely used. The kernel function is constructed using the multi-angle observation information of each pixel, and then the kernel function is used to estimate the BRDF information of each aggregated pixel at other angles. The kernel function of the algorithm used in this method is a combination of the Ross-Thick-Maignan kernel reflecting volume scattering and the Li-Sparse-Reciprocal kernel reflecting geometric optics scattering, namely the RTMLSR kernel function model. The specific expression is as follows:

[0022] Nuclear driven model:

[0023]

[0024] The f iso Isotropic scattering kernel coefficient, f vol is the volume scattering kernel coefficient, F1 is the volume scattering kernel function, f geo is the geometric optics kernel coefficient, and F2 is the geometric optics kernel function.

[0025] The formulas for the F1 and F2 kernel functions are as follows:

[0026]

[0027] cosξ=cosθ i cosθ v +sinθ i sinθ v cosφ (12)

[0028]

[0029] The θ i ,θ v 、 They are the solar zenith angle, observation zenith angle and the relative azimuth angle between the sun and the sensor. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] In order to more clearly illustrate the technical advantages of the present invention, the present invention and its application examples are briefly introduced below.

[0031] Figure 1 It is a schematic diagram of the process of the present invention;

[0032] Figure 2 These are four typical multi-angle observation modes of drones;

[0033] Figure 3 (a) Figure 3 (b) Comparison of the camera optical system before and after correction of the vignetting effect;

[0034] Figure 4 An example of an empirical line radiation correction model;

[0035] Figure 5 (a) Figure 5 (b) Digital elevation models and digital orthophotos reconstructed by computer vision software;

[0036] Figure 6 This is a visualization of the number of observation statistics for each grid when the grid size is 10m;

[0037] Figure 7 This is a schematic diagram for calculating the observation zenith angle and observation azimuth angle of a certain point on the ground;

[0038] Figure 8 A schematic diagram of the observation angles of all observation points on the main plane of a grid;

[0039] Figure 9 (a) Figure 9 (b) BRF of all observation points on the main plane of two grids and the BRF profile fitted by the kernel-driven BRDF model.

[0040] Specific implementation methods

[0041] The following will be combined with the example diagrams in the present invention to explain the technical solution in the present invention in more detail. Figure 1 As shown in Figure 1, a general method for obtaining and modeling BRDF information from multi-angle remote sensing of UAVs is implemented as follows:

[0042] Step 1: Multi-angle observation of drones can usually be achieved by controlling the rotation of the gimbal to achieve tilted observation of the camera, or by shooting images with a certain degree of overlap at different positions to achieve multi-angle observation within a relatively small angle, or by shooting with a wide-format camera with a large field of view to achieve multi-angle observation under certain assumptions. Figure 2 As shown, multi-angle observation modes such as "goniometer," "gyroscopic," "linear," "wide-format," and "tilted camera" can be used to acquire multi- and hyperspectral frame-format remote sensing imagery from drones. It's important to note that multi-angle observations from drones must be conducted under direct sunlight (clear sky or cloudless windows). The recommended observation time is a solar zenith angle of less than 35°, and a single flight mission is recommended to be less than 0.5 hours. During this time, the solar angle and radiation intensity typically vary little, so consistent lighting conditions can be assumed for each capture. Additionally, one or more standard Lambertian plates with different reflectances should be photographed before takeoff and landing.

[0043] Step 2: When performing radiometric correction on an image, since both direct sunlight and diffuse sky light exist in the field, the radiation received by the ground object includes both direct sunlight and diffuse sky radiation. Therefore, strictly speaking, the drone camera sensor obtains the hemispherical directional reflectance factor (HDRF). If the contribution of diffuse light is ignored at this time, it can be approximately considered as BRF. When performing radiometric correction, the sensor noise, lens vignetting effect, camera exposure time and gain must first be corrected to obtain a standardized DN value. (Formula 1); then the standardized DN value is converted into BRF (BRF) using the single-point correction method (using only one standard reflectance Lambertian plate, Formula 2) or the empirical line correction method (using two or more standard reflectance Lambertian plates, Formula 3) target ), in this process, implicit atmospheric correction is also achieved.

[0044]

[0045] Among them, DN raw is the DN value of the original image obtained, DN noise is the global noise of the camera sensor (or pixel-by-pixel noise value), t exposure is the exposure time of the camera when shooting, g iso is the sensitivity, r is the distance between the pixel point and the center pixel of the image (in pixels), k i BRF is the correction coefficient value of the polynomial model for the vignetting effect. ref is the known reflectivity value of the standard reflectivity Lambertian plate, Normalized DN values ​​of the Lambertian plate obtained for UAV imagery.

[0046] Step 3: Since the images observed from multiple angles have a high degree of overlap, they can be imported into computer vision software (such as Metashape or Pix4Dmapper) and used to perform 3D reconstruction using the structure-from-motion (SfM) algorithm to obtain dense point clouds and DSM ( Figure 5), and then obtain the orthophoto through orthorectification. Note that since the quality of 3D reconstruction is highly dependent on the image overlap, the quality of the generated dense point cloud and DSM is not high for those areas with low image overlap, so they need to be slightly cropped to retain only the point cloud and DSM in the high-quality 3D reconstruction area. To facilitate subsequent calculations, the geographic coordinate system of the DSM (i.e., longitude, latitude and elevation) needs to be converted to a projected coordinate system (the units of the xyz axes are all m). Note that it is necessary to specify the reference ellipsoid (generally WGS 1984), projection mode (generally UTM projection) and projection zone partition. This conversion can be directly implemented in the above software.

[0047] Step 4: Set up a grid for the DSM according to the specified size (the recommended grid size is 3m or more), or set up a grid according to the specified area of ​​interest, and calculate the three-dimensional coordinates of each corner point in the grid (in meters). Based on the projection transformation (or collinearity condition equation), the three-dimensional coordinates of the grid corner points are back-projected into the two-dimensional pixel coordinates of the original image. The intrinsic parameter matrix, extrinsic parameter matrix, and distortion coefficients required for the projection transformation can all be estimated by the computer vision software when executing the SfM algorithm. Note that since the geographical range of the DSM is usually larger than the geographical range covered by the original image, sometimes the image pixel coordinates after back-projection of some grid points exceed the range of the original image and therefore need to be discarded. Finally, the serial number of each grid, the three-dimensional coordinates of the corner points, the name of the image containing the grid, and the image pixel coordinates after back-projection of the grid in the image are saved ( Figure 6 ).

[0048]

[0049] Step 5: In order to calculate the observation geometry of all observation positions covered by each grid, the actual shooting point coordinates of all images containing the grid (not the waypoint coordinates, but the camera position coordinates in the external parameter matrix) are extracted, and then the spatial triangulation geometry relationship ( Figure 7 ), calculate the observation zenith angle θ (Formula 5) and azimuth angle of each grid's four corner points (Equation 6) The observation zenith angles of the four corner points are then averaged to obtain the grid's observation zenith angle, the observation azimuth angles of the four corner points are median to obtain the grid's observation azimuth angle, and the average BRF of all pixels within the grid is calculated to ultimately obtain the multi-angle BRF for each grid. Note that although the camera's principal optical axis observation angles for each waypoint have been predefined according to user design requirements during the multi-angle route planning phase of the drone, the actual camera's observation zenith angle and azimuth angle often deviate from the preset observation angles due to factors such as navigation positioning accuracy, wind force, self-vibration, and angle control accuracy during flight. Therefore, post-processing of the camera's observation geometry is necessary. Furthermore, the imaging sensor's field of view (FOV) effect causes the actual observation geometry of each pixel to be inconsistent. To meet the high-precision observation geometry requirements of certain applications, it is necessary to solve the observation geometry pixel by pixel or grid by grid.

[0050]

[0051]

[0052] Where x, y, and z are the easting and northing components and elevation (all in meters) in the UTM projected coordinate system. The subscripts p and c represent the ground point and camera point corresponding to a pixel in the image, respectively.

[0053] Step 6: Because multi-angle drone observations only capture BRFs for a limited number of directions, we use a kernel-driven BRDF model to perform a band-by-band fit to obtain BRFs for all directions within the hemisphere. During fitting, the number of observations and their angles may be concentrated within a very narrow angular range, so the number of observations and their angles must be carefully selected to ensure a relatively homogeneous distribution of observation angles. It should be noted that due to the camera's field of view, the BRF for a particular observation direction in each grid cell is not the value of a single observation angle, but rather the average of pixels from multiple different observation angles. This results in some offsetting, causing some discrepancies between the observed and actual BRF values. In particular, near-hotspot directions, the BRF values ​​are somewhat underestimated. After fitting, the accuracy of the fit is evaluated. If the fit accuracy is poor, the observation data to be fitted is readjusted until the fit accuracy meets a certain threshold.

[0054] The kernel-driven BRDF model used in this method is RTMLSR, which uses the Li-Sparse-Reciprocal kernel for geometric optics and the Ross-Thick-Maignan hotspot kernel for volume scattering. Based on selected multi-angle BRF observation data, the Levenberg-Marquardt algorithm is used to optimize the three coefficients of the RTLSR model. This allows for the prediction of the BRF in all directions of the hemisphere.

[0055] Kernel-driven model:

[0056]

[0057] where f iso is the isotropic kernel coefficient, f vol is the volume scattering kernel coefficient, F1 is the volume scattering kernel function, f geo is the geometric optics kernel coefficient, and F2 is the geometric optics kernel function.

[0058]

[0059] cosξ=cosθ i cosθ v +sinθ i sinθ v cosφ (12)

[0060]

[0061] The θ i ,θ v 、 They are the solar zenith angle, observation zenith angle and the relative azimuth angle between the sun and the sensor.

Claims

1. A general method for obtaining and modeling BRDF information from multi-angle remote sensing of unmanned aerial vehicles, comprising the following steps: (1) Obtain multi- / hyperspectral frame images from multi-angle remote sensing observations by UAVs; (2) Perform radiometric correction on the image to obtain the reflectivity factor of the bottom of the atmosphere; (3) Import the image into computer vision software and perform 3D reconstruction using structure-from-motion algorithm to obtain dense point cloud and DSM; (4) back-projecting the three-dimensional coordinates of the grid corner points set for the DSM into the two-dimensional pixel coordinates of the original image based on the projection transformation or the collinearity condition equation; (5) Calculate the observation geometry of the four corner points of each grid and the average BRF of all pixels in the grid, and finally obtain the multi-angle BRF of each grid; (6) A semi-empirical linear kernel driven model is used to construct a kernel function equation for each aggregated pixel, and the constructed kernel function equation is used to estimate the BRDF of the unknown angle. The combination of the kernel function is the Ross-Thick-Maignan kernel that characterizes volume scattering and the Li-Sparse-Reciprocal kernel that characterizes geometric optics scattering, namely the RTMLSR kernel function model; When performing radiometric correction on images, since both direct light and diffuse light from the sky exist in the wild environment, the radiation received by the ground object includes both direct solar radiation and diffuse radiation from the sky. Therefore, strictly speaking, the drone camera sensor obtains the hemispherical directional reflectance factor, HDRF. If the contribution of diffuse light is ignored at this time, it can be approximately considered as BRF. When performing radiometric correction, the sensor noise, lens vignetting effect, camera exposure time and gain must first be corrected to obtain a standardized DN value. See formula 1; then use the single-point correction method, that is, only one standard reflectivity Lambertian plate, see formula 2, or the empirical line correction method, that is, two or more standard reflectivity Lambertian plates, see formula 3, to convert the standardized DN value into BRF target ,In this process, implicit atmospheric correction is also achieved; Among them, DN raw is the DN value of the original image obtained, DN noise is the global noise or pixel-by-pixel noise value of the camera sensor, t exposure is the exposure time of the camera when shooting, g iso is the sensitivity, r is the distance between the pixel point and the center pixel of the image, the unit is pixel, k i is the correction coefficient value of the polynomial model for the vignetting effect, BRF ref is the known reflectivity value of the standard reflectivity Lambertian plate, Normalized DN values ​​of the Lambertian plate obtained for UAV imagery.

2. A general UAV multi-angle remote sensing BRDF information acquisition and modeling method according to claim 1, characterized in that: Since the images observed from multiple angles have a high degree of overlap, they can be imported into computer vision software and the structure-from-motion algorithm can be used for 3D reconstruction to obtain dense point clouds and DSM, and then orthophotos can be obtained through orthorectification. Note that the quality of 3D reconstruction is highly dependent on image overlap, so for those areas with low image overlap, the quality of the generated dense point clouds and DSM is not high, so they need to be slightly cropped to retain only the point clouds and DSMs in the high-quality 3D reconstruction areas. To facilitate subsequent calculations, the geographic coordinate system of the DSM, i.e., longitude, latitude and elevation, needs to be converted to a projected coordinate system. The units of the xyz axes are all in m. Note that the reference ellipsoid, projection mode, and projection zone partition need to be specified. This conversion can be implemented directly in the above software.

3. A general UAV multi-angle remote sensing BRDF information acquisition and modeling method according to claim 1, characterized in that: Based on the projection transformation or collinearity condition equation, the three-dimensional coordinates of the grid corner points set for the DSM are back-projected into the two-dimensional pixel coordinates of the original image, as shown in Equation 4.

4. A general UAV multi-angle remote sensing BRDF information acquisition and modeling method according to claim 1, characterized in that: The observation geometry needs to be solved pixel by pixel or grid by grid. In order to calculate the observation geometry at all observation positions covered by each grid, the actual shooting point coordinates of all images containing the grid are extracted. These are not the waypoint coordinates, but the camera position coordinates in the external parameter matrix. Then, according to the spatial triangulation relationship, the observation zenith angle θ of each of the four corner points of each grid is calculated, as shown in Equation 5, and the azimuth angle See Equation 6. Then, the observation zenith angle of the four corner points is averaged to obtain the observation zenith angle of the grid, the observation azimuth angle of the four corner points is median to obtain the observation azimuth angle of the grid, and the average BRF of all pixels in the grid is calculated to finally obtain the multi-angle BRF of each grid. Where x, y, and z are the Easting component, Northing component, and elevation in the UTM coordinate system, respectively. The units are all in meters. The subscripts p and c represent the ground point and camera point corresponding to a pixel in the image, respectively.

5. A general UAV multi-angle remote sensing BRDF information acquisition and modeling method according to claim 1, characterized in that: Based on the kernel-driven model method, the semi-empirical linear kernel-driven BRDF algorithm has strong fitting ability, can ensure inversion accuracy and speed, and has been widely used. It uses the multi-angle observation information of each pixel to construct a kernel function, and then uses this function to estimate the BRDF information of other angles of each aggregated pixel. The kernel function combination of the algorithm used in this method is the Ross-Thick-Maignan kernel reflecting volume scattering and the Li-Sparse-Reciprocal kernel reflecting geometric optical scattering, that is, the RTMLSR kernel function model. The specific expression is as follows: Kernel-driven model: The f iso Isotropic scattering kernel coefficient, f vol is the volume scattering kernel coefficient, F1 is the volume scattering kernel function, f geo is the geometric optics kernel coefficient, F2 is the geometric optics kernel function; The formulas for the F1 and F2 kernel functions are as follows: cosξ=cosθ i cosθ v +sinθ i sinθ v cosφ (12) The θ i ,θ v 、 They are the solar zenith angle, observation zenith angle and the relative azimuth angle between the sun and the sensor.

Images