Sea wave observation method based on binocular camera

Abstract

The invention provides a sea wave observation method based on a binocular camera, and belongs to the technical field of image processing, and the method comprises the steps: S1, enabling the binocular camera to shoot a sea wave image, and generating sea wave point cloud data at a current moment t; s2, converting the sea wave point cloud data in the S1 into sea wave point cloud data in a geodetic coordinate system; s3, performing quality screening and optimization on the sea wave data in the S2; the method comprises the steps of selecting a region of interest, performing quality screening, converting noise data points into null values, performing null value filling and performing Gaussian filtering. The quality of the sea wave data is divided into excellent, good and poor through quality screening; s4, calculating the wave height, wavelength and period of the sea wave at the moment based on the sea wave point cloud data with high quality and good quality; and S5, repeating the steps S1 to S4 to obtain multi-frame sea wave point cloud data at different moments, and calculating effective sea wave height, wavelength, period and sea wave spectrum. According to the method, the null value can be filled, so that the overall data is more reasonable, and the obtained elevation data is more accurate.

Description

Technical Field

[0001] The present invention belongs to the technical field of image processing, and in particular relates to an ocean wave observation method based on a binocular camera. Background Art

[0002] Ocean waves are a very common ocean phenomenon, and the observation of ocean wave data is of great significance to marine disaster prevention and mitigation.

[0003] Wave buoys are small buoy systems used to measure ocean waves. They are an in-situ measurement technology. They can remotely measure parameters such as wave period, direction, height, directional spectrum, and power spectrum. Common methods include GPS and gravity wave measurement. However, these methods have limitations, such as limited battery life and difficulty deploying them in distant waters.

[0004] Radar remote sensing wave detection involves transmitting radio waves to the sea surface, receiving the echoes, and using digital images to simulate sea surface conditions and calculate wave parameters. However, this method suffers from insufficient image data simulation, a wide sampling range, and limited resolution. It is generally suitable for medium- and high-scale waves, and has limited applicability to small-scale waves.

[0005] Satellite wave measurement technology, similar to radar, uses satellite equipment to remotely monitor ocean waves. Currently, two types of satellite remote sensing equipment are used to observe ocean surface wind and wave information: radar altimeters measure significant wave height; and synthetic aperture radars (SARs) measure and infer wave directional spectra and significant wave height. However, these methods have limitations, such as accuracy affected by long equipment distances, significant weather influences, discontinuous sampling, and resolution limitations.

[0006] In summary, there are currently numerous wave measurement methods, each with limited applicability. Buoy deployment is difficult in the open ocean, and radar and satellites are not suitable for measuring small-scale waves. Using binocular vision to observe waves could effectively address this problem. However, existing binocular vision technology remains in the research phase and suffers from low efficiency and inability to perform real-time observations and parameter inversion. To address this issue, improvements to existing technologies are urgently needed to address the shortcomings of traditional observation methods and existing binocular vision technology. Summary of the Invention

[0007] In order to solve the above technical problems, the present invention proposes a wave observation method based on a binocular camera.

[0008] To achieve the above objectives, the present invention adopts the following technical content: a wave observation method based on a binocular camera, comprising the following steps:

[0009] S1: The binocular camera captures the wave image and generates the wave point cloud data at the current time t [X ct ,Y ct ,Z ct ]; Wave point cloud data[X ct ,Yct ,Z ct ] is established in the left camera coordinate system;

[0010] S2: Based on the IMU inertial measurement unit, the wave point cloud data [X ct ,Y ct ,Z ct ]Convert to the ocean wave point cloud data in the geodetic coordinate system[X w ,Y w ,Z w ];

[0011] S3: The ocean wave data in the geodetic coordinate system in S2 [X w ,Y w ,Z w ] to perform quality screening and optimization; quality screening and optimization include: selecting regions of interest, quality screening, converting noise data points to null values, null value filling, and Gaussian filtering;

[0012] Quality screening will wave data [X w ,Y w ,Z w The quality of the data is divided into: excellent, good and poor; the poor quality wave data is eliminated [X w ,Y w ,Z w ], only retaining good quality and good wave data [X w ,Y w ,Z w ]Proceed to the next step;

[0013] S4: Based on high-quality and good wave point cloud data [X w ,Y w ,Z w ], calculate the wave height, wave wavelength and wave period at that moment;

[0014] S5: loop S1-S4, with each loop time advancing to obtain multiple frames of wave point cloud data at different times, and calculate the effective wave height, effective wave wavelength, effective wave period and wave spectrum based on the multiple frames of wave point cloud data at different times.

[0015] Furthermore, the selection of the region of interest and quality screening in step S3 includes the following steps:

[0016] S3.1: Select the region of interest. The formula is:

[0017] X1=int(x1*Height)

[0018] X2=int(x2*Height)

[0019] Y1=int(y1*Width)

[0020] Y2=int(y2*Width)

[0021] Among them, x1, x2, y1, y2 are the row and column ratios of the area of ​​interest, Height and Width are the resolutions of the images captured by the camera, that is, the number of pixels in rows and columns, X1, X2, Y1 and Y2 are the index values ​​of the selected wave point cloud matrix data, that is, the wave point cloud data of interest is intercepted according to [X1:X2,Y1:Y2];

[0022] S3.2: Count the number of noise and missing values ​​in the area of ​​interest; set quality screening rules and perform screening, and judge the quality of the wave point cloud data based on the screening results [X w ,Y w ,Z w ] quality, remove the wave point cloud data with poor quality, and retain the wave point cloud data with good quality for the following steps;

[0023] S3.2 specifically includes the following steps:

[0024] S3.2.1: Obtain wave elevation data for each data point in the area of ​​interest;

[0025] Assume that there are S data points in the area of ​​interest, including S1 data points with elevation data and S2 data points without elevation data; S = S1 + S2; the data points without elevation data are "null values".

[0026] The wave height data of the i-th data point is recorded as Z i , i∈[1,S1], calculate the noise value of the i-th data point, the formula is:

[0027]

[0028] Where μ and σ are the mean and standard deviation of the wave height data of S1 “data points” respectively; A i is the noise value of the i-th data point; if |A i |>3, then the i-th data point is recorded as a "noise data point"; traverse S data points and obtain the number of "noise data points" in S1 "data points", recorded as M, M∈[1,S1];

[0029] While traversing S data points, the number of null values ​​S2 mentioned above is also obtained;

[0030] S3.2.2: Calculate the ratio of the sum of the number of blank values ​​and the number of noise points to the total number of data points in the region of interest. This ratio is referred to as P. The formula is:

[0031] P=(M+S2) / S

[0032] S3.2.3: Set the quality screening rule Judge to judge the wave point cloud data [X w ,Y w ,Z w ] quality;

[0033] The judgment formula is:

[0034]

[0035] When Judge is -1, it indicates that the wave point cloud data of this frame [X w ,Y w ,Z w ]The quality is poor, and the wave point cloud data of this frame should be discarded;

[0036] When Judge is 1, it indicates that the wave point cloud data of this frame [X w ,Y w ,Z w ]Excellent quality;

[0037] When Judge is 0, it indicates that the wave point cloud data of this frame [X w ,Y w ,Z w ]Good quality; retain the wave point cloud data of excellent quality and good quality for subsequent calculations.

[0038] Furthermore, step S3 further includes the following steps:

[0039] S3.3: If the wave point cloud data [X w ,Y w ,Z w ] is judged as excellent or good quality, and the wave point cloud data is optimized [X w ,Y w ,Z w ] in the noise data points and null values;

[0040] S3.3 specifically includes the following steps:

[0041] S3.3.1: Use row and split dilution factors to convert the wave point cloud data [X w ,Y w ,Z w ] dilute the diluted wave point cloud data [X w ,Y w ,Z w ] are converted to null values;

[0042] Set the diluted [X w ,Y w,Z w ] has a total of r1 null values ​​and r2 non-null values; r1 includes the null values ​​in the original wave point cloud data and the total number of null values ​​converted from “noise data points”;

[0043] Let the coordinates of the bth null value be (x ib ,y jb );

[0044] S3.3.2: With null coordinates (x ib ,y jb ) is the center, and the side length R is the search range; if there is a valid value in the search range, the empty value needs to be filled; if there is no valid value in the search range, the empty value is not filled; the valid value is the data point with elevation data mentioned above;

[0045] The filling process of step S3.3.2 includes the following steps:

[0046] S3.3.2.1: Calculate the Euclidean distance from the null value to each non-null value;

[0047] The formula is:

[0048]

[0049] In the formula, (x f ,y f ) represents the coordinates of the fth non-null value; f∈[1,r2];

[0050] S3.3.2.2: Calculate the weight of the null value;

[0051] The formula is:

[0052]

[0053] In the formula, p is the set weight index;

[0054] S3.3.2.3: Fill in the gaps with the following formula:

[0055]

[0056] Where Z valid_values It is the selected valid non-empty value Z coordinate data. When calculating the X and Y axis data, just replace it with the corresponding axis data;

[0057] S3.3.2.4: Gaussian filtering;

[0058] Gaussian function:

[0059]

[0060] Construct a kernel matrix of size (2k+1)×(2k+1) and calculate the weight of each matrix position:

[0061]

[0062] Where i and j are the coordinates of each matrix data relative to the center of the kernel, i = -k, -k+1, ..., k-1, k, j = -k, -k+1, ..., k-1, k;

[0063] In order to retain the information of subsequent elevation data, the Gaussian kernel is normalized to ensure that the sum of the filter kernels is 1:

[0064]

[0065] The original wave elevation matrix Z w Expand k pixels symmetrically around the periphery to generate the expansion matrix Z wk ; Original wave elevation matrix Z w That is, the wave elevation data after gap filling;

[0066] For the wave elevation data matrix Z w Each original data point Z w (x,y), as the kernel center, calculate the smoothed value:

[0067]

[0068] That is, the Gaussian kernel is covered at the corresponding position of the extended matrix, multiplied point by point and summed up; finally, the extended boundary part is removed, and the smoothed result with the same size as the original data is retained for subsequent operations.

[0069] There is also a lack of Gaussian filtering here

[0070] Furthermore, the wave point cloud data [X w ,Y w ,Z w ] The wave height, wavelength and period of ocean waves include the following steps:

[0071] S4.1: Transform the wave point cloud data [X w ,Y w ,Z w ] is considered as a matrix (M t ,N t ,3); then the matrix (M t ,N t ) is further compressed to obtain (M t , 1) vector G t ; Use the zero-crossing method to extract the vector G t The maximum and minimum elevation values ​​are denoted as Hmax and H min；

[0072] S4.1: Calculate the wave height at the current time t;

[0073] The formula is:

[0074] H t =H max -H min ;

[0075] Where H t is the wave height at the current moment t;

[0076] S4.2: Calculate the wavelength of the ocean wave at the current time t;

[0077] Get H max and H min Then, according to the vector G t The index in can be used to query the corresponding H max and H min The data of the waves on the X and Y axes; for the uncompressed data matrix, there are M t data, there are corresponding N along the X direction t The number of data, assuming H max The corresponding data on the X-axis and Y-axis are [X max_1 、X max_2 …、X max_Nt ]、[Y max_1 、Y max_2 …、Y max_Nt ], H min The corresponding data on the X-axis and Y-axis are [X min_1 、X min_2 …、X min_Nt ]、[Y min_1 、Y min_2 …、Y min_Nt ]; thus the wavelength of the wave at the current moment t can be calculated; the formula is:

[0078]

[0079] Calculate the wave point cloud data[X w ,Y w ,Z w ] is regarded as the wavelength of the waves at the current time t; the formula is:

[0080]

[0081] Where, L t is the wavelength of the ocean wave at the current moment t;

[0082] S4.3: Calculate the period of the ocean wave at the current time t;

[0083] There are two situations:

[0084] Case 1: For deep water waves with a water depth greater than or equal to half the wavelength;

[0085] The calculation formula for the wave period is:

[0086]

[0087] Where c is the wave velocity; g is the acceleration due to gravity; L t is the wavelength of the ocean waves at the current time t; T t is the wave period at the current time t;

[0088] Case 2: For shallow water waves where the water depth is less than half of the wavelength;

[0089] The calculation formula for the wave period is:

[0090]

[0091] Where c is the wave velocity; g is the acceleration due to gravity; L t is the wavelength of the ocean wave at the current time t; h is the water depth; T t is the wave period at the current time t.

[0092] Furthermore, the calculation of significant wave height, significant wavelength, significant period and wave spectrum includes the following steps:

[0093] S5.1: Continue to cycle through S1-S4 multiple times, with each cycle time increasing. The binocular camera continuously captures images, and the computer connected to the binocular camera continuously generates wave point cloud data at different times. The computer then continuously calculates and obtains wave height, wave wavelength, and wave period data for the wave point cloud data at different times.

[0094] Assume that n sets of wave data are generated. Sort the wave heights in the n sets of data from large to small, and add the sorted data into set A in the above order. Set A is expressed as:

[0095] A={{H (1) 、L (1) , T (1)},{H (2) 、L (2) , T (2)},…,{H (n) 、L (n) , T (n)}}, where H (1) ≥H (2) ≥…≥H (n) ;

[0096] set up Select {{H (1) , L (1) , T (1)},{H (2) , L (2) , T (2)},…,{H (K) , L (K) , T (K)}};

[0097] S5.2: Calculate the effective wave height, effective wavelength, and effective period of ocean waves;

[0098] Effective wave height H s The formula is:

[0099]

[0100] Effective wavelength L s The formula is:

[0101]

[0102] Effective period T s The formula is:

[0103]

[0104] Furthermore, in step S5, calculating the wave spectrum includes the following steps:

[0105] S5.3.1: Set the loop S1-S4 multiple times to generate u-frame wave point cloud data;

[0106] S5.3.2: In the u-frame ocean wave point cloud data, select a fixed point whose position remains unchanged in each frame of ocean wave point cloud data. Obtain the elevation of the same position in each frame of ocean wave point cloud data, denoted as H x1 、H x2 ,……,H xu ; The generation time of each frame of wave point cloud data is set to: t b1 , t b2 ,……,t bu ;but:

[0107] t b1 to t bu The wave spectrum of the time series segment is:

[0108]

[0109] Where ω is the angular frequency of the wave. P(ω) is the power of the wave. τ is the time offset.

[0110] The formula for τ is:

[0111]

[0112] By employing the above-mentioned technical solution, the present invention achieves the following beneficial effects: It utilizes a binocular camera to collect wave elevation data, a continuous, contactless, and high-density observation method. This method enables continuous wave observation and capture of more detailed wave information without affecting the wave's motion. Furthermore, the use of an IMU (Inertial Measurement Unit) ensures that captured data remains reliable even when the camera is in motion or slightly offset, making it possible to mount the binocular camera observation platform on a moving platform. Furthermore, to address issues such as gaps and noisy data captured by the binocular camera, quality screening and optimization algorithms are employed to ensure and improve data quality, resulting in higher accuracy in subsequent calculations of wave spectra and statistical parameters. This allows the overall algorithm to achieve fully automatic and efficient data observation and parameter inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0113] Figure 1 is a flow chart of the present method;

[0114] Figure 2 It is the ocean wave matrix compression and index map. DETAILED DESCRIPTION

[0115] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of them. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0116] like Figure 1 As shown, a wave observation method based on a binocular camera specifically includes the following steps:

