CN113466854B - High-frequency ground wave radar inversion vector flow velocity method based on ocean power model - Google Patents

High-frequency ground wave radar inversion vector flow velocity method based on ocean power model Download PDF

Info

Publication number
CN113466854B
CN113466854B CN202110730290.5A CN202110730290A CN113466854B CN 113466854 B CN113466854 B CN 113466854B CN 202110730290 A CN202110730290 A CN 202110730290A CN 113466854 B CN113466854 B CN 113466854B
Authority
CN
China
Prior art keywords
flow velocity
grid
aver
sea area
steps
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
CN202110730290.5A
Other languages
Chinese (zh)
Other versions
CN113466854A (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.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
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 Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202110730290.5A priority Critical patent/CN113466854B/en
Publication of CN113466854A publication Critical patent/CN113466854A/en
Application granted granted Critical
Publication of CN113466854B publication Critical patent/CN113466854B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M10/00—Hydrodynamic testing; Arrangements in or on ship-testing tanks or water tunnels
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Abstract

A vector flow velocity inversion method for a high-frequency ground wave radar based on an ocean dynamic model belongs to the field of sea state inversion of the high-frequency ground wave radar. The method aims to solve the problem that the single-station high-frequency ground wave radar can only obtain the radial flow velocity during sea state inversion. The method comprises the steps of obtaining the water depth of an observation sea area, carrying out grid division on the sea surface according to the water depth and the radar beam position, determining a flow velocity expression mode corresponding to grid points, and mutually converting an x/y orthogonal basis and a group of radial/tangential orthogonal bases which express the flow velocity through a conversion relation; then obtaining the radial flow velocity and wave height of a selected sea area as an initial field, determining a driving condition for carrying out updating calculation on time advance on the initial field, then carrying out time step advance, carrying out wave height updating on the whole sea area, carrying out radial flow velocity updating on the whole sea area, and carrying out tangential flow velocity updating on the whole sea area; until all time steps are updated; and finally, synthesizing the radial flow velocity and the tangential flow velocity into a vector flow velocity. The method is mainly used for inverting the sea surface vector flow velocity.

Description

High-frequency ground wave radar inversion vector flow velocity method based on ocean power model
Technical Field
The invention belongs to the field of high-frequency ground wave radar sea state inversion, and relates to a vector flow velocity inversion method for a high-frequency ground wave radar.
Background
The high-frequency ground wave radar can invert sea state information such as ocean storm flow fields and the like by utilizing a echo spectrum received by the radar for sea detection. Sea state inversion is mature and accurate in the aspect of inverting flow fields, but for a single-station radar, the sea state inversion can only detect the radial velocity on a wave beam and cannot detect the vector flow velocity; in the field, wind field, wave field and flow field are usually extracted and inverted respectively, and each sea state parameter is relatively independent, so that the constraint relation among parameters in marine dynamics and hydrodynamics is not utilized in pure sea state inversion.
Disclosure of Invention
The method aims to solve the problem that only the radial flow velocity can be obtained when the high-frequency ground wave radar is used for sea state inversion.
The vector flow velocity inversion method for the high-frequency ground wave radar based on the ocean dynamic model comprises the following steps of:
the method comprises the following steps: acquiring the water depth of an observed sea area, and performing grid division on the sea surface according to the water depth and the radar beam position;
step two: determining a flow velocity representation mode corresponding to each grid point in a selected sea area, wherein the flow velocity is a two-dimensional vector and can be represented by a group of x/y orthogonal bases and a group of radial/tangential orthogonal bases, and two groups of orthogonal bases can be mutually converted through a conversion relation;
step three: acquiring the radial flow velocity and wave height of a selected sea area as an initial field;
step four: determining drive conditions for a time-marching up-to-date computation of the initial field: step of time
Figure BDA0003139042470000011
Wherein h is max The maximum water depth in the model domain, wherein delta x and delta y are the sizes of the grid points in the x direction and the y direction, and g is the gravity acceleration;
step five: advancing according to time steps, and updating the wave height of the whole sea area;
step six: advancing according to time steps, and updating the radial flow velocity of the whole sea area;
step seven: advancing according to time steps, and updating the tangential flow velocity of the whole sea area;
step eight: and (3) continuing to advance calculation according to time steps, and updating by adopting a staggered updating method: firstly, repeating the fifth step, the seventh step and the sixth step; repeating the fifth step, the sixth step and the seventh step when the next time step is advanced to calculate; namely, when the flow velocity of the adjacent time step is updated, the radial flow velocity and the tangential flow velocity are alternately updated until all the time steps are updated;
step nine: the radial flow velocity and the tangential flow velocity are combined into a vector flow velocity.
Further, the process of determining the flow rate representation corresponding to each grid point in the selected sea area includes the following steps:
firstly, calculating an angle theta of grid points of a grid in a polar coordinate position according to the relative positions of the grid and a radar station, wherein the theta is a negative direction included angle between a set beam and a parallel weft;
the radial flow velocity on the wave beam is Vr, and the direction far away from the radar station is the positive direction of the wave beam; the tangential flow velocity on the wave beam is Vs, and the direction which forms an acute angle with the positive direction of the parallel latitude line is the positive direction of the wave beam; u and v are flow velocities in the positive directions of the parallel weft and the warp respectively; thus, the following switching relationship exists between the two sets of orthogonal bases:
u＝-Vrcosθ+Vssinθ
v＝Vrsinθ+Vscosθ
Vr＝-ucosθ+vsinθ
Vs＝usinθ+vcosθ。
further, the θ is obtained by calculating an arc tangent of a ratio of a latitude difference between the latitude of the grid point position and the latitude of the radar station, and a longitude difference between the longitude of the radar station and the longitude of the grid point position.
Further, the process of obtaining the radial flow velocity and the wave height of the selected sea area as the initial field comprises the following steps of carrying out spatial interpolation on the wave height and the radial flow velocity on the radar wave beam by adopting a cubic spline interpolation method, and further obtaining wave height and radial flow velocity data of the whole sea area as the initial field of numerical calculation.
Further, in the process of determining the driving condition for performing the time advance update calculation on the initial field in step four, if the time interval corresponding to the two times of data acquired by the radar does not meet the condition, performing time smoothing linear interpolation processing on the two times of adjacent data.
Further, after the sea surface is gridded, each grid type is determined, which specifically comprises the following steps:
dividing grid points of the selected sea area into three categories, and marking the categories by type masks: if the land part type mask without water depth is set to be 0, the sea state parameter of the lattice point is represented; the default type mask of the sea with water depth is 1, which indicates that sea state parameters at the lattice points need to be calculated; the grid point type mask at the observation beam position is 2, which indicates that the water level data and the radial flow velocity at the position are obtained by direct observation.
Further, the flow rate mask of each grid is determined while the grid type is determined, ctrl and ctrl v are flow rate masks in x direction and y direction of the grid, respectively, ctrl and ctrl v need to be calculated as 1, and do not need to be calculated as 0.
Further, the specific process of the step five comprises the following steps:
firstly, converting the radial flow velocity Vr and the tangential flow velocity Vs into a flow velocity u and a flow velocity v according to the conversion relation of the step two;
traversing all grid points with the type mask being 1, and then determining the wave height of the next time step according to the following formula; wherein h is the water depth at the lattice point position;
Figure BDA0003139042470000031
Figure BDA0003139042470000032
Figure BDA0003139042470000033
Figure BDA0003139042470000034
Figure BDA0003139042470000035
Δ x- Δ y-2 π R × spatial resolution/360
Wherein i and j represent the corresponding values of grid points in the x direction and the y direction;
Figure BDA0003139042470000036
showing the wave height of grid points at the i and j positions at the n +1 moment,
Figure BDA0003139042470000037
the flow velocity u and the flow velocity v, h of grid points at the i and j positions at n moments r 、h l 、h u 、h d Are respectively intermediate variables, h j,i The water depth of grid points at the i and j positions is shown, and R is the earth radius.
Further, the specific process of the step six includes the following steps:
traversing all grid points with the type mask being 1, and determining the radial flow velocity of the next time step according to the following formula;
Figure BDA0003139042470000038
Figure BDA0003139042470000039
f vr viscosity term ＝-f u viscous item cosθ+f v item of viscosity sinθ
Figure BDA00031390424700000310
Figure BDA00031390424700000311
Vs_aver＝u_aver×sinθ+v_aver×cosθ
Figure BDA00031390424700000312
Figure BDA0003139042470000041
Figure BDA0003139042470000042
Wherein A is a fluid viscosity coefficient, K is a bottom friction coefficient, and a is a coefficient of a semi-implicit semi-explicit differential format;
Figure BDA0003139042470000043
vr of grid points at the i and j positions at n moments; vs aver 2 Is the square of Vs _ aver, and Vs _ aver, v _ aver and u _ aver are intermediate variables; f. of Vr viscosity term 、f Vr bottom friction 、f Vr advection item Is an intermediate variable representing the right side of the respective equation; ctrl v j,i 、ctrlu j,i And the i and j positions are respectively the ctrl v and ctrl u of the grid point.
Further, the specific process of the seventh step includes the following steps:
traversing all grid points with masks of different types not being 0, and determining the tangential flow velocity of the next time step according to the following formula;
f vs viscosity term ＝f u viscous item sinθ+f v item of viscosity cosθ
Figure BDA0003139042470000044
Figure BDA0003139042470000045
Vr_aver＝-u_aver×cosθ+v_aver×sinθ
Figure BDA0003139042470000046
Figure BDA0003139042470000047
Figure BDA0003139042470000048
Wherein the content of the first and second substances,
Figure BDA0003139042470000049
vs of grid points at the i and j positions at the n moments; vr _ aver 2 The square of Vr _ aver is represented, and the Vr _ aver is an intermediate variable; f. of Vs viscosity term 、f Vs bottom friction 、f Vs advection term Are intermediate variables that represent the right hand side of the respective equations.
Has the advantages that:
the method is based on the characteristics of ocean dynamics, utilizes a two-dimensional shallow water wave model to couple various sea state parameters, and provides a method for solving the ocean current vector flow velocity by using the coupling relation by a single-station radar.
Drawings
FIG. 1 is a schematic diagram of a vector flow velocity inversion method of a high-frequency ground wave radar based on two-dimensional shallow water waves according to the present invention;
FIG. 2 is a schematic diagram of a method for expressing two orthogonal bases of vector flow rate;
FIG. 3 is a schematic diagram of a radar beam for a validation experiment;
FIG. 4 is a graph of average wave height error as calculated by time step advance;
FIG. 5 is a graph of mean radial flow rate error as calculated by time step advance;
FIG. 6 is a graph of average tangential flow velocity error as calculated by time step advance;
FIG. 7 is a plot of radar raw inversion wave height data; wherein FIG. 7(a) is a perspective view and FIG. 7(b) is a plan view;
FIG. 8 is a graph of radar raw inversion radial flow velocity data;
FIG. 9 is a plot of the initial field wave height after interpolation of the radar original inversion wave height data; wherein FIG. 9(a) is a perspective view and FIG. 9(b) is a plan view;
FIG. 10 is a radial velocity plot of the initial field wave obtained after the radar original inversion radial velocity data is interpolated;
FIG. 11 is a subsequent time step height distribution graph calculated by the method; wherein FIG. 11(a) is a perspective view and FIG. 11(b) is a plan view;
FIG. 12 is a vector velocity distribution diagram at a subsequent time step calculated by the method;
fig. 13 is a schematic diagram of a time-stepping mode updating process.
Detailed Description
The first embodiment is as follows:
the embodiment is a vector flow velocity inversion method for a high-frequency ground wave radar based on an ocean power model, and aims to solve the problem that the high-frequency ground wave radar can only obtain radial flow velocity during sea state inversion. The method takes the radial velocity and wave height of the radar as input, and can solve the vector velocity.
Specifically, as shown in fig. 1, the method for inverting the vector flow velocity by using the high-frequency ground wave radar based on the ocean dynamic model according to the embodiment includes the following steps:
the method comprises the following steps: acquiring the water depth of an observation sea area, carrying out mesh division on the sea surface according to the water depth and the radar beam position, and determining the type of each mesh;
the process of determining the mesh type includes the steps of:
dividing grid points of the selected sea area into three categories, marking the categories through type masks (j, i): if the land part type mask without water depth is set to be 0, the sea state parameter of the lattice point is represented; the default type mask of the sea with water depth is 1, which indicates that sea state parameters at the lattice points need to be calculated; the grid point type mask at the observation beam position is 2, indicating that the water level data and the radial flow velocity at that position are obtained by direct observation and are not obtained by calculation. Thereby labeling all sea area lattice points as three categories.
Simultaneously generating a flow rate mask: the ctrl u and the ctrl v are flow rate masks in x and y directions of the grid, respectively, and the ctrl u and the ctrl v need to be calculated as 1 or 0. The unneeded computation refers to when an edge corresponding to land is encountered, the edge is a flow mask of 0.
if(mask(j，i)＝＝0||mask(j，min(i+1，numX))＝＝0)
ctrlU(j，i)＝0；
if(mask(j，i)＝＝0||mask(min(j+1，numY)，i)＝＝0)
ctrlV(j，i)＝0；
Step two: determining a flow speed representation mode corresponding to each grid point in the selected sea area; the flow velocity is a two-dimensional vector that can be represented by a set of x/y orthogonal bases and a set of radial/tangential orthogonal bases. The method comprises the following specific steps:
firstly, calculating an angle theta of grid points under the corresponding polar coordinate position of a grid according to the relative position of the grid and a radar station, wherein the theta is a negative direction included angle between a set wave beam and a parallel weft, and is shown in figure 2; theta is obtained by calculating the arctangent of the ratio of the latitude difference and the longitude difference, wherein the latitude difference is the difference between the latitude of the grid point position and the latitude of the radar station, and the longitude difference is the difference between the longitude of the radar station and the longitude of the grid point position.
The radial flow velocity on the wave beam is Vr, and the direction far away from the radar station is the positive direction of the wave beam; the tangential flow velocity on the wave beam is Vs, and the direction forming an acute angle with the positive direction of the parallel weft is the positive direction; u and v are flow velocities in positive directions of parallel weft and warp threads (or x direction and y direction) respectively; thus, the following switching relationship exists between the two sets of orthogonal bases:
u＝-Vrcosθ+Vssinθ
v＝Vrsinθ+Vscosθ
Vr＝-ucosθ+vsinθ
Vs＝usinθ+vcosθ
vflow is a two-dimensional vector representing horizontal flow velocity of ocean current, and is required to be represented by u, v or Vr, Vs;
step three: obtaining the radial flow velocity and wave height of a selected sea area as an initial field, and specifically comprising the following steps:
and (3) carrying out spatial interpolation on the wave height and the radial flow velocity on the radar wave beam by adopting a cubic spline interpolation method (cubic), and further obtaining wave height and radial flow velocity data of the whole sea area as an initial field of numerical calculation.
Step four: determining a driving condition for performing time-marching up-to-date calculation on the initial field, which comprises the following specific steps:
when calculating the difference equation, the difference in time needs to be matched with the spatial resolution, i.e. the corresponding time step needs to be selected to satisfy the CFL condition
Figure BDA0003139042470000061
Wherein h is max The maximum water depth in the model domain, wherein delta x and delta y are the sizes of the grid points in the x direction and the y direction, and g is the gravity acceleration;
if the time interval corresponding to the two times of data acquired by the radar does not meet the condition, performing time smoothing linear interpolation processing on the two times of adjacent data, namely: for data with larger time interval between two observations, a linear interpolation method is used to enable the time interval to meet the condition. And taking the wave height and the radial flow velocity on the corresponding radar beam under the time step scale as driving conditions of time step advancing updating calculation.
Step five: advancing according to time steps, and updating the wave height of the whole sea area, wherein the method comprises the following specific steps:
firstly, converting the radial flow velocity Vr and the tangential flow velocity Vs into a flow velocity u and a flow velocity v according to the conversion relation of the step two;
traversing all grid points with the type mask being 1, and then determining the wave height of the next time step according to the following formula; wherein h is the water depth at the lattice point position;
Figure BDA00031390424700000710
Figure BDA0003139042470000071
Figure BDA0003139042470000072
Figure BDA0003139042470000073
Figure BDA0003139042470000074
Δ x- Δ y-2 π R × spatial resolution/360
Wherein i and j represent the corresponding values of grid points in the x direction and the y direction;
Figure BDA0003139042470000075
showing the wave height of grid points at the i and j positions at the n +1 moment,
Figure BDA0003139042470000076
the flow velocity u and the flow velocity v, h of grid points at the i and j positions at n moments r 、h l 、h u 、h d Are respectively intermediate variables, h j,i Representing the water depth of grid points at the positions i and j, wherein R is the radius of the earth; the spatial resolution is the size of the divided unit space, and may be set as required, in some embodiments, 1/6 longitudes and latitudes may be set, and in some embodiments, 0.02 longitudes and latitudes may be set;
step six: advancing according to time steps, and updating the radial flow velocity of the whole sea area, wherein the method comprises the following specific steps:
traversing all grid points with the type mask being 1, and determining the radial flow velocity of the next time step according to the following formula;
Figure BDA0003139042470000077
Figure BDA0003139042470000078
f vr viscosity term ＝-f u viscous item cosθ+f v item of viscosity sinθ
Figure BDA0003139042470000079
Figure BDA0003139042470000081
Vs_aver＝u_aver×sinθ+v_aver×cosθ
Figure BDA0003139042470000082
Figure BDA0003139042470000083
Figure BDA0003139042470000084
Wherein A is a fluid viscosity coefficient, K is a bottom friction coefficient, and a is a coefficient of a semi-implicit semi-explicit differential format;
Figure BDA0003139042470000085
vr of grid points at the i and j positions at the n moments; vs aver 2 Is the square of Vs _ aver, and Vs _ aver, v _ aver and u _ aver are intermediate variables; f. of Vr viscosity term 、f Vr bottom friction 、f Vr advection item Can be directly thought of as an intermediate variable to the right of the respective equation; ctrl v j,i 、ctrlu j,i The i position grid point and the j position grid point are respectively ctrl v and ctrl u;
step seven: the method comprises the following steps of advancing according to time steps, and updating the tangential flow velocity of the whole sea area, wherein the method comprises the following specific steps:
traversing all grid points with masks of different types not being 0, and determining the tangential flow velocity of the next time step according to the following formula;
f vs viscosity term ＝f u viscous item sinθ+f v item of viscosity cosθ
Figure BDA0003139042470000086
Figure BDA0003139042470000087
Vr_aver＝-u_aver×cosθ+v_aver×sinθ
Figure BDA0003139042470000088
Figure BDA0003139042470000089
Figure BDA00031390424700000810
Wherein the content of the first and second substances,
Figure BDA0003139042470000091
vs of grid points at the i and j positions at the n time points; vr _ aver 2 Representing the square of Vr _ aver, wherein Vr _ aver is an intermediate variable; f. of Vs viscosity term 、f Vs bottom friction 、f Vs advection term Can be directly thought of as an intermediate variable to the right of the respective equation;
step eight: the calculation is continuously advanced according to the time step, the updating method adopts a staggered updating method, and the specific steps are as follows:
firstly, repeating the fifth step, the seventh step and the sixth step; repeating the fifth step, the sixth step and the seventh step when the next time step is advanced to calculate; i.e. updating the flow rate at the adjacent time step, an interleaved update of the radial flow rate and the tangential flow rate is performed. The time step advancing pattern is shown in fig. 13 until all time steps are updated;
step nine: according to the step two conversion method, the radial flow velocity and the tangential flow velocity are combined into a vector flow velocity Vflow.
Example 1
Now, the yellow sea area (122.83-125.83E, 34.17-37.67N) is selected, and the M2 tide division simulation result is taken as an example, because the tide process is a periodic process and has symmetry particularly about the middle time, the first 1-60 time steps are selected to carry out the propulsion experiment in the total 120 time steps. The spatial resolution is 1/6 latitudes and longitudes. The tidal cycle is 12 hours, divided into 120 time steps, so Δ t is 360 s. In this case, the CFL condition is satisfied, and temporal interpolation processing in the fourth step is not required.
A is the fluid viscosity coefficient, here taking the value of 1000. K is a bottom friction coefficient, the bottom friction coefficient K of the Bohai sea yellow sea is a value within [0.001,0.002], and the value is 0.0018. Alpha is a coefficient of a semi-implicit semi-explicit differential format, and is 0.5.
Given the initial state of the full range of M2 tidal radial flow velocities and wave heights as the initial conditions, the observations on the simulated radar beam at subsequent time steps are the driving conditions. The test was carried out in the sea area of the yellow sea, which was selected from (122.83 ° -125.83 ° E, 34.17-37.67 ° N). No land part in the sea area, namely all areas have sea state parameters.
Now, assume that there are five radar beams in the area, each beam of which is 0 °, 26 °, 45 °, 63 ° and 90 °, respectively, and the schematic diagram of the radar beams in the verification experiment is shown in fig. 3;
the projection of the vector flow velocity of the grid points on the beam in the radial direction is calculated by taking the M2 tidal process as a reference, and the value is set as the radial flow velocity on the beam observed by the radar in the radar application scene. The wave height on the beam is the exact wave height at the corresponding beam position during the tide.
The corresponding tangential flow velocity is calculated, thereby solving the problem from the radial flow velocity to the vector flow velocity. As shown in fig. 4, 5, and 6, the initial state sea area has accurate values of wave height and radial flow velocity, and as the time step advances, the accuracy of wave height and radial flow velocity in the sea area decreases, while the accuracy of tangential flow velocity, which cannot be sensed originally, increases. The wave height during this process is in the order of 0.5m and the radial and tangential flow velocities are in the order of 0.15 m/s. Therefore, at the expense of a small fraction of the radial flow velocity and wave height indicators, the accuracy of the tangential flow velocity is traded for improvement.
Example 2
In this embodiment, actual radar data of 03, 26 and 2019 are selected for analysis. A is the fluid viscosity coefficient, here taking the value of 1000. K is a bottom friction coefficient, the bottom friction coefficient K of the Bohai sea yellow sea is a value within [0.001,0.002], and the value is 0.0018 in the text. Alpha is a coefficient of a semi-implicit semi-explicit differential format, and is 0.5.
And selecting the space resolution of 0.02 longitude and latitude according to the distance resolution of the radar. The raw data are shown in fig. 7 and 8. The initial field obtained by spatial interpolation of the beam data using the cubic method is shown in fig. 9 and 10.
With a 24 second time step, 20 sets of data were linearly interpolated between two adjacent batches of data. After the observation data are used as the input of the two-dimensional shallow hydrodynamics model, the flow velocity and the wave height of a resolution unit between wave beams can be obtained through calculation. The calculated wave height and flow velocity effect of the resolution cell are shown in fig. 11 and 12. The first batch of data as input in this case has a high wave height around 121.3 ° longitude and 37.6 ° latitude. The higher wave height of the sheet area can be diffused to the lower wave height of the periphery to form ocean current flowing to the periphery; this phenomenon is consistent with the constraints on marine dynamics.
In conclusion, the method can effectively solve the problem that the tangential flow velocity cannot be sensed in the high-frequency ground wave radar sea state inversion, and is beneficial to offshore information monitoring.
The present invention is capable of other embodiments and its several details are capable of modifications in various obvious respects, all without departing from the spirit and scope of the present invention.

Claims (7)

1. The method for inverting the vector flow velocity by the high-frequency ground wave radar based on the ocean dynamic model is characterized by comprising the following steps of:
the method comprises the following steps: acquiring the water depth of an observation sea area, and meshing the sea surface according to the water depth and the radar beam position;
step two: determining a flow velocity representation mode corresponding to each grid point in a selected sea area, wherein the flow velocity is a two-dimensional vector and can be represented by a group of x/y orthogonal bases and a group of radial/tangential orthogonal bases, and two groups of orthogonal bases can be mutually converted through a conversion relation;
step three: acquiring the radial flow velocity and wave height of a selected sea area as an initial field;
step four: determining a drive condition for a time-marching update calculation on the initial field: step of time
Figure FDA0003781934590000011
Wherein h is max The maximum water depth in the model domain, wherein delta x and delta y are the sizes of the grid points in the x direction and the y direction, and g is the gravity acceleration;
step five: advancing according to time steps, and updating the wave height of the whole sea area; the specific process comprises the following steps:
firstly, converting the radial flow velocity Vr and the tangential flow velocity Vs into a flow velocity u and a flow velocity v according to the conversion relation of the step two;
traversing all grid points with the type mask being 1, and then determining the wave height of the next time step according to the following formula; wherein h is the water depth at the lattice point position;
Figure FDA0003781934590000012
Figure FDA0003781934590000013
Figure FDA0003781934590000014
Figure FDA0003781934590000015
Figure FDA0003781934590000016
Δ x- Δ y-2 π R × spatial resolution/360
Wherein i and j represent the corresponding values of grid points in the x direction and the y direction;
Figure FDA0003781934590000017
showing the wave height of grid points at the i and j positions at the n +1 moment,
Figure FDA0003781934590000018
the flow velocity u and the flow velocity v, h of grid points at the i and j positions at n moments r 、h l 、h u 、h d Are respectively intermediate variables, h j,i Representing the water depth of grid points at the positions i and j, wherein R is the radius of the earth;
step six: advancing according to time steps, and updating the radial flow velocity of the whole sea area; the specific process comprises the following steps:
traversing all grid points with the type mask being 1, and determining the radial flow velocity of the next time step according to the following formula;
Figure FDA0003781934590000019
Figure FDA0003781934590000021
f vr viscosity term ＝-f u viscous item cosθ+f v item of viscosity sinθ
Figure FDA0003781934590000022
Figure FDA0003781934590000023
Vs_aver＝u_aver×sinθ+v_aver×cosθ
Figure FDA0003781934590000024
Figure FDA0003781934590000025
Figure FDA0003781934590000026
Wherein A is a fluid viscosity coefficient, K is a bottom friction coefficient, and a is a coefficient of a semi-implicit semi-explicit differential format;
Figure FDA0003781934590000027
vr of grid points at the i and j positions at n moments; vs aver 2 Is the square of Vs _ aver, and Vs _ aver, v _ aver and u _ aver are intermediate variables; f. of Vr viscosity term 、f Vr bottom friction 、f Vr advection item Are intermediate variables representing the right of the respective equations; ctrl v j,i 、ctrlu j,i The i position grid point and the j position grid point are respectively ctrl v and ctrl u;
step seven: advancing according to time steps, and updating the tangential flow velocity of the whole sea area; the specific process comprises the following steps:
traversing all grid points with masks of different types not being 0, and determining the tangential flow velocity of the next time step according to the following formula;
f vs viscosity term ＝f u viscous item sinθ+f v item of viscosity cosθ
Figure FDA0003781934590000028
Figure FDA0003781934590000029
Vr_aver＝-u_aver×cosθ+v_aver×sinθ
Figure FDA0003781934590000031
Figure FDA0003781934590000032
Figure FDA0003781934590000033
Wherein, the first and the second end of the pipe are connected with each other,
Figure FDA0003781934590000034
vs of grid points at the i and j positions at the n time points; vr _ aver 2 Representing the square of Vr _ aver, wherein Vr _ aver is an intermediate variable; f. of Vs viscosity term 、f Vs bottom friction 、f Vs advection term Are intermediate variables representing the right of the respective equations;
step eight: and (3) continuing to advance calculation according to time steps, and updating by adopting a staggered updating method: firstly, repeating the fifth step, the seventh step and the sixth step; repeating the fifth step, the sixth step and the seventh step when the next time step is advanced to calculate; namely, when the flow velocity of the adjacent time step is updated, the radial flow velocity and the tangential flow velocity are alternately updated until all the time steps are updated;
step nine: the radial flow velocity and the tangential flow velocity are combined into a vector flow velocity.
2. The marine power model-based high-frequency ground wave radar inversion vector flow velocity method according to claim 1, wherein the process of determining the flow velocity representation mode corresponding to each grid point in the selected sea area comprises the following steps:
firstly, calculating an angle theta of grid points under the corresponding polar coordinate position of a grid according to the relative position of a grid and a radar station, wherein the theta is a negative direction included angle between a set wave beam and a parallel weft;
the radial flow velocity on the wave beam is Vr, and the direction far away from the radar station is the positive direction of the wave beam; the tangential flow velocity on the wave beam is Vs, and the direction forming an acute angle with the positive direction of the parallel weft is the positive direction; u and v are flow velocities in the positive directions of the parallel weft and the warp respectively; thus, the following switching relationship exists between the two sets of orthogonal bases:
u＝-Vrcosθ+Vssinθ
v＝Vrsinθ+Vscosθ
Vr＝-ucosθ+vsinθ
Vs＝usinθ+vcosθ。
3. the method for inverting vector flow velocity by using high-frequency ground wave radar based on an ocean dynamic model as claimed in claim 2, wherein θ is obtained by calculating the arctangent of the ratio of the difference in latitude between the grid point position and the radar station latitude, and the difference in longitude between the radar station longitude and the grid point position.
4. The method for inverting the vector flow velocity by the high-frequency ground wave radar based on the ocean dynamic model according to claim 3, wherein the process of obtaining the radial flow velocity and the wave height of the selected sea area as the initial field comprises the following steps of performing spatial interpolation on the wave height and the radial flow velocity on the radar wave beam by adopting a cubic spline interpolation method, and further obtaining wave height and radial flow velocity data of the whole sea area as the initial field of numerical calculation.
5. The method for inverting the vector flow velocity by the high-frequency ground wave radar based on the ocean power model as recited in claim 4, wherein in the process of determining the driving condition for performing the time-marching up-and-down updating calculation on the initial field in the step four, if the time interval corresponding to two times of data acquired by the radar does not meet the condition, the two times of adjacent data are subjected to time-smoothing linear interpolation processing.
6. The method for inverting vector flow velocity by using high-frequency ground wave radar based on the marine power model as claimed in one of claims 1 to 5, wherein each grid type is determined after the sea surface is subjected to grid division, and the method specifically comprises the following steps:
dividing grid points of the selected sea area into three categories, and marking the categories by type masks: the land part type mask without water depth is set to be 0, and the grid point sea-state-free parameter is represented; the default type mask of the sea with water depth is 1, which indicates that sea state parameters at the lattice points need to be calculated; the grid point type mask at the observation beam position is 2, which indicates that the water level data and the radial flow velocity at the position are obtained by direct observation.
7. The method of claim 6, wherein the grid type is determined and the flow mask of each grid is determined, the flow masks are respectively set to x-direction and y-direction of the grid, and the flow masks are set to 1 when the flow masks are needed to be calculated and are not needed to be calculated to 0 when the flow masks are needed to be calculated.
CN202110730290.5A 2021-06-29 2021-06-29 High-frequency ground wave radar inversion vector flow velocity method based on ocean power model Active CN113466854B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110730290.5A CN113466854B (en) 2021-06-29 2021-06-29 High-frequency ground wave radar inversion vector flow velocity method based on ocean power model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110730290.5A CN113466854B (en) 2021-06-29 2021-06-29 High-frequency ground wave radar inversion vector flow velocity method based on ocean power model

Publications (2)

Publication Number Publication Date
CN113466854A CN113466854A (en) 2021-10-01
CN113466854B true CN113466854B (en) 2022-09-30

Family

ID=77874017

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110730290.5A Active CN113466854B (en) 2021-06-29 2021-06-29 High-frequency ground wave radar inversion vector flow velocity method based on ocean power model

Country Status (1)

Country Link
CN (1) CN113466854B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114812710B (en) * 2022-06-28 2022-10-04 北京海兰信数据科技股份有限公司 Vector flow synthesis method and system for radar radial flow
CN114997250A (en) * 2022-08-01 2022-09-02 北京海兰信数据科技股份有限公司 Radar radial flow post-processing method and system
CN115877377B (en) * 2022-12-12 2023-10-03 中山大学 Radar networking vector flow field synthesis method, system, equipment and storage medium
CN117393062B (en) * 2023-12-13 2024-02-23 上海交通大学四川研究院 Simulation method for rigid chemical reaction flow rollback self-adaptive semi-hidden semi-explicit coupling time

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0412248A2 (en) * 1989-08-10 1991-02-13 STN ATLAS Elektronik GmbH Method for the passive determination of target data
CN105182317A (en) * 2015-08-20 2015-12-23 电子科技大学 Resource management method based on search pattern of centralized MIMO radar
CN107843895A (en) * 2017-10-20 2018-03-27 厦门市气象灾害防御技术中心(海峡气象开放实验室、厦门市避雷检测技术中心) A kind of Dual-Doppler weather radar dimensional wind inversion method

Family Cites Families (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4152895A (en) * 1978-02-21 1979-05-08 Lockheed Corporation Wave powered motor
CN101788683B (en) * 2009-12-29 2012-05-23 华东师范大学 Tsunami motion forecasting method based on multi-hierarchy interaction
DK3060866T3 (en) * 2013-10-24 2019-12-16 Spx Flow Tech Danmark A/S Gas distributor for a convection dryer with improved radial gas velocity control
CN103604944B (en) * 2013-12-11 2015-05-27 哈尔滨工业大学 Surface flow measurement method based on monostation shipborne high-frequency ground wave radar
US9122934B2 (en) * 2013-12-27 2015-09-01 Automotive Research & Testing Center Object detection method with a rising classifier effect and object detection device with the same
US11125866B2 (en) * 2015-06-04 2021-09-21 Chikayoshi Sumi Measurement and imaging instruments and beamforming method
CN105445730B (en) * 2015-11-27 2017-09-15 南京信息工程大学 A kind of Sea Current inverting Spaceborne SAR System and its method based on angle diversity
CN106093936B (en) * 2016-08-29 2018-08-17 中船重工鹏力（南京）大气海洋信息系统有限公司 Sweep the unrestrained stream information extracting method under pattern slowly based on coherent radar
CN109856628B (en) * 2019-01-11 2023-09-26 中国船舶集团有限公司第七二四研究所 Target three-dimensional acceleration motion model calculation method based on scanning radar
DE112019006380B4 (en) * 2019-01-24 2022-12-08 Mitsubishi Electric Corporation STATUS PREDICTION DEVICE AND STATUS PREDICTION METHOD
CN109884337B (en) * 2019-03-05 2021-01-19 哈尔滨工业大学 Method for detecting sea surface wind direction by using high-frequency ground wave radar
CN110455493A (en) * 2019-09-11 2019-11-15 中国船舶科学研究中心（中国船舶重工集团公司第七0二研究所） A kind of model test method of nearly islands and reefs target wave spectrum inverting input wave spectrum
CN112231969A (en) * 2020-09-15 2021-01-15 北京工业大学 System and method for demonstrating assembly sequence planning of swing angle milling head

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0412248A2 (en) * 1989-08-10 1991-02-13 STN ATLAS Elektronik GmbH Method for the passive determination of target data
CN105182317A (en) * 2015-08-20 2015-12-23 电子科技大学 Resource management method based on search pattern of centralized MIMO radar
CN107843895A (en) * 2017-10-20 2018-03-27 厦门市气象灾害防御技术中心(海峡气象开放实验室、厦门市避雷检测技术中心) A kind of Dual-Doppler weather radar dimensional wind inversion method

Also Published As

Publication number Publication date
CN113466854A (en) 2021-10-01

Similar Documents

Publication Publication Date Title
CN113466854B (en) High-frequency ground wave radar inversion vector flow velocity method based on ocean power model
CN111950211B (en) Seabed foundation local scouring depth determination method and system based on ocean current design flow velocity
CN102004856B (en) Rapid collective Kalman filtering assimilating method for real-time data of high-frequency observation data
CN107944608B (en) Sea surface drift and oil spill drift diffusion forecasting method based on satellite remote sensing
CN107341341B (en) A kind of river mouth point source sudden water pollution event source tracing method
CN106990404A (en) A kind of autoscale algorithm using X-band radar inverting sea wave height of navigating
CN110851790B (en) Ocean current force optimization forecasting method based on deep learning algorithm
CN112113545B (en) Inner wave amplitude inversion method based on multi-dimensional sea surface information
Winter Macro scale morphodynamics of the German North Sea coast
CN103389076A (en) Submarine topography change detection and analysis method based on mesh reconstruction
CN102253385A (en) Ocean internal wave forecast method based on synthetic aperture radar image and internal wave model
CN104040378B (en) Weather prognosis device and Predictive meteorological methods
CN104933291A (en) Method for the production of mean sea surface height products based on satellite altimeter data network function interpolation
CN102589528B (en) Multi-temporal imaging island shoreline surveying method
CN104091065A (en) Intermittent flow numerical simulation method for solving shallow water problem
CN111859748A (en) Ocean internal wave simulation method based on vertical hybrid coordinates
CN104977583A (en) Method for X-band radar wave retrieval based on empirical orthogonal decomposition
CN102508946B (en) Method for simulating spilled oil sea surface under finite water depth
CN111951204B (en) Sea surface wind speed inversion method for Tiangong No. two detection data based on deep learning
CN109785443A (en) A kind of three-dimensional model simplifying method for large ocean engineer equipment
CN102867184B (en) Extraction method for sea ice motion features in SAR (synthetic aperture radar) images
CN102930146B (en) Method for quantitatively evaluating fidelity precision of digital elevation model
CN103778282A (en) Effective channel width design method based on maximum transverse flow cumulative frequency
CN103745118A (en) Geomagnetic abnormal data meshing method based on magnetic dipole equivalent source method
CN106097408A (en) A kind of coastline key element continuous multi-scale expression 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