CN111504276B - Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method - Google Patents

Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method Download PDF

Info

Publication number
CN111504276B
CN111504276B CN202010368099.6A CN202010368099A CN111504276B CN 111504276 B CN111504276 B CN 111504276B CN 202010368099 A CN202010368099 A CN 202010368099A CN 111504276 B CN111504276 B CN 111504276B
Authority
CN
China
Prior art keywords
matrix
objective function
propeller
equation
scale factor
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010368099.6A
Other languages
Chinese (zh)
Other versions
CN111504276A (en
Inventor
李东博
张贵阳
方乃文
魏荣婕
金浩
侯伟
王珺玮
刘明鹏
于春洋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Bojue Technology Co ltd
Original Assignee
Harbin Bojue Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Bojue Technology Co ltd filed Critical Harbin Bojue Technology Co ltd
Priority to CN202010368099.6A priority Critical patent/CN111504276B/en
Publication of CN111504276A publication Critical patent/CN111504276A/en
Application granted granted Critical
Publication of CN111504276B publication Critical patent/CN111504276B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
    • G01C11/04—Interpretation of pictures
    • G01C11/06—Interpretation of pictures by comparison of two or more pictures of the same area
    • G01C11/12—Interpretation of pictures by comparison of two or more pictures of the same area the pictures being supported in the same relative position as when they were taken
    • G01C11/14—Interpretation of pictures by comparison of two or more pictures of the same area the pictures being supported in the same relative position as when they were taken with optical projection
    • G01C11/20—Interpretation of pictures by comparison of two or more pictures of the same area the pictures being supported in the same relative position as when they were taken with optical projection in separate planes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
    • G01C11/02—Picture taking arrangements specially adapted for photogrammetry or photographic surveying, e.g. controlling overlapping of pictures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Multimedia (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

The invention discloses a method for acquiring a joint target function multi-propeller attitude angle based on a visual projection scale factor set. Step 1: establishing a coordinate system of the propellers of the multi-aircraft; step 2: the projection scale factor set is γ ═ γ { (γ })1,γ2,…,γn}; and step 3: optimizing a projection scale factor set, and establishing a multi-propeller attitude angle to obtain a minimum objective function; and 4, step 4: carrying out anti-symmetric matrix deformation processing on the minimum objective function in the step 3, and solving through a multi-element polynomial equation; and 5: further processing the step 4 by using a Kronecker product and a matrix vectorization function; step 6: the minimum objective function can therefore be finally expressed as a polynomial. The area array CCD camera realizes high-speed acquisition of images of multiple propellers, avoids the limitation of the response speed of a sensor in traditional contact measurement, has high response speed, can quickly arrange high-speed cameras, is convenient to implement, and is not limited by space.

Description

Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method
Technical Field
The invention belongs to the technical field of aircraft ground simulation test; in particular to a method for acquiring a combined target function multi-propeller attitude angle based on a visual projection scale factor set.
Background
The precise control of the attitude angle of the propeller is the basis for realizing the stable and safe operation of the aircraft, and the precision of the acquisition of the attitude angle of the propeller is directly related to the accuracy and the reliability of the control performance of the aircraft; however, the propeller is affected by high-pressure load, assembly error, machining precision, structural deformation and other factors, so that the actual pose change of the propeller and control input data have certain deviation, and the precise measurement technology of the pose of the propeller is one of the main factors for restricting whether the aircraft engine can realize engineering application.
From the existing literature and patent research, the research results of the method for acquiring the attitude angle information of the propeller are not too many. The method for acquiring the attitude of the thruster in the angular contact manner has high cost, is easy to introduce large manual errors, and is difficult to ensure the accuracy and the repeatability of data acquisition, so that only a non-contact acquisition method needs to be compared. Related researchers in the early century propose a propeller pivot angle acquisition method combining a linear array CCD and a self-luminous reticle. The method comprises the steps of respectively calculating the swing angle components of a propeller in the directions of an x axis and a y axis and the displacement in the direction of a z axis through the displacement of an upper marked line, and then converting the swing angles of the propeller in the directions of the x axis and the y axis into swing angles of the propeller, wherein the swing angle measuring range is only 0 to +/-8 degrees; in addition, the center of gravity of the propeller is necessarily moved under the condition of actual high-pressure load, so that the method can obtain high measurement precision theoretically, but the actual measurement precision is difficult to evaluate; then, measuring the attitude angle by a combination method is proposed, and by fixing the laser on the propeller, when the propeller swings, the laser beam swings along with the propeller; the method has higher precision on displacement measurement, but because the method adopts the displacement measurement to replace the angle measurement and only adopts a single laser beam to carry out the displacement measurement, the method is essentially a one-dimensional measurement method, so the measurement precision problem caused by the problems of the swing center movement, the axis movement and the like is difficult to solve; researchers also adopt a non-contact measurement method for parameters such as propeller axis, swing angle and the like of infrared photoelectric detection and conical surface fitting, some infrared reflective balls are arranged at specific positions of the propeller, and feature points are captured by an infrared camera, but the system is set up more complexly, and the reliability of measurement cannot be ensured; there are also some documents in which a method based on the fitting of a conical surface of a propeller is adopted, and the measurement error is inevitably increased due to the processing error of the surface of the propeller and the like.
However, the above method has the following limitations: (1) the attitude angle acquisition precision is low, and the requirement for high-precision measurement is difficult to meet; (2) the measurement cost is high, the flexibility is poor, and the measurement is difficult to realize in the actual propeller measurement; (3) the acquisition of the attitude angle of a single propeller is tried; (4) the attitude angle space-time consistency in different coordinate axis directions is not realized.
Disclosure of Invention
The invention adopts a non-contact sensor, namely an area array CCD camera to realize the high-speed acquisition of images of the multiple propellers, avoids the limitation of the response speed of the sensor in the traditional contact measurement, has high response speed, can quickly arrange the high-speed camera, is convenient to implement, and is not limited by space; identifying characteristic points belonging to different propellers through characteristic point marks; considering the working environment of the propeller and noise interference factors, introducing a projection scale factor set, and establishing an optimal objective function for estimating the motion parameters of the multi-propeller under the actual condition; and carrying out parametric decomposition on the attitude angle matrix to obtain a multi-element linear space equation set, and solving an optimal equation to finish high-precision acquisition of the attitude angles of the multiple propellers.
The invention is realized by the following technical scheme:
a method for acquiring a combined objective function multi-propeller attitude angle based on a visual projection scale factor set specifically comprises the following steps:
step 1: establishing a coordinate system of the propellers of the multi-aircraft;
step 2: in the coordinate system of step 1, the projection scale factor set is γ ═ γ1,γ2,…,γnIn which γiIs of a size satisfying gammai＝||Rqi w+T||；
And step 3: establishing a multi-propeller attitude angle to obtain a minimum objective function through an optimized projection scale factor set;
and 4, step 4: performing staggered matrix deformation processing on the minimum objective function in the step 3 to realize the rotating matrix RτObtaining a form solved by a multivariate linear space equation;
and 5: processing the step 4 by a matrix vectorization mode;
step 6: obtaining the minimum objective function of the multiple linear space equations processed in the step 4 and the step 5 and finally expressing the minimum objective function as containing
Figure BDA0002477206270000021
A polynomial of (a);
and 7: solving the polynomial in step 6.
Further, in step 1, the coordinates of the points in the world coordinate system are recorded as
Figure BDA0002477206270000022
Image point coordinates are noted as piEstablishing a two-dimensional point mapping relation from the space point coordinates to the imaging plane under a world coordinate system:
Figure BDA0002477206270000023
in the above formulaiIs a single projection factor, can accurately reflect the optical axis direction information of a target, is an important parameter for realizing the conversion from two-dimensional image information to three-dimensional stereo coordinates,
Figure BDA0002477206270000024
is a 3 x 3-dimensional orthogonal matrix,
Figure BDA0002477206270000025
is a 3 x 1 dimensional matrix.
Further, the step 2 is to mark an image point of a certain point on the propeller on the image plane as eiThe unit direction vector corresponding to is denoted as eiBased on vector unitization
Figure BDA0002477206270000031
This gives:
Figure BDA0002477206270000032
γ＝{γ1,γ2,…,γnis a set of projection scale factors, gammaiIs of a size satisfying
Figure BDA0002477206270000033
Further, the establishing of the attitude angles of the multiple propellers in the step 3 to obtain the minimum objective function is specifically,
Figure BDA0002477206270000034
wherein n represents the number of mark points participating in operation on each propeller; τ denotes the number of propellers, R*,T*And respectively representing a rotation matrix and a translation matrix of the attitude angle of the multi-propeller.
Further, in the step 4, specifically, let the vector υ ═ υ1υ2υ3]Then R isτWriting into:
Figure BDA0002477206270000035
formula (4)) In (1),
Figure BDA0002477206270000036
G3-a 3 x 3 dimensional unit array;
Figure BDA0002477206270000037
then it is a staggered matrix of v,
the number of parameters to be solved in the formula (3) is (n +6) tau, the number of elements in the projection scale factor set gamma is n tau, and the matrix R is rotatedτAnd TτThe number of parameters in the translation vector is 6 tau, and the translation vector T passes throughτAnd projection scale factor set gamma and rotation matrix RτThe space transformation relationship between the two is used for further deforming the expansion of the objective function.
Further, said further deforming the expansion of the objective function is specifically,
step 4.1: expression of the formula (2) by matrix operation is as follows
Figure BDA0002477206270000038
Order to
Figure BDA0002477206270000039
Equation (5) is written as:
ξΥ＝Hβ (6)
step 4.2: two sides of the formula (6) are simultaneously multiplied by xiTThen xiTXi there is an inverse matrix, resulting in:
Υ＝(ξTξ)-1ξTHβ (7)
constraint condition
Figure BDA0002477206270000041
Substitution, xiTThe xi matrix is represented as:
Figure BDA0002477206270000042
by
Figure BDA0002477206270000043
Deducing (xi)Tξ)-1，
Figure BDA0002477206270000044
Writing equation (6) in block matrix form:
Figure BDA0002477206270000045
step 4.3: united (9), matrix W(3×3n)Expressed as:
Figure BDA0002477206270000046
according to W(3×3n)Then obtain V(n×3n)Expression (c):
Figure BDA0002477206270000047
by the formulae (2) and (10), TτIs expressed as the following formula, and a projection scale factor set gamma is established at the same timeτWith respect to the transformation matrix RτAnd TτThe relational equation of (1):
Figure BDA0002477206270000051
reducing the unknown parameters in the discovered objective function from (n +6) tau to 3 tau;
step 4.4: by substituting equation (13) for equation (3), the objective function is transformed into:
Figure BDA0002477206270000052
further, the step 5 of processing the above expression by using a matrix vectorization function specifically includes the following steps:
Figure BDA0002477206270000053
wherein the matrix RτIs expressed as Vec (R)τ)＝[r11 r12 r13 r21 r22 r23 r31 r32 r33]T，
Three-dimensional point and unit matrix G under world coordinate system(3×3)The Kronecker product operation is carried out,
then the vector T is translatedτThe expression is simplified as follows:
Figure BDA0002477206270000054
the combined type (15) and the formula (16) obtain:
Figure BDA0002477206270000055
equation (17) is substituted for equation (14), and the objective function in this case is simplified as follows:
Figure BDA0002477206270000056
wherein
Figure BDA0002477206270000057
Further, the step 6 specifically includes: r in the formula (4)τIs substituted into equation (18), and the same holds for the objective function by the parameter upsilon ═ upsilon1υ2υ3]TConversion into a multivariate linear space equation xi (upsilon)1,υ2,υ3)＝0。
Further, the step 7 specifically includes: solving a polynomial:
Figure BDA0002477206270000061
then, the deviation is calculated for the formula (19) to make
Figure BDA0002477206270000062
Vs. upsilon1,υ2,υ3Is zero, and the inclusion parameter v is determined1,υ2,υ3The ternary polynomial equation set is continuously updated by assigning upsilon ← upsilon-delta upsilon through a formula, and the value of delta upsilon can be updated according to the equation
Figure BDA0002477206270000063
Updating a variation law, adjusting a damping factor mu in an iteration process to enable f (upsilon-delta upsilon) to be less than or equal to f (upsilon), accelerating convergence speed, and solving a parameter upsilon1,υ2,υ3After the values of (3), the rotation matrix R can be obtained by substituting the values of the equations (4) and (6)τAnd translation vector Tτ。
Drawings
FIG. 1 is a schematic view of the measurement of attitude angles of multiple propellers.
FIG. 2 is a graph showing the variation trend of the attitude angle error of the propeller 1# along with the noise standard deviation, (b) a graph showing the variation trend of the attitude angle error of the propeller 2# along with the noise standard deviation, (c) a graph showing the variation trend of the attitude angle error of the propeller 3# along with the noise standard deviation, and (d) a graph showing the variation trend of the attitude angle error of the propeller 4# along with the noise standard deviation.
FIG. 3 is a schematic representation of the operating efficiency of the present invention.
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, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
A method for acquiring a combined objective function multi-propeller attitude angle based on a visual projection scale factor set specifically comprises the following steps:
step 1: establishing a coordinate system of the propellers of the multi-aircraft;
step 2: in the coordinate system of step 1, the projection scale factor set is γ ═ γ1,γ2,…,γnIn which γiIs of a size satisfying
Figure BDA0002477206270000064
And step 3: establishing a multi-propeller attitude angle to obtain a minimum objective function through an optimized projection scale factor set;
and 4, step 4: performing staggered matrix deformation processing on the minimum objective function in the step 3 to realize the rotating matrix RτObtaining a form solved by a multivariate linear space equation;
and 5: processing the step 4 by a matrix vectorization mode;
step 6: obtaining the minimum objective function of the multiple linear space equations processed in the step 4 and the step 5 and finally expressing the minimum objective function as containing
Figure BDA0002477206270000071
A polynomial of (a);
and 7: solving the polynomial in step 6.
According to the coordinate system of the multi-aircraft propeller established in fig. 1, further, in step 1, the coordinates of the lower point of the world coordinate system are recorded as
Figure BDA0002477206270000072
Image point coordinates are noted as piThe two-dimensional point mapping relation from the space point coordinates to the imaging plane under the world coordinate system can be established:
Figure BDA0002477206270000073
in the above formulaiIs a single projection factor, can accurately reflect the optical axis direction information of a target, is an important parameter for realizing the conversion from two-dimensional image information to three-dimensional stereo coordinates,
Figure BDA0002477206270000074
is a 3 x 3-dimensional orthogonal matrix,
Figure BDA0002477206270000075
is a 3 x 1 dimensional matrix.
Further, the step 2 is to mark an image point of a certain point on the propeller on the image plane as eiThe unit direction vector corresponding to is noted as
Figure BDA0002477206270000076
Based on vector unitization
Figure BDA0002477206270000077
Thus can obtain
Figure BDA0002477206270000078
γ＝{γ1,γ2,…,γnIs a set of projection scale factors, gammaiIs of a size satisfying
Figure BDA0002477206270000079
In some existing methods, a homography relation between a space three-dimensional coordinate and a plane two-dimensional coordinate is established directly under the constraint condition of a vision measurement model, however, in the actual working process of a propeller, due to the existence of serious noise interference (smoke, vibration and environmental factors of the propeller), a propeller attitude angle obtained according to a traditional method has a large deviation from a true value. The step 3 of establishing the attitude angles of the multiple propellers to obtain the minimum objective function is specifically that,
Figure BDA00024772062700000710
wherein n represents the number of mark points participating in operation on each propeller; τ denotes the number of propellers, R*Rotation matrix, T, representing attitude angles of multiple propellers*A translation matrix representing multi-propeller attitude angles.
If the Hessian matrix and the Jacobian matrix are directly solved for the formula (3), huge calculation amount is undoubtedly caused due to more involved characteristic points, the dimension of the matrix is very high, and the reliability of global solution is difficult to guarantee along with the increase of the number of the thrusters. Therefore, in order to improve the solving efficiency and reliability, the invention adopts the deformation processing of the staggered matrix to realize the rotating matrix RτThe form that can be solved by a multi-element linear space equation is obtained. The step 4 is specifically to make the vector upsilon ═ upsilon1υ2υ3]Then R isτWriting into:
Figure BDA0002477206270000081
in the formula (4), the reaction mixture is,
Figure BDA0002477206270000082
G3-a 3 x 3 dimensional unit array;
Figure BDA0002477206270000083
then is a staggered matrix of υ.
The number of parameters to be solved in the formula (3) is (n +6) tau, the number of elements in the projection scale factor set gamma is n tau, and the matrix R is rotatedτAnd TτThere are 6 τ total parameters in the translation vector. Will next pass through translation vector TτAnd projection scale factor set gamma and rotation matrix RτThe space transformation relationship between the two is used for further deforming the expansion of the objective function.
Further, said further deforming the expansion of the objective function is specifically,
step 4.1: expression of the formula (2) by matrix operation is as follows
Figure BDA0002477206270000084
Order to
Figure BDA0002477206270000085
Equation (5) can be written as:
ξΥ＝Hβ (6)
step 4.2: two sides of the formula (6) are simultaneously multiplied by xiTThen xiTXi there is an inverse matrix, where xiTA transposed matrix for ξ, resulting in:
Υ＝(ξTξ)-1ξTHβ (7)
where H is the rotation matrix R*Formed diagonal matrix, constraint condition
Figure BDA0002477206270000086
Substitution, xiTXi matrix can be expressed as
Figure BDA0002477206270000091
By
Figure BDA0002477206270000092
Deducing (xi)Tξ)-1
Figure BDA0002477206270000093
Wherein G is an n × n dimensional unit matrix, the formula (6) is written in the form of a block matrix,
Figure BDA0002477206270000094
step 4.3: united (9), matrix W(3×3n)Can be expressed as
Figure BDA0002477206270000095
According to W(3×3n)Then V can be obtained(n×3n)Expression (2)
Figure BDA0002477206270000096
Wherein V(n×3n)For the submatrix after Q matrix operation, by the formula (2) and the formula (10), TτCan be expressed as the following formula, and a projection scale factor set gamma is established at the same timeτWith respect to the transformation matrix RτAnd TτThe equation of the relationship (c) of (c),
Figure BDA0002477206270000097
the unknown parameters in the objective function which can be found at this time are reduced from (n +6) tau to 3 tau;
step 4.4: by substituting equation (13) for equation (3), the objective function is converted to
Figure BDA0002477206270000101
Wherein
Figure BDA0002477206270000102
In step 5, the matrix vectorization function is used to process the above formula, specifically, the processed form is as follows:
Figure BDA0002477206270000103
wherein the matrix RτIs expressed as Vec (R)τ)＝[r11 r12 r13 r21 r22 r23 r31 r32 r33]T，
Three-dimensional point and unit matrix G under world coordinate system(3×3)The Kronecker product operation is carried out,
Figure BDA0002477206270000104
then the vector T is translatedτThe expression can be simplified as:
Figure BDA0002477206270000105
wherein the combination of formula (15) and formula (17) results in:
Figure BDA0002477206270000106
step 4.5: then, formula (17) is substituted into formula (14), in which case the objective function can be simplified to
Figure BDA0002477206270000107
Wherein
Figure BDA0002477206270000108
R in the formula (4)τIs substituted into equation (18), and similarly the objective function may be represented by the parameter upsilon ═ upsilon1υ2υ3]TConversion to a spatially linear equation xi (upsilon)1,υ2,υ3) 0, the maximum power of the expanded formula is 4, as shown in formula (19)
Figure BDA0002477206270000109
Then, the deviation is calculated for the formula (19) to make
Figure BDA00024772062700001010
Vs. upsilon1,υ2,υ3The partial derivative of (a) is zero, namely the contained parameter upsilon can be obtained1,υ2,υ3The ternary space linear equation set is continuously updated by assigning upsilon ← upsilon-delta upsilon through a formula, and the value of delta upsilon can be updated according to the equation set
Figure BDA0002477206270000111
Updating a variation law, adjusting a damping factor mu in an iteration process to enable f (upsilon-delta upsilon) to be less than or equal to f (upsilon), accelerating convergence speed, and solving a parameter upsilon1,υ2,υ3After the values of (3), the rotation matrix R can be obtained by substituting the values of the equations (4) and (6)τAnd translation vector Tτ。
In the following, we shall verify the effectiveness of the method of the present invention through simulation experiments, respectively.
Example 2
Simulation experiment and analysis
The algorithm simulation program is executed on a Windows10 operating system with a CPU main frequency of 3.2GHz and an 8GB running memory, and is based on a Matlab 2017b software development platform. The feature points in the world coordinate system are randomly generated and distributed within a range of 4m × 4m × 4 m. The camera simulation parameters are set as follows: focal length f 20mm, radial distortion factor k1＝3.5×10-4，k2＝1.9×10-5Principal point coordinate is cx1280 and cy1240. The Euler angles at which the initial state is set are [ -10 °, -10 °,0 ° ] respectively]Relative translation matrix T between multiple propellers11＝[0,0,0]′，T21＝[0,300,0]′，T31＝[300,300,0]′，T41＝[300,0,0]', the distance units are millimeters. Firstly, the influence of noise on the measurement accuracy of the attitude angle of the propeller is analyzed. Gaussian noise with a mean value of zero and a standard deviation σ ranging from 0 to 5 was added to the simulation data. The experiment of 100 groups is repeated, and the attitude angle measurement errors of the four propellers are marked along with noiseThe trend of the increase in the tolerance is shown in fig. 2.
As can be seen from fig. 2, since the propellers themselves have no roll angle, but self-rotation occurs with the aircraft itself in actual flight, there is no need to separately consider the roll angle during the algorithmic simulation. Even if the roll angle is theoretically zero, when the rotational matrix is solved and euler angle decomposition is performed, the component of the roll angle exists, and as can be seen from the figure, the roll angle error is small and is relatively stable. In contrast, as the standard deviation of the added noise of the yaw angle and the pitch angle is increased, the measurement error of the attitude angle is gradually increased, but still can be stabilized within 0.15 degrees, and for the traditional method, when the mean square error of Gaussian noise reaches 5, the error range can reach 0.6-1.2 degrees, so that the method has better stability and robustness and is superior to the performance of the traditional method.
In order to test the operation efficiency of the method, the change of the operation time of the method along with the number of the participating operation points is recorded, as shown in fig. 3.
It can be seen from fig. 3 that when the number of feature points participating in the operation does not exceed 10 points, the operation time of the method of the present invention does not change significantly, when the number of feature points exceeds 40, the operation time starts to increase, and when the number of feature points reaches 200, the operation time approaches 105 milliseconds. It must be noted, however, that the attitude angles of four propellers are measured simultaneously when the method is operated, so that compared with the reference method 2 with the highest operation efficiency, if the operation time is multiplied by 4, the obtained time value is larger than that of the proposed method, and therefore, when the method is used for processing the information of four propellers simultaneously, the higher processing efficiency is still maintained. In addition, the reference method 1 has the largest operation error and relatively low precision and stability, so that the method has good comprehensive performance and popularization value in comprehensive evaluation of the performance.

Claims (8)

1. A method for acquiring a combined objective function multi-propeller attitude angle based on a visual projection scale factor set is characterized by specifically comprising the following steps:
step 1: establishing a coordinate system of the propellers of the multi-aircraft;
step 2: in the coordinate system of step 1, the projection scale factor set is γ ═ γ1,γ2,…,γnIn which γiIs of a size satisfying
Figure FDA0003543327510000011
Wherein R, T are the intermediate values of the rotation matrix and the translation vector in the iterative process,
Figure FDA0003543327510000012
three-dimensional coordinates of space points under a world coordinate system;
and step 3: establishing a multi-propeller attitude angle to obtain a minimum objective function through an optimized projection scale factor set;
and 4, step 4: performing staggered matrix deformation processing on the minimum objective function in the step 3 to realize the rotating matrix RτObtaining a form solved by a multivariate linear space equation;
and 5: processing the step 4 by a matrix vectorization mode;
step 6: obtaining the minimum objective function of the multiple linear space equations processed in the step 4 and the step 5 and finally expressing the minimum objective function as containing
Figure FDA0003543327510000013
A polynomial of (a);
2. the method according to claim 1, wherein the step 1 is to record coordinates of points in the world coordinate system as
Figure FDA0003543327510000014
Image point coordinates are noted as piEstablishing a two-dimensional point mapping relation from the space point coordinates to the imaging plane under a world coordinate system:
Figure FDA0003543327510000015
in the above formulaiIs a single projection factor, can accurately reflect the optical axis direction information of a target, is an important parameter for realizing the conversion from two-dimensional image information to three-dimensional stereo coordinates,
Figure FDA0003543327510000016
is a 3 x 3-dimensional orthogonal matrix,
Figure FDA0003543327510000017
is a 3 x 1 dimensional matrix.
3. The method according to claim 1, wherein step 2 is implemented by recording an image point of a point on the propeller on the image plane as eiThe unit direction vector corresponding to is noted as
Figure FDA0003543327510000018
Based on vector unitization
Figure FDA0003543327510000019
This gives:
Figure FDA00035433275100000110
γ＝{γ1,γ2,…,γnis a set of projection scale factors, gammaiIs of a size satisfying
Figure FDA00035433275100000111
Figure FDA00035433275100000112
Is the average of the coordinates of the points mapped to the image plane.
4. The method according to claim 1, wherein the establishing the minimum objective function for obtaining the attitude angles of the multiple propellers in the step 3 is specifically,
Figure FDA0003543327510000021
wherein n represents the number of mark points participating in operation on each propeller; τ denotes the number of propellers, R*,T*Respectively representing a rotation matrix and a translation matrix of the attitude angles of the multiple propellers;
Figure FDA0003543327510000022
is three-dimensional coordinate of space points under a plurality of propeller world coordinate systems,
Figure FDA0003543327510000023
mapping to an average, T, of image plane point coordinates for multiple thrustersτIs the translation vector of a plurality of propellers in the iterative process.
5. The method according to claim 3, wherein the step 4 is specifically to make the vector upsilon ═ upsilon1υ2υ3]Then R isτWriting into:
Figure FDA0003543327510000024
in the formula (4), the reaction mixture is,
Figure FDA0003543327510000025
G3-a 3 x 3 dimensional unit array;
Figure FDA0003543327510000026
then it is a staggered matrix of v,
the number of parameters to be solved in the formula (3) is (n +6) tau, the number of elements in the projection scale factor set gamma is n tau, and the matrix R is rotatedτAnd TτThe number of parameters in the translation vector is 6 tau, and the parameters are translatedVector TτAnd projection scale factor set gamma and rotation matrix RτThe space transformation relationship between the two is used for further deforming the expansion of the objective function.
6. The method according to claim 5, wherein said further deforming the expansion of the objective function is in particular,
step 4.1: expression of the formula (2) by matrix operation is as follows
Figure FDA0003543327510000027
Order to
Figure FDA0003543327510000028
Equation (5) is written as:
ξΥ＝Hβ (6)
step 4.2: two sides of the formula (6) are simultaneously multiplied by xiTThen xiTXi there is an inverse matrix, resulting in:
Υ＝(ξTξ)-1ξTHβ (7)
constraint condition
Figure FDA0003543327510000031
Substitution, xiTThe xi matrix is represented as:
Figure FDA0003543327510000032
by
Figure FDA0003543327510000033
Deducing (xi)Tξ)-1，
Figure FDA0003543327510000034
Writing equation (6) in block matrix form:
Figure FDA0003543327510000035
step 4.3: united (9), matrix W(3×3n)Expressed as:
Figure FDA0003543327510000036
according to W(3×3n)Then obtain V(n×3n)Expression (c):
Figure FDA0003543327510000037
by the formulae (2) and (10), TτIs expressed as the following formula, and a projection scale factor set gamma is established at the same timeτWith respect to the transformation matrix RτAnd TτThe relational equation of (1):
Figure FDA0003543327510000041
reducing the unknown parameters in the discovered objective function from (n +6) tau to 3 tau;
step 4.4: by substituting equation (13) for equation (3), the objective function is transformed into:
Figure FDA0003543327510000042
7. the method according to claim 6, wherein the step 5 processes the above formula by using a matrix vectorization function, specifically, the processed form is as follows:
Figure FDA0003543327510000043
wherein the matrix RτIs expressed as Vec (R)τ)＝[r11 r12 r13 r21 r22 r23 r31 r32 r33]T，
Three-dimensional point and unit matrix G under world coordinate system(3×3)The Kronecker product operation is carried out,
then the vector T is translatedτThe expression is simplified as follows:
Figure FDA0003543327510000044
the combined type (15) and the formula (16) obtain:
Figure FDA0003543327510000045
equation (17) is substituted for equation (14), and the objective function in this case is simplified as follows:
Figure FDA0003543327510000046
wherein
Figure FDA0003543327510000047
8. The method according to claim 7, wherein the step 6 is specifically: r in the formula (4)τIs substituted into equation (18), and the same holds for the objective function by the parameter upsilon ═ upsilon1υ2υ3]TConversion into a multivariate linear space equation xi (upsilon)1,υ2,υ3)＝0。
CN202010368099.6A 2020-04-30 2020-04-30 Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method Active CN111504276B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010368099.6A CN111504276B (en) 2020-04-30 2020-04-30 Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010368099.6A CN111504276B (en) 2020-04-30 2020-04-30 Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method

