IEEE TRANSACTIONS ON MOBILE COMPUTING,VOL.XX.NO.XX,2020 oPa (a)Translation (b)Rotation (c)Translation with rotation Fig.7.Micro-movement decomposition 4.4.3 Translation with Rotation track the rigid motion of the tagged object,the design of the According to the definition of the rigid transformation,any arbi- rigid motion is shown in Fig.2.Detailed steps are as follows. trary rigid body motion can be decomposed into the combination of 1)Data Preprocessing:With RE-signals from the tag array, rotation and translation.Suppose a rigid body is attached with we extract the series of phase values for each tag from the a tag array T,when the center of the rigid body translates two antennas and segment them into the phase value for from the position P:to the position Pe,the rigid body also each snapshot.Then,we estimate the initial state of the rotates around a local rotation center Pa by the angle of a, tagged objects,including rotation state and rough position. the local rotation center has the same translation as the rigid 2)Movement Tracking:We derive the tag movement ac- body as well.Without loss of generality,we can model the cording to the phase variations,respectively,according to process of the rigid body motion into two successive steps the situations of the linear region and vast region(includ- one after the other,i.e.,performing the rotation first and ing both the linear region and non-linear region).Further, then the translation.Specifically,the rigid body first rotates we estimate the rigid transformation of the tagged object, around Pa by the angle of o,then it translates from position including the rotation and translation. P,to position Pe.According to Eq.(⑦and Eq.(⑧)，let 3)Movement Calibration:We detect the outliers of the and T be the coordinates of tag T when phase values from the tag array,by comparing the estimated the rigid body starts moving and ends moving,respectively, movement of each single tag with the estimated movement let (xa,ya)be the coordinates of local rotation center Pa of the tag array.Then,we eliminate the outlier(s)and re- when the rigid body starts moving,then: estimate the rigid transformation of the tagged object with the remaining tags for calibration. -Za=RTi-Ta +S. (9) Yi.e Ya Yi,8-ya According to Eq.(9),the movement of tag Ti,i.e. 5.1 Data Preprocessing for Rigid Motion Tracking [Ari,Ayi]T,can be decomposed into the following com- 5.1.1 Data Segmentation ponents: Due to the issues such as the multi-path effect and ambi- [△ Ti.e-Ti,s ent noises,the measured phase values may contain some △ ,e-,s =(R-IE,。xa] +S,(10) yi,s-Va fluctuations in the waveforms.Hence,after receiving the RF-signals from the tag array,we extract the series of phase where I is an identity matrix.Fig.7(c)shows an example values for each tag from the two antennas,and calibrate of the translation with rotation,when the rigid body is the phase values first.Specifically,due to the operation of attached with a rectangle tag array.Such rigid body motion mod in Eq.(1),the measured phase values are discontin- is equivalent to first rotating around Pa with angle o,i.e., uous.Thus,we stitch the phase values and remove the from the green array to the blue one,then translating from periodicity among the phase values.Besides,due to the P to Pe,i.e.,from the blue array to the yellow one. diversity term in Eq.(1),each tag has its own phase offset, so we measure the diversity term among tags in advance 5 DESIGN OF RIGID MOTION TRACKING and eliminate the tag diversity by offsetting the diversity We use a set of tags to track the rigid transformation of term.Further,we utilize the Kalman Filter to filter the cor- the tagged object,including the translation and the rotation. responding noises in the phase values.After that,suppose To support the tracking of small objects,we choose the there are n tags in the tag array,we segment the phases of small tag AZ9629 as the movement tag.Note that,during n tags from the two antennas into m snapshots,denoted the movement,the polarization angle of the tag relative to Θ1,1Θ1,2…Θ1，m the antenna changes as well,which can bring in additional as⊙= whereθi,=(日z,j,fu,i,》 phase offsets apart from the change of position.Whereas, Θn.1θn,2…θn,m for the movement tracking,we take each movement tag as means the phase values of tag Ti in the jth snapshot,0r.i.j a particle,and only focus on the phase offset caused by the and 0y.i.j represent the phase values from antenna A and position.Thus,we should reduce the phase offset due to the Ay,respectively.The time interval At for each snapshot change of the polarization angle during the movement of is usually set to a small value,e.g,At 200ms in our the tagged object.Empirically,we fix the tag orientation as implementation. shown in Fig.20,which is insensitive to the rotation due to the movement on the table. 5.1.2 Initial State Estimation RF-Dial extracts phase variations of the tag array re- According to Eq.(10),to compute the rigid transformation ceived by the orthogonally deployed RFID antenna pair to R,S,it is essential to determine the initial state of the
IEEE TRANSACTIONS ON MOBILE COMPUTING, VOL. XX, NO. XX, 2020 6 ௬ݏ ܲ௦ ܲ ௫ݏ ܺ ܻ ܱ (a) Translation ܺ ܻ ܲ ܲ௦ ߙ ܲ ܱ (b) Rotation ܲ ௬ݏ ௫ݏ ܲ௦(ܲ) ߙ ܺ ܻ ܱ (c) Translation with rotation Fig. 7. Micro-movement decomposition 4.4.3 Translation with Rotation According to the definition of the rigid transformation, any arbitrary rigid body motion can be decomposed into the combination of rotation and translation. Suppose a rigid body is attached with a tag array T, when the center of the rigid body translates from the position Ps to the position Pe, the rigid body also rotates around a local rotation center Pa by the angle of α, the local rotation center has the same translation as the rigid body as well. Without loss of generality, we can model the process of the rigid body motion into two successive steps one after the other, i.e., performing the rotation first and then the translation. Specifically, the rigid body first rotates around Pa by the angle of α, then it translates from position Ps to position Pe. According to Eq. (7) and Eq. (8), let [xi,s, yi,s] T and [xi,e, yi,e] T be the coordinates of tag Ti when the rigid body starts moving and ends moving, respectively, let (xa, ya) be the coordinates of local rotation center Pa when the rigid body starts moving, then: xi,e − xa yi,e − ya = R xi,s − xa yi,s − ya + S. (9) According to Eq. (9), the movement of tag Ti , i.e. [∆xi , ∆yi ] T , can be decomposed into the following components: ∆xi ∆yi = xi,e − xi,s yi,e − yi,s = (R − I) xi,s − xa yi,s − ya + S, (10) where I is an identity matrix. Fig. 7(c) shows an example of the translation with rotation, when the rigid body is attached with a rectangle tag array. Such rigid body motion is equivalent to first rotating around Pa with angle α, i.e., from the green array to the blue one, then translating from Ps to Pe, i.e., from the blue array to the yellow one. 5 DESIGN OF RIGID MOTION TRACKING We use a set of tags to track the rigid transformation of the tagged object, including the translation and the rotation. To support the tracking of small objects, we choose the small tag AZ9629 as the movement tag. Note that, during the movement, the polarization angle of the tag relative to the antenna changes as well, which can bring in additional phase offsets apart from the change of position. Whereas, for the movement tracking, we take each movement tag as a particle, and only focus on the phase offset caused by the position. Thus, we should reduce the phase offset due to the change of the polarization angle during the movement of the tagged object. Empirically, we fix the tag orientation as shown in Fig. 20, which is insensitive to the rotation due to the movement on the table. RF-Dial extracts phase variations of the tag array received by the orthogonally deployed RFID antenna pair to track the rigid motion of the tagged object, the design of the rigid motion is shown in Fig. 2. Detailed steps are as follows. 1) Data Preprocessing: With RF-signals from the tag array, we extract the series of phase values for each tag from the two antennas and segment them into the phase value for each snapshot. Then, we estimate the initial state of the tagged objects, including rotation state and rough position. 2) Movement Tracking: We derive the tag movement according to the phase variations, respectively, according to the situations of the linear region and vast region (including both the linear region and non-linear region). Further, we estimate the rigid transformation of the tagged object, including the rotation and translation. 3) Movement Calibration: We detect the outliers of the phase values from the tag array, by comparing the estimated movement of each single tag with the estimated movement of the tag array. Then, we eliminate the outlier(s) and reestimate the rigid transformation of the tagged object with the remaining tags for calibration. 5.1 Data Preprocessing for Rigid Motion Tracking 5.1.1 Data Segmentation Due to the issues such as the multi-path effect and ambient noises, the measured phase values may contain some fluctuations in the waveforms. Hence, after receiving the RF-signals from the tag array, we extract the series of phase values for each tag from the two antennas, and calibrate the phase values first. Specifically, due to the operation of mod in Eq. (1), the measured phase values are discontinuous. Thus, we stitch the phase values and remove the periodicity among the phase values. Besides, due to the diversity term in Eq. (1), each tag has its own phase offset, so we measure the diversity term among tags in advance and eliminate the tag diversity by offsetting the diversity term. Further, we utilize the Kalman Filter to filter the corresponding noises in the phase values. After that, suppose there are n tags in the tag array, we segment the phases of n tags from the two antennas into m snapshots, denoted as Θ = Θ1,1 Θ1,2 · · · Θ1,m · · · · · · · · · · · · Θn,1 Θn,2 · · · Θn,m , where Θi,j = hθx,i,j , θy,i,j i means the phase values of tag Ti in the jth snapshot, θx,i,j and θy,i,j represent the phase values from antenna Ax and Ay, respectively. The time interval ∆t for each snapshot is usually set to a small value, e.g, ∆t = 200ms in our implementation. 5.1.2 Initial State Estimation According to Eq. (10), to compute the rigid transformation R, S , it is essential to determine the initial state of the
IEEE TRANSACTIONS ON MOBILE COMPUTING,VOL.XX,NO.XX,2020 and T;collected from antenna Ar,the estimation of Ari.j, i.e.,Aij can be computed as follows: (0z,i-0z,)×/(4r), l9z,i-6x<元 △i,= (0x,i-8x-2r)×/(4π)，8z,i-0z.>π (0x,i-9x,+2r)×/(4π).8x,i-0z,j<-π Similarly,we can compute the estimated distance between Fig.8.Estimate the initial rotation angle of the object Ti and Ti along the Y-axis,i.e.,Ati.j.Meanwhile,according to Eq.(11),given the rotation angle B,we can also compute rigid body first.E.g.,the matrix T= Ti,s-Ta the theoretical value for Ari.j =i-xj and Ayi.j=yi-yj. yi,s-Ya in Eg.(10) Therefore,we leverage the Minimum Mean Square Error depends on the relative positions of tags in the tag array,i.e., (MMSE)estimator to estimate the rotation angle B,we the topology and rotation state of the tag array.Hence,since define the mean squared error e by enumerating the squared the rotation state of the tag array depends on both initial errors for all pairs of tags,that is, and subsequent rotation angles of the tag array,we need to estimate the initial rotation angle first. e(B)= ∑∑（△r()-△立）2+（△(）-△P) Without loss of generality,we use the rectangle tag array i=1j=1 as an example to illustrate our method to compute the initial We then compute the optimal value of B that achieves the rotation angle for the tag array.As shown in Fig.8,we set minimal mean squared error for e(B): the rotation center Pa at the center of the rectangle,and set two orthogonal polar axes l and ly according to the X B*=arg min e(B). and Y-axis in the global coordinate system.Each tag can be regarded as a particle point in the coordinate system,e.g., Vast Region:The vast region usually covers a large area T~T4.The initial rotation angle of the tag array,i.e.,the of 4 x 4m2 in the intersection of the scanning areas for angle between Ir and PaT1,is B. two orthogonally deployed antennas.As the vast region Topology Matrix T.According to the rotation angle B, includes both the linear region and non-linear region,we Ti,s-Ta thus consider the solution for the non-linear region,as it is we can depict the matrix T for each tag Yi.8-Ya a more generalized solution for both regions.Suppose the Ti.E.g.,for the rectangle tag array in Fig.8,let TiTall tag array is located at the non-linear region,the distance T2T3ll h,ITT2ll IT3Tll w,PaTill L, between the tags are not linear to the phase differences LT1PaT2 =n,these parameters can be regarded as con- between the tags,which depends on the exact positions of stants,as long as the topology of the tag array is fixed.Then, the corresponding tags.Therefore,it is essential to estimate according to the relative positions of the tags in the rectangle both the position and rotation state of the tag array in the tag array, non-linear region.For the tag array,we set the position of the rotation center as (x,y),and set the rotation angle as [1-za:y-ya]T=[l cos B,Isin B]T, B.Then,different from the solution in the linear region,in [2-a:y2-ya]T=[l cos(8+n),Isin (8+n)], regard to the non-linear region,for any tag Ti in the tag array,we estimate the Euclidean distance di between tag Ti [s-za:U3-ya]T=[-lcos(B),-Isin (B)]T, (11) and the corresponding antenna(A-or Ay),and compute [z4-xa,4-aT=【-lcos(B+n),-lsin(B+nj】T the difference of the Euclidean distance Adi.j=di-dj for any two arbitrary tags Ti and Ti.Specifically,for antenna As aforementioned in Section 4.2,the distances between Ar,the difference of the Euclidean distance Adr.i.j can be tag pairs are linear/non-linear to the phase differences computed as follows: between tag pairs in the linear region/non-linear region, (0z,i-0z,）×入/(4π)： 10zi-0z.jl<T respectively.In the following,we provide solutions to figure out the rotation angle B for the linear region and the vast △dzij (0x,i-0z.-2r）×λ/（4r),0x.i-0z,j>T region (including both the linear region and non-linear (0z,i-0x,1+2x）×λ/(4r).0z,i-8x,<-元 region),respectively. Similarly,we can also compute the difference of the Eu- Linear Region:The linear region usually has an area clidean distance Ady.i.j for antenna Ay.Meanwhile,ac- of 0.6 x 0.6m2 in the intersection of the central beam cording to the position (x,y)and rotation angle of the tag regions for two orthogonally deployed antennas.If we set array B,we can also compute the theoretical value for the the operation area in this linear region,e.g.,a tabletop with difference of the Euclidean distance Adz.i.j and Ady.i.j. the size smaller than 0.6 x 0.6m2,then all the movement Hence,by leveraging the MMSE estimator,we can estimate can be controlled in this region.Suppose the tag array is the,,B}by finding the optimal set of parameters that located in the linear region,we first consider the relationship minimizes the squared errors for all pairs of tags. between the distance and the phase difference among tags for antenna Ar along the X-axis.For any two arbitrary tags (,y,B)"=arg min e'(,y,B), Ti and Tj in the tag array,if the distance between Ti and x,y,8 Tj along the X-axis,i.e.,Axi.j ri-xj,is smaller than half-wavelength A/2,according to the phase values of Ti e=∑∑（ad-△dP+(ad,-△d,P) 11
IEEE TRANSACTIONS ON MOBILE COMPUTING, VOL. XX, NO. XX, 2020 7 ܶଵ ȟݔଵଶ ܶଶ ܶଷ ܶସ ȟݕଵଶ ߚ ଷସݔȟ ߚ ȟݕଷସ ܶଵ ܶଶ ܶଷ ܶସ ߚ ݈ ݓ ݄ ߟ ܱ(ܲ) ܠܔ ܡܔ Fig. 8. Estimate the initial rotation angle of the object rigid body first. E.g., the matrix T = xi,s − xa yi,s − ya in Eq. (10) depends on the relative positions of tags in the tag array, i.e., the topology and rotation state of the tag array. Hence, since the rotation state of the tag array depends on both initial and subsequent rotation angles of the tag array, we need to estimate the initial rotation angle first. Without loss of generality, we use the rectangle tag array as an example to illustrate our method to compute the initial rotation angle for the tag array. As shown in Fig. 8, we set the rotation center Pa at the center of the rectangle, and set two orthogonal polar axes lx and ly according to the X and Y -axis in the global coordinate system. Each tag can be regarded as a particle point in the coordinate system, e.g., T1 ∼ T4. The initial rotation angle of the tag array, i.e., the angle between lx and PaT1, is β. Topology Matrix T. According to the rotation angle β, we can depict the matrix T = xi,s − xa yi,s − ya for each tag Ti . E.g., for the rectangle tag array in Fig.8, let kT1T4k = kT2T3k = h, kT1T2k = kT3T4k = w, kPaTik = l, 6 T1PaT2 = η, these parameters can be regarded as constants, as long as the topology of the tag array is fixed. Then, according to the relative positions of the tags in the rectangle tag array, [x1 − xa, y1 − ya] T = [l cos β, lsin β] T , [x2 − xa, y2 − ya] T = [l cos (β + η), lsin (β + η)]T , [x3 − xa, y3 − ya] T = [−l cos (β), −lsin (β)]T , [x4 − xa, y4 − ya] T = [−l cos (β + η), −lsin (β + η)]T . (11) As aforementioned in Section 4.2, the distances between tag pairs are linear/non-linear to the phase differences between tag pairs in the linear region/non-linear region, respectively. In the following, we provide solutions to figure out the rotation angle β for the linear region and the vast region (including both the linear region and non-linear region), respectively. Linear Region: The linear region usually has an area of 0.6 × 0.6m2 in the intersection of the central beam regions for two orthogonally deployed antennas. If we set the operation area in this linear region, e.g., a tabletop with the size smaller than 0.6 × 0.6m2 , then all the movement can be controlled in this region. Suppose the tag array is located in the linear region, we first consider the relationship between the distance and the phase difference among tags for antenna Ax along the X-axis. For any two arbitrary tags Ti and Tj in the tag array, if the distance between Ti and Tj along the X-axis, i.e., ∆xi,j = xi − xj , is smaller than half-wavelength λ/2, according to the phase values of Ti and Tj collected from antenna Ax, the estimation of ∆xi,j , i.e., ∆xbi,j can be computed as follows: ∆xbi,j = (θx,i − θx,j ) × λ/(4π), |θx,i − θx,j | < π (θx,i − θx,j − 2π) × λ/(4π), θx,i − θx,j > π (θx,i − θx,j + 2π) × λ/(4π). θx,i − θx,j < −π Similarly, we can compute the estimated distance between Ti and Tj along the Y -axis, i.e., ∆ybi,j . Meanwhile, according to Eq. (11), given the rotation angle β, we can also compute the theoretical value for ∆xi,j = xi−xj and ∆yi,j = yi−yj . Therefore, we leverage the Minimum Mean Square Error (MMSE) estimator to estimate the rotation angle β, we define the mean squared error e by enumerating the squared errors for all pairs of tags, that is, e(β) = Xn i=1 Xn j=1 (∆xi,j (β) − ∆xbi,j ) 2 + (∆yi,j (β) − ∆ybi,j ) 2 . We then compute the optimal value of β that achieves the minimal mean squared error for e(β): β ∗ = arg min β e(β). Vast Region: The vast region usually covers a large area of 4 × 4m2 in the intersection of the scanning areas for two orthogonally deployed antennas. As the vast region includes both the linear region and non-linear region, we thus consider the solution for the non-linear region, as it is a more generalized solution for both regions. Suppose the tag array is located at the non-linear region, the distance between the tags are not linear to the phase differences between the tags, which depends on the exact positions of the corresponding tags. Therefore, it is essential to estimate both the position and rotation state of the tag array in the non-linear region. For the tag array, we set the position of the rotation center as (x, y), and set the rotation angle as β. Then, different from the solution in the linear region, in regard to the non-linear region, for any tag Ti in the tag array, we estimate the Euclidean distance dbi between tag Ti and the corresponding antenna (Ax or Ay), and compute the difference of the Euclidean distance ∆dbi,j = dbi − dbj for any two arbitrary tags Ti and Tj . Specifically, for antenna Ax, the difference of the Euclidean distance ∆dbx,i,j can be computed as follows: ∆dbx,i,j = (θx,i − θx,j ) × λ/(4π), |θx,i − θx,j | < π (θx,i − θx,j − 2π) × λ/(4π), θx,i − θx,j > π (θx,i − θx,j + 2π) × λ/(4π). θx,i − θx,j < −π Similarly, we can also compute the difference of the Euclidean distance ∆dby,i,j for antenna Ay. Meanwhile, according to the position (x, y) and rotation angle of the tag array β, we can also compute the theoretical value for the difference of the Euclidean distance ∆dx,i,j and ∆dy,i,j . Hence, by leveraging the MMSE estimator, we can estimate the {x, y, β} by finding the optimal set of parameters that minimizes the squared errors for all pairs of tags. (x, y, β) ∗ = arg min x,y,β e 0 (x, y, β), e 0 = Xn i=1 Xn j=1 (∆dx,i,j − ∆dbx,i,j ) 2 + (∆dy,i,j − ∆dby,i,j ) 2