2150 IEEE/ACM TRANSACTIONS ON NETWORKING,VOL.26.NO.5.OCTOBER 2018 Z axis "地 Z axis constant f System Bi Y axis 阳 gravitational Y axis acceleration g Fig.5.fi(t)consists of consistent acceleration f(t)and inconsistent acceleration f(t).the projection of fe(t)on the r.y and 2 axes are Fig.6.Deriving the global coordinate. fe(t)cos ai,fe(t)cos Bi,and fe(t)cos Yi.respectively. Theorem 1:For a fixed coordinate system relative to the A vector quantity v defined in the local coordinate system earth coordinate system,for any inconsistent acceleration f.(t) is equivalent to the vector v=C.v defined in the global caused by intra-body movements,the expected value of f(t) coordinate.The inverse transformation is v=C.v as the is approximately equal to 0 within a sufficiently large time inverse of a rotation matrix is equal to its transpose. interval [ts,te].That is, For the local coordinate of a device,assuming that we can extract a constant acceleration as a vector fe from the E':（训=te-t。J 1 E(t)dt≈0. (2) acceleration measurements and extract a constant gravitational acceleration as a vector g from the low pass filter (such as the Proof:Without loss of generality,we prove this theorem Butterworth filter [271).we can build the global coordinate and for z-axis.Our reasoning is applicable for other axes.Let compute the rotation matrix as follows.After we obtain the Avi.denote the change of the relative velocity between gravity vector g,we derive its opposite value and normalize device Di and the body after the intra-body movements this vector as we then set this vector to represent occurred during time interval [ts,te].According to the rela- the global Z-axis as it is in the opposite direction of the tion between velocity and acceleration,we have Avi.r= gravitational acceleration.We set the vector zo to represent f(t)dt.As body parts move back and forth,the ampli- to be perpendicular to the horizontal plane.After computing tudes of the intra-body movement are usually small and do the cross product y =g x fc,we obtain a vector y that is not change significantly in any direction.Thus,Avi.should perpendicular to the plane determined by the two distinct but be smaller than a constant threshold C.In general situations, intersecting lines corresponding to g and fe.We normalize this C<0.4m/s due to the back-and-forth moving property. vector as yo=y.Since the vector yo is on the horizontal Thus,let At =te-ta,for E[f'i.(t)],i.e.,the expected plane,we set this vector to represent the global Y-axis.After value of the acceleration f(t)during [ts,te],we have that,by computing the cross product x =g x y,we obtain a vector x that is orthogonal to the plane determined by the two EU(. 1 fiz(t)dt △i< C (3) distinct but intersecting lines corresponding to g and y.We △t △t normalize this vector as x=.Since the vector xo is on If the time interval At =te-ts is sufficiently larger than the horizontal plane,and it is orthogonal to the global Y-axis the value of C,then E[fi,z(t)]≈0，e.g,△t=2 Os and and Z-axis,we set it to represent the global X-axis.After we C-0.4ms,△t>C,then Ef(d≤%ms2=2× obtain the vectors xo,yo,and zo in the local coordinate,which 10-2m/s2≈0. correspond to the global X,Y,and Z-axes,we can derive the rotation matrix CT.Fig.6 illustrates the above process. Theorem 1 implies that,during a sufficiently large time interval,e.g.,10~20s,the expected value of the inconsistent acceleration can be negligible,as the inconsistent accelerations B.Maintain Space Synchronization cancel each other out during the back and forth movements As the orientation of a device may be continuously chang- This explains our observations that the mean value of the ing,we need to track and update the rotation matrix C over inconsistent acceleration is in the order of 10-5m/s2. time.For example,if the rotation matrix at time t is given by Ct,then the rotation matrix C+6t at time t+6t can be IV.SPACE SYNCHRONIZATION MODEL computed as the product of two matrices A.Achieve Space Synchronization C1+8t=CiAt.t+8t (4) The objective of space synchronization for mobile devices where At.t+6t is the rotation matrix relating the local coordi- is to align their local coordinate to their global coordinate nate at time t to the one at time t+ot.To estimate the matrix We use the direction cosine representation to quantify the orientation difference between the local and global coordi- At.t+6t,we use gyroscope measurements.According to the nates of each device.In the direction cosine representation, small angle approximation [12],if we use w(t),w(t),and the orientation of the local coordinate relative to the global w(t)to represent the rotation rate of small rotations about its coordinate system is specified by a 3 x 3 rotation matrix C. x,y and z axes between time t and t+ot,then where each column is a unit vector along one axis in the local rt+6t coordinate specified in terms of the global coordinate axes. At.t+6t exp( (t)dt) (5)2150 IEEE/ACM TRANSACTIONS ON NETWORKING, VOL. 26, NO. 5, OCTOBER 2018 Fig. 5. fi(t) consists of consistent acceleration fc(t) and inconsistent acceleration f i(t), the projection of fc(t) on the x, y and z axes are fc(t) cos αi, fc(t) cos βi, and fc(t) cos γi, respectively. Theorem 1: For a fixed coordinate system relative to the earth coordinate system, for any inconsistent acceleration f i(t) caused by intra-body movements, the expected value of f i(t) is approximately equal to 0 within a sufficiently large time interval [ts, te]. That is, E[f’i(t)] = 1 te − ts  te ts f i(t)dt ≈ 0. (2) Proof: Without loss of generality, we prove this theorem for x-axis. Our reasoning is applicable for other axes. Let Δvi,x denote the change of the relative velocity between device Di and the body after the intra-body movements occurred during time interval [ts, te]. According to the rela￾tion between velocity and acceleration, we have Δvi,x =  te ts f i,x(t)dt. As body parts move back and forth, the ampli￾tudes of the intra-body movement are usually small and do not change significantly in any direction. Thus, Δvi,x should be smaller than a constant threshold C. In general situations, C ≤ 0.4m/s due to the back-and-forth moving property. Thus, let Δt = te − ts, for E[f i,x(t)], i.e., the expected value of the acceleration f i,x(t) during [ts, te], we have E[f i,x(t)] = 1 Δt  te ts f i,x(t)dt = Δvi,x Δt < C Δt . (3) If the time interval Δt = te − ts is sufficiently larger than the value of C, then E[f i,x(t)] ≈ 0, e.g., Δt = 20s and C = 0.4m/s, Δt  C, then E[f i,x(t)] ≤ 0.4 20 m/s2 = 2 × 10−2m/s2 ≈ 0. Theorem 1 implies that, during a sufficiently large time interval, e.g., 10∼20s, the expected value of the inconsistent acceleration can be negligible, as the inconsistent accelerations cancel each other out during the back and forth movements. This explains our observations that the mean value of the inconsistent acceleration is in the order of 10−5m/s2. IV. SPACE SYNCHRONIZATION MODEL A. Achieve Space Synchronization The objective of space synchronization for mobile devices is to align their local coordinate to their global coordinate. We use the direction cosine representation to quantify the orientation difference between the local and global coordi￾nates of each device. In the direction cosine representation, the orientation of the local coordinate relative to the global coordinate system is specified by a 3 × 3 rotation matrix C, where each column is a unit vector along one axis in the local coordinate specified in terms of the global coordinate axes. Fig. 6. Deriving the global coordinate. A vector quantity vl defined in the local coordinate system is equivalent to the vector vg = C · vl defined in the global coordinate. The inverse transformation is vl = CT · vg as the inverse of a rotation matrix is equal to its transpose. For the local coordinate of a device, assuming that we can extract a constant acceleration as a vector fc from the acceleration measurements and extract a constant gravitational acceleration as a vector g from the low pass filter (such as the Butterworth filter [27]), we can build the global coordinate and compute the rotation matrix as follows. After we obtain the gravity vector g, we derive its opposite value and normalize this vector as zo = −g g , we then set this vector to represent the global Z-axis as it is in the opposite direction of the gravitational acceleration. We set the vector zo to represent to be perpendicular to the horizontal plane. After computing the cross product y = g × fc, we obtain a vector y that is perpendicular to the plane determined by the two distinct but intersecting lines corresponding to g and fc. We normalize this vector as yo = y y . Since the vector yo is on the horizontal plane, we set this vector to represent the global Y -axis. After that, by computing the cross product x = g × y, we obtain a vector x that is orthogonal to the plane determined by the two distinct but intersecting lines corresponding to g and y. We normalize this vector as xo = x x . Since the vector xo is on the horizontal plane, and it is orthogonal to the global Y -axis and Z-axis, we set it to represent the global X-axis. After we obtain the vectors xo, yo, and zo in the local coordinate, which correspond to the global X, Y , and Z-axes, we can derive the rotation matrix CT . Fig. 6 illustrates the above process. B. Maintain Space Synchronization As the orientation of a device may be continuously chang￾ing, we need to track and update the rotation matrix C over time. For example, if the rotation matrix at time t is given by Ct, then the rotation matrix Ct+δt at time t + δt can be computed as the product of two matrices Ct+δt = CtAt,t+δt (4) where At,t+δt is the rotation matrix relating the local coordi￾nate at time t to the one at time t+δt. To estimate the matrix At,t+δt, we use gyroscope measurements. According to the small angle approximation [12], if we use ωx(t), ωy(t), and ωz(t) to represent the rotation rate of small rotations about its x, y and z axes between time t and t + δt, then At,t+δt = exp( t+δt t Ω(t)dt) (5)