0 Introduction

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.

1 Positioning Data Preprocessing

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.

$ \left\{ \begin{array}{l} {P_t}.x = {\rm{sort}}\left( {{P_{t - 2n}}.x, \cdots ,{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_t}.x} \right).{\rm{middle}}\\ {P_t}.y = {\rm{sort}}\left( {{P_{t - 2n}}.y, \cdots ,{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_t}.y} \right).{\rm{middle}}\\ {P_t}.z = {\rm{sort}}\left( {{P_{t - 2n}}.z, \cdots ,{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_t}.z} \right).{\rm{middle}} \end{array} \right. $

(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, the n-th value is output as the median. Independent median filtering is performed in three dimensions, and the coordinates obtained are constrained as

$ {P_{t + 1}} = \left\{ \begin{array}{l} {P_t} + \frac{{{V_{{\rm{limit}}}}\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_{t + 1}} - {P_t}} \right)}}{{\left| {\frac{{{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_{t + 1}} - {P_t}}}{{\Delta t}}} \right|}}\;\;\;\;\;{V_{{\rm{limit}}}} > \left| {\frac{{{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_{t + 1}} - {P_t}}}{{\Delta t}}} \right|\\ {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_{t + 1}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{V_{{\rm{limit}}}} \le \left| {\frac{{{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} }_{t + 1}} - {P_t}}}{{\Delta t}}} \right| \end{array} \right. $

(2)

where P_t and P_t+1 denote the coordinates after processing at time t and t+1, ${{\overset\frown{P}}_{t}}$ and ${{\overset\frown{P}}_{t+1}}$ the original coordinates at time t and t+1, respectively. Δt is the time interval between two refreshes of the positioning coordinate value and V_limit the maximum threshold for moving speed.

2 HMM-Based Target Tracking Model

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.

2.1 Hidden Markov model

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 variables x and observation variables y. x={x₁, x₂, …, x_n}, where x_i∈X represents the system state at time i. It is assumed that the state variables are state and unobservable. y={y₁, y₂, …, y_n}, where y_i∈Y represents the observational value at time i. 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, y_t is determined by x_t. At the same time, the state x_t at time t is only dependent of the state x_t-1 at time t-1 rather than other states.

2.2 Model of target motion

(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

$ \dot x\left( t \right) = \mathit{\boldsymbol{F}}x\left( t \right) + \mathit{\boldsymbol{G}}a\left( t \right) $

(3)

where

$ x\left( t \right) = \left\{ \begin{array}{l} {\rm{target}}\;{\rm{position}}\;{\rm{at}}\;{\rm{time}}\;t\\ {\rm{target}}\;{\rm{speed}}\;{\rm{at}}\;{\rm{time}}\;t \end{array} \right. $

$ a\left( t \right) = {\rm{target}}\;{\rm{acceleration}}\;{\rm{at}}\;{\rm{time}}\;t $

$ \mathit{\boldsymbol{F}} = \left[ {\begin{array}{*{20}{c}} 0&1\\ 0&0 \end{array}} \right],\mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} 0\\ 1 \end{array}} \right]. $

The acceleration a(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 is P₀ when 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 [-a_lm, a_lm] as "no acceleration", the target undergoes no acceleration with a probability P₀.

(2) The target can accelerate at a maximum rate a_max(-a_max) and will do each with a probability P_max.

(3) In other cases, the acceleration is averagely distributed within (-a_max, -a_lm) and (a_lm, a_max).

(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.

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 the z axis. k is 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 radius a_lm as "no acceleration", the target undergoes no acceleration with a probability P₀.

(2) The target can accelerate at a maximum rate a_max with a probability P_max.

(3) In other cases, the acceleration is averagely distributed within (a_lm, a_max), the corresponding acceleration probability P_t is related to the angle α and the acceleration intensity k.

The cross-section of 3D acceleration distribution sphere is taken, as shown in Fig. 4.

The corresponding formulas can be obtained as

$ \left\{ \begin{array}{l} {h_0} = r\sin \alpha \\ {h_1} = r\sin \left( {\alpha + \gamma } \right)\\ {l_1} = r\cos \left( {\alpha + \gamma } \right)\\ {l_2} = r\cos \alpha \\ v = \frac{1}{3}{\rm{ \mathsf{ π} }}\left( {l_2^2{h_0} - l_1^2{h_1}} \right) + \int_{{h_0}}^{{h_1}} {{\rm{ \mathsf{ π} }}\left( {{r^2} - {z^2}} \right){\rm{d}}z} \end{array} \right. $

(4)

where v is the volume of the space body formed by 2π rotation of AOB around the z axis. Without considering P₀ and P_max, the probability density of acceleration distribution in interval (a_lm, a_max) could be calculated.

$ f = \frac{{3\left( {1 - {P_0} - {P_{\max }}} \right)}}{{4{\rm{ \mathsf{ π} }}\left( {a_{\max }^3 - a_{\min }^3} \right)}} $

(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

$ {P_{{\rm{block}}}} = \frac{{2{\rm{ \mathsf{ π} }}f\left( {{k^3} - a_{lm}^3} \right)}}{3}\left[ {\sin \left( {\alpha + \gamma } \right) - \sin \alpha } \right] $

(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.

$ P = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} {P_0}\\ {P_{{\rm{block}}}}\\ {P_{\max }} \end{array}&\begin{array}{l} k \in \left[ {0,{a_{lm}}} \right]\\ k \in \left( {{a_{lm}},{a_{\max }}} \right)\\ k \in \left( {{a_{\max }}, + \infty } \right) \end{array} \end{array}} \right. $

(7)

The hidden state is discretized 3D accelerationa. It corresponds to the discretized block in the 3D acceleration sphere.

2.3 HMM parameters

A HMM could be concisely represented by a three tuple λ=[A, B, Π], where A represents the hidden state transition probability matrix (HSTPM), B the observation state transition probability matrix (confusion matrix, CM), and Π the initial state probability matrix (ISPM). Let N denote the number of hidden states and M the number of observable states. The details are as follows.

HSTPM describes the transition probability between different hidden states in HMM, which is expressed as A =(A_ij)_N×N, A_ij=P(S_j|S_i), 1≤i, j≤N. It represents the probability that the state at time t+1 is S_j when the state at time t is S_i. 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 as B =(B_ij)_M×N, B_ij=P(O_j|S_i), 1≤i≤M, 1≤j≤N. It represents the probability that the observed value is O_j when the state at time t is S_i. 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 time t=1. It is expressed as Π =[P_S1 P_S2 … P_SN], where P_Si(1≤i≤N) is the probability of that the initial hidden state is S_i. The initial state is subject to the equal probability distribution.

3 MH-K Based Positioning Optimization

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.

3.1 Kalman filter

(1) State prediction stage

From the kinematic formula, which is

$ \left[ {\begin{array}{*{20}{c}} {{P_t}}\\ {{V_t}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&{\Delta t}\\ 0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{P_{t - 1}}}\\ {{V_{t - 1}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\Delta {t^2}/2}\\ {\Delta t} \end{array}} \right]{a_t} $

(8)

the acceleration a_t can be obtained by HMM. And the first formula of Kalman filter^[18] can be obtained as

$ {\mathit{\boldsymbol{X}}_{t\left| {t - 1} \right.}} = \mathit{\boldsymbol{F}}{\mathit{\boldsymbol{X}}_{t - 1\left| {t - 1} \right.}} + \mathit{\boldsymbol{B}}{a_t} $

(9)

where F is the state transition matrix, B the control matrix, and X the state vector.

$ \mathit{\boldsymbol{F}} = \left[ {\begin{array}{*{20}{c}} 1&{}&{}&{\Delta t}&{}&{}\\ {}&1&{}&{}&{\Delta t}&{}\\ {}&{}&1&{}&{}&{\Delta t}\\ {}&{}&{}&1&{}&{}\\ {}&{}&{}&{}&1&{}\\ {}&{}&{}&{}&{}&1 \end{array}} \right] $

(10)

$ \mathit{\boldsymbol{X}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_x}}&{{\mathit{\boldsymbol{P}}_y}}&{{\mathit{\boldsymbol{P}}_z}}&{{\mathit{\boldsymbol{V}}_x}}&{{\mathit{\boldsymbol{V}}_y}}&{{\mathit{\boldsymbol{V}}_z}} \end{array}} \right]^{\rm{T}}} $

(11)

$ \mathit{\boldsymbol{B}} = \Delta t\left[ {\begin{array}{*{20}{c}} {\Delta t/2}&{}&{}\\ {}&{\Delta t/2}&{}\\ {}&{}&{\Delta t/2}\\ 1&{}&{}\\ {}&1&{}\\ {}&{}&1 \end{array}} \right] $

(12)

After updating the coordinates through the motion model, the Kalman gain can be calculated as

$ {\mathit{\boldsymbol{P}}_{t\left| {t - 1} \right.}} = {\mathit{\boldsymbol{F}}_{t - 1\left| {t - 1} \right.}}{\mathit{\boldsymbol{F}}^{\rm{T}}} + {\mathit{\boldsymbol{Q}}_t} $

