CN116500649B - Inversion method and device for boundary layer height - Google Patents

Inversion method and device for boundary layer height Download PDF

Info

Publication number
CN116500649B
CN116500649B CN202310778764.2A CN202310778764A CN116500649B CN 116500649 B CN116500649 B CN 116500649B CN 202310778764 A CN202310778764 A CN 202310778764A CN 116500649 B CN116500649 B CN 116500649B
Authority
CN
China
Prior art keywords
path unit
distance
boundary layer
obstacle
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310778764.2A
Other languages
Chinese (zh)
Other versions
CN116500649A (en
Inventor
苏练
夏海云
袁金龙
王璐
李蔓一
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Science and Technology of China USTC
Original Assignee
University of Science and Technology of China USTC
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Science and Technology of China USTC filed Critical University of Science and Technology of China USTC
Priority to CN202310778764.2A priority Critical patent/CN116500649B/en
Publication of CN116500649A publication Critical patent/CN116500649A/en
Application granted granted Critical
Publication of CN116500649B publication Critical patent/CN116500649B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/95—Lidar systems specially adapted for specific applications for meteorological use
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/148—Wavelet transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Software Systems (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Environmental & Geological Engineering (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Algebra (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Electromagnetism (AREA)
  • Ecology (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Environmental Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Atmospheric Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)
  • Processing Or Creating Images (AREA)

Abstract

The invention provides a boundary layer height inversion method and device, which are characterized in that coordinates of a starting path unit are determined in a search map, convolution operation is carried out by utilizing a Gaussian filter convolution function and the coordinates of the starting path unit, a center path unit and a target path unit are determined, a non-obstacle path unit with the maximum weight distribution value in the target path unit is used as a reference path unit, the moment corresponding to the reference path unit is determined as inversion moment, and the product of a distance library number corresponding to the reference path unit and a distance resolution is used as the boundary layer height at the inversion moment. In the scheme, a non-obstacle path unit corresponding to the height change range of the boundary layer is obtained based on convolution operation simulation, a distance library number corresponding to the non-obstacle path unit with the largest weight distribution value is determined, and the product of the distance library number and the distance resolution is used as the height of the boundary layer, so that the purpose of avoiding the problems of uncertainty in wavelet scale selection and inaccurate results caused by severe jump of the height of the inversion boundary layer is achieved.

Description

Inversion method and device for boundary layer height
Technical Field
The invention relates to the technical field of atmosphere detection, in particular to a boundary layer height inversion method and device.
Background
The atmospheric boundary layer is positioned at the bottom of the troposphere, is a gas layer directly influenced by the earth surface and directly acts on the earth surface, is an important meteorological variable in the vertical direction, influences the longitudinal diffusion degree of pollutants, is widely applied to various meteorological models, and has important auxiliary effects on monitoring the air pollution condition of various places and formulating pollution control policies and the like due to the fact that the atmospheric boundary layer is suitable.
The traditional wavelet transformation method is to confirm the boundary layer height by the abrupt change position of the gradient profile of the air backscattering signal of the laser radar in space.
Because the uncertainty of wavelet scale selection and the severe jump of the inversion boundary layer height lead to inaccurate results, how to avoid the problem of inaccurate results caused by the uncertainty of wavelet scale selection and the severe jump of the inversion boundary layer height is a problem which needs to be solved at present.
Disclosure of Invention
In view of the above, the embodiment of the invention provides a boundary layer height inversion method and device, so as to achieve the purpose of avoiding the problems of uncertainty in wavelet scale selection and inaccurate results caused by severe jump of the inversion boundary layer height.
In order to achieve the above object, the embodiment of the present invention provides the following technical solutions:
the first aspect of the embodiment of the invention discloses a boundary layer height inversion method, which comprises the following steps:
determining a target time for starting boundary layer height inversion from a pre-established search map, wherein the search map is pre-constructed based on a history distance correction signal, the horizontal axis of the search map represents n times of preset intervals, the vertical axis represents k distance library numbers, the search map comprises n multiplied by k path units, each path unit corresponds to a different distance library, each distance library is divided by the history distance correction signal according to the same distance resolution, and each path unit stores the noise ratio of the corresponding distance library, the vertical speed value of the corresponding distance library and the history distance correction signal, and n and k are positive integers;
determining a start path unit and coordinates of the start path unit in a search map based on the target time, the historical distance correction signal and a wavelet covariance transformation method;
performing convolution operation by using a Gaussian filter convolution function and coordinates of the initial path unit, and determining a central path unit and a plurality of target path units from each path unit based on the convolution operation result;
Setting a target path unit with the noise ratio of the corresponding distance library larger than a threshold value as a non-obstacle path unit;
setting weights for the non-obstacle path units based on the distances from the non-obstacle path units to the central path unit, calculating weight distribution values corresponding to the non-obstacle path units based on a distribution value formula and corresponding weights, and taking the non-obstacle path unit with the largest weight distribution value as a reference path unit;
and determining the moment corresponding to the reference path unit as inversion moment, determining a distance library number corresponding to the reference path unit, and taking the product of the distance library number and the distance resolution as the boundary layer height of the inversion moment.
Preferably, the determining the coordinates of the start path unit and the start path unit in the search map based on the target time, the historical distance correction signal and the wavelet covariance transformation method includes:
performing wavelet covariance transformation based on the target time and the historical distance correction signal to obtain a covariance function; wherein the covariance function is:，Rthe height is represented by the height of the plate,bthe distance library corresponding to the smoothed center position of the wavelet transform is numbered, aAs the scale factor of the wavelet transform,R t representing the top height of the historical distance correction signal,R b representing the bottom height of the historical distance correction signal, X%R) Correcting signals for said historical distanceProfile, ->Is Harr function;
calculating covariance function values of all path units corresponding to the target time by using the covariance function and coordinates of all path units corresponding to the target time;
and taking the path unit with the maximum covariance function value as a starting path unit, and determining the coordinates of the starting path unit in the search map.
Preferably, the convolution operation is performed by using a gaussian filter convolution function and coordinates of the initial path unit, and the central path unit and the plurality of target path units are determined from the path units based on a result of the convolution operation, including:
selecting a preset number of adjacent path units as convolution windows in the search map by taking the coordinates of the initial path units as the center;
performing convolution operation based on the vertical velocity value stored in the path unit in the convolution window and a Gaussian filter convolution function to obtain the vertical change velocity of the boundary layer height corresponding to the target moment;
Determining a time interval between the target time and a time previous to or next to the target time;
calculating the product of the time interval and the vertical change speed to obtain a vertical movement distance;
dividing the vertical movement distance by the distance resolution to obtain the number of distance libraries in the vertical movement distance;
taking the sum of the ordinate of the initial path unit and the number of the distance libraries as the ordinate, and determining a central path unit from each path unit corresponding to the last time or the next time of the target time;
obtaining a target ordinate range based on the coordinates of the central path unit, the number of distance libraries in the vertical movement distance and a preset minimum search depth; the target ordinate range is: h c -b t -b min To H c +b t +b min Wherein H is C To search the ordinate of the depth center, b t B is the number of distance bins within the vertical travel distance min The method comprises the steps of setting a preset minimum search depth;
and taking a path unit with the abscissa as the previous time or the next time of the target time and the ordinate in the range of the target ordinate as a target path unit.
Preferably, the setting weights for the respective non-obstacle path units based on the distances from the respective non-obstacle path units to the central path unit, calculating weight allocation values corresponding to the respective non-obstacle path units based on an allocation value formula and the corresponding weights, and taking the non-obstacle path unit with the largest weight allocation value as the reference path unit includes:
Obtaining covariance function values of each non-obstacle path unit by using a wavelet covariance transformation method;
calculating to obtain the distance from each non-obstacle path unit to the central path unit in the search map;
calculating the weight of each non-obstacle path unit by using a Gaussian function and the distance from each non-obstacle path unit to the central path unit;
calculating to obtain weight distribution values of all the non-obstacle path units by using a distribution value formula, covariance function values and weights of all the non-obstacle path units; the distribution value formula is as follows:F(i)=D(i)×W(i) Wherein, the method comprises the steps of, wherein,Fa value is assigned to the weight and,ithe distance library is numbered in order to be able to,Das the weight of the material to be weighed,Wis a covariance function value;
and taking the non-obstacle path unit with the maximum weight distribution value as a reference path unit.
Preferably, the method further comprises:
determining a reference moment at which the height of a boundary layer shows a linear rising trend and has no cloud influence after sunrise;
selecting a preset verification time based on the reference time and a preset time interval;
calculating the boundary layer height corresponding to each verification moment by using a wavelet covariance transformation method;
selecting a median from boundary layer heights corresponding to all verification moments as a verification height;
Inverting to obtain the boundary layer height at the verification moment;
and when the difference value between the boundary layer height at the verification time obtained by inversion and the verification height is in a preset range, determining that the boundary layer height at the inversion time obtained by inversion is accurate.
The second aspect of the embodiment of the invention discloses a boundary layer height inversion device, which comprises:
the determining unit is used for determining target time for starting boundary layer height inversion from a pre-established search map, the search map is pre-constructed based on a history distance correction signal, the horizontal axis of the search map represents time of n preset intervals, the vertical axis of the search map represents k distance library numbers, the search map comprises n multiplied by k path units, each path unit corresponds to a different distance library, each distance library is divided by the history distance correction signal according to the same distance resolution, the noise ratio of the corresponding distance library, the vertical speed value of the corresponding distance library and the history distance correction signal are stored in each path unit, and n and k are positive integers;
a transformation unit for determining a start path unit and coordinates of the start path unit in a search map based on the target time, the history distance correction signal, and a wavelet covariance transformation method;
The convolution unit is used for carrying out convolution operation by utilizing a Gaussian filter convolution function and the coordinates of the initial path unit, and determining a central path unit and a plurality of target path units from the path units corresponding to the target time based on the convolution operation result; setting a target path unit with the noise ratio of the corresponding distance library larger than a threshold value as a non-obstacle path unit;
the distribution unit is used for setting weights for the non-obstacle path units based on the distances from the non-obstacle path units to the central path unit, calculating weight distribution values corresponding to the non-obstacle path units based on a distribution value formula and corresponding weights, and taking the non-obstacle path unit with the largest weight distribution value as a reference path unit;
and the calculating unit is used for determining the moment corresponding to the reference path unit as the inversion moment, determining the distance library number corresponding to the reference path unit, and taking the product of the distance library number and the distance resolution as the boundary layer height of the inversion moment.
Preferably, the transformation unit includes:
a transformation subunit, configured to perform wavelet covariance transformation based on the target time and the historical distance correction signal, to obtain a covariance function; wherein the covariance function is: ，RThe height is represented by the height of the plate,bthe distance library corresponding to the smoothed center position of the wavelet transform is numbered,aas the scale factor of the wavelet transform,R t representing the top height of the historical distance correction signal,R b representing the bottom height of the historical distance correction signal, X%R) Profile for said history distance correction signal, < >>Is Harr function;
a first calculating subunit, configured to calculate, using the covariance function and coordinates of each path unit corresponding to the target time, a covariance function value of each path unit corresponding to the target time;
and the determining subunit is used for taking the path unit with the maximum covariance function value as a starting path unit and determining the coordinates of the starting path unit in the search map.
Preferably, the convolution unit includes:
a selecting subunit, configured to select a preset number of adjacent path units as a convolution window in the search map with the coordinates of the initial path unit as a center;
the convolution subunit is used for carrying out convolution operation based on the vertical speed value stored by the path unit in the convolution window and the Gaussian filter convolution function to obtain the vertical change speed of the boundary layer height corresponding to the target moment;
A second calculating subunit, configured to determine a time interval between the target time and a time previous to or next to the target time; calculating the product of the time interval and the vertical change speed to obtain a vertical movement distance; dividing the vertical movement distance by the distance resolution to obtain the number of distance libraries in the vertical movement distance; taking the sum of the ordinate of the initial path unit and the number of the distance libraries as the ordinate, and determining a central path unit from each path unit corresponding to the last time or the next time of the target time; obtaining a target ordinate range based on the coordinates of the central path unit, the number of distance libraries in the vertical movement distance and a preset minimum search depth; the target ordinate range is: h c -b t -b min To H c +b t +b min Wherein H is C To search the ordinate of the depth center, b t B is the number of distance bins within the vertical travel distance min The method comprises the steps of setting a preset minimum search depth; and taking a path unit with the abscissa as the previous time or the next time of the target time and the ordinate in the range of the target ordinate as a target path unit.
Preferably, the dispensing unit is specifically configured to:
Obtaining covariance function values of each non-obstacle path unit by using a wavelet covariance transformation method; calculating to obtain the distance from each non-obstacle path unit to the central path unit in the search map; calculating the weight of each non-obstacle path unit by using a Gaussian function and the distance from each non-obstacle path unit to the central path unit; calculating to obtain each non-obstacle path unit by using the distribution value formula, covariance function value and weightWeight assignment values for the non-obstacle path units; the distribution value formula is as follows:F(i)=D(i)×W(i) Wherein, the method comprises the steps of, wherein,Fa value is assigned to the weight and,ithe distance library is numbered in order to be able to,Das the weight of the material to be weighed,Wis a covariance function value; and taking the non-obstacle path unit with the maximum weight distribution value as a reference path unit.
Preferably, the method further comprises:
the verification unit is used for determining a reference moment which shows a linear rising trend of the boundary layer height after sunrise and has no cloud influence; selecting a preset verification time based on the reference time and a preset time interval; calculating the boundary layer height corresponding to each verification moment by using a wavelet covariance transformation method; selecting a median from boundary layer heights corresponding to all verification moments as a verification height; inverting to obtain the boundary layer height at the verification moment; and when the difference value between the boundary layer height at the verification time obtained by inversion and the verification height is in a preset range, determining that the boundary layer height at the inversion time obtained by inversion is accurate.
Based on the inversion method and the inversion device for the boundary layer height provided by the embodiment of the invention, the target moment for starting to invert the boundary layer height is determined in a pre-established search map, the search map is pre-constructed based on a history distance correction signal, the horizontal axis of the search map represents the moment of n preset intervals, the vertical axis represents k distance library numbers, each path unit corresponds to a different distance library respectively, each distance library is divided by the history distance correction signal according to the same distance resolution, the noise ratio of the corresponding distance library, the vertical speed value of the corresponding distance library and the history distance correction signal are stored in each path unit, and n and k are positive integers; determining a start path unit and coordinates of the start path unit in a search map based on the target time, the historical distance correction signal and a wavelet covariance transformation method; performing convolution operation by using a Gaussian filter convolution function and coordinates of the initial path unit, and determining a central path unit and a plurality of target path units from each path unit based on the convolution operation result; setting a target path unit with the noise ratio of the corresponding distance library larger than a threshold value as a non-obstacle path unit; setting weights for the non-obstacle path units based on the distances from the non-obstacle path units to the central path unit, calculating weight distribution values corresponding to the non-obstacle path units based on a distribution value formula and corresponding weights, and taking the non-obstacle path unit with the largest weight distribution value as a reference path unit; and determining the moment corresponding to the reference path unit as inversion moment, determining a distance library number corresponding to the reference path unit, and taking the product of the distance library number and the distance resolution as the boundary layer height of the inversion moment. In the scheme, coordinates of a target moment and an initial path unit are determined in a search map comprising a plurality of path units, convolution operation is carried out based on the coordinates of the initial path unit, a non-obstacle path unit corresponding to the boundary layer height is obtained through simulation, weights are distributed to the non-obstacle path unit, a distance library number corresponding to the non-obstacle path unit with the largest weight distribution value is determined, and the product of the distance library number and the distance resolution is used as the boundary layer height, so that the purpose of avoiding the problems of uncertainty in wavelet scale selection and inaccurate results caused by severe jump of the inversion boundary layer height is achieved.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings that are required to be used in the embodiments or the description of the prior art will be briefly described below, and it is obvious that the drawings in the following description are only embodiments of the present invention, and that other drawings can be obtained according to the provided drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a method for inverting boundary layer height according to an embodiment of the present invention;
FIG. 2 is an exemplary diagram of a search map according to an embodiment of the present invention;
FIG. 3 is a schematic diagram of a convolution window according to an embodiment of the present disclosure;
FIG. 4 is a schematic diagram of a target path unit distribution according to an embodiment of the present invention;
FIG. 5 is a schematic diagram of an obstacle path unit according to an embodiment of the invention;
FIG. 6 is a schematic diagram of inverting boundary layer heights at various moments in time according to an embodiment of the present invention;
FIG. 7 is a schematic diagram of a boundary layer height inversion result according to an embodiment of the present invention;
FIG. 8 is a flowchart of a method for verifying boundary layer height inversion results according to an embodiment of the present invention;
Fig. 9 is a block diagram of an inversion apparatus for boundary layer height according to an embodiment of the present application.
Detailed Description
The following description of the embodiments of the present application will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present application, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the application without making any inventive effort, are intended to be within the scope of the application.
In the present disclosure, the terms "comprises," "comprising," or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising one … …" does not exclude the presence of other like elements in a process, method, article, or apparatus that comprises the element.
According to the background technology, in the process of inversion of the boundary layer height, the problem of inaccurate results caused by uncertainty of wavelet scale selection and severe jump of the inversion boundary layer height is solved.
In the scheme, the coordinates of a target moment and an initial path unit are determined in a search map comprising a plurality of path units, convolution operation is performed on the coordinates of the initial path unit, a non-obstacle path unit corresponding to the boundary layer height is obtained through simulation, weights are distributed to the non-obstacle path unit, a distance library number corresponding to the non-obstacle path unit with the largest weight distribution value is determined, and the product of the distance library number and the distance resolution is used as the boundary layer height, so that the purpose of avoiding the problems of uncertainty in wavelet scale selection and inaccuracy in the result caused by severe jump of the inversion boundary layer height is achieved.
Referring to fig. 1, a flow chart of a boundary layer height inversion method according to an embodiment of the present invention is disclosed, and the method mainly includes the following steps:
step S101: and determining the target moment for starting the boundary layer height inversion from a pre-established search map.
In step S101, a search map is constructed in advance based on a history distance correction signal acquired by a coherent wind lidar.
Fig. 2 is a diagram showing an example of searching a map according to an embodiment of the present invention.
Specifically, the horizontal axis of the search map represents the time of n preset intervals, the vertical axis represents k distance library numbers, the search map contains n×k path units, and n and k are positive integers.
Each path unit corresponds to a different distance library, and each distance library is divided according to the same distance resolution based on the historical distance correction signal.
The range bin is a small range unit divided into ranges along the ray direction in radar echo signal processing.
In one embodiment, a history file acquired by the coherent wind lidar is read to obtain a history distance correction signal (Range Corrected Signals, RCS), and a noise ratio (Carrier to Noise Ratio, CNR) corresponding to each distance bin and a vertical velocity value corresponding to each distance bin.
In a specific implementation, because the distance resolutions of the partitioned distance libraries have differences, the distance resolution normalization process is required, for example, the distance resolution of the 1 st to 100 th distance libraries is 30m, the distance resolution of the 101 st to 150 th distance libraries is 60m, the distance resolution of the rest of the remaining distance libraries is 150m, and the distance resolution of each distance library is normalized to 30m.
In the embodiment of the invention, the preferred distance resolution is 30m, the minimum boundary layer height threshold is set to 240m, the corresponding distance bin number is 240/30=8, the minimum boundary layer height threshold is set to 4000m, the corresponding distance bin number is 4000/30=133, that is, the inversion of the boundary layer height is not performed by searching the path units of the part with the distance bin number greater than 133 or less than 8 in the map.
The noise ratio of the corresponding distance library, the vertical velocity value of the corresponding distance library and the historical distance correction signal are stored in each path unit.
In the specific implementation process of step S101, a target time to be subjected to boundary layer height inversion is determined from n times included in a pre-established search map.
Step S102: the coordinates of the start path unit and the start path unit are determined in the search map based on the target time, the historical distance correction signal, and the wavelet covariance transformation method.
In the specific implementation process of step S102, the profile of the historical distance correction signal is corrected based on the target timeX(R) And carrying out wavelet covariance transformation to obtain a covariance function.
Wherein the covariance function is:
（1）
in the formula (1),Rthe height is represented by the height of the plate, bThe distance library corresponding to the smoothed center position of the wavelet transform is numbered,aas the scale factor of the wavelet transform,R t representing the top height of the historical distance correction signal,R b representing the bottom height of the historical distance correction signal, X%R) For the profile of the historical distance correction signal,as a Harr function.
In particular, the method comprises the steps of,
（2）
and calculating covariance function values of the path units corresponding to the target time by using the equation (1) and the ordinate of the path units corresponding to the target time, namely the distance library numbers corresponding to the path units.
And taking the path unit with the maximum covariance function value as a starting path unit, and determining the coordinates of the starting path unit in the search map.
Specifically, as shown in fig. 2, a path unit with a darkened color and marked as T in the figure is a start path unit, and its coordinates are [ ]t 4, b 4 )。
Wherein,,t 4 for the moment of time of the object,b 4 and numbering the distance library corresponding to the initial path unit.
It should be noted that, the product of the distance library number corresponding to the start path unit and the distance resolution is the boundary layer height at the target time.
Step S103: and carrying out convolution operation by utilizing the Gaussian filter convolution function and the coordinates of the initial path unit, and determining a central path unit and a plurality of target path units from the path units corresponding to the target time based on the convolution operation result.
In step S103, the gaussian filter convolution function is:
（3）
wherein, the method comprises the following steps oft,b) Representing the coordinates in the search map,σrepresenting standard deviation.
In the specific implementation process of step S103, a preset number of adjacent path units are selected as convolution windows in the search map with the coordinates of the initial path units as the center.
Preferably, the size of the convolution window is chosen to be 3 x 3.
Fig. 3 is a schematic diagram of a convolution window according to an embodiment of the present disclosure.
Wherein,,Vrepresenting the vertical velocity value stored in the path element,brepresenting the number of the distance library,tas a representative point in time of day,kis the number of the sequence number,V(b,k) For the vertical velocity value stored in the initial path #b,t k ) Is the coordinates of the starting path element.
The standard deviation in the formula (3)σThe calculation rule of (2) is: when calculating the boundary layer height at the next time of the target time, the standard deviation is calculated using 3×3 pairs of right two columns of vertical velocity values in the convolution window, as shown in the left-hand partial graph (a) of fig. 3. When calculating the boundary layer height at the time immediately preceding the target time, the standard deviation is calculated using the two left-hand columns of vertical velocity values in the convolution window corresponding to 3×3, as shown in the right-hand division diagram (b) of fig. 3.
And (3) performing convolution operation based on the vertical speed value stored by the path unit in the convolution window and the formula (3) to obtain the vertical change speed of the boundary layer height at the corresponding target moment.
Specifically, for a convolution window of size (2k+1) × (2k+1) in equation (3), in the Gaussian kerneli,j) The calculation formula of the vertical velocity value at the position is:
（4）
based on the formula (3) and the formula (4), a calculation formula of the vertical change speed is obtained:
（5）
wherein,,Velis of vertical change speedDegree, the degree is%t,b) Representing coordinates in the search map.
And determining a time interval between the target time and the time at the last time or the next time of the target time, and determining a central path unit and a plurality of target path units from the path units based on the vertical change speed and the time interval.
Fig. 4 is a schematic diagram of a target path unit distribution according to an embodiment of the present invention.
The dotted line rectangular frame on the left side of the initial path unit is a target path unit corresponding to the last time of the target time, and the dotted line rectangular frame on the right side of the initial path unit is a target path unit corresponding to the next time of the target time.
Specifically, the product of the time interval and the vertical change speed is calculated to obtain the vertical movement distance, and the vertical movement distance is divided by the distance resolution to obtain the number of distance libraries in the vertical movement distance.
And determining a central path unit with the ordinate as the sum value from each path unit corresponding to the previous time or the next time of the target time by taking the sum value of the ordinate of the initial path unit and the number of the distance libraries as the ordinate.
And obtaining a target ordinate range based on the coordinates of the central path unit, the number of distance libraries in the vertical movement distance and the preset minimum search depth.
Specifically, the target ordinate range is: h c -b t -b min To H c +b t +b min Wherein H is C To search for the ordinate of the depth center, i.e. the ordinate of the central path element, b t B is the number of distance bins within the vertical travel distance min The minimum boundary layer height threshold value specified in the embodiment of the present application corresponds to the distance library number 8.
And taking a path unit with the abscissa as the previous time or the next time of the target time and the ordinate in the range of the target ordinate as a target path unit.
The range occupied by each target path element in the search map, i.e., the search depth range of the boundary layer height, is denoted by H in the present application.
Considering that when the vertical wind speed is too small, the search depth range is small even 0 distance bins, the present application specifies the minimum search depth range to be 8 distance bins, and confirms the movement depth of the boundary layer height in consideration of the vertical movement speed and the elapsed time.
Step S104: and setting a target path unit with the noise ratio of the corresponding distance library being larger than a threshold value as a non-obstacle path unit.
In step S104, the noise ratio of its corresponding distance library is stored in each target path unit.
In the specific implementation process of step S104, a path unit with a noise ratio of nan or less than a threshold is used as an obstacle path unit, B is used for identification in the search map, and a path unit with a noise ratio of not nan and greater than the threshold is used as a non-obstacle path unit.
Wherein the obstacle unit does not participate in subsequent calculations and determination of the boundary layer height arrival path unit. The preferred noise ratio threshold is set to-19 dB.
Fig. 5 is a schematic diagram of an obstacle path unit according to an embodiment of the present invention.
The dark path units in the search map are obstacle units, the search map is marked by B, and the rest path units are non-obstacle path units.
Step S105: weights are set for each of the non-obstacle path units based on the distances of each of the non-obstacle path units from the center path unit.
In the specific implementation of step S105, the covariance function value of each non-obstacle path unit is obtained by using a wavelet covariance transformation method.
And calculating covariance function values of the non-obstacle path units by using the formula (1), the historical distance correction signals and coordinates of the non-obstacle path units in the search map.
The present invention uses wavelet covariance transformation as a reference judgment basis, which can greatly reduce the calculation amount of the reference judgment basis, and the reference judgment basis can be replaced by a threshold method, a gradient method, a standard deviation method, etc.
And calculating the distance from each non-obstacle path unit to the central path unit in the search map.
And calculating the weight of each non-obstacle path unit by using the Gaussian function and the distance from each non-obstacle path unit to the central path unit.
Specifically, the gaussian function is:
（6）
Das the weight of the material to be weighed,H R in order to search for a radius,ito search depth rangesHThe distance library numbers corresponding to the non-obstacle path units in the network,dto search depth center H c Distance to each non-obstacle path element within the search depth range.
Step S106: and calculating a weight distribution value corresponding to each non-obstacle path unit based on the distribution value formula and the corresponding weight, and taking the non-obstacle path unit with the maximum weight distribution value as a reference path unit.
In step S106, the assigned value formula is:
F(i)=D(i)×W(i)（7）
wherein,,Fa value is assigned to the weight and,Das the weight in step S105,Wthe reference discrimination in the implementation of step S105 is based on the calculated discrimination value, i.e. the covariance function value of each non-obstacle path element, iAnd numbering a distance library corresponding to each non-obstacle path unit in the search depth range H.
In the specific implementation process of step S106, the weight distribution value corresponding to each non-obstacle path unit is calculated based on the distribution value formula and the corresponding weight, and the weight distribution value (i.e.FSequence) corresponding to the maximum value in the sequence) of the distance library numberi t As to boundary layer heightI.e. the reference path element, the path is clear.
For example, the coordinates in FIG. 5 are [ ]t 3 ,b 3 ) And%t 6 ,b 3 ) Is a reference path element.
Step S107: and determining the time corresponding to the reference path unit as the inversion time, determining the distance library number corresponding to the reference path unit, and taking the product of the distance library number and the distance resolution as the boundary layer height of the inversion time.
In step S107, the distance resolution is 30m.
For example, in FIG. 5t 3 Andt 6 the time corresponding to the reference path unit is the inversion time,b 3 namely the distance library number corresponding to the reference path unit.
Specifically, the inversion timet 3 Andt 6 the corresponding boundary layer heights are all 30 mXb 3 。
In one embodiment, the period of time to be inverted is presett 0 To the point oft n Wherein each time in the search map contains t 0 To t n The target time in the applicationtIs thatt 0 To the point oft n In (a), i.et 0 ≦t≦t n 。
Fig. 6 is a schematic diagram of inversion of boundary layer heights at various moments according to an embodiment of the present application.
Wherein the broken line connection line represents the boundary layer height, and each path unit through which the connection passes is a path unit up to the boundary layer height, i.e., the reference path unit mentioned in the above embodiment.
Specifically, starting with the target moment, int 0 To the point oft n Traversing the data forward and backward respectively in the time interval, repeating the steps S103 to S107, and finally calculatingt 0 To the point oft n Boundary layer height values corresponding to respective moments in time.
The method for inverting the height of the convenient layer disclosed by the embodiment of the application is shown in fig. 7, and is a schematic diagram of the inversion result of the height of the boundary layer disclosed by the embodiment of the application.
Wherein, the inversion time period corresponding to the boundary layer height inversion result is 0 to 24 hours, which corresponds to the time period to be inverted in the embodiment of the applicationt 0 To the point oft n 。
Based on the inversion method of the boundary layer height disclosed by the embodiment of the application, the coordinates of the initial path unit and the initial path unit are determined in the search map, convolution operation is carried out by utilizing a Gaussian filter convolution function and the coordinates of the initial path unit, a central path unit and a plurality of target path units are determined from all path units based on the convolution operation result, non-obstacle path units are set, weights are set for all non-obstacle path units, weight distribution values corresponding to all non-obstacle path units are calculated, the non-obstacle path unit with the largest weight distribution value is taken as a reference path unit, the moment corresponding to the reference path unit is determined as the inversion moment, the distance library number corresponding to the reference path unit is determined, and the product of the distance library number and the distance resolution is taken as the boundary layer height at the inversion moment. In the scheme, convolution operation is carried out based on the coordinates of the initial path unit, a non-obstacle path unit corresponding to the height change range of the boundary layer is obtained through simulation, weights are distributed to the non-obstacle path unit, the distance library number corresponding to the non-obstacle path unit with the largest weight distribution value is determined, and the product of the distance library number and the distance resolution is used as the height of the boundary layer, so that the purpose of avoiding the problem of inaccurate results caused by uncertainty of wavelet scale selection and severe jump of the height of the inversion boundary layer is achieved.
The method for inverting the boundary layer height disclosed based on the embodiment of the invention is shown in fig. 8, and is a flowchart of a method for verifying the boundary layer height inversion result disclosed by the embodiment of the invention.
The verification method of the boundary layer height inversion result disclosed by the embodiment of the invention is to verify the boundary layer height obtained by inversion in the corresponding embodiment of fig. 1, and mainly comprises the following steps:
step S801: and determining a reference moment which shows a linear rising trend of the boundary layer height after sunrise and has no cloud influence.
In step S801, with the enhancement of solar radiation after sunrise, the underlying surface starts to heat, the heat flux increases, the boundary layer is developed, and the change of the boundary layer height appears to have a good linear trend, so that a reference time at which the boundary layer height shows a linear rising trend and has no cloud influence after sunrise is determined.
It should be noted that the reference time is not equal to the period to be inverted in the corresponding embodiment of fig. 1t 0 To the point oft n Any one of the times.
Step S802: and selecting a preset verification time based on the reference time and a preset time interval.
In the specific implementation process of step S802, 10 verification times are selected at 1 minute intervals within 10 minutes adjacent to the reference time.
Step S803: and calculating the boundary layer height corresponding to each verification moment by using a wavelet covariance transformation method.
In the specific implementation process of step S803, the boundary layer height corresponding to each verification time may also be calculated by using a threshold method, a gradient method and a standard deviation method.
Step S804: and selecting a median from boundary layer heights corresponding to all verification moments as a verification height.
In the specific implementation process of step S804, the maximum value and the minimum value in the boundary layer heights corresponding to 10 verification moments are removed, and the median of the remaining boundary layer heights is taken as the final verification height.
Step S805: inversion is performed to obtain the boundary layer height at the verification time.
In the specific implementation process of step S805, the boundary layer height at the verification time is inverted by using the boundary layer height inversion method disclosed in the corresponding embodiment of fig. 1, so as to obtain the boundary layer height at the verification time.
For example, the step S102 is executed back to obtain the boundary layer height at the verification time, taking the last time of the verification time in the search map as the target time.
Step S806: when the difference value between the boundary layer height and the inspection height of the verification moment obtained by inversion is in a preset range, determining that the boundary layer height of the inversion moment obtained by inversion is accurate.
In step S806, since the minimum distance resolution of the original data is 30m and the maximum distance resolution is 150m, in order to avoid the effect of the resolution of the single distance library, a height of 210m is selected as the basis for effectiveness discrimination, i.e. when the difference between the boundary layer height at the verification time obtained by inversion and the verification height is less than 210m, the boundary layer height inverted in the corresponding embodiment of fig. 1 is considered to be accurate.
According to the verification method for the inversion result of the boundary layer height, which is disclosed by the embodiment of the invention, the preset verification time is selected, the boundary layer height corresponding to each verification time is calculated by utilizing a wavelet covariance transformation method, the median is selected from the boundary layer heights corresponding to each verification time as the verification height, the boundary layer height of the verification time is obtained through inversion, when the difference value between the boundary layer height of the verification time and the verification height obtained through inversion is within the preset range, the boundary layer height of the inversion time obtained through inversion is determined to be accurate, in the scheme, the boundary layer height of the verification time is calculated through the wavelet covariance transformation method, then the boundary layer height of the verification time is obtained through inversion, and when the difference value between the boundary layer height of the verification time and the verification height obtained through inversion is within the preset range, the boundary layer height of the inversion time obtained through inversion is determined to be accurate, so that the purpose of checking the correctness of the inversion result is achieved.
Corresponding to the inversion method of the boundary layer height disclosed in the embodiment of the present invention, as shown in fig. 9, the apparatus is a structure diagram of an inversion apparatus of the boundary layer height disclosed in the embodiment of the present invention, and the apparatus includes: a determination unit 901, a transformation unit 902, a convolution unit 903, an allocation unit 904, and a calculation unit 905.
Wherein, the determining unit 901 is configured to determine a target time for starting to perform boundary layer height inversion from a pre-established search map.
It should be noted that, the search map is constructed in advance based on the history distance correction signal, the horizontal axis of the search map represents the time of n preset intervals, the vertical axis represents k distance library numbers, the search map includes n×k path units, each path unit corresponds to a different distance library, each distance library is divided by the history distance correction signal according to the same distance resolution, the noise ratio of the corresponding distance library, the vertical velocity value of the corresponding distance library and the history distance correction signal are stored in each path unit, and n and k are positive integers.
A transformation unit 902 for determining the coordinates of the start path unit and the start path unit in the search map based on the target time, the history distance correction signal and the wavelet covariance transformation method.
In an embodiment, the transforming unit 902 comprises:
and the transformation subunit is used for carrying out wavelet covariance transformation based on the target moment and the historical distance correction signal to obtain a covariance function.
Wherein the covariance function is the formula (1) in the embodiment of the present invention.
And the first calculating subunit is used for calculating the covariance function value of each path unit corresponding to the target moment by using the covariance function and the coordinates of each path unit corresponding to the target moment.
And the determining subunit is used for taking the path unit with the maximum covariance function value as a starting path unit and determining the coordinates of the starting path unit in the searching map.
The convolution unit 903 is configured to perform convolution operation using a gaussian filter convolution function and coordinates of a start path unit, determine a center path unit and a plurality of target path units from path units corresponding to target moments based on a result of the convolution operation, and set a target path unit with a noise ratio greater than a threshold value in a corresponding distance library as a non-obstacle path unit.
In one embodiment, convolution unit 903 includes:
and the selecting subunit is used for selecting preset adjacent path units in the search map by taking the coordinates of the initial path unit as the center and taking the preset adjacent path units as convolution windows.
The convolution subunit is used for carrying out convolution operation based on the vertical velocity value stored in the path unit in the convolution window and the Gaussian filter convolution function to obtain the vertical change speed of the boundary layer height at the corresponding target moment.
The second calculating subunit is configured to determine a time interval between the target time and a time immediately before or immediately after the target time, calculate a product of the time interval and a vertical change speed, obtain a vertical movement distance, divide the vertical movement distance by a distance resolution, obtain a number of distance banks within the vertical movement distance, determine a central path unit from each path unit corresponding to the time immediately before or immediately after the target time by using a sum of an ordinate of a start path unit and the number of distance banks as an ordinate, and obtain a target ordinate range based on a coordinate of the central path unit, the number of distance banks within the vertical movement distance, and a preset minimum search depth, where the target ordinate range is: h c -b t -b min To H c +b t +b min 。
Wherein H is C To search the ordinate of the depth center, b t B is the number of distance bins within the vertical travel distance min Is a preset minimum search depth.
And taking a path unit with the abscissa as the previous time or the next time of the target time and the ordinate in the range of the target ordinate as a target path unit.
An allocation unit 904, configured to set weights for each of the non-obstacle path units based on distances from each of the non-obstacle path units to the central path unit, calculate weight allocation values corresponding to each of the non-obstacle path units based on an allocation value formula and corresponding weights, and use the non-obstacle path unit with the largest weight allocation value as a reference path unit.
Specifically, the allocation unit 904 is specifically configured to: the covariance conversion method is utilized to obtain covariance function values of all non-obstacle path units, distances from all the non-obstacle path units to a central path unit in a search map are calculated, weights of all the non-obstacle path units are calculated by utilizing a Gaussian function and distances from all the non-obstacle path units to the central path unit, and weight distribution values of all the non-obstacle path units are calculated by utilizing a distribution value formula, covariance function values of all the non-obstacle path units and the weights.
The calculating unit 905 is configured to determine that the time corresponding to the reference path unit is the inversion time, determine the distance library number corresponding to the reference path unit, and use the product of the distance library number and the distance resolution as the boundary layer height at the inversion time.
In one embodiment, the demarcation layer height inversion apparatus further comprises:
the verification unit is used for determining a reference time which shows a linear rising trend and has no cloud influence on the boundary layer height after sunrise, selecting preset verification times based on the reference time and preset time intervals, calculating the boundary layer height corresponding to each verification time by utilizing a wavelet covariance transformation method, selecting a median from the boundary layer heights corresponding to each verification time as a verification height, inverting to obtain the boundary layer height of the verification time, and determining that the boundary layer height of the inversion time obtained by inversion is accurate when the difference value between the boundary layer height of the inversion obtained verification time and the verification height is within a preset range.
The formulas mentioned in the embodiments of the present invention are explained in the above embodiments, and are not repeated here.
Based on the inversion device of boundary layer height disclosed by the embodiment of the invention, the coordinates of the initial path unit and the initial path unit are determined in the search map, convolution operation is carried out by utilizing a Gaussian filter convolution function and the coordinates of the initial path unit, a central path unit and a plurality of target path units are determined from all path units based on the convolution operation result, non-obstacle path units are set, weights are set for all non-obstacle path units, weight distribution values corresponding to all non-obstacle path units are calculated, the non-obstacle path unit with the largest weight distribution value is taken as a reference path unit, the moment corresponding to the reference path unit is determined as the inversion moment, the distance library number corresponding to the reference path unit is determined, and the product of the distance library number and the distance resolution is taken as the boundary layer height at the inversion moment. In the scheme, convolution operation is carried out based on the coordinates of the initial path unit, a non-obstacle path unit corresponding to the height change range of the boundary layer is obtained through simulation, weights are distributed to the non-obstacle path unit, the distance library number corresponding to the non-obstacle path unit with the largest weight distribution value is determined, and the product of the distance library number and the distance resolution is used as the height of the boundary layer, so that the purpose of avoiding the problem of inaccurate results caused by uncertainty of wavelet scale selection and severe jump of the height of the inversion boundary layer is achieved.
In this specification, each embodiment is described in a progressive manner, and identical and similar parts of each embodiment are all referred to each other, and each embodiment mainly describes differences from other embodiments. In particular, for a system or system embodiment, since it is substantially similar to a method embodiment, the description is relatively simple, with reference to the description of the method embodiment being made in part. The systems and system embodiments described above are merely illustrative, wherein the elements illustrated as separate elements may or may not be physically separate, and the elements shown as elements may or may not be physical elements, may be located in one place, or may be distributed over a plurality of network elements. Some or all of the modules may be selected according to actual needs to achieve the purpose of the solution of this embodiment. Those of ordinary skill in the art will understand and implement the present invention without undue burden.
Those of skill would further appreciate that the various illustrative elements and algorithm steps described in connection with the embodiments disclosed herein may be implemented as electronic hardware, computer software, or combinations of both, and that the various illustrative elements and steps are described above generally in terms of functionality in order to clearly illustrate the interchangeability of hardware and software. Whether such functionality is implemented as hardware or software depends upon the particular application and design constraints imposed on the solution. Skilled artisans may implement the described functionality in varying ways for each particular application, but such implementation decisions should not be interpreted as causing a departure from the scope of the present invention.
The previous description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims (10)

