A method and device for measuring carbon dioxide flux in forest areas based on drones

Abstract

The present invention provides a method and device for measuring carbon dioxide flux in forest areas based on drones. The device comprises: (1) a wind speed correction module; (2) a profile data processing module; (3) a horizontal data processing module; and (4) a carbon dioxide flux calculation module. The present invention constructs a three-dimensional carbon dioxide flux algorithm model and system for forest ecosystems based on drones, and uses the drone system to measure horizontal and vertical carbon dioxide fluxes above the forest canopy. This overcomes the problem that flux towers are limited by their height and cannot measure high-level fluxes or cannot be moved for measurement. For complex terrain such as forest ecosystems, the device takes into account the horizontal flux variation in non-horizontal, uniform, and flat terrain areas, making the calculated carbon dioxide flux results closer to the actual situation.

Description

Technical Field

[0001] The present invention relates to the field of environmental measurement technology, and in particular to a method and device for measuring carbon dioxide flux in forest areas based on drones. Background Art

[0002] As the impact of global climate change becomes increasingly serious, accurately predicting and controlling greenhouse gas emissions has become a top priority. Among them, carbon dioxide is one of the main greenhouse gases, and its source strength and sink weakness are of great significance for understanding and mitigating climate change. However, traditional carbon dioxide monitoring methods, such as ground-based observation stations and flux towers, are often restricted by geographical location and have difficulty covering large areas of forest ecosystems. Flux towers are limited by their height and cannot measure high-level fluxes and cannot be moved for measurement. For complex terrain such as forest ecosystems, there is no integrated algorithm research and equipment system for complex and comprehensive terrain such as non-horizontal, uniform, and flat terrain areas. Therefore, how to efficiently and accurately measure carbon dioxide flux in forest areas has become a hot topic of current research. Summary of the Invention

[0003] The object of the present invention is to provide a method for measuring forest carbon dioxide flux based on an unmanned aerial vehicle, the method comprising:

[0004] (1) Wind speed correction module: corrects the wind speed changes caused by the UAV propeller rotation, movement and attitude (pitch, roll and yaw);

[0005] (2) Profile data processing module: Obtain the vertical profile information of the UAV and perform Gaussian weighted average fitting on the wind profile, carbon dioxide concentration, and water vapor concentration; calculate the stability correction function based on the above wind profile using the empirical constant judgment formula; calculate the Obukhov length and turbulence characteristic quantity based on the above Richardson number; calculate the turbulence diffusion coefficient based on the above turbulence characteristic quantity;

[0006] (3) Horizontal data processing module: obtains horizontal data information of the UAV at different altitudes, and performs Gaussian weighted average fitting on the horizontal wind speed, carbon dioxide concentration, and water vapor concentration; calculates the normal wind speed in the longitude and latitude directions according to the flight trajectory based on the above horizontal wind speed; performs the carbon dioxide flux line integral based on the above longitude and latitude directions, that is, multiplying the normal wind speed by the carbon dioxide concentration gradient to obtain the inflow and outflow carbon dioxide flux in the longitude and latitude directions respectively;

[0007] (4) CO2 flux calculation module: For profile data, based on the flux gradient method, the turbulent diffusion coefficient at different heights is multiplied by the CO2 concentration gradient to obtain the vertical CO2 flux; for horizontal data, based on the mass balance method, the CO2 flux outflow in the longitude and latitude directions is subtracted from the CO2 flux inflow to obtain the horizontal CO2 flux. Furthermore, if the temperature of the two layers is known, and Two-layer wind speed and The overall Richardson number can be calculated, and the expression for the overall Richardson number is as follows: Where g is the acceleration due to gravity; The average temperature of the two layers; z is the height of the lower layer wind speed. Furthermore, different wind speed profile and temperature profile stability correction functions are selected according to different overall Richardson numbers. and

[0008] The stability correction function can generally be written as an empirical function: Among them, the determination of the empirical constants α and a, b is based on the different Richardson numbers. Further, if the average wind speed at different heights z1, z2 and z3 is known and According to the overall Richardson number R iB The relationship with the stability parameter z / L is: in,

