River flow measurement method based on satellite-borne SAR along-track interference

Abstract

The invention discloses a river flow measurement method based on satellite-borne SAR along-track interference, and relates to the technical field of river flow measurement. Comprising the following steps: S1, selecting a straight and along-track vertical river reach of a to-be-measured river as a flow measurement area, and obtaining a satellite-borne SAR along-track interference image of the flow measurement area, the cross section terrain of the flow measurement area and the water level; s2, calculating an interference phase of a river measurement area according to the satellite-borne SAR along-track interference image in the S1, and performing inversion calculation by using the interference phase of the river measurement area to obtain a river surface flow velocity of a flow measurement area; and S3, determining flow velocities of different areas on the cross section of the flow measurement area according to a cross section vertical flow velocity distribution rule on the basis of the cross section terrain and water level of the flow measurement area obtained in S1 and the surface flow velocity of the flow measurement area obtained in S2, and calculating the cross section flow of the flow measurement area so as to obtain the flow of the to-be-measured river. According to the invention, large-range, long-period and all-weather river flow measurement can be realized.

Description

Technical Field

[0001] The present invention relates to the technical field of river flow measurement, and particularly to a river flow measurement method based on spaceborne SAR along-track interferometry. Background Art

[0002] Spaceborne SAR has been applied to ocean scenarios such as sea wave and wind field inversion, oil spill monitoring, ship detection, and sea ice monitoring. However, velocity information about the observed object cannot be obtained from conventional SAR images. Therefore, to measure the ocean surface flow field, the along-track interferometry (ATI) technique was introduced. The ATI technique was first proposed by Goldstein and Zebker in 1987, and the sea surface velocity in the radar line-of-sight direction was obtained using the ATI-SAR phase image. The main principle is that the interference phase generated by ATI-SAR is proportional to the radial velocity of the target in the corresponding observed scene. In 2001, an airborne ATI-SAR experiment was conducted near Sylt Island, Germany, and the water depth in the observed sea area was measured to be between 0 and 30 m, with a maximum flow velocity reaching 2 m / s. In 2003, the dual-beam interferometric SAR system (Umass-DBI) developed by the Massachusetts Institute of Technology achieved two-dimensional sea surface flow field measurement. The TerraSAR-X satellite launched by the German Aerospace Center in 2007 and the TanDEM-X satellite launched in 2010 have taken the measurement of the sea surface flow field by spaceborne high-resolution ATI-SAR to a higher level. Based on this, a series of studies have been carried out: In 2010, the first along-track interferometric measurement result of the Elbe River by TerraSAR-X was successfully obtained, and the result was consistent with the numerical theoretical model UnTRIM, with the average flow velocity error between 0.11 and 0.08 m / s. The along-track interferometry technique can achieve the measurement of the ocean surface flow field. However, no one has yet proposed using spaceborne SAR along-track interferometry to measure the river flow field on a river, and combining river hydrodynamics to calculate the river flow, and finally realizing a river flow measurement method based on spaceborne SAR along-track interferometry.

[0003] Therefore, proposing a river flow measurement method based on spaceborne SAR along-track interferometry to solve the difficulties existing in the prior art is an urgent problem that needs to be solved by those skilled in the art. Summary of the Invention

[0004] In view of this, the present invention provides a river flow measurement method based on spaceborne SAR along-track interferometry, which can achieve large-scale, long-term, and all-weather river flow measurement.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] A river flow measurement method based on spaceborne SAR along-track interferometry includes the following steps:

[0007] S1. Select a straight and cross - track - perpendicular river section of the river to be measured as the flow measurement area, and obtain the along - track interferometric SAR image of the flow measurement area, the cross - sectional topography of the flow measurement area, and the water level.

[0008] S2. According to the along - track interferometric SAR image in S1, calculate the interferometric phase of the river measurement area, and inversely calculate the surface velocity of the river in the flow measurement area using the interferometric phase of the river measurement area.

[0009] S3. Based on the cross - sectional topography and water level of the flow measurement area obtained in S1 and the surface velocity of the flow measurement area obtained in S2, determine the velocities of different regions on the cross - section of the flow measurement area according to the vertical velocity distribution law of the cross - section, and calculate the cross - sectional flow of the flow measurement area, thereby obtaining the flow of the river to be measured.

[0010] For the above - mentioned method, optionally, the specific content of S1 is as follows:

[0011] Select a straight and cross - track - perpendicular river section of the river to be measured as the flow measurement area, obtain the cross - sectional topography and water level of the flow measurement area, send a shooting task of the flow measurement area to the spaceborne SAR or select the SAR image of the flow measurement area from existing data, and obtain the along - track interferometric SAR image of the river measurement area, so that the direction of the surface velocity of the flow measurement area is consistent with the moving direction of the spaceborne SAR.

[0012] For the above - mentioned method, optionally, the specific content of S2 is as follows:

[0013] S21: Based on the main and secondary along - track interferometric SAR images of the flow measurement area, perform image registration on the two images of the flow measurement area to obtain the registered image of the flow measurement area, including pixel - level registration and sub - pixel - level registration.

[0014] S22: Based on the registered image of the flow measurement area in S21, perform multi - looking processing on the auxiliary image of the flow measurement area.

[0015] S23: Based on the registered image of the flow measurement area in S21 and the multi - looked - processed auxiliary image of the flow measurement area in S2, perform conjugate multiplication to obtain the interferometric phase diagram of the flow measurement area in complex form.

[0016] S24: Perform flat - earth phase removal, phase filtering, and phase unwrapping on the interferometric phase diagram of the flow measurement area obtained in S23, and inversely calculate the surface velocity of the river in the flow measurement area.

[0017] For the above - mentioned method, optionally, the specific content of S3 is as follows:

[0018] S31: Divide the river cross - section into N segments in terms of width Br.

[0019] S32: Based on the surface velocity of the river in the flow measurement area obtained in S2, calculate the average surface velocity V of each section i , on the i-th section of the river cross-section, where i = 1, 2... N;

[0020] S33: Calculate the average velocity in the vertical direction of the i-th section according to the vertical velocity distribution law of the cross-section and the water level The flow rate of the i-th section of the river cross-section where, Br i is the width of the i-th section, and the cross-section flow rate of the flow measurement area is the sum of the flow rates of N sections

[0021] Compared with the prior art, the present invention provides a method for measuring river flow based on spaceborne SAR along-track interferometry. The beneficial effects achieved by the present invention are as follows: Based on spaceborne SAR and using along-track interferometry technology, the present invention realizes the measurement of the surface flow field of the river. Combining the cross-section shape and water level of the river measurement area, according to the vertical distribution law of the channel velocity in hydrodynamics, the relationship between the surface flow field of the channel and the average velocity of the channel cross-section is deduced, and the cross-section flow rate of the river is calculated; it can realize large-scale, long-term, and all-weather river flow measurement, and is an important means to realize the river flow measurement of wide-area river basins. Brief Description of the Drawings

[0022] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only the embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained according to the provided drawings.

[0023] Figure 1 is the method for measuring river flow based on spaceborne SAR along-track interferometry disclosed by the present invention;

[0024] Figure 2 is the schematic diagram of measuring the surface velocity of the flow measurement area by spaceborne SAR along-track interferometry disclosed by the present invention;

[0025] Figure 3 is the schematic diagram of the division and velocity of the river cross-section of the present invention. Detailed Embodiment

[0026] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0027] In this application, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. The terms "comprising", "including" or any other variant thereof are intended to cover non-exclusive inclusion, such that a process, method, article or device comprising a series of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article or device. Without further limitation, an element defined by the statement "comprising an..." does not exclude the presence of additional identical elements in the process, method, article or device comprising the element.

[0028] Referring to Figure 1 As shown, the present invention discloses a method for measuring river flow based on along-track interferometry of spaceborne SAR, comprising the following steps:

[0029] S1. Select a straight and along-track perpendicular section of the river to be measured as the flow measurement area, and obtain the spaceborne SAR along-track interferogram of the flow measurement area, the cross-sectional topography of the flow measurement area and the water level;

[0030] S2. According to the spaceborne SAR along-track interferogram in S1, calculate the interference phase of the river measurement area, and inversely calculate the surface velocity of the river in the flow measurement area by using the interference phase of the river measurement area;

[0031] S3. Based on the cross-sectional topography and water level of the flow measurement area obtained in S1 and the surface velocity of the flow measurement area obtained in S2, determine the velocities of different regions on the cross-section of the flow measurement area according to the vertical velocity distribution law of the cross-section, and calculate the cross-sectional flow of the flow measurement area, so as to obtain the flow of the river to be measured.

[0032] Further, the specific content of S1 is:

[0033] Select a straight and along-track perpendicular section of the river to be measured as the flow measurement area, obtain the cross-sectional topography and water level of the flow measurement area, send a shooting task of the flow measurement area to the spaceborne SAR or select the SAR image of the flow measurement area from existing data, and obtain the spaceborne SAR along-track interferogram of the river measurement area, so that the direction of the surface velocity of the flow measurement area is consistent with the moving direction of the spaceborne SAR.

[0034] Further, referring to Figure 2 As shown, the specific content of S2 is:

[0035] S21: Based on the master and slave images of the flow measurement area in the along-track interferometry of spaceborne SAR, perform image registration on the two flow measurement area images to obtain the registered flow measurement area images, including pixel-level registration and sub-pixel-level registration;

[0036] S22: Based on the registered flow measurement area images in S21, perform multi-look processing on the auxiliary image of the flow measurement area;

[0037] S23: Based on the registered flow measurement area images in S21 and the multi-looked auxiliary image of the flow measurement area in S2, perform conjugate multiplication to obtain the interferometric phase image of the flow measurement area in complex form;

[0038] S24: Perform flat-earth phase removal, phase filtering, and phase unwrapping on the interferometric phase image of the flow measurement area obtained in S23, and perform inversion calculation to obtain the surface velocity of the river in the flow measurement area.

[0039] Further, as shown in Figure 3 The specific content of S3 is as follows:

[0040] S31: Divide the river cross-section into N segments in the width Br;

[0041] S32: Based on the surface velocity of the river in the flow measurement area obtained in S2, statistically calculate the average surface velocity V i of each segment. On the i-th segment of the river cross-section, where i = 1, 2... N;

[0042] S33: Calculate the average velocity in the vertical direction of the i-th segment according to the vertical velocity distribution law of the cross-section and the water level The flow rate of the i-th segment of the river cross-section where Br i is the width of the i-th segment, and the cross-section flow rate of the flow measurement area is the sum of the flow rates of the N segments

[0043] In a specific embodiment, the specific implementation steps are as follows:

[0044] Step 1: Select a straight and along-track perpendicular section of the river to be measured as the flow measurement area ( Figure 1 the gray area of the river in the figure), and obtain the spaceborne SAR along-track interferometry image, the river cross-section topography, and the water level of the measurement area.

[0045] Step 2: As shown in Figure 2 The measurement of the surface velocity of the river to be measured by spaceborne SAR along-track interferometry includes:

[0046] Step 2A: Perform registration on the along-track interferometry image of the measurement area.

[0047] In an embodiment, the obtained images will be pixel-level registered using the cross-correlation function algorithm, a classic image registration method. In the main image, first select the first control point, determine an MxM matching window centered on the control point, and determine an NxN search window centered on the corresponding control point in the auxiliary image. Determine an M×M window within the search window, calculate the correlation coefficient between this window and the matching window, and use symbols s1 and s2 to represent this window and the matching window respectively. Take the modulus of s1 and s2 to obtain their corresponding amplitude images. The process of moving the maximum cross-correlation function of the window by Fourier transform is transformed into the process of locating the maximum value of the cross-correlation function of the two image windows. The window at the maximum is the window when it reaches rough registration with the matching window. Record the center positions of the matching window and the corresponding window at this time to find a pair of corresponding control points. The calculation method of the Fourier transform of the cross-correlation function of the two image windows is as follows:

[0048] ρ r =|FFTSHIFT(FFT2 -1 (FFT2(y1)..*FFT2(y2) * ))|

[0049] In the formula, FFT2(·) and FFT2 -1 (·) represent the two-dimensional forward Fourier transform and the inverse transform respectively,.* represents matrix dot multiplication, FFTSHIFT(·) represents shifting the zero-frequency component of the image spectrum to the middle of the spectrum, and y1 and y2 are the data of the two image windows.

[0050] The registration offset of the window center point can be obtained by the position where ρ r obtains the maximum value:

[0051]

[0052]

[0053] where Δx and Δy are the registration offsets in the azimuth direction and the range direction respectively, M and N are the azimuth and range dimensions of the data window respectively, and max_x and max_y are the azimuth and range position coordinates when ρ r obtains the maximum value.

[0054] The window at this time is the window when it reaches rough registration with the matching window. Record the center positions of the matching window and the corresponding window at this time to find a pair of corresponding control points.

[0055] When the maximum value of the window matching correlation coefficient is less than a certain threshold, abandon this control point, and finally calculate the rough registration control points and output the coordinate positions in the main and auxiliary images one by one.

[0056] The above work completes pixel-level registration, using the entire pixel as the offset, so its accuracy can only reach the pixel level. The next step is to start sub-pixel registration.

[0057] Sub-pixel image registration generally also uses window-based automatic registration technology. In this embodiment, the interpolation resampling method is used to improve the registration accuracy, and the Lancos interpolation method is adopted.

[0058] In interpolation resampling, generally, the interpolation method is first used to perform oversampling processing on the master and slave images. The interpolation is generally performed at 0.1×0.1 pixels. For example, for a 1×1 image, after interpolation, it becomes a 10×10 image. Then, similar to pixel-level registration, control points are reselected, the matching window and search window are set, the matching index is calculated for reliability detection, unreliable control points are removed, and the offset of the control points is determined.

[0059] Lanczos interpolation realizes image enlargement or reduction through the following steps:

[0060] (1) For each pixel point to be interpolated, find several original pixel points closest to the point to be interpolated in its surrounding neighborhood;

[0061] (2) Use the Lanczos window function to weight the pixel values in the neighborhood of the point to be interpolated to obtain a set of weighting coefficients;

[0062] (3) Use this set of weighting coefficients and the corresponding neighborhood pixel values to calculate the pixel value of the pixel point to be interpolated through the Lanczos interpolation formula.

[0063] The one-dimensional Lanczos interpolation formula is as follows:

[0064]

[0065] where s i is the pixel value of the i-th pixel in the original image, and x is the interpolation position.

[0066] The above method can perform rough registration and fine registration on the master and slave SAR amplitude images in sequence, and then calculate the coherence coefficient and perform interpolation resampling on the slave image based on the registration result. And calculate the azimuth and range position coordinates when the offset and the master image obtain the maximum value to construct a quadratic registration polynomial, check the polynomial fitting accuracy, and use the polynomial for sub-pixel registration of the image.

[0067] Step 2B: Perform multi-look processing on the registered image for image denoising. In this embodiment, the multi-look processing method is to perform phase neighborhood averaging on the phase diagram.

[0068] Step 2C: Calculate the interferometric phase diagram for the image produced in Step 2B. Multiply the SAR image pairs conjugated after Step 3 processing to obtain an interferogram in complex form. The phase of the interferogram is the wrapped form of the actual interferometric phase, which differs from the actual phase by an integer multiple of 2π.

[0069] The interferometric phase diagram can be calculated by multiplying the conjugate complex numbers of the main image s1 and the auxiliary image s2'.

[0070]

[0071] In the formula, arg(·) represents finding the phase value of a complex number.

[0072] Step 2D: Evaluate the interferometric phase diagram generated in Step 2C. If the result of the interferometric phase diagram is not satisfactory, phase filtering and phase unwrapping can be performed. In this embodiment, mean filtering is used for phase filtering: select an image window of size N×N, and use the average value of all pixels within the window as the target value. The window slides along with the image until it covers the entire image, thereby reducing image noise; minimum norm method is used for phase unwrapping, and this operation can also be skipped when the phase fringes are relatively sparse. As shown in the figure, the phase filtering adopts the mean filtering method: select an N×N-sized image window, and use the average value of all pixels within the window as the target value. The window slides along with the image until it covers the entire image, so as to reduce the image noise; the phase unwrapping adopts the minimum norm method, and this operation can also be not performed when the phase fringes are relatively sparse.

[0073] Step 2E: Perform surface flow velocity inversion calculation on the interferometric phase diagram generated in Step 2C. Obtain the flow velocity u r (y,t) of each pixel point in the measurement area, and the calculation method is as follows:

[0074]

[0075] In the above formula, u r is the surface flow velocity, the interferometric phase, λ is the radar wavelength, V is the platform moving speed, B is the effective baseline length, θ is the incident angle, as shown in Figure 2 the figure.

[0076] Step 3: Based on the cross-sectional topography, water level of the flow measurement area and the surface flow velocity of the flow measurement area, determine the flow velocities in different areas of the cross-section of the river measurement area according to the vertical flow velocity distribution law of the cross-section, and calculate the cross-sectional flow of the flow measurement area, so as to obtain the river flow, as shown in Figure 3 the figure, specifically including:

[0077] Step 3A: After obtaining the cross-sectional topography and water level of the measurement area, divide the cross-section of the measurement area into N segments in the river width (Br) direction (the number of divided segments N and the width Br of each segment i can be determined according to the degree of topographic change. Generally, Br i is an integer multiple n of the spatial resolution of the interferometric phase diagram), as shown in Figure 3 the figure;

[0078] Step 3B: Correlate the surface velocity of the river in the measured area obtained in Step 2 with each river reach i in Figure 3 and calculate the average surface velocity of each river reach i where is the velocity of the kth grid in river reach i

[0079] Step 3C: Determine the vertical velocity distribution law of the cross-section according to the water level and the cross-sectional shape of the river channel. In this embodiment, the exponential velocity distribution is used to describe the vertical distribution of the cross-sectional velocity, as shown in the following formula:

[0080]

[0081] In the above formula, u is the velocity at a distance y from the river bottom, y is the distance from the river bottom, h is the water depth, and V avr is the average velocity of the vertical line; m is the exponent, which takes 1 / 6 - 1 / 7 in clear water flow;

[0082] Step 3D: According to the vertical velocity distribution law and the surface velocity V i of the ith river reach, use integration to calculate the vertical average velocity of each river reach i as shown in the following formula:

[0083]

[0084] Step 3E: Calculate the flow rate Q i of each river reach i, as shown in the following formula:

[0085]

[0086] In the above formula, Br i is the width of the ith section, is the vertical average velocity of the ith river reach.

[0087] Step 3F: Aggregate the flow rates of all river reaches in the cross-section to obtain the river flow rate in the flow measurement area, as shown in the following formula:

[0088]

[0089] In the above formula: Q i is the flow rate of the ith river reach, and N is the number of sections divided in the river width (Br) direction of the cross-section.

[0090] Thus, the flow rate Q is the measurement result of the river flow rate measurement method based on spaceborne SAR along-track interferometry of the present invention.

[0091] In this embodiment, the registration method adopted in step 2A is only one of the available registration methods, the multi-view processing method adopted in step 2B is only one of the available multi-view processing methods, the phase filtering method adopted in step 2D is only one of the available methods, the method for retrieving flow velocity from interference phase adopted in step 2E is only one of the available methods, the calculation and statistical method for the average surface velocity of the river reach in step 3B is only one of them, and the vertical flow velocity distribution law and the calculation of the vertical average velocity in step 3D are only one of them. These are some embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the scope of protection of the present invention.

[0092] The embodiments in this specification are described in a progressive manner. For a system or a system embodiment, since it is basically similar to the method embodiment, the description is relatively simple. For the relevant parts, reference can be made to the partial description of the method embodiment. The systems and system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed to multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment. Those of ordinary skill in the art can understand and implement it without creative efforts.

[0093] To clearly illustrate the interchangeability of hardware and software, the composition and steps of each example have been generally described according to functions in the above description. Whether these functions are executed in a hardware or software manner depends on the specific application and design constraints of the technical solution. Professional technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered to exceed the scope of the present invention.

[0094] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but will be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for measuring river flow based on along-track interferometry of spaceborne SAR, characterized in that, It includes the following steps: S1. Select a straight and along-track perpendicular river section of the river to be measured as the flow measurement area, and obtain the along-track interferometric SAR image of the flow measurement area, the cross-sectional topography and water level of the flow measurement area; S2. According to the along-track interferometric SAR image in S1, calculate the interferometric phase of the river measurement area, and inversely calculate the surface velocity of the river in the flow measurement area by using the interferometric phase of the river measurement area; S3. Based on the cross-sectional topography and water level of the flow measurement area obtained in S1 and the surface velocity of the flow measurement area obtained in S2, determine the velocities of different regions on the cross-section of the flow measurement area according to the vertical velocity distribution law of the cross-section, and calculate the cross-sectional flow of the flow measurement area, so as to obtain the flow of the river to be measured.

2. The river flow measurement method based on along-track interferometry of spaceborne SAR according to claim 1, characterized in that The specific content of S1 is: Select a straight and along-track perpendicular river section of the river to be measured as the flow measurement area, obtain the cross-sectional topography and water level of the flow measurement area, send a shooting task of the flow measurement area to the spaceborne SAR or select the SAR image of the flow measurement area from existing data, and obtain the along-track interferometric SAR image of the river measurement area, so that the direction of the surface velocity of the flow measurement area is consistent with the moving direction of the spaceborne SAR.

3. The river flow measurement method based on along-track interferometry of spaceborne SAR according to claim 1, characterized in that The specific content of S2 is: S21: Based on the main and secondary along-track interferometric SAR images of the flow measurement area, perform image registration on the two flow measurement area images to obtain the registered flow measurement area images, including pixel-level registration and sub-pixel-level registration; S22: Based on the registered flow measurement area images in S21, perform multi-look processing on the auxiliary flow measurement area image; S23: Based on the registered flow measurement area images in S21 and the multi-look processed auxiliary flow measurement area image in S2, perform conjugate multiplication to obtain the interferometric phase diagram of the flow measurement area in complex form; S24: Perform flat-earth phase removal, phase filtering and phase unwrapping processing on the interferometric phase diagram of the flow measurement area obtained in S23, and inversely calculate the surface velocity of the river in the flow measurement area.

4. The river flow measurement method based on along-track interferometry of spaceborne SAR according to claim 1, characterized in that The specific content of S3 is: S31: Divide the river cross-section into N segments in width Br; S32: Based on the surface velocity of the river in the flow measurement area obtained in S2, calculate the average surface velocity V of each section i , on the i-th section of the river cross-section, where i = 1, 2... N; S33: Calculate the average velocity in the vertical direction of the $i$-th section based on the vertical velocity distribution law of the cross-section and the water level Discharge of the $i$-th section of the river cross-section where $Br$ i is the width of the $i$-th section, and the cross-section discharge of the flow measurement area is the sum of the discharges of $N$ sections

Images