1. A method of inverting boundary layer height, comprising:
determining a target time for starting boundary layer height inversion from a pre-established search map, wherein the search map is pre-constructed based on a history distance correction signal, the horizontal axis of the search map represents n times of preset intervals, the vertical axis represents k distance library numbers, the search map comprises n multiplied by k path units, each path unit corresponds to a different distance library, each distance library is divided by the history distance correction signal according to the same distance resolution, and each path unit stores the noise ratio of the corresponding distance library, the vertical speed value of the corresponding distance library and the history distance correction signal, and n and k are positive integers;
Determining a start path unit and coordinates of the start path unit in a search map based on the target time, the historical distance correction signal and a wavelet covariance transformation method;
performing convolution operation by using a Gaussian filter convolution function and coordinates of the initial path unit, and determining a central path unit and a plurality of target path units from each path unit based on the convolution operation result;
setting a target path unit with the noise ratio of the corresponding distance library larger than a threshold value as a non-obstacle path unit;
setting weights for the non-obstacle path units based on the distances from the non-obstacle path units to the central path unit, calculating weight distribution values corresponding to the non-obstacle path units based on a distribution value formula and corresponding weights, and taking the non-obstacle path unit with the largest weight distribution value as a reference path unit;
and determining the moment corresponding to the reference path unit as inversion moment, determining a distance library number corresponding to the reference path unit, and taking the product of the distance library number and the distance resolution as the boundary layer height of the inversion moment.
2. The method of claim 1, wherein the determining the starting path unit and coordinates of the starting path unit in a search map based on the target time, the historical distance correction signal, and a wavelet covariance transformation method comprises:
Performing wavelet covariance transformation based on the target time and the historical distance correction signal to obtain a covariance function; wherein the covariance function is:，Rthe height is represented by the height of the plate,bthe distance library corresponding to the smoothed center position of the wavelet transform is numbered,aas the scale factor of the wavelet transform,R t representing the top height of the historical distance correction signal,R b representing the bottom height of the historical distance correction signal, X%R) Profile for said history distance correction signal, < >>Is Harr function;
calculating covariance function values of all path units corresponding to the target time by using the covariance function and coordinates of all path units corresponding to the target time;
and taking the path unit with the maximum covariance function value as a starting path unit, and determining the coordinates of the starting path unit in the search map.
3. The method of claim 1, wherein the convolving with the gaussian filter convolution function and the coordinates of the starting path unit, determining a center path unit and a plurality of target path units from the path units based on the convolving result, comprises:
selecting a preset number of adjacent path units as convolution windows in the search map by taking the coordinates of the initial path units as the center;
Performing convolution operation based on the vertical velocity value stored in the path unit in the convolution window and a Gaussian filter convolution function to obtain the vertical change velocity of the boundary layer height corresponding to the target moment;
determining a time interval between the target time and a time previous to or next to the target time;
calculating the product of the time interval and the vertical change speed to obtain a vertical movement distance;
dividing the vertical movement distance by the distance resolution to obtain the number of distance libraries in the vertical movement distance;
taking the sum of the ordinate of the initial path unit and the number of the distance libraries as the ordinate, and determining a central path unit from each path unit corresponding to the last time or the next time of the target time;
obtaining a target ordinate range based on the coordinates of the central path unit, the number of distance libraries in the vertical movement distance and a preset minimum search depth; the target ordinate range is: h c -b t -b min To H c +b t +b min Wherein H is C To search the ordinate of the depth center, b t B is the number of distance bins within the vertical travel distance min The method comprises the steps of setting a preset minimum search depth;
and taking a path unit with the abscissa as the previous time or the next time of the target time and the ordinate in the range of the target ordinate as a target path unit.
4. The method according to claim 1, wherein the setting weights for the respective non-obstacle path units based on the distances from the respective non-obstacle path units to the center path unit, calculating weight allocation values corresponding to the respective non-obstacle path units based on allocation value formulas and the respective weights, and taking the non-obstacle path unit with the largest weight allocation value as a reference path unit, comprises:
obtaining covariance function values of each non-obstacle path unit by using a wavelet covariance transformation method;
calculating to obtain the distance from each non-obstacle path unit to the central path unit in the search map;
calculating the weight of each non-obstacle path unit by using a Gaussian function and the distance from each non-obstacle path unit to the central path unit;
calculating to obtain weight distribution values of all the non-obstacle path units by using a distribution value formula, covariance function values and weights of all the non-obstacle path units; the distribution value formula is as follows:F(i)=D(i)×W(i) Wherein, the method comprises the steps of, wherein,Fa value is assigned to the weight and,ithe distance library is numbered in order to be able to,Das the weight of the material to be weighed,Wis a covariance function value;
and taking the non-obstacle path unit with the maximum weight distribution value as a reference path unit.
5. The method according to any one of claims 1 to 4, further comprising:
determining a reference moment at which the height of a boundary layer shows a linear rising trend and has no cloud influence after sunrise;
selecting a preset verification time based on the reference time and a preset time interval;
calculating the boundary layer height corresponding to each verification moment by using a wavelet covariance transformation method;
selecting a median from boundary layer heights corresponding to all verification moments as a verification height;
inverting to obtain the boundary layer height at the verification moment;
and when the difference value between the boundary layer height at the verification time obtained by inversion and the verification height is in a preset range, determining that the boundary layer height at the inversion time obtained by inversion is accurate.
6. An inversion apparatus for boundary layer height, comprising:
the determining unit is used for determining target time for starting boundary layer height inversion from a pre-established search map, the search map is pre-constructed based on a history distance correction signal, the horizontal axis of the search map represents time of n preset intervals, the vertical axis of the search map represents k distance library numbers, the search map comprises n multiplied by k path units, each path unit corresponds to a different distance library, each distance library is divided by the history distance correction signal according to the same distance resolution, the noise ratio of the corresponding distance library, the vertical speed value of the corresponding distance library and the history distance correction signal are stored in each path unit, and n and k are positive integers;
A transformation unit for determining a start path unit and coordinates of the start path unit in a search map based on the target time, the history distance correction signal, and a wavelet covariance transformation method;
the convolution unit is used for carrying out convolution operation by utilizing a Gaussian filter convolution function and the coordinates of the initial path unit, and determining a central path unit and a plurality of target path units from the path units corresponding to the target time based on the convolution operation result; setting a target path unit with the noise ratio of the corresponding distance library larger than a threshold value as a non-obstacle path unit;
the distribution unit is used for setting weights for the non-obstacle path units based on the distances from the non-obstacle path units to the central path unit, calculating weight distribution values corresponding to the non-obstacle path units based on a distribution value formula and corresponding weights, and taking the non-obstacle path unit with the largest weight distribution value as a reference path unit;
and the calculating unit is used for determining the moment corresponding to the reference path unit as the inversion moment, determining the distance library number corresponding to the reference path unit, and taking the product of the distance library number and the distance resolution as the boundary layer height of the inversion moment.
7. The apparatus of claim 6, wherein the transformation unit comprises:
a transformation subunit, configured to perform wavelet covariance transformation based on the target time and the historical distance correction signal, to obtain a covariance function; wherein the covariance function is:，Rthe height is represented by the height of the plate,bthe distance library corresponding to the smoothed center position of the wavelet transform is numbered,aas the scale factor of the wavelet transform,R t representing the top height of the historical distance correction signal,R b representing the bottom height of the historical distance correction signal, X%R) Profile for said history distance correction signal, < >>Is Harr function;
a first calculating subunit, configured to calculate, using the covariance function and coordinates of each path unit corresponding to the target time, a covariance function value of each path unit corresponding to the target time;
and the determining subunit is used for taking the path unit with the maximum covariance function value as a starting path unit and determining the coordinates of the starting path unit in the search map.
8. The apparatus of claim 6, wherein the convolution unit comprises:
a selecting subunit, configured to select a preset number of adjacent path units as a convolution window in the search map with the coordinates of the initial path unit as a center;
The convolution subunit is used for carrying out convolution operation based on the vertical speed value stored by the path unit in the convolution window and the Gaussian filter convolution function to obtain the vertical change speed of the boundary layer height corresponding to the target moment;
a second calculating subunit, configured to determine a time interval between the target time and a time previous to or next to the target time; calculating the product of the time interval and the vertical change speed to obtain a vertical movement distance; dividing the vertical movement distance by the distance resolution to obtain the number of distance libraries in the vertical movement distance; taking the sum of the ordinate of the initial path unit and the number of the distance libraries as the ordinate, and determining a central path unit from each path unit corresponding to the last time or the next time of the target time; obtaining a target ordinate range based on the coordinates of the central path unit, the number of distance libraries in the vertical movement distance and a preset minimum search depth; the target ordinate range is: h c -b t -b min To H c +b t +b min Wherein H is C To search the ordinate of the depth center, b t B is the number of distance bins within the vertical travel distance min The method comprises the steps of setting a preset minimum search depth; and taking a path unit with the abscissa as the previous time or the next time of the target time and the ordinate in the range of the target ordinate as a target path unit.
9. The device according to claim 6, characterized in that said distribution unit is in particular adapted to:
obtaining covariance function values of each non-obstacle path unit by using a wavelet covariance transformation method; calculating to obtain the distance from each non-obstacle path unit to the central path unit in the search map; calculating the weight of each non-obstacle path unit by using a Gaussian function and the distance from each non-obstacle path unit to the central path unit; by dividingCalculating a weight distribution value of each non-obstacle path unit according to a distribution value formula, covariance function values and weights of each non-obstacle path unit; the distribution value formula is as follows:F(i)=D(i)×W(i) Wherein, the method comprises the steps of, wherein,Fa value is assigned to the weight and,ithe distance library is numbered in order to be able to,Das the weight of the material to be weighed,Wis a covariance function value; and taking the non-obstacle path unit with the maximum weight distribution value as a reference path unit.
10. The apparatus according to any one of claims 6 to 9, further comprising:
The verification unit is used for determining a reference moment which shows a linear rising trend of the boundary layer height after sunrise and has no cloud influence; selecting a preset verification time based on the reference time and a preset time interval; calculating the boundary layer height corresponding to each verification moment by using a wavelet covariance transformation method; selecting a median from boundary layer heights corresponding to all verification moments as a verification height; inverting to obtain the boundary layer height at the verification moment; and when the difference value between the boundary layer height at the verification time obtained by inversion and the verification height is in a preset range, determining that the boundary layer height at the inversion time obtained by inversion is accurate.
CN202310778764.2A 2023-06-28 2023-06-28 Inversion method and device for boundary layer height Active CN116500649B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310778764.2A CN116500649B (en) 2023-06-28 2023-06-28 Inversion method and device for boundary layer height

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310778764.2A CN116500649B (en) 2023-06-28 2023-06-28 Inversion method and device for boundary layer height