[0009] are the integral forms of the stability correction functions for the wind speed profile and temperature profile, respectively, ξ = z / L, ξ0 = z0 / L; z0 is the surface roughness, and L is the Monin-Obukhov length. For steep, non-streamlined rough surfaces, the fitting curve expression is based on the zero plane displacement (d0) and the positive area coefficient (λ): d0 = hβλ γ Where h is the canopy height; β and γ are general coefficients and are positive area coefficients. In addition, based on the relationship between zero plane displacement and surface roughness (z0), an empirical formula is given: d0 = Cz 0。 C is a general coefficient. Substituting the coefficient C into the above relationship, the surface roughness z0 is calculated, and then according to the overall Richardson number R iB The relationship between the Monin-Obukhov length and the stability parameter z / L is used to calculate the Monin-Obukhov length. Furthermore, based on the functional relationship between the average meteorological element profile near the surface layer and the turbulent flux, the surface roughness z0 is uniformly taken as the lower boundary condition to obtain the integral expressions of the wind speed profile and temperature profile:

[0010] Where κ is the vonKarman constant; u * is the friction speed; θ * is the potential temperature characteristic scale; is the average potential temperature near the ground. According to the integral expression of the wind speed profile and temperature profile, the turbulence characteristic scale θ can be calculated * Furthermore, the turbulent viscosity rate (K m ) and turbulent thermal conductivity (K h) is a physical quantity that characterizes the strength of turbulence development in the near-surface layer. It is expressed by the following flux profile relationship: Furthermore, the calculation formula of the vertical carbon dioxide flux is as follows: Where ρ is the air density; K is the turbulent diffusion coefficient; and c is the carbon dioxide mixing ratio. Furthermore, horizontal data from the drone at different altitudes was obtained, and a Gaussian weighted average fit was performed on the horizontal wind speed, carbon dioxide concentration, and water vapor concentration. The normal wind speed in the longitude and latitude directions was calculated based on the flight trajectory. The carbon dioxide flux line integral was performed along these longitude and latitude directions. The normal wind speed was multiplied by the carbon dioxide concentration gradient to obtain the inflow and outflow carbon dioxide flux in the longitude and latitude directions, respectively. Using the mass balance method, the outflow carbon dioxide flux in the longitude and latitude directions was subtracted from the inflow carbon dioxide flux to obtain the horizontal carbon dioxide flux. The calculation formula is as follows: ΔF NIt is the change in carbon dioxide flux corresponding to each data point in the longitude or latitude direction. In addition, the present invention also provides a device for measuring forest carbon dioxide flux based on a drone, which includes: a drone system, a three-dimensional ultrasonic wind speed system, a data transmission system, a gas sampling system and an RTK positioning system. Among them, the gas sampling system is mounted on the bottom of the drone, and the air inlet of the sampling tube is fixed at the same height as the sensor of the three-dimensional ultrasonic wind speed system to ensure that the data collected by the gas sampling system and the three-dimensional ultrasonic wind speed system have the same spatiality and avoid the influence caused by the downwash airflow of the drone; the ultrasonic anemometer is fixed above the propeller of the drone; the data transmission system is located under the rectangular base under the drone; the RTK system is to install two GPS systems behind the nose of the drone, and add a high-precision GNSS receiver on the ground to ensure continuous tracking and observation of satellites and the quality of satellite signals. Furthermore, the drone system is a multi-purpose medium-sized four-rotor or six-rotor drone (such as Shenzhen DJI M300, Tianjin Feiyan F-10L, F-30L, etc.). The drone system is equipped with a three-dimensional ultrasonic analyzer system, a data transmission system, a gas sampling system, etc. for flight, and ensures that the single flight time is more than half an hour. Furthermore, the three-dimensional ultrasonic wind speed system is small in size and low in power consumption. It is mounted above the drone and measures meteorological parameters such as wind speed, wind direction, and temperature during the flight of the drone. Furthermore, the gas sampling system samples the surrounding environmental gases during the flight of the drone and outputs the carbon dioxide and water vapor concentration values. Since it is mounted on the drone, it needs to be small in size and light in weight. Furthermore, the RTK positioning system, the real-time dynamic measurement system (RTK, Real Time Kinematic) (RTK, Haixingda, CN) adopts a multi-star multi-frequency GNSS unit, supporting BDS, GPS, GLONASS and other systems for navigation and positioning; the UBase receiver is installed on a tripod, fixed on relatively flat ground, to receive satellite signals, and form a three-point positioning with the two GPS positioning systems on the drone, to achieve high-precision measurement of the forest ecosystem. Since the canopy height data needs to be measured for forest measurement, the RTK positioning system can accurately measure the flight altitude of the drone to prevent the drone from colliding with the trees when measuring the forest canopy. Furthermore, the data transmission system mainly uses two data radios for data transmission, namely the drone's data transmission radio and the second data transmission radio, which transmit the drone's flight parameters and the measurement data of the meteorological sensor to the ground computer at the same time and save them locally for further data processing and analysis.The drone's data radio communicates with the ground receiver via a serial port, allowing the ground receiver to receive decoded GPS information and flight attitude data from the drone. This data can then be sent to the drone via the ground station software, enabling mission configuration and route planning. The second data radio, through communication between the drone's transmitter and the ground receiver, can also transmit various meteorological sensor measurement data to a local computer via a serial port.