Publications (2)

Publication Number Publication Date
CN111504276A CN111504276A (en) 2020-08-07
CN111504276B true CN111504276B (en) 2022-04-19

Family

ID=71871748

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010368099.6A Active CN111504276B (en) 2020-04-30 2020-04-30 Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method

Country Status (1)

Country Link
CN (1) CN111504276B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117485595B (en) * 2023-11-21 2024-04-05 北京易动宇航科技有限公司 Energy management method for electric propulsion system

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5388059A (en) * 1992-12-30 1995-02-07 University Of Maryland Computer vision system for accurate monitoring of object pose
CN103630120A (en) * 2013-07-16 2014-03-12 中国人民解放军信息工程大学 Mars surface linear array image epipolar ray resampling method based on strict geometric model
CN107680159A (en) * 2017-10-16 2018-02-09 西北工业大学 A kind of space non-cooperative target three-dimensional rebuilding method based on projection matrix
CN107871327A (en) * 2017-10-23 2018-04-03 武汉大学 The monocular camera pose estimation of feature based dotted line and optimization method and system
CN110375738A (en) * 2019-06-21 2019-10-25 西安电子科技大学 A kind of monocular merging Inertial Measurement Unit is synchronous to be positioned and builds figure pose calculation method
CN110500954A (en) * 2019-07-30 2019-11-26 中国地质大学（武汉） A kind of aircraft pose measuring method based on circle feature and P3P algorithm

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012044216A1 (en) * 2010-10-01 2012-04-05 Saab Ab Method and apparatus for solving position and orientation from correlated point features in images
US10126115B2 (en) * 2014-09-18 2018-11-13 Nec Corporation Triangulation device, triangulation method, and recording medium recording program therefor

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5388059A (en) * 1992-12-30 1995-02-07 University Of Maryland Computer vision system for accurate monitoring of object pose
CN103630120A (en) * 2013-07-16 2014-03-12 中国人民解放军信息工程大学 Mars surface linear array image epipolar ray resampling method based on strict geometric model
CN107680159A (en) * 2017-10-16 2018-02-09 西北工业大学 A kind of space non-cooperative target three-dimensional rebuilding method based on projection matrix
CN107871327A (en) * 2017-10-23 2018-04-03 武汉大学 The monocular camera pose estimation of feature based dotted line and optimization method and system
CN110375738A (en) * 2019-06-21 2019-10-25 西安电子科技大学 A kind of monocular merging Inertial Measurement Unit is synchronous to be positioned and builds figure pose calculation method
CN110500954A (en) * 2019-07-30 2019-11-26 中国地质大学（武汉） A kind of aircraft pose measuring method based on circle feature and P3P algorithm

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于直线特征的空间非合作目标位姿视觉测量方法研究;张跃强;《中国优秀博硕士学位论文全文数据库（博士） 工程科技Ⅱ辑》;20171115;第18-150页 *
基于视觉的空间目标位置姿态测量方法研究;尚洋;《中国优秀博硕士学位论文全文数据库（博士） 工程科技Ⅱ辑》;20080615;第40-66页 *