[0117] S1: The binocular camera captures ocean wave images and generates ocean wave point cloud data;

[0118] Aim the binocular camera at the wave area, avoid direct sunlight, and conduct observations in a well-lit environment.

[0119] After the binocular camera is turned on, the left and right cameras of the binocular camera simultaneously capture images of the waves. Through existing means such as image acquisition, image filtering, automatic correction, stereo matching, and three-dimensional reconstruction, the computer generates wave point cloud data.

[0120] S2: Based on the IMU inertial measurement unit, the wave point cloud data generated in the left camera coordinate system is converted to the wave point cloud data in the geodetic coordinate system.

[0121] The wave point cloud data generated in S1 is the wave point cloud data in the binocular camera coordinate system. Different binocular cameras have different coordinate systems, so for the same wave, the wave point cloud data generated by different binocular cameras will also be different. Therefore, it is necessary to convert the wave point cloud data in the binocular camera coordinate system into the unified geodetic coordinate system to facilitate measurement and evaluation under the same standard.

[0122] S2 specifically includes the following steps:

[0123] S2.1: Determine the rotation matrix R t .

[0124] The coordinate system of the left camera in the binocular camera has its optical center as the origin, the optical axis direction as the Z axis, the direction pointing to the optical center of the right camera as the X axis, and the Y axis is determined according to the right-hand rule.

[0125] The IMU built into the left camera will initialize a geodetic coordinate system with the optical center of the left camera as the origin and the Z axis parallel to and opposite to gravity. The IMU inertial measurement unit measures the deflection angles of the left camera coordinate system to the geodetic coordinate system along the Z, X, and Y axes, which are denoted as θ, ω. The IMU inertial measurement unit can also measure the translation vector The left camera rotates around the Z, X, and Y axes by an angle θ. ω, the corresponding rotation matrices R1, R2 and R3 will be obtained; the calculation formula is:

[0126]

[0127] Thus, the rotation matrix R of the left camera is obtained.

[0128] The formula is:

[0129] R=R1R2R3

[0130] S2.2: The wave point cloud data generated in S1 is converted to the geodetic coordinate system.

[0131] Set the current time t, at which a frame of wave point cloud data of the left camera coordinate system is generated [X ct ,Y ct ,Z ct ], convert it into point cloud data in the geodetic coordinate system [X w ,Y w ,Z w ].

[0132] The formula is:

[0133]

[0134] In formula (1), [X ct ,Y ct ,Z ct ] is the wave point cloud data in the left camera coordinate system; [X w ,Y w ,Z w ] is the wave point cloud data of the corresponding frame in the geodetic coordinate system, where Z w It represents the real elevation data of the waves.

[0135] S3: The wave point cloud data [X w ,Y w ,Z w ] to perform quality screening and optimization. Quality screening and optimization include: selecting regions of interest, quality screening, null value filling, and Gaussian filtering.

[0136] S3 specifically includes the following steps:

[0137] S3.1: The wave point cloud data of the frame at the current time t [X w ,Y w ,Z w ]Select the area of ​​interest.

[0138] X1=int(x1*Height)

[0139] X2=int(x2*Height)

[0140] Y1=int(y1*Width)

[0141] Y2=int(y2*Width)

[0142] Among them, x1, x2, y1, y2 (<1) are the row and column ratios of the area of ​​interest, Height and Width are the resolutions of the images captured by the camera, that is, the number of row and column pixels, X1, X2, Y1 and Y2 are the index values ​​of the selected wave point cloud matrix data, that is, the wave point cloud data of the area of ​​interest in the rectangle is intercepted according to [X1:X2,Y1:Y2].

[0143] Selecting the region of interest can reduce the amount of data for subsequent operations and facilitate the rapid determination of the wave point cloud data [X w ,Y w ,Z w The region of interest is selected manually.

[0144] S3.2: Count the number of noise and missing values ​​in the region of interest; set quality screening rules and perform screening, and judge the quality of the wave point cloud data of this frame based on the screening results. w ,Y w ,Z w ] quality.

[0145] S3.2 specifically includes the following steps:

[0146] S3.2.1: Obtain wave elevation data for each “data point” in the area of ​​interest.

[0147] Assume there are S data points in the region of interest. These S data points include S1 data points with elevation data and S2 data points without elevation data. S = S1 + S2. Data points without elevation data are called "null values."

[0148] The wave height data of the i-th (i∈[1,S1]) data point is recorded as Z i , calculate the noise value of the i-th data point, the formula is:

[0149]

[0150] In formula (2), μ and σ are the mean and standard deviation of the wave height data of S1 “data points”, respectively. i is the noise value of the i-th data point, if |A i |>3, then the i-th data point is recorded as a "noise data point"; through the above threshold comparison method, traverse S data points and obtain the number of noise data in S1 "data points", which is set to M, M∈[1,S1].

[0151] While traversing S data points, the number S2 of the aforementioned null values ​​is also obtained.

[0152] S3.2.2: Calculate the ratio of the sum of the number of blank values ​​and the number of noise points to the total number of data points in the region of interest. This ratio is recorded as P. The formula is:

[0153] P=(M+S2) / S

[0154] S3.2.3: Set the quality screening rule Judge to judge the wave point cloud data [X w ,Y w ,Z w ] quality.

[0155] The formula for the judgment rule is:

[0156]

[0157] When Judge is -1, it indicates that the wave point cloud data of this frame [X w ,Y w ,Z w ] is of poor quality, the frame of ocean wave point cloud data should be discarded; that is, if the current time t of the frame of ocean wave point cloud data [X w ,Y w ,Z w ] If the quality is poor, it will not be used in the subsequent calculation of effective wave height, effective wavelength and effective period;

[0158] When Judge is 1, it indicates that the wave point cloud data of this frame [X w ,Y w ,Z w ]Excellent quality;

[0159] When Judge is 0, it indicates that the wave point cloud data of this frame [X w ,Y w ,Z w ]Good quality.

[0160] This means that if the wave point cloud data of the frame at the current time t [X w ,Y w ,Z w ]Excellent and good quality will be used for subsequent calculation of effective wave height, effective wavelength and effective period;

[0161] If the wave point cloud data [X w ,Y w ,Z w ] is judged as excellent or good quality, and the wave point cloud data is optimized [X w ,Y w ,Z w ] in the noise data points and null values.

[0162] S3.3 specifically includes the following steps:

[0163] S3.3.1: Use row and column dilution factors to convert the wave point cloud data [X w ,Y w ,Z w ] to dilute the diluted wave point cloud data [X w ,Y w ,Z w ] are converted to null values;

[0164] Set diluted wave point cloud data [X w ,Y w ,Z wThere are r1 null values ​​and r2 non-null values. Here r1 includes the null values ​​in the original wave point cloud data and the total number of null values ​​converted from “noise data points”.

[0165] Let the coordinates of the bth null value be (x ib ,y jb ), b∈[1,r1].

[0166] The row and column dilution factors are existing methods in image processing. For example, if the row and column dilution factors are both 0.5, it means that the wave point cloud data [X w ,Y w ,Z w ]Discarding half of the data points in both rows and columns can reduce the overall data size.

[0167] S3.3.2: With null coordinates (x ib ,y jb ) is the center, and the side length R is the search range; if there is a valid value in the search range, the empty value needs to be filled; if there is no valid value in the search range, the empty value is not filled; the valid value is the data point with elevation data mentioned above;

[0168] The process of filling in blanks in S3.3.2 includes the following steps:

[0169] S3.3.2.1: Compute the Euclidean distance from the null value to each non-null value.

[0170] The formula is:

[0171]

[0172] In the formula, (x f ,y f ) represents the coordinates of the fth non-null value; f∈[1,r2].

[0173] S3.3.2.2: Calculate the weight of the null value to each non-null value.

[0174]

[0175] Where p is the set weight index.

[0176] S3.3.2.3: Fill in the gaps with the following formula:

[0177]

[0178] Where Z valid_values It is the selected valid non-empty value of the Z coordinate data. When calculating the X and Y axis data, just replace it with the corresponding axis data.

[0179] S3.4: Use a Gaussian filter to filter the wave point cloud data filled with null values ​​in S5.

[0180] The wave point cloud data in S3.4 is diluted in rows and columns in the aforementioned S3.3.

[0181] The wave height data after gap filling is smoothed using Gaussian smoothing filter; the Gaussian function is defined in two-dimensional space as:

[0182]

[0183] Where σ is the standard deviation of the Gaussian distribution; x, y are the offsets of the filter center;

[0184] Therefore,

[0185] Construct a kernel matrix of size (2k+1)×(2k+1) and calculate the weight of each matrix position:

[0186]

[0187] Where i and j are the coordinates of each matrix data relative to the center of the kernel, i = -k, -k+1, ..., k-1, k, j = -k, -k+1, ..., k-1, k.

[0188] In order to retain the information of subsequent elevation data, the Gaussian kernel is normalized to ensure that the sum of the filter kernels is 1:

[0189]

[0190] To prevent edge data loss, the original wave elevation matrix Z w Expand k pixels symmetrically around the periphery to generate the expansion matrix Z wk The original wave elevation matrix Z w That is, the wave elevation data after the gaps are filled.

[0191] For the wave elevation data matrix Z w Each original data point Z w (x,y), as the kernel center, calculate the smoothed value:

[0192]

[0193] This involves overlaying the Gaussian kernel over the corresponding position of the expanded matrix, multiplying the points point by point, and summing them. Finally, the extended boundary is removed, leaving the smoothed result with the same size as the original data for subsequent operations.

[0194] S4: Since the data obtained by the binocular camera has noise and vacancies, which will affect the feasibility and accuracy of subsequent parameter calculations, the wave point cloud data [Xw ,Y w ,Z w ], calculate the wave height H at the current time t t 、wave wavelength L t and the wave period T t。

[0195] S4 specifically includes the following steps:

[0196] S4.1: Since waves propagate in one direction along their propagation direction, the characteristics of waves in a certain direction do not change with the change of the position in that direction. That is, in the wave point cloud data, the waves propagate along the Y axis and their elevation changes along the Z axis. The waves are considered to be uniform in the X axis direction. Based on this, it can be assumed that the wave point cloud data [X w ,Y w ,Z w ] is a wave (M t ,N t ,3) matrix. M t Indicates the data length of the wave along the Y direction at the current time t, N t Indicates the data length of the wave along the X axis at the current time t, such as Figure 2 shown.

[0197] Then the matrix (M t ,N t ) is further compressed to obtain (M t , 1) vector G t ; G t That is the change curve along the Y-axis direction.

[0198] S4.2: Use the upper zero point method to identify the vector G t The maximum and minimum elevation values ​​are denoted as H max and H min , thus, we can calculate the wave point cloud data [X w ,Y w ,Z w The wave height of the waves in ] is H t =H max -H min .

[0199] S4.3: Calculate the wavelength of the ocean wave at the current time t.

[0200] Get H max and H min Then, according to the vector G t The index in can be used to query the corresponding H max and H minThe data of the waves on the X and Y axes. Figure 2 As shown, for the uncompressed data matrix, there are M t data, there are corresponding N along the X direction t The number of data, assuming H max The corresponding data on the X-axis and Y-axis are [X max_1 、X max_2 …、X max_Nt ]、[Y max_1 、Y max_2 …、Y max_Nt ], H min The corresponding data on the X-axis and Y-axis are [X min_1 、X min_2 …、X min_Nt ]、[Y min_1 、Y min_2 …、Y min_Nt ]. Thus, the wavelength of the wave at the current moment t can be calculated. The formula is:

[0201]

[0202] Calculate the wave point cloud data[X w ,Y w ,Z w ] is regarded as the wavelength of the waves at the current time t; the formula is:

[0203]

[0204] S4.4: Calculate the period of the ocean wave at the current time t.

[0205] There are two situations for discussion:

[0206] Case 1: For deep water waves, that is, the water depth is greater than or equal to half the wavelength.

[0207] The calculation formula for the wave period is:

[0208]

[0209] Where c is the wave velocity; g is the acceleration due to gravity (9.81 m / s 2 );L t is the wavelength of the ocean waves at the current time t; T t is the wave period at the current time t.

[0210] Case 2: For shallow water waves, that is, the depth of the water is less than half of the wavelength.

[0211] The calculation formula for the wave period is:

[0212]

[0213] Where c is the wave velocity; g is the acceleration due to gravity (9.81 m / s 2 );L t is the wavelength of the waves at the current time t; H is the water depth.

[0214] S5: Calculate the significant wave height, effective wavelength and effective period; calculate the wave spectrum.

[0215] In the above S1-S4, the wave height H in the point cloud data at the current time t is calculated. t 、wave wavelength L t and the wave period T t If only the above data H is used t 、L t and T t There will be large errors in expressing the effective wave height, effective wavelength and effective period; therefore, it is necessary to use the traditional "one-third large wave method" to calculate the effective wave height, effective wavelength and effective period.

[0216] The specific steps include:

[0217] S5.1: Continue to cycle S1-S4 multiple times, with the time of each cycle advancing. The binocular camera continuously captures images, and the computer connected to the binocular camera continuously generates wave point cloud data at different times, and continuously obtains wave height, wave wavelength and wave period data of the wave point cloud data at different times through calculation.

[0218] For example: loop through the operations S1-S4 once to obtain the wave height H at the next moment t+1 t+1 、wave wavelength L t+1 and the wave period T t+1 ;

[0219] Repeat the operation of S1-S4 to obtain the wave height H at time t+2 t+2 、wave wavelength L t+2 and the wave period T t+2 ;

[0220] And so on...

[0221] Assume that the operations S1-S4 are cycled n-1 times; combined with the original wave data at time t, there are n sets of data in total; each set of data includes wave height, wave wavelength and wave period.

[0222] Sort the wave heights in n sets of data from large to small, and add the sorted data into set A in the above order.

[0223] The set A is represented as:

[0224] A={{H (1) 、L (1) , T (1)},{H (2) 、L (2) , T (2)},…,{H (n) 、L (n) , T (n)}}, where H (1) ≥H (2) ≥…≥H (n)

[0225] Select the first third (round down) of the data in the set A; assuming the number of selected data is K, Where, Indicates rounding down. That is, select {{H (1) 、L (1) , T (1)},{H (2) 、L (2) , T (2)},…,{H (K) 、L (K) , T (K)}}

[0226] S5.2: Calculate the significant wave height, effective wavelength, and effective period of ocean waves.

[0227] Effective wave height H s The formula is:

[0228]

[0229] Effective wavelength L s The formula is:

[0230]

[0231] Effective period T s The formula is:

[0232]

[0233] Effective wave height H s , effective wavelength L s and effective period T s This is the wave data that needs to be observed, and the purpose of the plan has been completed.

[0234] S5.3: Calculate the ocean wave spectrum.

[0235] include:

[0236] S5.3.1: Assume that after loop S1-S4 u times, u frames of ocean wave point cloud data are generated. The generation time of u frames of ocean wave point cloud data is t b1 , t b2 ,……,t bu The u-frame ocean wave point cloud data in step S5.3 may be the same as or different from the n-frame ocean wave point cloud data in step S5.2.

[0237] S5.3.2: In the u-frame ocean wave point cloud data, select a fixed point whose position remains unchanged in each frame of ocean wave point cloud data. Obtain the elevation of the same position (fixed point) in each frame of ocean wave point cloud data and set it to: H x1 、H x2 ,……,H xu .

[0238] At this point, we can calculate t b1 to t bu The wave spectrum of the time series segment is:

[0239]

[0240] Where ω is the angular frequency of the wave, P(ω) is the power of the wave, and τ is the time offset.

[0241] The formula for τ is:

[0242]

[0243] The wave spectrum is also the wave data that needs to be observed, and this has also completed the purpose of the plan.

[0244] With the above-described preferred embodiments of the present invention as a guide, and with reference to the above description, relevant personnel are fully capable of making various changes and modifications without departing from the technical scope of this invention. The technical scope of this invention is not limited to the contents of the specification and must be determined according to the scope of the claims.

Claims

1. A method for observing ocean waves based on a binocular camera, characterized in that: The following steps are involved: S1: The binocular camera captures the wave image and generates the wave point cloud data at the current time t [X ct , Y ct , Z ct ]; Wave point cloud data[X ct , Y ct , Z ct ] is established in the left camera coordinate system; S2: Based on the IMU inertial measurement unit, the wave point cloud data [X ct , Y ct , Z ct ]Convert to the ocean wave point cloud data in the geodetic coordinate system[X w , Y w , Z w ]; S3: The ocean wave data in the geodetic coordinate system in S2 [X w , Y w , Z w ] for quality screening and optimization; Quality screening and optimization include: selecting regions of interest, quality screening, converting noisy data points to null values, null value filling, and Gaussian filtering; Quality screening will wave data [X w , Y w , Z w The quality of the data is divided into: excellent, good and poor; the poor quality wave data is eliminated [X w , Y w , Z w ], only retaining good quality and good wave data [X w , Y w , Z w ]Proceed to the next step; S4: Based on high-quality and good wave point cloud data [X w , Y w , Z w ], calculate the wave height, wave wavelength and wave period at that moment; S5: loop S1-S4, with each loop time advancing to obtain multiple frames of wave point cloud data at different times, and calculate the effective wave height, effective wave wavelength, effective wave period and wave spectrum based on the multiple frames of wave point cloud data at different times.

2. The method for observing ocean waves using a binocular camera according to claim 1, wherein: The selection of the region of interest and quality screening in step S3 include the following steps: S3.1: Select the region of interest. The formula is: X1=int(x1*Height) X2=int(x2*Height) Y1=int(y1*Width) Y2=int(y2·Width) Among them, x1, x2, y1, y2 are the row and column ratios of the area of ​​interest, Height and Width are the resolutions of the images captured by the camera, that is, the number of row and column pixels, X1, X2, Y1, Y2 are the index values ​​of the selected wave point cloud matrix data, that is, the wave point cloud data of interest is intercepted according to [X1:X2, Y1:Y2]; S3.2: Count the number of noise and missing values ​​in the area of ​​interest; set quality screening rules and perform screening, and judge the quality of the wave point cloud data based on the screening results [X w , Y w , Z w ] quality, remove the wave point cloud data with poor quality, and retain the wave point cloud data with good quality for the following steps; S3.2 specifically includes the following steps: S3.2.1: Obtain wave elevation data for each data point in the region of interest. Assume there are S data points in the region of interest, including S1 data points with elevation data and S2 data points without elevation data. S = S1 + S2. Data points without elevation data are considered null values. The wave height data of the i-th data point is recorded as Z i , i∈[1, S1], calculate the noise value of the i-th data point, the formula is: Where μ and σ are the mean and standard deviation of the wave height data of S1 "data points" respectively; A i is the noise value of the i-th data point; if |A i |>3, then the i-th data point is recorded as a "noise data point"; traverse S data points and obtain the number of "noise data points" in S1 "data points", recorded as M, M∈[1, S1]; While traversing S data points, the number of null values ​​S2 mentioned above is also obtained; S3.2.2: Calculate the ratio of the sum of the number of blank values ​​and the number of noise points to the total number of data points in the region of interest. This ratio is referred to as P. The formula is: P=(M+S2) / S S3.2.3: Set the quality screening rule Judge to judge the wave point cloud data [X w , Y w , Z w ] quality; The judgment formula is: When Judge is -1, it indicates that the wave point cloud data of this frame [X w , Y w , Z w ]The quality is poor, and the wave point cloud data of this frame should be discarded; When Judge is 1, it indicates that the wave point cloud data of this frame [X w , Y w , Z w ]Excellent quality; When Judge is 0, it indicates that the wave point cloud data of this frame [X w , Y w , Z w ]Good quality; retain the wave point cloud data of excellent quality and good quality for subsequent calculations.

3. The method for observing ocean waves using a binocular camera according to claim 2, wherein: Step S3 further includes the following steps: S3.3: If the wave point cloud data [X w , Y w , Z w ] is judged as excellent or good quality, and the wave point cloud data is optimized [X w , Y w , Z w ] in the noise data points and null values; S3.3 specifically includes the following steps: S3.3.1: Use row and column dilution factors to convert the wave point cloud data [X w , Y w , Z w ] dilute the diluted wave point cloud data [X w , Y w , Z w ] are converted to null values; Set the diluted [X w , Y w , Z w ] has a total of r1 null values ​​and r2 non-null values; r1 includes the null values ​​in the original wave point cloud data and the total number of null values ​​converted from "noise data points"; Let the coordinates of the bth null value be (x ib ,y jb ); S3.3.2: With null coordinates (x ib ,y jb ) is the center, and the side length R is the search range; if there is a valid value within the search range, the empty value needs to be filled; if there is no valid value within the search range, the empty value is not filled; the valid value is the data point with "elevation data" mentioned above; The filling process of step S3.3.2 includes the following steps: S3.3.2.1: Compute the Euclidean distance from the null value to each non-null value. The formula is: In the formula, (x f ,y f ) represents the coordinates of the fth non-null value; f∈[1, r2]; S3.3.2.2: Calculate the weight of the null value to each non-null value; The formula is: In the formula, p is the set weight index; S3.3.2.3: Fill in the gaps with the following formula: Where Z valid_values It is the selected valid non-empty value Z coordinate data. When calculating the X and Y axis data, just replace it with the corresponding axis data; S3.3.2.4: Gaussian filtering; Gaussian function: Construct a kernel matrix of size (2k+1)×(2k+1) and calculate the weight of each matrix position: Where i and j are the coordinates of each matrix data relative to the center of the kernel, i = -k, -k+1, ..., k-1, k, j = -k, -k+1, ..., k-1, k; In order to retain the information of subsequent elevation data, the Gaussian kernel is normalized to ensure that the sum of the filter kernels is 1: The original wave elevation matrix Z w Expand k pixels symmetrically around the periphery to generate the expansion matrix Z wk ; Original wave elevation matrix Z w That is, the wave elevation data after gap filling; For the wave elevation data matrix Z w Each original data point Z w (x, y), as the kernel center, calculate the smoothed value: That is, the Gaussian kernel is covered at the corresponding position of the extended matrix, multiplied point by point and summed up; finally, the extended boundary part is removed, and the smoothed result with the same size as the original data is retained for subsequent operations.

4. The method for observing ocean waves using a binocular camera according to claim 3, wherein: In step S4, the wave point cloud data [X w , Y w , Z w ] The wave height, wavelength and period of ocean waves include the following steps: S4.1: Transform the wave point cloud data [X w , Y w , Z w ] is considered as a matrix (M t , N t ,3); then the matrix (M t , N t ) is further compressed to obtain (M t , 1) vector G t ; Use the zero-crossing method to extract the vector G t The maximum and minimum elevation values ​​are denoted as H max and H min； S4.1: Calculate the wave height at the current time t; The formula is: H t =H max -H min ; Where H t is the wave height at the current moment t; S4.2: Calculate the wavelength of the ocean wave at the current time t; Get H max and H min Then, according to the vector G t The index in can be used to query the corresponding H max and H min The data of the waves on the X and Y axes; for the uncompressed data matrix, there are M t data, there are corresponding N along the X direction t The number of data, assuming H max The corresponding data on the X-axis and Y-axis are [X max_1 、X max_2 …、X max_Nt ]、[Y max_1 、Y max_2 …、Y max_Nt ], H min The corresponding data on the X-axis and Y-axis are [X min_1 、X min_2 …、X min_Nt ]、[Y min_1 、Y min_2 …、Y min_Nt ]; thus the wavelength of the wave at the current moment t can be calculated; the formula is: Calculate the wave point cloud data[X w , Y w , Z w ] is regarded as the wavelength of the waves at the current time t; the formula is: S4.3: Calculate the period of the ocean wave at the current time t; There are two situations: Case 1: For deep water waves with a water depth greater than or equal to half the wavelength; The calculation formula for the wave period is: Where c is the wave velocity; g is the acceleration due to gravity; L t is the wavelength of the ocean waves at the current time t; T t is the wave period at the current time t; Case 2: For shallow water waves where the water depth is less than half of the wavelength; The calculation formula for the wave period is: Where c is the wave velocity; g is the acceleration due to gravity; L t is the wavelength of the ocean wave at the current time t; h is the water depth; T t is the wave period at the current time t.

5. The method for observing ocean waves using a binocular camera according to claim 4, wherein: In step S5, calculating the significant wave height, effective wavelength, effective period and wave spectrum includes the following steps: S5.1: Continue to cycle through S1-S4 multiple times, with each cycle time increasing. The binocular camera continuously captures images, and the computer connected to the binocular camera continuously generates wave point cloud data at different times. The computer then continuously calculates and obtains wave height, wave wavelength, and wave period data for the wave point cloud data at different times. Assume that n sets of wave data are generated. Sort the wave heights in the n sets of data from large to small, and add the sorted data into set A in the above order. Set A is expressed as: A = {{H (1) , L (1) , T (1)}, {H (2) , L (2) , T (2)}, …, {H (n) , L (n) , T (n)}}, where H (1) ≥ H (2) ≥ … ≥ H (n) ; set up Select {{H (1) , L (1) , T (1) }, {H (2) , L (2) , T (2) },…,{H (K) , L (K) , T (K) }}; S5.2: Calculate the effective wave height, effective wavelength, and effective period of ocean waves; Effective wave height H s The formula is: Effective wavelength L s The formula is: Effective period T s The formula is:

6. The method for observing ocean waves using a binocular camera according to claim 5, wherein: In step S5, calculating the wave spectrum includes the following steps: S5.3.1: Set the loop S1-S4 multiple times to generate u-frame wave point cloud data; S5.3.2: In the u-frame ocean wave point cloud data, select a fixed point whose position remains unchanged in each frame of ocean wave point cloud data. Obtain the elevation of the same position in each frame of ocean wave point cloud data, denoted as H x1 、H x2 ,……,H xu ; The generation time of each frame of wave point cloud data is set to: t b1 , t b2 ,……,t bu ;but: t b1 to t bu The wave spectrum of the time series segment is: Where ω is the angular frequency of the wave. P(ω) is the power of the wave. τ is the time offset. The formula for τ is:

Images