[0011] This invention has the following beneficial effects: It establishes a three-dimensional carbon dioxide flux algorithm model for forest ecosystems based on drones. Using the drone system, it measures horizontal and vertical carbon dioxide fluxes in the forest canopy, thus overcoming the limitations of traditional flux towers, which are unable to observe high-level fluxes and are difficult to measure using mobile devices. For forest ecosystems in complex, non-horizontal, uniform, and flat terrain, the model also accounts for horizontal flux variations, making the calculated carbon dioxide flux results more consistent with actual environmental conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 Measurement calculation flow chart.

[0013] Specific implementation: The following is a further detailed description of the implementation of the present invention in conjunction with the accompanying drawings and examples. The following examples are used to illustrate the present invention, but are not intended to limit the scope of the present invention.

[0014] Example 1, as Figure 1 As shown,

[0015] A method for measuring forest carbon dioxide flux based on drones, the method comprising:

[0016] (S1) Wind speed correction module: corrects the wind speed changes caused by the UAV propeller rotation, movement and attitude (pitch, roll and yaw);

[0017] (S2) Profile data processing module: obtains the vertical profile information of the UAV, performs Gaussian weighted average fitting on the wind profile, carbon dioxide concentration, and water vapor concentration; calculates the stability correction function based on the above wind profile using the empirical constant judgment formula; calculates the Obukhov length and turbulence characteristic quantity based on the above Richardson number; calculates the turbulence diffusion coefficient based on the above turbulence characteristic quantity;

[0018] (S3) Horizontal data processing module: obtains horizontal data information of the UAV at different altitudes, performs Gaussian weighted average fitting on horizontal wind speed, carbon dioxide concentration, and water vapor concentration; calculates the normal wind speed in the longitude and latitude directions according to the flight trajectory based on the horizontal wind speed; performs a linear integral of the carbon dioxide flux based on the longitude and latitude directions, that is, multiplies the normal wind speed by the carbon dioxide concentration gradient to obtain the inflow and outflow carbon dioxide flux in the longitude and latitude directions, respectively;

[0019] (S4) CO2 flux calculation module: For profile data, based on the flux gradient method, the turbulent diffusion coefficient at different heights is multiplied by the CO2 concentration gradient to obtain the vertical CO2 flux; for horizontal data, based on the mass balance method, the CO2 flux outflow in the longitudinal and latitudinal directions is subtracted from the CO2 flux inflow to obtain the horizontal CO2 flux. and Two-layer wind speed and The overall Richardson number can be calculated, and the expression for the overall Richardson number is as follows: Where g is the acceleration due to gravity; The average temperature of the two layers; z is the height of the lower layer wind speed. Different wind speed profiles and temperature profile stability correction functions are selected according to different overall Richardson numbers. and The stability correction function can generally be written as an empirical function: The empirical constants α, a, and b are determined based on the different Richardson numbers. If the average wind speed at different heights z1, z2, and z3 is known, and According to the overall Richardson number R iB The relationship with the stability parameter z / L is:

[0020] in,

[0021] are the integral forms of the stability correction functions for the wind speed profile and temperature profile, respectively, ξ = z / L, ξ0 = z0 / L; z0 is the surface roughness, and L is the Monin-Obukhov length. For steep, non-streamlined rough surfaces, the fitting curve expression is based on the zero plane displacement (d0) and the positive area coefficient (λ): d0 = hβλ γ Where h is the canopy height; β and γ are general coefficients, which are positive area coefficients. Furthermore, the empirical formula d0 = Cz0 is given based on the relationship between zero plane displacement and surface roughness (z0). C is a general coefficient. Substituting coefficient C into the above relationship, the surface roughness z0 is calculated, and then the surface roughness z0 is calculated based on the overall Richardson number R. iBThe relationship between the Monin-Obukhov length and the stability parameter z / L is calculated. According to the functional relationship between the average meteorological element profile of the near-surface layer and the turbulent flux, the surface roughness z0 is uniformly taken as the lower boundary condition, and the integral expressions of the wind speed profile and temperature profile are obtained:

[0022] in,

[0023] κ is the von Karman constant; u * is the friction speed; θ * is the potential temperature characteristic scale; is the average potential temperature near the ground. According to the integral expression of the wind speed profile and temperature profile, the turbulence characteristic scale θ can be calculated * . Turbulent viscosity rate (K m ) and turbulent thermal conductivity (K h ) is a physical quantity that characterizes the strength of turbulence development in the near-surface layer. It is expressed by the following flux profile relationship: Furthermore, the calculation formula of the vertical carbon dioxide flux is as follows: Where ρ is the air density; K is the turbulent diffusion coefficient; and c is the carbon dioxide mixing ratio. Horizontal data from the drone at different altitudes was obtained, and a Gaussian weighted average fit was performed on the horizontal wind speed, carbon dioxide concentration, and water vapor concentration. The normal wind speed in the longitude and latitude directions was calculated based on the flight trajectory. The carbon dioxide flux was linearly integrated along these longitude and latitude directions. The normal wind speed was multiplied by the carbon dioxide concentration gradient to obtain the inflow and outflow carbon dioxide fluxes in the longitude and latitude directions, respectively. Using the mass balance method, the outflow carbon dioxide flux in the longitude and latitude directions was subtracted from the inflow carbon dioxide flux to obtain the horizontal carbon dioxide flux. The calculation formula is as follows:

[0024] Example 2

[0025] In addition, the technical solution of the present invention also includes a drone carbon dioxide flux measurement device, which mainly includes the following components: a drone system, a three-dimensional ultrasonic wind speed monitoring system, a data transmission system, a gas sampling system, and an RTK positioning system. The gas sampling system is designed as a detachable module on the bottom of the drone. The opening position of the measured gas inlet is set at the same height as the three-dimensional ultrasonic wind speed system sensor to ensure that the data obtained by the gas sampling system has the same spatial properties and avoid the impact of the drone's descent. The ultrasonic anemometer is installed on the top of the drone's propeller. The data transmission system is placed under the rectangular base of the drone's base. The RTK system consists of two GPS systems, which are installed behind the drone's head. At the same time, a high-precision GNSS receiver is added on the ground to achieve continuous satellite tracking and improve satellite signal quality. For the drone system, a multi-functional quad-rotor or hexacopter drone (such as the DJI M300 in Shenzhen, Tianjin Feiyan F-10L, F-30L, etc.) was selected. The drone system is equipped with a three-dimensional ultrasonic wind speed system, a data transmission system, a gas sampling system, and a flight time of more than half an hour each time. The 3D ultrasonic wind velocity system is a compact and energy-efficient device installed on a drone. It measures meteorological information such as wind speed, wind direction, and temperature in real time during flight. During flight, the gas sampling system collects ambient air samples and outputs carbon dioxide and water vapor concentration values. This design reduces the weight and size of the device, making it easy to carry on the drone. The RTK positioning system utilizes a real-time kinematic measurement system (RTK, such as the Haixingda, a Chinese product). It utilizes a multi-constellation, multi-frequency GNSS unit and supports navigation and positioning with multiple systems, including BDS, GPS, and GLONASS. This unit is installed behind the drone's nose, and a high-precision GNSS receiver is added on the ground to ensure continuous satellite tracking and improve satellite signal quality. Finally, the data transmission system primarily transmits data via two data radios, which can be used for the drone's own data transmission or as a secondary data radio. The data transmission system's primary function is to simultaneously transmit the drone's flight parameters and meteorological sensor measurements to ground-based computers for processing and analysis. Specifically, the drone's data transmission radio communicates with a ground receiver, allowing the receiver to receive decoded GPS information and flight attitude data from the drone. This control data is then transmitted to the drone via ground station software, enabling mission configuration and route planning. A second data radio also transmits all measurement data from the meteorological sensors to a ground computer via a serial port. A drone-based three-dimensional carbon dioxide flux algorithm model and hardware system for forest ecosystems have been invented and designed.With the help of the drone system, the model can accurately measure the lateral and longitudinal carbon dioxide fluxes in the forest canopy. This can solve the problem that the traditional flux tower cannot observe the top flux and is difficult to measure migration due to height limitations. In addition, for those non-horizontal, uniform and flat terrain areas, the model also takes into account the influence of lateral flux changes, so that the calculated carbon dioxide flux results can be closer to the real environment. The embodiments of the present invention are given for the purpose of illustration and description, and are not exhaustive or limit the invention to the disclosed form. Many modifications and changes are obvious to those of ordinary skill in the art. The embodiments are selected and described in order to better illustrate the principles and practical applications of the present invention, and to enable those of ordinary skill in the art to understand the invention and design various embodiments with various modifications suitable for specific purposes.

Claims

1. A method for calculating forest carbon dioxide flux based on drone measurements, characterized in that: The method comprises: (1) Wind speed correction: Correct the wind speed changes caused by the rotation, movement and attitude of the UAV propeller; (2) Profile data processing: Obtain the vertical profile information of the UAV and perform Gaussian weighted average fitting on the wind profile, carbon dioxide concentration, and water vapor concentration; (3) Horizontal data processing: Obtain horizontal data information of the drone at different altitudes and perform Gaussian weighted average fitting on horizontal wind speed, carbon dioxide concentration, and water vapor concentration; (4) Calculation of carbon dioxide flux: For the profile data, based on the flux gradient method, the turbulent diffusion coefficient at different heights is multiplied by the carbon dioxide concentration gradient to obtain the carbon dioxide vertical flux; the calculation formula of the carbon dioxide vertical flux is as follows: Where ρ is the air density; K is the turbulent diffusion coefficient; c is the carbon dioxide mixing ratio; for horizontal data, based on the mass balance method, the carbon dioxide flux is obtained by subtracting the carbon dioxide flux inflow from the carbon dioxide flux outflow in the longitude and latitude directions respectively. The calculation formula is as follows: ΔF N is the change in carbon dioxide flux corresponding to each data point in the longitude or latitude direction.

2. The method according to claim 1, characterized in that Known temperature and and wind speed and The overall Richardson number can be calculated, and the expression for the overall Richardson number is as follows: Where g is the acceleration due to gravity; is the average temperature of the two layers; z is the height where the wind speed of the lower layer is located.

3. The method according to claim 2, characterized in that Select different stability correction functions for wind speed profiles and temperature profiles according to different overall Richardson numbers and The stability correction function can generally be written as an empirical function: Among them, the empirical constants α and a, b are determined based on the different Richardson numbers and are selected accordingly, ξ=z / L, ξ0=z0 / L, z0 is the surface roughness, and L is the Monin-Obukhov length.

4. The method according to claim 3, characterized in that The average wind speed at different heights z1, z2 and z3 is known and According to the overall Richardson number R iB The relationship with the stability parameter z / L is: in, are the integral forms of the stability correction functions of the wind speed profile and temperature profile, respectively.

5. The method according to claim 4, characterized in that: Turbulent viscosity K m and turbulent thermal conductivity K h It is a physical quantity that characterizes the strength of turbulence development in the near-surface layer. According to the following flux profile relationship expression, it can be obtained: ； u * is the friction speed; θ * is the characteristic scale of turbulence.

6. The method according to claim 5, characterized in that According to the functional relationship between the average meteorological element profile near the surface layer and the turbulent flux, the surface roughness z0 is uniformly taken as the lower boundary condition, and the integral expressions of the wind speed profile and temperature profile are obtained: ； is the average potential temperature near the ground, d0 is the zero plane displacement, and is calculated by the formula d0 = Cz0, where C is the general coefficient. According to the integral expression of the wind speed profile and temperature profile, the turbulence characteristic scale θ can be calculated. * .

Images