Also Published As

Publication number Publication date
CN111504276A (en) 2020-08-07

Similar Documents

Publication Publication Date Title
CN113532311B (en) Point cloud splicing method, device, equipment and storage equipment
CN110189399B (en) Indoor three-dimensional layout reconstruction method and system
CN110497373B (en) Joint calibration method between three-dimensional laser radar and mechanical arm of mobile robot
JP5627325B2 (en) Position / orientation measuring apparatus, position / orientation measuring method, and program
Wang et al. A screw axis identification method for serial robot calibration based on the POE model
CN112907735B (en) Flexible cable identification and three-dimensional reconstruction method based on point cloud
CN111983620B (en) Target positioning method for underwater robot searching and exploring
CN111754464B (en) Part accurate alignment method combining PD-like algorithm with ICP algorithm
Wang et al. A vision-based fully-automatic calibration method for hand-eye serial robot
Hesch et al. Mirror-based extrinsic camera calibration
Mashita et al. Calibration method for misaligned catadioptric camera
CN111504276B (en) Visual projection scale factor set-based joint target function multi-propeller attitude angle acquisition method
Ding et al. An efficient solution to the relative pose estimation with a common direction
Sun et al. An orthogonal iteration pose estimation algorithm based on an incident ray tracking model
CN114529703A (en) Entropy increase optimization-based point cloud global matching method for large complex components
KR102372298B1 (en) Method for acquiring distance to at least one object located in omni-direction of vehicle and vision device using the same
Grudziński et al. Stereovision tracking system for monitoring loader crane tip position
CN113916128A (en) Method for improving precision based on optical pen type vision measurement system
Jiang et al. 3D SLAM based on NDT matching and ground constraints for ground robots in complex environments
CN111553954A (en) Direct method monocular SLAM-based online luminosity calibration method
CN115619877A (en) Method for calibrating position relation between monocular laser sensor and two-axis machine tool system
CN115493616A (en) Method for evaluating field accuracy of laser tracking attitude angle
CN114266824A (en) Non-cooperative target relative pose measurement method and system based on deep learning
Liang et al. An integrated camera parameters calibration approach for robotic monocular vision guidance
Ming et al. A real-time monocular visual SLAM based on the bundle adjustment with adaptive robust kernel

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant