Detailed Description
The invention provides a sea surface fluctuation target SAR image refocusing-based sea surface parameter acquisition method, which realizes accurate acquisition of sea surface parameters by adopting an algorithm based on inverse transformation and ISAR motion compensation and refocusing.
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings in conjunction with specific embodiments.
In one embodiment of the invention, a sea wave parameter acquisition method based on sea surface fluctuation target SAR image refocusing is provided. Fig. 1 is a technical flowchart of a sea wave parameter acquisition method based on sea surface fluctuation target SAR image refocusing. Referring to fig. 1, the sea wave parameter acquisition method based on sea surface fluctuation target SAR image refocusing comprises the following steps:
step A: establishing a model of wave fluctuation motion, and quantitatively analyzing an inclined distance error caused by the model, thereby obtaining a phase error expression caused by the inclined distance error;
the heave motion expression:
slope distance error:
phase error:
whereinA vector representing the heave motion of the target point, Δ R (η) representing the slope error,representing the modelled derived phase error, AvRepresenting the amplitude of the heave motion, λ representing the pulse wavelength, γ representing the sum of the downward viewing angle and the earth's center angle, ω representing the period of the heave motion, η representing any time instant corresponding to the heave motion,indicates the initial phase angle.
A small buoy is placed in an observation sea area to serve as a target point for the fluctuation motion of sea waves. Fig. 2 is a schematic view of observation target points-satellite-earth geography, as shown in fig. 2, T represents a stationary target point, S represents a satellite position of η at any azimuth, and S' represents a satellite subsatellite point at the time. Assuming that T has undulations in the form of sinusoidal vibrations along the radial direction of the earth, P represents the instantaneous position of the target point undulation motion.
Based on the above, the derivation process of the 3 formulas in step a is as follows:
vector representing the undulation motion of the target point:
the instantaneous slope distance vector between the static target and the satellite at the moment isTherefore, the instantaneous slope distance of the undulating target at the moment can be obtained as follows:
by applying the geographical relation and the vector knowledge, the approximate instantaneous slope distance of the fluctuating target can be obtained as follows:
then it can be deduced that the instantaneous slope error caused by the rolling object motion is:
the phase error caused by the instantaneous slope error is:
as can be seen from the phase error expression of equation (3): the phase error caused by the heave motion is in a simple harmonic vibration mode, and the phase error caused by the heave motion and the heave motion have the same period as seen by referring to a target point heave motion vector expression, and the amplitude of the phase error is determined by the heave motion amplitude, the downward viewing angle, the geocentric angle between the target and the satellite and the pulse wavelength.
After the sea wave motion target model is established in the step A to derive the phase error expression, the phase error expression can be directly called in the subsequent step D execution process without repeated derivation. Of course, step a may be omitted if the model is known and the phase error expression is also known.
Particularly, when analyzing the wave information of different areas, the wave fluctuation motion model established in the step a is established once, and modeling is not required to be carried out again each time, and the deduced error expression can be directly called when analyzing the wave information of different areas subsequently.
And B: acquiring an SAR image of a target area, and applying a refocusing algorithm based on inverse transformation and ISAR motion compensation to data containing a moving target sub-image in the SAR image to obtain a focused SAR image and a phase error compensation curve corresponding to the image;
fig. 3 is a flowchart of a refocusing processing algorithm based on Inverse transform and ISAR motion compensation, and as shown in fig. 3, the step B of applying an Inverse transform and ISAR (Inverse SAR/Inverse synthetic aperture radar) based motion compensation refocusing algorithm to the target data of the heave motion further includes:
substep B1: extracting a subimage containing a fluctuating moving target from an SAR image obtained by acquiring information in a target area;
substep B2: performing azimuth Fourier transform on the sub-image containing the fluctuating moving target to obtain equivalent original data before azimuth compression;
the data can be regarded as azimuth unfocused data after completely compensating the motion between the stationary target and the satellite, and can be used as raw data of input ISAR processing focusing, and the distance frequency domain can be equivalent to the following formula:
wherein f is0Representing the carrier frequency, B representing the transmission pulse bandwidth, TsRepresenting the imaging time, R (η) represents the instantaneous slope distance between the target and the satellite due to the heave motion of the target:
R(η)=R0+ΔR(η) (7)
substep B3: performing refocusing processing on the equivalent original data obtained in the step B2 by adopting a distance alignment algorithm to correct a distance migration error caused by the fluctuation motion to obtain a distance migration error correction curve;
fig. 4 is a processing flow chart of a distance alignment algorithm in the ISAR algorithm, and as shown in fig. 4, the distance alignment algorithm adopts a global minimum entropy distance alignment algorithm, performs 32-fold interpolation on distance direction data to obtain distance alignment accuracy better than 32-fold distance resolution, and can completely correct distance migration errors caused by rolling motion.
The minimum entropy distance alignment algorithm applied to the orientation uncompressed data in fig. 4 achieves the distance-to-envelope alignment of the target, and eliminates the dislocation of adjacent echoes generated by motion in the distance direction without providing a priori knowledge. Establishing a model of the optimization problem, wherein the objective function is entropy values of all one-dimensional range image sums, and can be expressed as:
pave(r) represents the sum of all target one-dimensional range images:
where r is the range-oriented sample point and Δ r (n) represents the range offset of the nth echo.
paveThe degree of sharpening of (r) may be used to measure the degree of alignment of the echo envelopes, and when the degree of sharpening reaches the highest, the echoes are substantially aligned. And the Shannon entropy can be used to measure pave(r) degree of sharpening, the higher the waveform degree of sharpening, the smaller the entropy value, and therefore the optimization criterion is paveThe entropy value of (r) is minimum, and a global optimal solution of delta r (n) can be obtained through iteration, and the specific algorithm steps are described as follows:
step S401: initializing an offset Δ r (n);
step S402: calculate the current pave(r);
Step S403: calculating pave(r), if the entropy value is not reduced, obtaining a final estimated value of the offset delta r (n), and exiting the loop algorithm, otherwise, executing the step S404;
step S404: computingIs calculated according to the following formula:
step S405: using FFT to convertCorrelating with each pulse;
step S406: taking the delta r (n) when the absolute value of the cross-correlation function reaches the maximum as the distance offset required by the iteration to obtain new pave(r), return to step 402.
Substep B4: further performing phase correction on the data of the range migration error correction curve obtained in the step B3, and correcting a phase error caused by the fluctuating motion by adopting a phase compensation algorithm to obtain a phase error compensation curve;
fig. 5 is a flowchart of distance phase correction processing in the ISAR algorithm, and as shown in fig. 5, the phase compensation algorithm adopts a minimum entropy phase correction algorithm based on a variable step gradient descent method to obtain optimal operation efficiency and result. Establishing an optimization problem model for solving the compensation phase to ensure that the image focusing effect is the best, wherein an objective function is an entropy value of an image and is expressed as follows:
where n is the azimuth sample point, g (r, n) is the phase compensated image, which can be expressed as:
wherein f (r, m) represents data of the azimuth frequency domain after the distance alignment,indicating the amount of phase compensation for the mth echo.
The degree of focus of an image can be measured by the entropy of the image, and the higher the degree of focus, the smaller the entropy of the image and the sharper the image. Therefore, the optimization criterion is that the value of entropy reaches the minimum, under the condition of giving an initial value of a compensation phase, a numerical value iterative algorithm is utilized to gradually approximate to finally solve a target function to obtainIn addition, the experiment searches for the global optimal solution by a gradient descent method with variable step lengthThe value of (2) further improves the iterative computation speed, and simultaneously ensures the final image quality.
The minimum entropy phase correction algorithm based on the variable step gradient descent method comprises the following specific implementation steps:
step S501: initializing a compensation phase to be zero;
step S502: calculating a gradient vector grad of the entropy function at the current position;
step S503: determining magnitude of negative direction change along gradientNamely, a Step length parameter Step needs to be determined;
the specific principle of the method is that the step length is halved each time to make the algorithm converge towards the minimum entropy quickly, and when the convergence is close to the minimum entropy, the optimization result is obtained by continuously fine-tuning the step length.
Step S504: will obtainAndand is added up, as a new acquisition value,and compensating the original image by using the new phase acquisition value to obtain a new image, and circulating to the Step 502 until Step is equal to 0.
Substep B5: and B4, performing azimuth focusing processing on the data of the phase error compensation curve obtained in the step B4 to finally obtain a focused SAR image.
And the azimuth focusing is realized by adopting Fourier transform or inverse Fourier transform.
In the algorithm execution process, the final output result is an SAR image with good focusing effect and a phase error compensation curve, and the SAR image and the phase error compensation curve have a one-to-one correspondence relationship.
And C: fitting a sine function to the phase error compensation curve obtained in the step B4 in combination with the first-order fitting term in the range migration correction curve obtained in the step B3 to obtain a phase error compensation fitting curve in the form of a sine function;
step D: c, correspondingly analyzing the phase error compensation fitting curve in the form of the sine function obtained in the step C in combination with the phase error expression in the sea wave fluctuating motion model established in the step A to obtain the motion amplitude and the motion period of the fluctuating motion target;
the phase error compensation curve represents the phase of each pulse corresponding to the compensationAnd the pulse number represents the azimuth sampling point, i.e. equivalent to the azimuth time:
the phase correction of the defocused target is performed to obtain an SAR image with a good focusing effect, and it can be considered that the phase error is completely compensated, that is, the compensation phase in the obtained phase error compensation curve is the error phase obtained by the analysis in the first step:
the period T of the function can be obtained by sinusoidal fitting of the compensation phase0And amplitude A0Then, the period T and the amplitude a of the rolling motion can be obtained by equation (14) as:
step E: and deducing parameters such as wave height, wave period and the like of the sea waves according to the motion amplitude and the motion period of the fluctuating motion target.
Specifically, the derivation process of obtaining the sea wave parameters according to the motion parameters of the undulating motion target is as follows:
assuming a target mass m, buoyancy and gravity are balanced at rest
Wave height of sea wave is h0Then the buoyancy force applied to the object is increased to:
Fs=ρgS(h+h0) (17)
the object is acted by upward buoyancy to cause the fluctuation, the movement from the rest position to the highest point is a quarter of a cycle, and the height which can be reached is the amplitude of the fluctuation of the object.
Where A is the amplitude of the heave motion of the object and a is the acceleration.
It should be noted that the above formula (19) is a rough approximation, and the increase of the buoyancy of the object caused by the actual sea wave should be smaller than the above calculation, so a should also be smaller than the result in the above formula, and the ratio is about one fourth to one half, that is:
it can be seen that A is inversely proportional to h, and h reflects the draught of the object and is directly proportional to the weight of the object itself, so that the fluctuation range of the object is inversely proportional to the weight of the object itself, which is consistent with the common sense.
It is to be noted that, in the attached drawings or in the description, the implementation modes not shown or described are all the modes known by the ordinary skilled person in the field of technology, and are not described in detail. Furthermore, the above definition of the model for establishing the wave motion is not limited to the various specific expressions and specific observation models mentioned in the embodiments, and those skilled in the art can easily modify or replace them, for example:
(1) the SAR image can be acquired by using other forms except a satellite, such as a communication form carrying radar transmitting and receiving devices, such as an airplane, a hot air balloon, a microwave shore-based radar and the like, and then a corresponding heave motion model is established;
(2) the wave motion model can be replaced by other functions describing the wave motion besides simple harmonic motion.
By using the technical scheme of the embodiment, the wave parameters of the coast of a certain island in the sea area near the Chinese voyage are obtained. FIG. 6 is an SAR image of the offshore area of a certain island in China voyage under a Terra-SAR sliding bunching high-resolution mode; fig. 7 shows an optical image corresponding to fig. 6 for google earth. The data used for processing comes from a certain island coast of the sea area around the China voyage in the Terra-SAR satellite sliding bunching mode. Table 1 shows the parameters associated with the corresponding image of fig. 6. Combining fig. 6 and 7, it can be seen that the line of mustache defocusing of the island edge is tied to the floating barrel on the bank, because of the up-and-down fluctuating movement of the sea wave, a certain bright line is taken for experiment, a comparison graph before and after refocusing treatment is obtained, as shown in fig. 8, then a phase curve for compensating the line is obtained, as shown in fig. 9, by performing sine fitting on fig. 9 and combining a distance migration correction curve, considering that a certain error still exists, the result of obtaining parameters of the fluctuating movement is obtained as follows:
Ae=0.14~0.25m
Te=2.8~3.5s
wherein A iseRepresenting the magnitude of the acquisition, TeThe obtained period is represented, and considering that the floating bucket has weight, the wave height of the fluctuation motion of the floating bucket is about ten times of the obtained amplitude, and the period is equivalent to the period of the fluctuation motion, so the obtained sea wave parameter result is as follows:
As=1.4~2.5m
Ts=2.8~3.5s
wherein A issIndicating the wave height, TsRepresenting the wave period.
According to the data, the Chinese vogue belongs to the Bohai sea area, the average wave height of the sea wave is within the range of 0.5-2 m, the period is within 3-5 s, and the obtained wave parameters are within the range, so that the correctness of the method is demonstrated, and meanwhile, the accuracy is certain.
TABLE 1 Terra-SAR Tanshutan data corresponding parameters
Radar parameter |
Numerical value |
Radar parameter |
Numerical value |
λ |
0.03112m |
PRF |
42300 |
Bandwidth of |
300MHz |
Downward viewing angle |
38.799° |
Sampling rate |
329.658MHz |
Earth's center angle |
3.13° |
Short distance |
613.981km |
Distance resolution |
0.5m |
Synthetic pore size time |
5.39s |
Azimuthal resolution |
0.25m |
Thus, the embodiment of the present invention is described.
In summary, in the embodiment, an algorithm based on inverse transformation and ISAR motion compensation and refocusing is used for accurately obtaining the phase error of the SAR image data containing the moving target, and the obtained phase error is compared with an error expression in the established sea wave motion model, so that the obtained results of parameters such as the wave height and the period of the sea wave can be obtained.
Of course, according to actual needs, the sea wave parameter acquisition method based on sea surface fluctuation target SAR image refocusing provided by the invention further comprises other common algorithms and steps, and is not repeated herein because the method is irrelevant to the innovation of the invention.
The above-mentioned embodiments are intended to illustrate the objects, technical solutions and advantages of the present invention in further detail, and it should be understood that the above-mentioned embodiments are only illustrative of the present invention and are not intended to limit the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the protection scope of the present invention.