Disclosure of Invention
The invention solves the problems: the method overcomes the defects of the prior art, provides a horizontal attitude determination method based on underwater Snell window edge identification, and can realize the full-autonomous determination of the horizontal attitude angle without accumulated errors in an underwater environment.
The invention provides the pitch angle and the roll angle for the carrier independently by acquiring the external optical information. The Snell window is a high-brightness circular area formed at the water-air interface due to refraction when underwater air is imaged. In the imaging system, the observation vector of the circle center is vertical to the water surface upwards, namely the sky direction under the navigation coordinate system. Because the light inside the Snell window is formed by refracting sky light and has higher brightness relative to the outside, the Snell window in the null imaging has more obvious edges under the condition of a flat still water surface, and the window edges are conveniently extracted by using an image processing method so as to obtain the circle center of the Snell window. By integrating the two points, the invention provides a method for calculating the horizontal attitude of the carrier by extracting the circle center of a Snell window in the empty imaging.
The circle detection in the image is usually performed by using the hough circle transform, however, due to the distortion of the lens in the imaging system, when the underwater carrier is placed obliquely, the snell window in the image is not in a regular circle shape, so the conventional hough circle transform is not applicable any more. Aiming at the problem, the method is improved on the basis of Hough circle transformation, and a camera lens model is added, so that a Snell window and the circle center under the distortion of a camera lens can be accurately extracted.
The technical scheme adopted by the invention for solving the technical problems is as follows: a horizontal attitude determination method based on underwater Snell window edge identification comprises the following implementation steps:
step (1), carrying out noise filtering and binarization processing on the acquired image, obtaining an edge image according to light intensity information in the image, and extracting coordinates [ m, n ] of a non-zero pixel p in the edge image in a pixel coordinate system]. Wherein any point is represented by non-zero pixelsIs piWith the coordinate of [ mi,ni]Wherein i is 1,2, …, s, s is the number of nonzero pixels;
step (2), calculating a value set of the radius r of the included angle, wherein r is a possible value of the radius of the view angle of a Snell window, establishing a three-dimensional Hough space H (m, n, r), and setting all elements in the space to zero;
step (3) utilizing the edge image subjected to the binarization processing in the step (1) to convert each non-zero pixel p into
iThe observation vector under the world coordinate system, namely the w system, corresponding to the camera lens model theta is inverted through the camera lens model theta, and the observation vector is used as the z axis to establish the pixel local coordinate system, namely the l system
iSolving the non-zero pixel p
iCorresponding to l
iCoordinate transformation matrix between system and w system
Step (4) with a non-zero pixel p
iDetermining a cone set by using the included angle radius r, the redundancy radius delta alpha and the radius stepping delta calculated in the step (2) and taking all generatrix vectors g of the cone set as an axis
iReverting to the pixel coordinate system, at each generatrix vector g
iqVoting in a three-dimensional Hough space H (m, n, r) of corresponding pixel coordinates, wherein
To obtain non-zero pixel point p
iA three-dimensional Hough space as a circle center;
and (5) performing the operations of the step (3) and the step (4) on all non-zero pixels, wherein the pixel point with the highest ticket number in the three-dimensional Hough space H (m, n, r) is the circle center of a Snell window, calculating an observation vector of the pixel coordinate of the circle center under a world coordinate system, and inverting the horizontal posture.
Further, the radius r of the included angle calculated in the step (2) is a value set, where r is a possible value of the radius of the view angle of the snell window, and a three-dimensional hough space H (m, n, r) is established, and all elements in the space are set to zero, which is specifically realized as follows:
according to Snell refractionDetermining an angle between the edge and the center of a window by a law and a lens model, taking the angle as an included angle radius in Hough circle transformation, and determining a redundant range of the included angle radius and radius stepping; the refractive indexes of atmosphere and water are respectively naAnd nwThen the angle α between the edge and the center of the snell window satisfies:
considering the errors of refractive indexes of a camera and a medium in practical application, the imaging size of a Snell window can not strictly accord with the theoretical model of the formula, so that the value range of the included angle radius r is [ alpha-delta alpha, alpha + delta alpha ] when the redundant radius delta alpha is set according to the practical water body environment and the lens model. Then, setting the included angle radius step as delta, and setting the value set of the included angle radius r as:
r∈{α-Δα,α-Δα+δ,α-Δα+2δ,…,α-Δα+kδ}(α-Δα+kδ≤α+Δα,k∈N+)
and establishing a three-dimensional Hough space H (m, n, r), wherein (m, n) is pixel coordinates, and initial values of the three-dimensional Hough space are all set to be 0.
Further, the step (3) uses the edge image binarized in the step (1) to convert each non-zero pixel p into
iThe observation vector under the world coordinate system, namely the w system, corresponding to the camera lens model theta is inverted through the camera lens model theta, and the observation vector is used as the z axis to establish the pixel local coordinate system, namely the l system
iSolving the non-zero pixel p
iCorresponding to l
iCoordinate transformation matrix between system and w system
The concrete implementation is as follows:
in the pixel coordinate system, an arbitrary point is a nonzero pixel piHas the coordinate of [ mi,ni]Through camera lens calibration, a camera lens model theta (p) can be calibrated, and a pixel point piCorresponding incident ray propagation vector kappa under system wiExpressed as:
κi=Θ(pi) (16)
κ
iis a unit vector, an observation vector v in the direction
iIs represented by in the w system
Order to
Then there are:
[ai bi ci]T=-Θ(pi) (17)
with v
iEstablishing a local coordinate system of picture elements, i.e./for the z-axis
iIs, definition w is and
iin the process of conversion between systems, the rotation angles around the x-axis and the y-axis are respectively omega
xi,
The rotation angle around the z-axis is 0, yielding l
iThe coordinate transformation matrix from system to w system is:
due to the observation vector v
iIn l
iIs represented as
Thus, the following steps are obtained:
combined formula (17), omegaxiAnd omegayiThe trigonometric function form of (a) can be expressed as follows:
thus, the method can be obtained by using any non-zero pixel point p on the imageiCorresponding observation vector viPicture element part of z-axisCoordinate system liThere is a coordinate transformation matrix with the world coordinate system of the camera as:
further, the step (4) is to use a non-zero pixel p
iDetermining a cone set by using the included angle radius r, the redundancy radius delta alpha and the radius stepping delta calculated in the step (2) and taking all generatrix vectors g of the cone set as an axis
iReverting to the pixel coordinate system, at each generatrix vector g
iqVoting in a three-dimensional Hough space H (m, n, r) of corresponding pixel coordinates, wherein
To obtain non-zero pixel point p
iThe three-dimensional Hough space with the circle center is specifically realized as follows:
establishing random non-zero pixel point piCorresponding observation vector viTaking the included angle radius r as a cone angle half-angle cone as an axis, the circumferential step is tau, and one generatrix vector g of the coneiqIn liIs represented as follows:
wherein
Then the generatrix vector g
iqCan be calculated from the following formula under the w system:
then, the model of the camera lens is used to calculate
Imaging of directionally incident light on an imagePixel coordinates are as follows:
[miq,niq]=round(Θ-1(-gi w)) (24)
where round () denotes rounding all the elements of () when point (m) is in the corresponding three-dimensional hough space H (m, n, r)iq,niqAnd r) voting once, namely accumulating the point value by 1, traversing q and r to obtain a non-zero pixel point piA three-dimensional Hough space as a circle center.
Further, in the step (5), the operations in the steps (3) and (4) are performed on all non-zero pixels, a pixel point with the highest ticket number in the three-dimensional hough space H (m, n, r) is the center of a snell window, an observation vector of a pixel coordinate of the center of the circle under a world coordinate system is solved, and a horizontal attitude is inverted, and the method is specifically realized as follows:
performing the operation of the step (3) and the operation of the step (4) on all non-zero pixel points in the binary edge image, voting and accumulating in the same three-dimensional Hough space H (m, n, r), and selecting a point (m, n, r) corresponding to the maximum value from the three-dimensional Hough space H (m, n, r)0,n0,r0) I.e. (m)0,n0,r0) argmaxH (m, n, r), then [ m0,n0]Is the coordinate of the center of the Snell window in the pixel coordinate system, r0The radius of the included angle of the Snell window, namely the actual included angle between the edge of the Snell and the center of the window;
zeta of observation vector of w system corresponding to central pixel of Snell windowwThe following are calculated by a camera lens model:
ζw=Θ([m0,n0]) (25)
the horizontal attitude of the carrier is determined by a roll angle gamma and a pitch angle theta, and under the condition of not considering a course angle, a carrier coordinate system, namely a b system, and a navigation coordinate system, namely an n system, and a coordinate conversion matrix are as follows:
because the world coordinate system (w system) of the camera and the carrier coordinate system (b system) are completely overlappedThus ζw=ζbSince the center of the Snell window is vertically upward under the n series, ζ isn=[0 0 1]TAccording to the coordinate transformation relationship, the following are provided:
the following formula is solved: zetaw=[-cosθsinγ sinθ cosθcosγ]TFrom this, the carrier pitch angle θ and roll angle γ can be calculated as:
therein, ζw(. indicates) vector ζwThe central element. Thus, the horizontal posture is determined.
Compared with the prior art, the invention has the following advantages: the existing underwater carrier attitude determination is mostly based on a combination mode of inertial navigation and other navigation modes or other navigation modes, and has the defects of error accumulation or non-autonomy and the like. The horizontal attitude determination method based on the Snell window is a method for solving space information through optical information in nature, can acquire absolute information of the horizontal attitude of a carrier in real time in a fully autonomous manner, has no error accumulation, and provides a new idea for underwater autonomous navigation.
Detailed Description
The technical solutions in the embodiments of the present invention will be described clearly and completely with reference to the accompanying drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, rather than all embodiments, and all other embodiments obtained by a person skilled in the art based on the embodiments of the present invention belong to the protection scope of the present invention without creative efforts.
As shown in fig. 1, the horizontal attitude determination method based on underwater snell window edge identification of the present invention specifically includes the following steps:
step 1, carrying out noise filtering and binarization processing on the acquired image, obtaining an edge image according to light intensity information in the image, and extracting coordinates [ m, n ] of a non-zero pixel p in the edge image in a pixel coordinate system]. Where any non-zero pixel at any point is denoted as piWith the coordinate of [ mi,ni]And i is 1,2, …, s, s is the number of nonzero pixels.
And 2, calculating a value set of the radius r of the included angle, wherein r is a possible value of the radius of the view angle of a Snell window, establishing a three-dimensional Hough space H (m, n, r), and setting all elements in the space to zero. The initialization of the included angle radius r is as follows:
as shown in fig. 2, the view angle of the whole sky in the atmosphere is approximately a celestial sphere, the image in the celestial sphere is compressed into an inverted cone view angle with an underwater observer as a vertex under the action of refraction of the water surface, and the bottom of the cone is a circular area on the water surface. In an empty observation image of an underwater observer, the circular area has higher brightness relative to the area outside the circle. And determining the angle between the edge and the center of the window according to the Snell's law of refraction and a lens model, taking the angle as the radius of an included angle in Hough circle transformation, and determining the redundant range of the radius of the included angle and radius stepping. The refractive indexes of atmosphere and water are respectively naAnd nwThen the angle α between the edge and the center of the snell window satisfies:
considering the errors of refractive indexes of a camera and a medium in practical application, the imaging size of a Snell window can not strictly accord with the theoretical model of the formula, so that the value range of the included angle radius r is [ alpha-delta alpha, alpha + delta alpha ] when the redundant radius delta alpha is set according to the practical water body environment and the lens model. Then, setting the included angle radius step as delta, and setting the value set of the included angle radius r as:
r∈{α-Δα,α-Δα+δ,α-Δα+2δ,…,α-Δα+kδ}(α-Δα+kδ≤α+Δα,k∈N+)
and establishing a three-dimensional Hough space H (m, n, r), wherein (m, n) is pixel coordinates, and initial values of the three-dimensional Hough space are all set to be 0.
Step 3, utilizing the edge image subjected to the binarization processing in the step 1 to divide each nonzero pixel p into a plurality of non-zero pixels
iThe observation vector under the world coordinate system, namely the w system, corresponding to the camera lens model theta is inverted through the camera lens model theta, and the observation vector is used as the z axis to establish the pixel local coordinate system, namely the l system
iSolving the non-zero pixel p
iCorresponding to l
iCoordinate transformation matrix between system and w system
The concrete implementation is as follows:
as shown in FIG. 3, in the pixel coordinate system, a non-zero pixel point piHas the coordinate of [ mi,ni]. Through camera lens calibration, a camera lens model Θ (p) can be calibrated. Then pixel point piCorresponding incident ray propagation vector kappa under system wiCan be expressed as:
κi=Θ(pi) (30)
κ
iis a unit vector, an observation vector v in the direction
iIs represented by in the w system
Order to
Then there are:
[ai bi ci]T=-Θ(pi) (31)
with v
iEstablishing a local coordinate system of picture elements, i.e./for the z-axis
iIs, definition w is and
ithe rotation angles around the x-axis and the y-axis during the conversion between the systemsIs omega
xi,
The rotation angle around the z-axis is 0, then w is l corresponding to the pixel
iCoordinate transformation matrix between systems
Can be expressed as the angle of 3-1-2 euler:
then l can be obtainediThe coordinate transformation matrix from system to w system is:
due to the observation vector v
iIn l
iIs represented as
Thus, the following steps are obtained:
combined formula (31), omegaxiAnd omegayiThe trigonometric function form of (a) can be expressed as follows:
thus, the method can be obtained by using any non-zero pixel point p on the imageiCorresponding observation vector viLocal coordinate system l of picture element as z-axisiThere is a coordinate transformation matrix with the world coordinate system of the camera as:
step 4, a non-zero pixel p is used
iDetermining a cone set by using the included angle radius r, the redundancy radius delta alpha and the radius stepping delta calculated in the step (2) and taking all generatrix vectors g of the cone set as an axis
iReverting to the pixel coordinate system, at each generatrix vector g
iqVoting in a three-dimensional Hough space H (m, n, r) of corresponding pixel coordinates, wherein
To obtain non-zero pixel point p
iThe three-dimensional Hough space with the circle center is specifically realized as follows:
establishing a non-zero pixel point piCorresponding observation vector viThe included angle radius r is used as a cone of a cone angle half angle. The circumference is stepped by τ. One of the generatrix vectors g of the coneiqIn liIs represented as follows:
wherein
Then the generatrix vector g
iqCan be calculated from the following formula under the w system:
then, the model of the camera lens is used to calculate
Imaging pixel coordinates on the image of the incident ray in the direction:
[miq,niq]=round(Θ-1(-gi w)) (39)
where round () denotes rounding all the elements of () to the whole. At this time at corresponding threeMiddle point (m) in dimension Hough space H (m, n, r)iq,niqAnd r) voting once, namely accumulating 1 by the point value. Traversing q and r to obtain non-zero pixel point piA three-dimensional Hough space as a circle center.
And 5, performing the operations of the step 3 and the step 4 on all non-zero pixels, taking the pixel point with the highest ticket number in the three-dimensional Hough space H (m, n, r) as the circle center of a Snell window, calculating an observation vector of the pixel coordinate of the circle center under a world coordinate system, and inverting the horizontal posture. The concrete implementation is as follows:
and (4) performing the operations of the step (3) and the step (4) on all non-zero pixel points in the binary edge image, and voting and accumulating in the same three-dimensional Hough space H (m, n, r). Then the circles with the non-zero pixel points at the edge of the snell window as the center of the circle intersect at the center of the snell window, i.e. the center coordinate of the snell window gets the most vote in the three-dimensional hough space. Therefore, the point (m) corresponding to the maximum value is selected from the three-dimensional Hough space H (m, n, r)0,n0,r0) I.e. (m)0,n0,r0) argmaxH (m, n, r). Then [ m0,n0]Is the coordinate of the center of the Snell window in the pixel coordinate system, r0Is the radius of the included angle of the Snell window, namely the actual included angle between the edge of the Snell and the center of the window.
Zeta of observation vector of w system corresponding to central pixel of Snell windowwCan be calculated by a camera lens model to obtain:
ζw=Θ([m0,n0]) (40)
the horizontal attitude of the carrier is determined by the roll angle gamma and the pitch angle theta. In the case of not considering the heading angle, the carrier coordinate system (b system) and the navigation coordinate system (n system) have the coordinate transformation matrix as follows:
since the world coordinate system (w system) of the camera and the carrier coordinate system (b system) completely coincide with each other, ζ is obtainedw=ζb. Since the center of the Snell window is vertically upward under the n series, ζ isn=[0 0 1]TAccording to the coordinate transformation relationship, the following are provided:
the following formula is solved: zetaw=[-cosθsinγ sinθ cosθcosγ]TThe carrier pitch angle θ and roll angle γ can thus be calculated as:
therein, ζw(. indicates) vector ζwThe central element. Thus, the horizontal posture is determined.
Although illustrative embodiments of the present invention have been described above to facilitate the understanding of the present invention by those skilled in the art, it should be understood that the present invention is not limited to the scope of the embodiments, but various changes may be apparent to those skilled in the art, and it is intended that all inventive concepts utilizing the inventive concepts set forth herein be protected without departing from the spirit and scope of the present invention as defined and limited by the appended claims.