(13)

where P_t-1|t-1 is the posteriori covariance error matrix based on motion model prediction at time t-1, and Q_t the process noise error which accords with the Gauss distribution.

(2) State update stage

The coordinate updating in the measurement system needs to be considered.

$ {\mathit{\boldsymbol{Z}}_t} = \mathit{\boldsymbol{H}}{\mathit{\boldsymbol{X}}_{t\left| {t - 1} \right.}} + {\mathit{\boldsymbol{R}}_t} $

(14)

$ {\mathit{\boldsymbol{S}}_t} = \mathit{\boldsymbol{H}}{\mathit{\boldsymbol{P}}_{t\left| {t - 1} \right.}}{\mathit{\boldsymbol{H}}^{\rm{T}}} + {\mathit{\boldsymbol{R}}_t} $

(15)

where Z_t is the observation result at time t, S_t the measurement covariance at time t, R_t the measurement noise at time t, and H the observation vector.

$ \mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} 1&{}&{}&0&{}&{}\\ {}&1&{}&{}&0&{}\\ {}&{}&1&{}&{}&0 \end{array}} \right] $

(16)

The noise formula is as follows.

$ {\mathit{\boldsymbol{R}}_t} = \left[ {\begin{array}{*{20}{c}} {{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_x},{\mathit{\boldsymbol{P}}_x}} \right)}&{{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_x},{\mathit{\boldsymbol{P}}_y}} \right)}&{{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_x},{\mathit{\boldsymbol{P}}_z}} \right)}\\ {{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_y},{\mathit{\boldsymbol{P}}_x}} \right)}&{{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_y},{\mathit{\boldsymbol{P}}_y}} \right)}&{{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_y},{\mathit{\boldsymbol{P}}_z}} \right)}\\ {{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_x},{\mathit{\boldsymbol{P}}_z}} \right)}&{{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_y},{\mathit{\boldsymbol{P}}_z}} \right)}&{{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{P}}_z},{\mathit{\boldsymbol{P}}_z}} \right)} \end{array}} \right] $

(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.

$ {\mathit{\boldsymbol{K}}_t} = {\mathit{\boldsymbol{P}}_{t\left| {t - 1} \right.}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{S}}_t^{ - 1} $

(18)

$ {\mathit{\boldsymbol{X}}_{t\left| t \right.}} = {\mathit{\boldsymbol{X}}_{t\left| {t - 1} \right.}} + {\mathit{\boldsymbol{K}}_t}\left( {{\mathit{\boldsymbol{Z}}_t} - \mathit{\boldsymbol{H}}{\mathit{\boldsymbol{X}}_{t\left| {t - 1} \right.}}} \right) $

(19)

where K_t is 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.

$ {\mathit{\boldsymbol{P}}_{t\left| t \right.}} = {\mathit{\boldsymbol{P}}_{t\left| {t - 1} \right.}} - {\mathit{\boldsymbol{K}}_t}\mathit{\boldsymbol{H}}{\mathit{\boldsymbol{P}}_{t\left| {t - 1} \right.}} $

(20)

3.2 Implementation process of MH-K algorithm

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.

Step 1  UWB 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.

Step 2  The 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 {tag_id}.

Step 3  When the update event from the UWB positioning platform is obtained, the object storage space corresponding to the tag_id is 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 tag_id value is instantiated and registered to OMC.

Step 4  When 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".

Step 5  Data 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 (Stream_Median) 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 Stream_Median and transfer to Steps 6, 7 respectively.

Step 6  SSTP in Stream_Median flows into HMM, which outputs the displacement prediction at time t+1 relative to time t. The output is packaged into data stream (Stream_HMM) and transferred to the next step.

Step 7  The corresponding SSTP objects in Stream_Median and Stream_HMM are extracted. Kalman filtering is used to fuse the two results.

Step 8  The result is the output. The state value at time t is saved and RC is updated.

4 Experiment and Analysis

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 m²/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.

$ {\mathit{\boldsymbol{R}}_t} = \left[ {\begin{array}{*{20}{c}} {8.409\;964}&{3.509\;666}&{ - 1.230\;04}\\ {3.509\;666}&{10.344\;68}&{ - 0.768\;9}\\ { - 1.230\;04}&{ - 0.768\;9}&{0.583\;671} \end{array}} \right] $

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(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 cm². 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 cm², 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.

5 Conclusions

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.

Acknowledgements

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).