, *, , ,
1.College of Mechanical and Electrical Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, P.R. China;2.Jiangxi Hongdu Aviation Industry Group Co. Ltd., Nanchang 330001, P.R. China
Abstract: Due to the complexity of strong metal interference and multiple occlusions in aircraft assembly workshop, the random “drift” phenomenon often happens in the ultra wide band (UWB) based positioning system. To solve this, a fusion positioning optimization algorithm is proposed based on median filtering, hidden Markov model (HMM) and Kalman filtering. Firstly, based on the three-dimensional (3D) median filtering, a queue optimization method with weights is introduced to smooth the measurement data and eliminate the abnormal value. Secondly, taking Singer model as a reference, a single-dimension acceleration distribution model is designed. In order to further consider the spatial motion characteristics of objects in workshop, the distribution is extended from 1D to 3D, and discretized into the state quantity of HMM. Subsequently, the data obtained by the two methods are fused by taking Kalman filter as an iterator, and then the optimized location solution is obtained based on dynamic weights. Finally, an experiment is conducted in an aircraft assembly workshop. Results show that 99.3% of dynamic positioning errors are less than 15 cm after using the proposed algorithm. Even in the situation with large signal-fluctuation, there are 71.6% of positioning data whose errors are reduced. The random “drift” is greatly decreased.
Key words: ultra wide band (UWB); aircraft assembly workshop; hidden Markov model(HMM); median filtering; Kalman filtering
In “industry 4.0” plan, intelligent manufacturing has become a key development trend in today’s manufacturing industry. As a significant component, real time positioning technology plays an important role in constructing intelligent and transparent manufacturing workshop[1]. In aircraft assembly workshop, real time positioning can effectively solve the common problems that tools are easy to forget, material management is chaotic, and delivery vehicles are difficult to track. Real time accurate location and status information are key factors in promoting management level and assembly capacity[2]. The location information of manufacturing resources could provide data support for resource schedule[3], material distribution[4], and other activities. The manufacturing resources that need to be located include work in process, materials, tools, vehicles, workers, etc. By their accurate positioning, it is convenient for real time searching, monitoring, and tracking.
At present, UWB[5], radio frequency indentification (RFID)[6], Bluetooth[7]and Lidar[8]are the main positioning methods in manufacturing workshops. Due to its advantages of high security, great capacity, strong anti-interference ability, low power dissipation, and low cost, UWB technology has attracted more attention in recent years[9]. Furthermore, UWB adopts the active positioning method, which has high data transmission rate and strong ability to track the target. Its positioning accuracy could achieve the centimeter level, and the positioning requirement of the manufacturing resources could be satisfied. Therefore, many researchers have applied UWB to manufacturing workshop. Huang et al.[2]constructed a holographic workshop map to monitor the whole workshop site in a visual way based on UWB and RFID. Jiang et al.[5]combined automatic guided vehicle (AGV) with UWB-based real time location platform to solve the problem of poor flexibility and monitoring ability in traditional material delivery process. Zhou et al.[4]developed an UWB-based location system to monitor and navigate the forklifts in real time. The accurate material distribution of digital manufacturing workshop was realized. He et al.[10]and Leune et al.[11]applied UWB to AGV navigation and tracking, respectively. Baboli et al.[12]optimized and configured manufacturing systems by locating manufacturing resources in real time. Zhang et al.[13]improved the management efficiency of manufacturing resources through an UWB-based location system.
The UWB-based positioning of static objects in aircraft assembly workshop has high stability. However, for AGV, material and many other moving objects, the location data in workshop always have the phenomenon of irregular “drift”. This is mainly caused by the following reasons. On the one hand, the metal equipment and materials are densely distributed in workshop. These metals have strong ability to absorb electromagnetic waves. UWB signal is a microwave with shorter wavelengths, leading to its diffraction ability weak. When the microwave is injected into the metal, induction current could be produced on the metal surface, which leads to the signal attenuation. On the other hand, there are many large production and transportation equipments in the workshop. Because of this complicated environment, an object is often accompanied by non-line of sight and multipath propagation when it moves. This causes the UWB positioning system to “misjudge” the positioning data, so that it is seriously deviated from the real position.
Considering the above problems, this paper proposes an UWB-based target tracking model and a fusion positioning optimization algorithm (median filter and HMM based Kalman filter, MH-K). Median filter is introduced to smooth and denoise abnormal data. Then the Singer model[14]is extended from a single dimension to a discrete 3D model. Taking HMM[15]as a speculating method, the predicted value is output based on the displacement and acceleration state. The smoothed denoising value is combined with the model output value by the Kalman filter. The solution with the highest confidence degree is output. An experiment is conducted to verify that the proposed method could effectively solve the random “drift” problem.
Strong metal interference, frequent position change of assembly resources and high dynamicity of assembly process are typical characteristics in aircraft assembly workshop, which have great influence on the stability and accuracy of the UWB positioning results. Positioning data preprocessing is designed to reduce the interference of these complex factors as much as possible. By filtering, the positioning data is smoothed to reduce the noise of the elements’ movement path and the burr amplitude. For tackling such problems, median filtering has the advantages of low complexity, small hardware requirements and fast processing speed.
Median filtering is a sort of nonlinear signal processing technology, which can effectively suppress noise based on sorting statistics theory. In a number sequence, the value of a point is replaced by the median of all points in its one neighborhood. Median filtering is widely used in image processing. For real time location data, there is only preordering data without postorder data. In order to achieve similar results, a slide window is used as the data input mode. The purpose of sliding window is to record the location information of the same label ID, and reduce the computation cost.
The single dimension median filtering introduces a sliding window to add or reject data. The size of the sliding window is set as an odd number (2n+1). The coordinates in three dimensions of manufacturing resource need to be acquired simultaneously in UWB-based positioning system. The single dimension median filtering cannot be used directly. Therefore, median filtering processes of three dimensions are combined. On this basis, the priority queues (PQ) based data structure is used to improve operation efficiency.
(1)
As shown in Eq.(1), heap removes the most advanced data while adding new data. Then all data is sorted. Because the size of the window is 2n+1, then-th value is output as the median. Independent median filtering is performed in three dimensions, and the coordinates obtained are constrained as
(2)
UWB-based positioning platform could realize accurate positioning with 15 cm for static positioning in the manufacturing workshop[5]. It could meet the real time positioning requirements in workshop. However, the stability of the dynamic positioning error drops sharply when positioning the moving objects in the aircraft assembly workshop. A section of trajectory of an AGV is recorded and analyzed. The maximum deviation from the real trajectory is almost 70 cm, and the local area has the overall distortion of the trajectory. Obviously, the positioning data directly provided by UWB positioning platform are difficult to meet the requirements of locating moving objects accurately in aircraft assembly workshop. In this regard, a HMM-based target tracking method is proposed to predict the location of manufacturing resources in aircraft assembly workshop.
Target tracking model describes the changing process of the target with time, and it determines the effect of the target tracking to a great extent. In industry, the precise system model and the statistical characteristics of noise are not easy to be obtained, so it is difficult to describe them by mathematics. In this paper, a target tracking method based on HMM is used to predict the location of manufacturing resources in assembly workshop.
HMM is the simplest dynamic Bayesian network, which has been widely applied in time series data modeling, speech recognition and natural language processing. HMM focuses on three basic issues: (1) For the given model, how to effectively calculate the probability of producing the observation sequence; (2) For the given model and observation sequence, how to find the most matched state sequence with the observation sequence; (3) For the given observation sequence, how to adjust the model parameters to make the sequence most likely to appear.
Fig.1 shows a typical Markov process. Variables in HMM include state variablesxand observation variablesy.x={x1,x2,…,xn}, wherexi∈Xrepresents the system state at timei. It is assumed that the state variables are state and unobservable.y={y1,y2,…,yn}, whereyi∈Yrepresents the observational value at timei. The arrows in Fig.1 indicate the dependence between variables. At any time, the values of the observed variables depend only on the state variables, and are independent of other state and observation variables. That is,ytis determined byxt. At the same time, the statextat timetis only dependent of the statext-1at timet-1 rather than other states.
Fig.1 Schematic diagram of HMM structure
(1) Improved single-dimensional model of acceleration distribution
In a single physical dimension, the motion equation of mobile objects represented by AGV is represented as
(3)
where
a(t)=target acceleration at timet
The accelerationa(t) accounts for the target deviations from a straight line trajectory, so it is termed the target maneuver variable. In the Singer model[14], the probability isP0when the acceleration equals zero. However, the moving speed of the manufacturing resources is relatively low and is constrained by the measurement method. The acceleration is almost impossible to be “zero”. Therefore, combined with the characteristics of moving objects in manufacturing workshop, the Singer model is altered. A single-dimensional acceleration distribution model is given as shown in Fig.2, and the following assumptions are made.
(1) Regarding the interval [-alm,alm] as “no acceleration”, the target undergoes no acceleration with a probabilityP0.
(2) The target can accelerate at a maximum rateamax(-amax) and will do each with a probabilityPmax.
(3) In other cases, the acceleration is averagely distributed within (-amax,-alm) and (alm,amax).
Fig.2 Model of single-dimensional acceleration distribution
(2) 3D model of acceleration distribution
In ordinary indoor environment, most targets move in the horizontal plane. However, materials and tools sometimes move in the vertical direction in aircraft assembly workshop. Considering the spatial motion characteristics of objects in aircraft assembly workshop, the single-dimensional distribution is extended to three dimensions. In order to apply HMM, the 3D acceleration distribution sphere is divided into a great many blocks with the thickness ofδ, as shown in Fig.3.
Fig.3 3D acceleration distribution
Taking one of the blocks, the original point and the center point of the block is connected to be a connecting line.αis the coangle of the included angle between the connecting line and thezaxis.kis the length of the connecting line, that is, the intensity of acceleration.αis segmented withγas the basic unit. Based on the above description, the following contents are stipulated.
(1) Regarding the space range with radiusalmas “no acceleration”, the target undergoes no acceleration with a probabilityP0.
(2) The target can accelerate at a maximum rateamaxwith a probabilityPmax.
(3) In other cases, the acceleration is averagely distributed within (alm,amax), the corresponding acceleration probabilityPtis related to the angleαand the acceleration intensityk.
The cross-section of 3D acceleration distribution sphere is taken, as shown in Fig.4.
Fig.4 Cross-section of sphere
The corresponding formulas can be obtained as
(4)
wherevis the volume of the space body formed by 2π rotation ofAOBaround thezaxis. Without consideringP0andPmax, the probability density of acceleration distribution in interval (alm,amax)could be calculated.
(5)
The circumference is divided into 36 equal parts (γ=π/18). Each part is divided into 10 layers in the direction of the center according to the equal volume. The probability of a block can be calculated as
(6)
By such a discretization method, the 3D acceleration sphere could be discretized into a lot of blocks. Thus, the complete probability of any point/block in the 3D acceleration distribution sphere could be obtained after discretization.
(7)
The hidden state is discretized 3D accelerationa. It corresponds to the discretized block in the 3D acceleration sphere.
A HMM could be concisely represented by a three tupleλ=[A,B,Π], whereArepresents the hidden state transition probability matrix (HSTPM),Bthe observation state transition probability matrix (confusion matrix, CM), andΠthe initial state probability matrix (ISPM). LetNdenote the number of hidden states andMthe number of observable states. The details are as follows.
HSTPM describes the transition probability between different hidden states in HMM, which is expressed asA=(Aij)N×N,Aij=P(Sj|Si), 1≤i,j≤N. It represents the probability that the state at timet+1 isSjwhen the state at timetisSi. Here, HSTPM is subject to the 3D acceleration distribution after discretization.
CM describes the transition probability between different observation states in HMM, which is expressed asB=(Bij)M×N,Bij=P(Oj|Si), 1≤i≤M, 1≤j≤N. It represents the probability that the observed value isOjwhen the state at timetisSi. Here, displacement is taken as the observation state, and can be obtained by direct observation.
ISPM represents the probability matrix of the hidden state at timet=1. It is expressed asΠ=[PS1PS2…PSN], wherePSi(1≤i≤N) is the probability of that the initial hidden state isSi. The initial state is subject to the equal probability distribution.
Kalman filter is utilized to fuse the pose of moving objects acquired by HMM and the measured data after median filtering, making the output result more consistent with the true position trajectory[16]. Kalman filter belongs to a posteriori method[17], which is mainly composed of state prediction stage and state update stage.
(1) State prediction stage
From the kinematic formula, which is
(8)
the accelerationatcan be obtained by HMM. And the first formula of Kalman filter[18]can be obtained as
Xt|t-1=FXt-1|t-1+Bat
(9)
whereFis the state transition matrix,Bthe control matrix, andXthe state vector.
(10)
(11)
(12)
After updating the coordinates through the motion model, the Kalman gain can be calculated as
Pt|t-1=FPt-1|t-1FT+Qt
(13)
wherePt-1|t-1is the posteriori covariance error matrix based on motion model prediction at timet-1, andQtthe process noise error which accords with the Gauss distribution.
(2) State update stage
The coordinate updating in the measurement system needs to be considered.
Zt=HXt|t-1+Rt
(14)
St=HPt|t-1HT+Rt
(15)
whereZtis the observation result at timet,Stthe measurement covariance at timet,Rtthe measurement noise at timet, andHthe observation vector.
(16)
The noise formula is as follows.
(17)
At this point, the state observation results, the respective error matrix of the prediction model, and the measurement system are obtained. Then, the data of two independent systems are fused.
(18)
Xt|t=Xt|t-1+Kt(Zt-HXt|t-1)
(19)
whereKtis the error gain, and also a “balance coefficient” between the predicting values and the measuring values. The larger its value is, the lower the confidence level is, and vice versa.
Finally, the minimum mean square error matrix of Kalman filter is modified to prepare for the next iteration.
Pt|t=Pt|t-1-KtHPt|t-1
(20)
The proposed MH-K algorithm combines HMM, median filtering and Kalman filtering. In detail, it is mainly composed of the following parts.
(1) The object management container (OMC) assigns label data to the corresponding running container (RC).
(2) 3D-median algorithm is used to stabilize the fluctuation of data.
(3) HMM gives localization conjecture.
(4) Kalman filter is used to fuse the measured data processed by 3D-median filter and the data predicted by HMM.
The implementation process of MH-K algorithm is shown in Fig.5. The detailed operations are described as follows.
Step1UWB positioning platform is started from the service application programming interface (API). The downloading, burning and running of the binary function blocks of location base station are completed.
Step2The OMC in Utag (location tag in UWB positioning platform) is initialized. Key and value mapping based on the HashMap structure is used to register, address mapping, and life cycle control of an instantiated object under a set of Utag objects {tagid}.
Step3When the update event from the UWB positioning platform is obtained, the object storage space corresponding to the tagidis first tried to obtain from the OMC. If the address is successfully implemented, the data is directly passed into the RC. If it does not exist, the RC object that binds tagidvalue is instantiated and registered to OMC.
Fig.5 Implementation process of MH-K algorithm
Step4When RC accepts incoming data, it formats it into a standard space-time packet (SSTP), which is intended to be transferred to other component instances in the form of a “middleware”.
Step5Data enters 3D-median components. 3D-median tries to discard an outdated SSTP object from the smallest heap of PQ. At the initialization stage, the number of PQ objects needs to enter the climbing phase first, and the output will seriously deviate from the real value, so the data stream (StreamMedian) will stop at this point. During the normal operation, 3D-median components perform Sort and Correction operations after enqueue and dequeue operations. Then the output is packaged into StreamMedianand transfer to Steps 6, 7 respectively.
Step6SSTP in StreamMedianflows into HMM, which outputs the displacement prediction at timet+1 relative to timet. The output is packaged into data stream (StreamHMM) and transferred to the next step.
Step7The corresponding SSTP objects in StreamMedianand StreamHMMare extracted. Kalman filtering is used to fuse the two results.
Step8The result is the output. The state value at timetis saved and RC is updated.
In order to verify the proposed target tracking model and location algorithm, an UWB-based positioning platform (Ubisense7000) is constructed in an assembly workshop. A real time positioning data processing system is developed by using Intellij IDEA 2017.1 platform. Source data and processing results are stored in text documents intput.txt and output.txt, respectively.
The experiment is carried out by an AGV with active UWB tags. The AGV walks around a workstation in the workshop according to the designed route. The sampling target is a total of 2 000 records. The sampling interval is 0.1 s, the maximum speed is 2 m/s, and the maximum acceleration is 2 m2/s. In order to facilitate the inspection, the AGV runs as straight as possible at a constant speed except for turning.
By sampling 2 000 records of UWB tags in static state, the noise of measurement system is obtained as follows.
Classifing the collected data and extracting the corresponding label location data of ID number “020-000-192-253”, data processing is implemented by 3D-median and MH-K algorithm respectively, whose trajectories are shown in Fig.6.
Fig.6 AGV trajectories based on UWB positioning
Fig.6(a) is the global graph of sampling trajectory. The enlarged part of the upper left corner is shown in Fig. 6(b). Through the observation and analysis of the above experimental results, it is found that
(1) After using the proposed algorithm, UWB positioning platform greatly reduces the location instability caused by the complex environment of the workshop, and the location error is less than 20 cm.
(2) From the left side of the partial enlargement graph, it can be seen that the electromagnetic interference has an effect on the ultra wide band signal transmitted by the UWB tag. The trajectory before processing is obviously drifting, and the maximum error in the sampling range is produced. MH-K algorithm and 3D-median all can effectively alleviate this phenomenon. And the trajectory processed by MH-K is closer to the actual uniform linear motion trajectory.
(3) 3D-median could play a good role in smoothing the trajectory, and has strong processing power for instantaneous “drift”. However, it is affected by the window size. The larger the window is, the better the smoothing effect is, but the less sensitive it is for the direction change. Through multiple adjustment attempts, the size of the window is finally determined to be 9, which can better balance trajectory smoothing and directional response sensitivity.
In order to further analyze the influence of proposed algorithm on location accuracy, the experimental data are statistically analyzed, and the location error histogram is plotted, as shown in Fig.7. The following conclusions are drawn.
(1) Under the condition of great fluctuation, the average location error of UWB positioning platform without MH-K algorithm is 8.82 cm, the maximum is 57.89 cm, and the variance is 76.79 cm2. After using MH-K algorithm, the location error is 4.42 cm, which is 49.89% lower than that of before processing. The maximum is 16 cm, reduced by 72.36%. The variance is 10.5 cm2, reduced by 86.32 percent.
(2) After using MH-K algorithm, 71.6% of positioning error is less than that of before processing, 98.3% of positioning error is less than 10 cm, and 99.3% of positioning error is less than 15 cm. It can be seen that the proposed algorithm greatly improves the location stability of UWB in aircraft assembly workshop, and effectively solves the “drift” problem. In conclusion, it is able to meet the positioning requirements in workshop.
Fig.7 Histogram of positioning error
UWB could effectively improve the intelligent control and management ability for manufacturing resources in aircraft assembly workshop. However, the complex environment brings the problem of location instability, resulting in the “drift” phenomenon. So this paper proposes an optimization method which combines median filtering, HMM and Kalman filtering. The 3D median filtering is used to smooth positioning data and stabilize fluctuations. A discrete 3D acceleration distribution model is constructed based on HMM for location prediction. Data provided by the first two methods are fused by Kalman filter, and the final position coordinates are outputted. The positioning data obtained in the real manufacturing environment show that the proposed algorithm can effectively reduce the positioning error and greatly alleviate the positioning “drift” problem.
This work was supported in part by the National Natural Science Foundation of China (No.51575274), the National Defense Basic Scientific Research of China (No.JCKY2016605B006) and the Jiangxi Provincial Key R & D Program Project (No.20161ACE50004).
Transactions of Nanjing University of Aeronautics and Astronautics2018年6期