Publications (2)

Publication Number Publication Date
CN116500649A CN116500649A (en) 2023-07-28
CN116500649B true CN116500649B (en) 2023-10-20

Family

ID=87325338

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310778764.2A Active CN116500649B (en) 2023-06-28 2023-06-28 Inversion method and device for boundary layer height

Country Status (1)

Country Link
CN (1) CN116500649B (en)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20180072584A (en) * 2016-12-21 2018-06-29 목원대학교 산학협력단 Altitude determination method using laser radar
CN108256546A (en) * 2017-11-30 2018-07-06 中国人民解放军国防科技大学 Method and system for detecting height of atmospheric boundary layer under non-precipitation condition
CN110031868A (en) * 2019-05-15 2019-07-19 国耀量子雷达科技有限公司 A method of based on coherent wind laser radar carrier-to-noise ratio inversion boundary layer height
KR102200282B1 (en) * 2019-12-12 2021-01-07 한국외국어대학교 연구산학협력단 System for atmospheric boundary layer heights estimation using ceilometer-derived backscattering coefficient
AU2021106025A4 (en) * 2021-08-19 2021-10-28 PLA Naval Engineering University Comprehensive verification system for modified atmospheric refractivity situation distribution
CN115144871A (en) * 2022-06-27 2022-10-04 兰州大学 Atmospheric boundary layer height identification method and device based on laser radar
CN115980697A (en) * 2022-12-29 2023-04-18 无锡中科光电技术有限公司 Method for inverting boundary layer height by using laser radar under different weather conditions
CN116125494A (en) * 2023-01-09 2023-05-16 武昌首义学院 Boundary layer height inversion method based on dual-wavelength laser radar

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20180072584A (en) * 2016-12-21 2018-06-29 목원대학교 산학협력단 Altitude determination method using laser radar
CN108256546A (en) * 2017-11-30 2018-07-06 中国人民解放军国防科技大学 Method and system for detecting height of atmospheric boundary layer under non-precipitation condition
CN110031868A (en) * 2019-05-15 2019-07-19 国耀量子雷达科技有限公司 A method of based on coherent wind laser radar carrier-to-noise ratio inversion boundary layer height
KR102200282B1 (en) * 2019-12-12 2021-01-07 한국외국어대학교 연구산학협력단 System for atmospheric boundary layer heights estimation using ceilometer-derived backscattering coefficient
AU2021106025A4 (en) * 2021-08-19 2021-10-28 PLA Naval Engineering University Comprehensive verification system for modified atmospheric refractivity situation distribution
CN115144871A (en) * 2022-06-27 2022-10-04 兰州大学 Atmospheric boundary layer height identification method and device based on laser radar
CN115980697A (en) * 2022-12-29 2023-04-18 无锡中科光电技术有限公司 Method for inverting boundary layer height by using laser radar under different weather conditions
CN116125494A (en) * 2023-01-09 2023-05-16 武昌首义学院 Boundary layer height inversion method based on dual-wavelength laser radar

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Relationship analysis of PM2:5 and boundary layer height using an aerosol and turbulence detection lidar;Chong Wang 等;《Atmospheric Measurement Techniques》;第12卷(第06期);3303-3315 *
小波变换在风廓线雷达探测大气边界层高度中的应用研究;魏浩 等;《热带气象学报》;第31卷(第06期);811-820 *
激光雷达反演大气边界层高度的优化方法;于思琪 等;《光学学报》;第41卷(第07期);159-167 *

Also Published As

Publication number Publication date
CN116500649A (en) 2023-07-28

Similar Documents

Publication Publication Date Title
CN106054194B (en) A kind of spaceborne radar and ground-based radar reflectivity factor data three-dimensional fusion method
CN108414998B (en) Satellite laser altimeter echo waveform analog simulation method and device
Zhao et al. Theoretical analysis and numerical experiments of variational adjoint approach for refractivity estimation
CN112505651B (en) Automatic processing method for atmospheric detection laser radar
CN103698305A (en) Method and system for observing atmospheric transmittance in real time
CN116297068B (en) Aerosol optical thickness inversion method and system based on earth surface reflectivity optimization
CN105044039B (en) A kind of method according to laser radar data automatic inversion horizontal visibility
CN112801413A (en) Photovoltaic power station generated power prediction method and device
CN109191408B (en) Rapid circulation ground weather fusion method and device and server
CN111487621B (en) Sea surface flow field inversion method based on radar image and electronic equipment
CN109856608A (en) A kind of high confidence rate detection method of radar target based on recurrence clutter map
CN113740934A (en) Rainfall estimation method based on S-band dual-polarization weather radar
Ehmele et al. Flood-related extreme precipitation in southwestern Germany: development of a two-dimensional stochastic precipitation model
CN116500649B (en) Inversion method and device for boundary layer height
US20210223039A1 (en) Method and device for determining height, electronic device, and computer-readable storage medium
CN106546958B (en) A kind of radar data assimilation method of optimization
CN117848501A (en) Aviation multispectral image radiation correction method, device, equipment and medium
CN113593006A (en) Meteorological data spatial interpolation refining method and system based on deep learning
CN111859255A (en) Radar detection range calculation method under influence of terrain shielding
CN110888142A (en) Spacecraft hidden target point measuring method based on MEMS laser radar measuring technology
CN115015842B (en) Atmospheric refraction correction method for double-station intersection of theodolite without ranging
CN113064130B (en) Method, device, storage medium and program product for determining particle spectral distribution
CN112987010B (en) System and method for multi-radar mapping of robot
CN115755012A (en) Method, device and equipment for training attitude angle prediction model and predicting attitude angle
CN110726982B (en) Signal smoothing method and system

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant