WO2017000563A1 - 智能设备的实时定位方法及系统、手机运动姿态的判定方法 - Google Patents

智能设备的实时定位方法及系统、手机运动姿态的判定方法 Download PDF

Info

Publication number
WO2017000563A1
WO2017000563A1 PCT/CN2016/073895 CN2016073895W WO2017000563A1 WO 2017000563 A1 WO2017000563 A1 WO 2017000563A1 CN 2016073895 W CN2016073895 W CN 2016073895W WO 2017000563 A1 WO2017000563 A1 WO 2017000563A1
Authority
WO
WIPO (PCT)
Prior art keywords
acceleration
smart device
axis
real
time
Prior art date
Application number
PCT/CN2016/073895
Other languages
English (en)
French (fr)
Inventor
赵然
张弢
王周红
高民
Original Assignee
广州市香港科大霍英东研究院
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
Priority claimed from CN201510376966.XA external-priority patent/CN104964685B/zh
Priority claimed from CN201610051156.1A external-priority patent/CN105758404B/zh
Application filed by 广州市香港科大霍英东研究院 filed Critical 广州市香港科大霍英东研究院
Priority to US15/504,025 priority Critical patent/US10641608B2/en
Publication of WO2017000563A1 publication Critical patent/WO2017000563A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/183—Compensation of inertial measurements, e.g. for temperature effects
    • G01C21/188—Compensation of inertial measurements, e.g. for temperature effects for accumulated errors, e.g. by coupling inertial systems with absolute positioning systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C22/00—Measuring distance traversed on the ground by vehicles, persons, animals or other moving solid bodies, e.g. using odometers, using pedometers
    • G01C22/006—Pedometers
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P15/00—Measuring acceleration; Measuring deceleration; Measuring shock, i.e. sudden change of acceleration
    • G01P15/18—Measuring acceleration; Measuring deceleration; Measuring shock, i.e. sudden change of acceleration in two or more dimensions

Definitions

  • the invention relates to the field of real-time positioning, in particular to a real-time positioning method and system for a smart device, and a method for determining a mobile gesture.
  • Portable smart devices are becoming more and more popular and indispensable in the daily lives of modern people.
  • various MEMS sensors are integrated into portable smart devices. Therefore, many techniques for using pedestrians for navigation and tracking have been developed.
  • One method combines a dead reckoning algorithm with a zero speed correction technique to achieve this. It obtains the current position of the device by integrating the acceleration data to obtain the motion speed and motion direction of the smart device.
  • the zero-speed correction technique is combined with the estimation process to reduce the cumulative error generated during the integration process.
  • the integration of the acceleration sensor data produces a very large cumulative error.
  • the zero-speed correction technique is used to reduce the cumulative error, the accuracy of the data is very limited due to the high randomness and complexity of human motion.
  • the direction of motion of the smart device is determined by the orientation during its movement.
  • the method requires the user to keep a certain coordinate axis of the device pointing in the moving direction during the movement, and calculate the moving direction of the device according to the gravity acceleration and the magnetic field direction.
  • the displacement data of the smart device is obtained, and the real-time positioning result is displayed.
  • this method imposes certain requirements on the user's usage method. If the user does not maintain the specified orientation of the smart device during the motion, the positioning result may cause a large error.
  • the third method detects the deceleration phase in each step period by analyzing the acceleration waveform generated by the user while walking, and obtains an instantaneous acceleration in the deceleration phase according to the waveform feature, and projects the acceleration on the horizontal plane to determine the motion of the device. direction.
  • the purpose of the embodiments of the present invention is to provide a real-time positioning method and system for an intelligent device, which can effectively improve the accuracy of positioning.
  • an embodiment of the present invention provides a real-time positioning method for a smart device, which includes the following steps:
  • S1 Calculating the X, Y, and Z axis accelerations in the world coordinate system when the smart device moves according to the data collected by the acceleration sensor, the gravity acceleration sensor, and the magnetic field sensor in the smart device;
  • the X, Y, and Z axis acceleration sequences in the world coordinate system are sequentially intercepted and obtained according to a preset intercepting period.
  • the intercepting period includes 6 to 8 step periods, and each of the step periods refers to carrying The average step size of the user of the smart device;
  • the current Z-axis taken period taken acceleration sequence fast Fourier transform main frequency f step, using passband [f step -0.5Hz, f step + 0.5Hz] bandpass filter, respectively Filtering the X, Y, and Z axis acceleration sequences intercepted during the current interception period;
  • S is the displacement of the smart device in the current interception period
  • N is the number of user steps of the current interception period
  • V move is the moving direction of the smart device in the current interception period
  • d is the average step size of the user.
  • the method further comprises the steps of:
  • step S7 Connect the current position obtained in step S6 to the position recorded in the last interception period, thereby obtaining the motion track of the smart device and displaying it.
  • the method further comprises the steps of:
  • the preset minimum step interval is 0.5 seconds.
  • the preset step threshold is 0.5 m/s 2 .
  • the preset intercept period is 4 seconds; and/or, the acceleration sensor has an acquisition frequency of 50 Hz.
  • the step S1 specifically includes:
  • the rotation matrix is:
  • R is the transformation matrix
  • the azimuth, elevation, and roll angles are respectively represented and satisfy the formula (2):
  • G X , G Y and G Z respectively represent the X, Y and Z axis gravitational accelerations in the intelligent device coordinate system acquired by the gravity acceleration sensor
  • B X , B Y and B Z respectively represent the smart devices collected by the magnetic field sensor in real time.
  • the step S5 specifically includes:
  • the filtered X-axis acceleration sequence and the Y-axis acceleration sequence form a horizontal acceleration matrix
  • the step S52 specifically includes:
  • the feature vector corresponding to the calculated maximum feature value is used as the first principal component vector.
  • the step S53 specifically includes:
  • the horizontal acceleration matrix is multiplied by the first principal component vector to obtain the horizontal acceleration sequence after dimensionality reduction.
  • the step S54 specifically includes:
  • the embodiment of the invention provides a real-time positioning system for a smart device, including:
  • the triaxial acceleration calculation module is configured to calculate the X, Y, and Z axis accelerations in the world coordinate system when the smart device moves according to the real-time data collected by the acceleration sensor, the gravity acceleration sensor, and the magnetic field sensor in the smart device;
  • the triaxial acceleration sequence intercepting module is configured to sequentially intercept and obtain the X, Y, and Z axis acceleration sequences in the world coordinate system according to a preset intercepting period; wherein the intercepting period includes 6 to 8 step periods, each The step period refers to an average step size of a user who carries the smart device;
  • a Z-axis taken the current cycle sequence of acceleration taken fast Fourier transform main frequency f step, using passband [f step -0.5Hz, f step + 0.5Hz] bandpass
  • the filter filters the X, Y, and Z axis acceleration sequences intercepted in the current interception period respectively;
  • a step counting module configured to determine each peak in the time domain waveform diagram formed by the filtered Z-axis acceleration sequence that is greater than a preset step threshold as a step, thereby obtaining a user step number of the current intercepting period;
  • the motion direction calculation module calculates the filtered X-axis acceleration sequence and the Y-axis acceleration sequence by using a principal component analysis algorithm to obtain a motion direction of the smart device in the current interception period;
  • the position calculation and recording module calculates the displacement of the smart device in the current interception period by using formula (1), and obtains and records the current position in combination with the position recorded in the previous interception period:
  • S is the displacement of the smart device in the current interception period
  • N is the number of user steps of the current interception period
  • V move is the moving direction of the smart device in the current interception period
  • d is the average step size of the user.
  • the system further includes:
  • the motion track display module is configured to connect the current position obtained by the position calculation and recording module with the position recorded by the last interception period, thereby obtaining a motion track of the smart device and displaying the motion track.
  • the system further includes:
  • the step quantity correction module is configured to: in the time domain waveform diagram formed by the filtered Z-axis acceleration sequence, if there is a time interval between two adjacent peaks of the M pair greater than the preset step threshold, the time interval between the two adjacent peaks is less than the preset minimum step At intervals, the number of user steps is reduced by M, where M ⁇ 1.
  • the preset minimum step interval is 0.5 seconds.
  • the preset step threshold is 0.5 m/s 2 .
  • the preset intercept period is 4 seconds; and/or, the acceleration sensor has an acquisition frequency of 50 Hz.
  • the three-axis acceleration calculation module specifically includes:
  • the rotation matrix calculation unit is configured to calculate a rotation matrix between the smart device coordinate system and the world coordinate system according to the data in the intelligent device coordinate system acquired by the gravity acceleration sensor and the magnetic field sensor in real time;
  • the three-axis acceleration calculation unit is configured to multiply the X, Y, and Z-axis accelerations in the smart device coordinate system acquired by the acceleration sensor in real time, and obtain the X, Y in the world coordinate system when the smart device moves. , Z-axis acceleration.
  • the rotation matrix is:
  • R is the transformation matrix
  • the azimuth, elevation, and roll angles are respectively represented and satisfy the formula (2):
  • G X , G Y and G Z respectively represent the X, Y and Z axis gravitational accelerations in the intelligent device coordinate system acquired by the gravity acceleration sensor
  • B X , B Y and B Z respectively represent the smart devices collected by the magnetic field sensor in real time.
  • the motion direction calculation module specifically includes:
  • a horizontal acceleration matrix forming unit configured to form the filtered X-axis acceleration sequence and the Y-axis acceleration sequence into a horizontal acceleration matrix
  • a first principal component vector calculation unit obtains a first principal component vector of the horizontal acceleration matrix by using a principal component analysis algorithm
  • a horizontal acceleration sequence calculation unit configured to perform dimensionality reduction on the first principal component vector to obtain a horizontal acceleration sequence
  • a motion direction correcting unit configured to correct the first principal component vector according to a phase difference between the horizontal acceleration sequence and the filtered Z-axis acceleration sequence, thereby obtaining a motion direction of the smart device.
  • the first principal component vector calculation unit specifically includes:
  • a covariance matrix calculation unit for calculating a covariance matrix of the horizontal acceleration matrix according to formula (3):
  • a maximum eigenvector calculation unit for calculating the eigenvector v of the covariance matrix C and its corresponding eigenvalue ⁇ according to formula (4):
  • the feature vector corresponding to the calculated maximum feature value is used as the first principal component vector.
  • the horizontal acceleration sequence calculation unit is specifically configured to:
  • the horizontal acceleration matrix is multiplied by the first principal component vector to obtain the horizontal acceleration sequence after dimensionality reduction.
  • the motion direction correcting unit is specifically configured to:
  • the real-time positioning method and system of the smart device disclosed in the present invention sequentially collects and calculates the world coordinates from the acceleration sensor, the gravity acceleration sensor and the magnetic field sensor according to a preset intercepting period (time window).
  • the X, Y, and Z-axis accelerations are intercepted to obtain the X, Y, and Z-axis acceleration sequences under the current interception period, and then the fast Fourier transform is performed based on the Z-axis acceleration sequence intercepted by the current interception period to obtain the main frequency f step, using passband [f step -0.5Hz, f step + 0.5Hz] bandpass filter, respectively, taken in the current period taken X, Y, Z axis acceleration sequence filter, the filtered
  • Each of the peaks in the time domain waveform diagram formed by the Z-axis acceleration sequence that is greater than the preset step threshold is determined as a step, thereby obtaining the number of user steps of the current interception period, and using the principal component analysis algorithm to filter the X-
  • the Y-axis acceleration sequence is calculated to obtain the motion direction of the smart device in the current interception period, and then combined with the number of user steps and the smart device
  • the direction of motion is calculated to obtain the displacement and current position of the smart device during the current interception period.
  • the preset intercept period time window
  • the set intercept period makes it more accurate to determine the main frequency of the Z-axis acceleration sequence by using the fast Fourier transform in the subsequent steps. And ensure the real-time performance of the step.
  • the present invention does not employ a bandpass filter with a fixed passband, but obtains a main frequency by performing a fast Fourier transform on the Z-axis acceleration sequence intercepted for the current intercept period (because the acceleration data has the largest amplitude on the Z-axis).
  • f step and according to the main frequency f step , the passband filter with [f step -0.5Hz, f step +0.5Hz] is designed to respectively perform the X, Y and Z axis acceleration sequences intercepted in the current interception period. Filtering, therefore, the passband of the bandpass filter used in the present invention is automatically adjusted according to the acquired Z-axis acceleration data, and the interference data can be filtered more effectively than the existing fixed passband bandpass filter. .
  • the present invention also adds the condition that the step determination threshold is Each peak in the time domain waveform diagram formed by the filtered Z-axis acceleration sequence is greater than the preset step threshold, and is determined as a step, thereby improving the accuracy of the user's step calculation, thereby improving the accuracy of the real-time positioning of the smart device. .
  • an object of the present invention is to provide a method for determining a mobile phone's motion posture.
  • the technical solution adopted by the present invention to solve the technical problem thereof is:
  • a method for determining a mobile gesture of a mobile phone includes:
  • the current position of the mobile phone is obtained by combining the step count calculation.
  • step S1 includes:
  • step S2 includes:
  • the preset time sliding window refers to a time window in which the time length is one step period.
  • step S22 is specifically as follows:
  • the covariance matrix of the matrix is calculated, and then the eigenvectors of the covariance matrix and their corresponding eigenvalues are calculated, and finally Obtain the feature vector corresponding to the maximum eigenvalue of the covariance matrix as the walking direction vector of the mobile phone at the current time.
  • step S3 includes:
  • step S3 the step of performing dimension reduction processing on the acceleration in the horizontal direction in the vector direction is specifically:
  • the acceleration data of the X-axis and the Y-axis in the horizontal direction in one step period is formed into a sampling matrix, and the sampling matrix is multiplied by the forward direction vector, and the obtained result is taken as the acceleration in the horizontal direction after the dimensionality reduction.
  • step S4 includes:
  • step S41 is specifically:
  • Step S42 which is specifically:
  • s current represents the current position of the mobile phone
  • s initial represents the initial position of the mobile phone
  • N represents the current total value of the step count
  • i is a natural number
  • s i represents the displacement of the mobile phone at the i-th step
  • L i represents the i-th a pace length
  • v gi represents the i-th walking pace heading vector.
  • the conversion matrix in the step S11 is:
  • R represents the transformation matrix
  • Azimuth, elevation, and roll angle are respectively indicated, and the following formula is satisfied:
  • g represents the gravitational acceleration
  • G x , G y , and G z represent the real-time data of the X-axis, the Y-axis, and the Z-axis respectively in the cell phone coordinate system acquired by the gravitational accelerometer
  • B x , B y , and B z respectively Represents real-time data of the X-axis, Y-axis, and Z-axis in the cell phone coordinate system acquired by the magnetic field sensor.
  • the method for determining the movement posture of the mobile phone comprises: S1, real-time data collected by the accelerometer, the gravity accelerometer and the magnetic field sensor, and calculating the real-time in the world coordinate system when the mobile phone is in motion.
  • S3, calculating the current The forward direction vector of the step, and the horizontal direction acceleration is subjected to dimensionality reduction processing in the vector direction; S4, after correcting the forward direction vector of the current step, the current position of the mobile phone is obtained by combining the step count calculation.
  • the method comprehensively utilizes the data collected by the accelerometer, the gravity accelerometer and the magnetic field sensor built in the smart phone to determine the movement posture of the mobile phone under the user walking state, and processes and compares the acceleration data of the mobile phone in the horizontal and vertical directions.
  • the motion posture of the mobile phone in the walking state of the user is calculated, the cumulative error is effectively avoided, the accuracy is high, the measurement stability is strong, the user is not required to hold the handshake in a fixed posture, and has good applicability for different users.
  • FIG. 1 is a schematic flowchart of a real-time positioning method of a smart device according to Embodiment 1 of the present invention.
  • FIG. 2 is a schematic flow chart of step S1 in FIG. 1.
  • Embodiment 3 is a three-dimensional rotation diagram from the smart device coordinate system to the world coordinate system in Embodiment 1 of the present invention.
  • FIG. 4 is a frequency domain diagram of a Z-axis acceleration sequence taken in the current interception period in Embodiment 1 of the present invention.
  • Fig. 5 is a time-domain waveform diagram of the filtered Z-axis acceleration sequence in the first embodiment of the present invention, and shows the set step determination threshold.
  • FIG. 6 is a schematic flow chart of step S5 in FIG. 1.
  • Fig. 7 is a scatter distribution diagram of the filtered X-axis acceleration sequence and the Y-axis acceleration sequence in the horizontal plane in the first embodiment of the present invention.
  • Figure 8 is a comparison of time domain waveforms of a horizontal acceleration sequence and a filtered Z-axis acceleration sequence in Embodiment 1 of the present invention.
  • FIG. 9 is a schematic flowchart of a real-time positioning method of a smart device according to Embodiment 2 of the present invention.
  • FIG. 10 is a schematic flowchart of a real-time positioning method of a smart device according to Embodiment 3 of the present invention.
  • FIG. 11 is a structural block diagram of a real-time positioning system of a smart device according to Embodiment 4 of the present invention.
  • FIG. 12 is a block diagram showing the structure of the three-axis acceleration calculation module of FIG.
  • Figure 13 is a block diagram showing the structure of the motion direction calculation module of Figure 11.
  • FIG. 14 is a structural block diagram of a real-time positioning system of a smart device according to Embodiment 5 of the present invention.
  • FIG. 15 is a structural block diagram of a real-time positioning system of a smart device according to Embodiment 6 of the present invention.
  • FIG. 1 is a schematic flowchart diagram of a real-time positioning method of a smart device according to an embodiment of the present invention. The method includes steps S1 to S6:
  • S1 Calculate the X, Y, and Z axis accelerations in the world coordinate system when the smart device moves according to the data collected by the acceleration sensor, the gravity acceleration sensor, and the magnetic field sensor in the smart device.
  • step S1 specifically includes steps S11 to S12:
  • FIG. 3 is a schematic diagram of three-dimensional rotation from a smart device coordinate system to a world coordinate system.
  • the coordinate axis definition of the smart device coordinate system referenced in the built-in acceleration sensor of the smart device is as shown in FIG. 3: the X axis is parallel to the screen and points to the right side of the screen, the Y axis is parallel to the screen and points to the top of the screen, and the Z axis is perpendicular to the screen and Point to the screen.
  • the gravity acceleration sensor is first used to determine the direction perpendicular to the horizontal direction, and then the magnetic field sensor is used to determine the magnetic north direction. Then calculate the value of the triaxial acceleration in the world coordinate system.
  • the three-axis definition of the world coordinate system is shown in Figure 1. 102: The three axes of X, Y, and Z point to the east, north, and sky directions, respectively.
  • the direction of the vertical horizontal direction is determined by the gravitational acceleration
  • the direction of the magnetic north pole is determined by the magnetic flux.
  • the azimuth can be calculated
  • the pitch angle ⁇ , the roll angle ⁇ , and the three-dimensional rotation matrix R according to the following formula:
  • R is the transformation matrix
  • the azimuth, elevation, and roll angles are respectively represented and satisfy the formula (2):
  • G X , G Y and G Z respectively represent the X, Y and Z axis gravitational accelerations in the intelligent device coordinate system acquired by the gravity acceleration sensor
  • B X , B Y and B Z respectively represent the smart devices collected by the magnetic field sensor in real time.
  • the X, Y and Z axis accelerations of the intelligent device coordinate system acquired by the acceleration sensor in real time are:
  • the X, Y, and Z axis accelerations of the obtained smart device in the world coordinate system are:
  • ACC x represents the X-axis acceleration in the world coordinate system
  • ACC y represents the Y-axis acceleration in the world coordinate system, that is, ACC x and ACC y are the horizontal accelerations in the world coordinate system
  • ACC x represents the world coordinate system.
  • the lower Z-axis acceleration is the acceleration in the vertical direction in the world coordinate system.
  • the X, Y, and Z axis acceleration sequences in the world coordinate system are sequentially intercepted and obtained according to a preset intercepting period.
  • the intercepting period includes 6 to 8 step periods, and each of the step periods refers to carrying The average step size of the user of the smart device.
  • the intercept period (i.e., the time sliding window) is set to 4 seconds.
  • the acquisition frequency of the acceleration sensor of the present invention is 50 Hz
  • the length L of the time sliding window is 200.
  • the time window of this length may include about 6-8 step periods, and each of the step periods refers to an average step size of the user carrying the smart device, that is, each average step time.
  • the acquisition of the length L of the time sliding window can improve the accuracy when determining the main frequency by using the fast Fourier transform in the subsequent steps, and the real-time performance of the step counting is ensured.
  • an empty array of length 600 (3-axis acceleration, 200-time sliding window length) can be designed.
  • the acceleration sensor fills the array in the order of xyz axis every time data is acquired and rotated, and the empty array is filled. Subsequent processing and restarting an empty array at the same time for data recording of the next segment.
  • the current interception is required.
  • the periodic X, Y, and Z axis acceleration sequences are processed for subsequent processing. First, the X, Y, and Z axis acceleration sequences of the current interception period are filtered.
  • the Z-axis acceleration sequence of the current interception period is extracted from the above empty array and subjected to fast Fourier transform, and the amplitude of the obtained frequency response map at the step frequency is obvious. Higher than other frequency bands, as shown in Figure 4.
  • the frequency amplitude at 1.855 Hz is significantly higher than other frequencies, and the step frequency (ie, the main frequency f step ) of the set of data acquisition is determined to be 1.855 Hz.
  • the obtained step frequency amplitude will be used for the filter design.
  • a pass band design [f step -0.5Hz, f step + 0.5Hz] of the band-pass filter, respectively, taken the current cycle, Y, Z axis acceleration sequence filter.
  • a lower delay IIR filter is selected, and a more stable FIR can also be used. Since the filtering will produce a phase delay, filtering the X, Y, and Z-axis acceleration sequences of the current interception period by the same filter filter can ensure the consistency of the filtered data.
  • the bandpass filter designed by the present invention filters the X, Y, and Z axis acceleration sequences of the current interception period, and obtains the X, Y, and Z axis acceleration sequences filtered by the current intercept period:
  • fACC x represents the X axis after the current taken periodic filtered acceleration sequence
  • fACC y represents the Y axis after the current taken periodic filtered acceleration sequences, i.e. fACC x and fACC y acceleration sequence currently interception period filtering in the horizontal direction
  • fACC x represents the Z-axis acceleration sequence filtered by the current intercept period, that is, the acceleration sequence (vertical acceleration sequence) in the vertical direction after the current intercept period is filtered.
  • S4 Determine each peak of the time domain waveform diagram formed by the filtered Z-axis acceleration sequence that is greater than a preset step threshold as a step, thereby obtaining the number of user steps of the current interception period.
  • the filtered Z-axis acceleration sequence fACC z which is the vertical acceleration sequence ACC V , determines the threshold according to the set pace, performs step detection and records the number of user steps.
  • a band pass filter is designed according to the step frequency to filter the environmental noise in the original acceleration data.
  • the waveform of the filtered vertical acceleration sequence ACC V exhibits a waveform similar to a sine wave, judging its periodicity and referring to the step threshold to detect the number of user steps.
  • the waveform 301 of the filtered vertical acceleration sequence is relatively smooth, which greatly reduces the probability of false positives.
  • the slight fluctuation of the smart device also generates the fluctuation of the acceleration data and the corresponding peak and valley.
  • the embodiment adds the condition of the step determination threshold.
  • each peak above the step determination threshold 302 is determined to be a step, thereby obtaining the number N of user steps of the current intercept period.
  • the preset step threshold is 0.5 m/s 2 .
  • step S5 specifically includes steps S51-S54:
  • the filtered X-axis acceleration sequence and the Y-axis acceleration sequence form a horizontal acceleration matrix
  • the intercepted and filtered X-axis acceleration sequence fACC x and the Y-axis acceleration sequence fACC y are all sequences of length L. These two sequences are combined into a matrix of L*2, which is recorded as a horizontal acceleration matrix X.
  • the step S52 specifically includes steps S521 to S522:
  • the feature vector v max corresponding to the calculated maximum eigenvalue ⁇ max is used as the first principal component vector.
  • FIG. 7 is a scatter distribution diagram of the filtered X-axis acceleration sequence and the Y-axis acceleration sequence in the horizontal plane in Embodiment 1 of the present invention.
  • the acceleration sensor collects raw data collected during a westward walk.
  • the filtered horizontal acceleration sequence constitutes a set of two-dimensional scatter values 401 distributed near the X-axis.
  • the principal component analysis is performed on the horizontal acceleration sequence, and the obtained first principal component vector 402 points to the maximum direction of the variance of the original array.
  • the horizontal acceleration matrix X is multiplied by the first principal component vector v max to obtain the horizontal acceleration sequence ACC H after dimensionality reduction.
  • FIG. 8 is a time-domain waveform comparison diagram of the horizontal acceleration sequence and the filtered Z-axis acceleration sequence (ie, the vertical acceleration sequence) in the present embodiment.
  • the horizontal acceleration sequence 501 and the vertical acceleration sequence 502 have the same frequency and different phases.
  • the vector direction of the first principal component is corrected according to the phase difference between the two sets of acceleration sequence waveforms, and the motion direction of the smart device is obtained.
  • S is the displacement of the smart device in the current interception period
  • N is the number of user steps of the current interception period
  • V move is the moving direction of the smart device in the current interception period
  • d is the average step size of the user.
  • FIG. 9 is a schematic flowchart of a real-time positioning method of a smart device according to Embodiment 2 of the present invention.
  • the method is based on Embodiment 1, and further includes the steps of:
  • the step of correcting the calculated number of user steps is also added on the basis of the first embodiment.
  • the acceleration data is filtered by a bandpass filter, but the data still has irregular waveforms, such as two small peaks in one cycle.
  • a shortest step interval needs to be set to reduce the false positive. If the interval between the two peaks is less than the shortest step interval, it is only recorded as one step. Therefore, in order to effectively avoid the influence of the irregular waveform on the step counting result, in the time domain waveform diagram formed by the filtered Z-axis acceleration sequence, the time interval between each pair of adjacent peaks will be detected.
  • the time interval is shorter than the preset minimum step interval (eg, 0.5 seconds)
  • the pair of peaks is only recorded as one step.
  • N correct alternative to N Example Formula 1 (1) is implemented to Calculate the displacement S of the smart device during the current interception period, thereby improving the accuracy of the calculation result.
  • FIG. 10 is a schematic flowchart of a real-time positioning method of a smart device according to Embodiment 3 of the present invention.
  • the method is based on Embodiment 1, and further includes the steps of:
  • step S7 Connect the current position obtained in step S6 to the position recorded in the last interception period, thereby obtaining the motion track of the smart device and displaying it.
  • the motion track of the smart device in addition to being able to locate the current location of the smart device in real time, can also be Displayed to the user in real time.
  • the current position of each intercepting period can be sequentially obtained by step S6, and the current position of each intercepting period is connected with the position recorded in the last intercepting period, so that the motion track of the smart device can be obtained and passed through the display screen. Display to the user.
  • FIG. 11 is a structural block diagram of a real-time positioning system of a smart device according to Embodiment 4 of the present invention.
  • the system includes modules 41-46:
  • the triaxial acceleration calculation module 41 is configured to calculate, according to the real-time data collected by the acceleration sensor, the gravity acceleration sensor and the magnetic field sensor in the smart device, the X, Y, and Z axis accelerations in the world coordinate system when the smart device moves;
  • the triaxial acceleration sequence intercepting module 42 is configured to sequentially intercept and obtain the X, Y, and Z axis acceleration sequences in the world coordinate system according to a preset intercepting period; wherein the intercepting period includes 6 to 8 step periods, each The step period refers to an average step size of a user carrying the smart device;
  • the interception period (ie, the time sliding window) in the present invention is set to 4 seconds. If the acquisition frequency of the acceleration sensor of the present invention is 50 Hz, the length L of the time sliding window is 200.
  • the step counting module 44 is configured to determine each peak in the time domain waveform diagram formed by the filtered Z-axis acceleration sequence that is greater than a preset step threshold as a step, thereby obtaining the number of user steps of the current intercept period;
  • the motion direction calculation module 45 calculates the filtered X-axis acceleration sequence and the Y-axis acceleration sequence by using a principal component analysis algorithm to obtain a motion direction of the smart device in the current interception period;
  • the position calculation and recording module 46 calculates the displacement of the smart device in the current interception period by using formula (1), and obtains and records the current position according to the position recorded in the last interception period:
  • S is the displacement of the smart device in the current interception period
  • N is the number of user steps of the current interception period
  • V move is the moving direction of the smart device in the current interception period
  • d is the average step size of the user.
  • FIG. 12 is a structural block diagram of the three-axis acceleration calculation module in FIG.
  • the three-axis acceleration calculation module 41 specifically includes units 411 to 412:
  • the rotation matrix calculation unit 411 is configured to calculate a rotation matrix between the smart device coordinate system and the world coordinate system according to the data in the smart device coordinate system acquired by the gravity acceleration sensor and the magnetic field sensor in real time;
  • the triaxial acceleration calculation unit 412 is configured to multiply the X, Y, and Z axis accelerations in the smart device coordinate system acquired by the acceleration sensor in real time, and obtain the X in the world coordinate system when the smart device moves. Y, Z axis acceleration.
  • FIG. 13 is a structural block diagram of the motion direction calculation module in FIG.
  • the motion direction calculation module 45 specifically includes units 451 to 454:
  • the horizontal acceleration matrix constituting unit 451 is configured to form the filtered X-axis acceleration sequence and the Y-axis acceleration sequence into a horizontal acceleration matrix;
  • the first principal component vector calculation unit 452 obtains a first principal component vector of the horizontal acceleration matrix by using a principal component analysis algorithm
  • a horizontal acceleration sequence calculation unit 453, configured to perform dimensionality reduction on the first principal component vector to obtain a horizontal acceleration sequence
  • the horizontal acceleration matrix is multiplied by the first principal component vector to obtain the horizontal acceleration sequence after dimension reduction.
  • the motion direction correction unit 454 is configured to correct the first principal component vector according to a phase difference between the horizontal acceleration sequence and the filtered Z-axis acceleration sequence, thereby obtaining a motion direction of the smart device.
  • the moving direction of the smart device and the first main The component vector direction is the same; if the phase of the horizontal acceleration sequence is ahead of the filtered Z-axis acceleration sequence, then The direction of motion of the smart device is opposite to the direction of the first principal component vector.
  • the first principal component vector calculation unit 452 specifically includes:
  • a covariance matrix calculation unit for calculating a covariance matrix of the horizontal acceleration matrix according to formula (3):
  • a maximum eigenvector calculation unit for calculating the eigenvector v of the covariance matrix C and its corresponding eigenvalue ⁇ according to formula (4):
  • the feature vector corresponding to the calculated maximum feature value is used as the first principal component vector.
  • FIG. 14 is a structural block diagram of a real-time positioning system of a smart device according to Embodiment 5 of the present invention. Based on Embodiment 4, the system further includes:
  • the step quantity correction module 441 is configured to: in the time domain waveform diagram formed by the filtered Z-axis acceleration sequence, if there is a time interval between two adjacent peaks of the M pair being greater than the preset step threshold, the time interval between the two adjacent peaks is less than a preset minimum At the step interval, the number of user steps is reduced by M, where M ⁇ 1.
  • FIG. 15 is a structural block diagram of a real-time positioning system of a smart device according to Embodiment 6 of the present invention. Based on Embodiment 4, the system further includes:
  • the motion track display module 47 is configured to connect the current position obtained by the position calculation and recording module with the position recorded in the last interception period, thereby obtaining a motion track of the smart device and displaying the motion track.
  • the real-time positioning method and system of the smart device disclosed in the present invention sequentially collect and calculate the real-time coordinate system, the gravity acceleration sensor and the magnetic field sensor according to a preset intercepting period (time window).
  • the X, Y, and Z-axis accelerations are intercepted to obtain the X, Y, and Z-axis acceleration sequences in the current interception period, and then the fast-Fourier transform is performed based on the Z-axis acceleration sequence intercepted by the current interception period to obtain the main frequency f.
  • Step and use the bandpass filter with [f step -0.5Hz, f step +0.5Hz] to filter the X, Y and Z axis acceleration sequences intercepted in the current interception period respectively, and filter the Z
  • Each peak in the time domain waveform diagram formed by the axial acceleration sequence that is greater than the preset step threshold is determined as a step, thereby obtaining the number of user steps of the current interception period, and using the principal component analysis algorithm to filter the X-axis acceleration sequence and Y
  • the axis acceleration sequence is calculated to obtain the motion direction of the smart device in the current interception period, and then combined with the number of user steps and the motion of the smart device. Obtain the current displacement cycles taken intelligent devices and computing the current position.
  • the preset intercept period includes 6-8 step periods, and the set intercept period makes it more accurate to determine the main frequency of the Z-axis acceleration sequence by using the fast Fourier transform in the subsequent steps. And ensure the real-time performance of the step.
  • the present invention does not employ a bandpass filter with a fixed passband, but obtains a main frequency by performing a fast Fourier transform on the Z-axis acceleration sequence intercepted for the current intercept period (because the acceleration data has the largest amplitude on the Z-axis).
  • the passband filter with [f step -0.5Hz, f step +0.5Hz] is designed to respectively perform the X, Y and Z axis acceleration sequences intercepted in the current interception period. Filtering, therefore, the passband of the bandpass filter used in the present invention is automatically adjusted according to the acquired Z-axis acceleration data, and the interference data can be filtered more effectively than the existing fixed passband bandpass filter. .
  • the present invention also adds the condition that the step determination threshold is Each peak in the time domain waveform diagram formed by the filtered Z-axis acceleration sequence is greater than the preset step threshold, and is determined as a step, thereby improving the accuracy of the user's step calculation, thereby improving the accuracy of the real-time positioning of the smart device. .
  • the invention provides a method for determining a mobile gesture of a mobile phone, comprising:
  • the current position of the mobile phone is obtained by combining the step count calculation.
  • the step S1 includes:
  • the step S2 includes:
  • the preset time sliding window refers to a time window in which the time length is one step period.
  • step S22 is specifically:
  • the covariance matrix of the matrix is calculated, and then the eigenvectors of the covariance matrix and their corresponding eigenvalues are calculated, and finally Obtain the feature vector corresponding to the maximum eigenvalue of the covariance matrix as the walking direction vector of the mobile phone at the current time.
  • the step S3 includes:
  • the step of performing the dimension reduction processing on the acceleration in the horizontal direction in the vector direction in the step S3 is specifically:
  • the acceleration data of the X-axis and the Y-axis in the horizontal direction in one step period is formed into a sampling matrix, and the sampling matrix is multiplied by the forward direction vector, and the obtained result is taken as the acceleration in the horizontal direction after the dimensionality reduction.
  • the step S4 includes:
  • step S41 is specifically as follows:
  • Step S42 which is specifically:
  • s current represents the current position of the mobile phone
  • s initial represents the initial position of the mobile phone
  • N represents the current total value of the step count
  • i is a natural number
  • s i represents the displacement of the mobile phone at the i-th step
  • L i represents the i-th
  • v gi represents the walking forward direction vector of the i-th step.
  • the conversion matrix in the step S11 is:
  • R represents the transformation matrix
  • Azimuth, elevation, and roll angle are respectively indicated, and the following formula is satisfied:
  • g represents the gravitational acceleration
  • G x , G y , and G z represent the real-time data of the X-axis, the Y-axis, and the Z-axis respectively in the cell phone coordinate system acquired by the gravitational accelerometer
  • B x , B y , and B z respectively Represents real-time data of the X-axis, Y-axis, and Z-axis in the cell phone coordinate system acquired by the magnetic field sensor.
  • a method for determining a mobile gesture of a mobile phone includes:
  • S1 Calculate the real-time triaxial acceleration in the world coordinate system when the mobile phone moves according to the real-time data collected by the accelerometer, the gravity accelerometer and the magnetic field sensor, and specifically includes steps S11 and S12:
  • the three axes X W , Y W , and Z W of the world coordinate system are defined as the east, north, and sky directions, respectively, and the three axes X P , Y P , and Z P of the mobile coordinate system are respectively along the screen to the right and along the screen.
  • the vertical screen is outward
  • the azimuth, pitch angle, and roll angle are respectively represented, and the conversion matrix of the world coordinate system to the mobile phone coordinate system is:
  • the real-time data of the X-axis, Y-axis and Z-axis in the mobile phone coordinate system acquired by the gravity accelerometer is [G x G y G z ] T
  • the real-time data of the X-axis, Y-axis and Z-axis in the mobile phone coordinate system acquired by the magnetic field sensor is [B x B y B z ] T.
  • the real-time triaxial acceleration in the mobile phone coordinate system is:
  • the real-time triaxial acceleration in the world coordinate system is: ACCx and ACCy represent the acceleration in the horizontal direction, and ACCz represents the acceleration in the vertical direction.
  • step counting specifically including steps S21-S23:
  • fACCy constitutes an N*2 horizontal acceleration matrix X, where N is the number of samples in the preset time sliding window, and the value of N is determined by the pace period and the acceleration sampling frequency.
  • the preset time sliding window refers to a time window in which the time length is one step period.
  • the step of performing dimensionality reduction processing on the acceleration in the horizontal direction in the vector direction is specifically as follows:
  • the acceleration data fACCx and fACCy of the X-axis and the Y-axis in the horizontal direction in one step period constitute a sampling matrix P of M*2, M is the number of samples obtained in the period, and the sampling matrix P is compared with the forward direction vector. After multiplication, the obtained result is taken as the acceleration ACCh in the horizontal direction after dimensionality reduction;
  • the current position of the mobile phone is obtained by combining the step count calculation, including steps S41 and S42:
  • s current represents the current position of the mobile phone
  • s initial represents the initial position of the mobile phone
  • N represents the current total value of the step count
  • i is a natural number
  • s i represents the displacement of the mobile phone at the i-th step
  • L i represents the i-th
  • v gi represents the walking forward direction vector of the i-th step.
  • the method comprehensively utilizes the data collected by the accelerometer, the gravity accelerometer and the magnetic field sensor built in the smart phone to determine the movement posture of the mobile phone under the user walking state.
  • the method calculates the motion posture of the mobile phone under the walking state of the user, effectively avoids the accumulated error, has high precision, and has strong measurement stability.
  • the user holds the handshake in a fixed posture and has good applicability to different users, and can be used in various applications such as indoor positioning and walking navigation.
  • the computer program is programmed to execute the associated hardware, and the program can be stored in a computer readable storage medium, which, when executed, can include the flow of an embodiment of the methods described above.
  • the storage medium may be a magnetic disk, an optical disk, a read-only memory (ROM), or a random access memory (RAM).

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Telephone Function (AREA)
  • Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
  • User Interface Of Digital Computer (AREA)

Abstract

一种智能设备的实时定位方法，包括：根据预设截取周期依次截取并获得世界坐标系下的X、Y、Z轴加速度序列（S2）；对当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率f step，利用通带为[f step-0.5Hz,f step+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波（S3）；将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，得到当前截取周期的用户步伐数量（S4）；对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向（S5）；根据用户步伐数量和智能设备的运动方向计算得到智能设备的当前位置（S6）。

Description

智能设备的实时定位方法及系统、手机运动姿态的判定方法 技术领域
本发明涉及实时定位领域，尤其是涉及一种智能设备的实时定位方法及系统、一种手机运动姿态的判定方法。
背景技术
便携式智能设备在现代人的日常生活中越来越普及和不可或缺。随着硬件的制造技术不断提高，各种MEMS传感器被集成到便携式智能设备中。因此，很多利用智能设备进行步行导航和追踪的技术被开发出来。
一种方法结合了航位推算法和零速校正技术来实现这一目的。其通过对加速度数据的积分处理来获得智能设备的运动速度和运动方向，来进一步推算设备的当前位置。推算过程中结合零速校正技术来降低积分过程中产生的累计误差。然而由于智能设备中的MEMS传感器的精度相对较低，对加速度传感器数据的积分会产生很严重的累计误差。虽然零速校正技术被运用来降低累计误差，但是由于人运动的随机性和复杂性都很高，对数据精度的提升十分有限。
另一种方法中，智能设备的运动方向由其运动过程中的朝向来确定。该方法要求用户在运动过程中需要保持设备的某一坐标轴指向运动方向，并根据重力加速度和磁场方向来计算出设备的运动方向。结合计步算法，得到智能设备的位移数据，并显示其实时的定位结果。显然，该方法对用户的使用方法做出了一定的要求，如果用户在运动过程中没有保持智能设备的指定朝向，定位结果会产生很大的误差。
第三种方法通过分析用户步行时产生的加速度波形，检测出每个步伐周期中的减速相，并根据波形特征获得减速相中的一个瞬时加速度，将该加速度在水平面上投影来确定设备的运动方向。该方法虽然避免了累计误差，但是由于人运动的随机性和复杂性并不普遍适用。
发明内容
本发明实施例的目的是提供一种智能设备的实时定位方法及系统，能有效提高定位的精确度。
为实现上述目的，本发明实施例提供了一种智能设备的实时定位方法，其包括步骤：
S1、根据智能设备中的加速度传感器、重力加速度传感器和磁场传感器实时采集的数据，计算获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度；
S2、根据预设截取周期依次截取并获得所述世界坐标系下的X、Y、Z轴加速度序列；其中，所述截取周期包括6～8个步伐周期，每一所述步伐周期是指携带所述智能设备的用户的平均步长；
S3、对当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波；
S4、将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量；
S5、利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向；
S6、利用公式(1)计算当前截取周期内智能设备的位移，并结合上个截取周期所记录的位置得到并记录当前位置：
S＝N*d*Vmove；  公式(1)
其中，S为当前截取周期内智能设备的位移，N为当前截取周期的用户步伐数量，Vmove为当前截取周期内智能设备的运动方向，d为用户的平均步长。
作为上述方案的改进，所述方法还包括步骤：
S7、将步骤S6得到的当前位置与上个截取周期所记录的位置连接，从而得到智能设备的运动轨迹并显示。
作为上述方案的改进，所述方法还包括步骤：
S41、在滤波后的Z轴加速度序列所构成的时域波形图中，若存在M对大于预设步伐阈值的相邻两个波峰之间的时间间隔小于预设的最短步伐间隔，则将所述用户步伐数量减M，其中，M≥1。
作为上述方案的改进，所述预设的最短步伐间隔为0.5秒。
作为上述方案的改进，所述预设步伐阈值为0.5m/s2。
作为上述方案的改进，所述预设截取周期为4秒；和/或，所述加速度传感器的采集频率为50Hz。
作为上述方案的改进，所述步骤S1具体包括：
S11、根据重力加速度传感器和磁场传感器实时采集的智能设备坐标系下的数据，计算得到智能设备坐标系与世界坐标系之间的旋转矩阵；
S12、将加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴加速度与所述旋转矩阵相乘，从而获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度。
作为上述方案的改进，所述旋转矩阵为:
Figure PCTCN2016073895-appb-000001
其中，R为转换矩阵，
Figure PCTCN2016073895-appb-000002
分别表示方位角、俯仰角、滚动角，且满足公式(2)：
Figure PCTCN2016073895-appb-000003
  公式(2)
其中，GX、GY、GZ分别表示重力加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴重力加速度，BX、BY、BZ分别表示磁场传感器实时采集的智能设备坐标系下的X、Y、Z轴磁通量。
作为上述方案的改进，所述步骤S5具体包括：
S51、将滤波后的X轴加速度序列和Y轴加速度序列构成一个水平加速度矩阵；
S52、利用主成分分析算法得到所述水平加速度矩阵的第一主成分向量；
S53、将所述水平加速度矩阵在所述第一主成分向量上进行降维，得到水平加速度序列；
S54、根据所述水平加速度序列和滤波后的Z轴加速度序列之间的相位差对所述第一主成分向量进行校正，从而得到智能设备的运动方向。
作为上述方案的改进，所述步骤S52具体包括：
S521、根据公式(3)计算出所述水平加速度矩阵的协方差矩阵:
C＝E{(X-E[X])(X-E[X])T}  公式(3)
其中，X为水平加速度矩阵，C为协方差矩阵；
S522、根据公式(4)计算得到协方差矩阵C的特征向量v及其相应的特征值λ：
Cv＝λv  公式(4)
其中，将计算得到的最大特征值所对应的特征向量作为所述第一主成分向量。
作为上述方案的改进，所述步骤S53具体包括：
将所述水平加速度矩阵与所述第一主成分向量相乘，从而得到降维后的所述水平加速度序列。
作为上述方案的改进，所述步骤S54具体包括：
将所述水平加速度序列和滤波后的Z轴加速度序列进行对比，若滤波后的Z轴加速度序列的相位领先于所述水平加速度序列，则智能设备的运动方向与所述第一主成分向量方向相同；若所述水平加速度序列的相位领先于滤波后的Z轴加速度序列，则智能设备的运动方向与所述第一主成分向量方向相反。
本发明实施例对应提供了一种智能设备的实时定位系统，包括：
三轴加速度计算模块，用于根据智能设备中的加速度传感器、重力加速度传感器和磁场传感器实时采集的数据，计算获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度；
三轴加速度序列截取模块，用于根据预设截取周期依次截取并获得所述世界坐标系下的X、Y、Z轴加速度序列；其中，所述截取周期包括6～8个步伐周期，每一所述步伐周期是指携带所述智能设备的用户的平均步长；
滤波模块，用于对当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波；
步伐计数模块，用于将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量；
运动方向计算模块，利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向；
位置计算及记录模块，利用公式(1)计算当前截取周期内智能设备的位移，并结合上个截取周期所记录的位置得到并记录当前位置：
S＝N*d*Vmove；公式(1)
其中，S为当前截取周期内智能设备的位移，N为当前截取周期的用户步伐数量，Vmove为当前截取周期内智能设备的运动方向，d为用户的平均步长。
作为上述方案的改进，所述系统还包括：
运动轨迹显示模块，用于将位置计算及记录模块得到的当前位置与上个截取周期所记录的位置连接，从而得到智能设备的运动轨迹并显示。
作为上述方案的改进，所述系统还包括：
步伐数量纠正模块，用于在滤波后的Z轴加速度序列所构成的时域波形图中，若存在M对大于预设步伐阈值的相邻两个波峰之间的时间间隔小于预设的最短步伐间隔，则将所述用户步伐数量减M，其中，M≥1。
作为上述方案的改进，所述预设的最短步伐间隔为0.5秒。
作为上述方案的改进，所述预设步伐阈值为0.5m/s2。
作为上述方案的改进，所述预设截取周期为4秒；和/或，所述加速度传感器的采集频率为50Hz。
作为上述方案的改进，所述三轴加速度计算模块具体包括：
旋转矩阵计算单元，用于根据重力加速度传感器和磁场传感器实时采集的智能设备坐标系下的数据，计算得到智能设备坐标系与世界坐标系之间的旋转矩阵；
三轴加速度计算单元，用于将加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴加速度与所述旋转矩阵相乘，从而获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度。
作为上述方案的改进，所述旋转矩阵为:
Figure PCTCN2016073895-appb-000004
其中，R为转换矩阵，
Figure PCTCN2016073895-appb-000005
分别表示方位角、俯仰角、滚动角，且满足公式(2)：
Figure PCTCN2016073895-appb-000006
  公式(2)
其中，GX、GY、GZ分别表示重力加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴重力加速度，BX、BY、BZ分别表示磁场传感器实时采集的智能设备坐标系下的X、Y、Z轴磁通量。
作为上述方案的改进，所述运动方向计算模块具体包括：
水平加速度矩阵构成单元，用于将滤波后的X轴加速度序列和Y轴加速度序列构成一个水平加速度矩阵；
第一主成分向量计算单元，利用主成分分析算法得到所述水平加速度矩阵的第一主成分向量；
水平加速度序列计算单元，用于将所述水平加速度矩阵在所述第一主成分向量上进行降维，得到水平加速度序列；
运动方向校正单元，用于根据所述水平加速度序列和滤波后的Z轴加速度序列之间的相位差对所述第一主成分向量进行校正，从而得到智能设备的运动方向。
作为上述方案的改进，所述第一主成分向量计算单元具体包括：
协方差矩阵计算单元，用于根据公式(3)计算出所述水平加速度矩阵的协方差矩阵:
C＝E{(X-E[X])(X-E[X])T}  公式(3)
其中，X为水平加速度矩阵，C为协方差矩阵；
最大特征向量计算单元，用于根据公式(4)计算得到协方差矩阵C的特征向量v及其相应的特征值λ：
Cv＝λv  公式(4)
其中，将计算得到的最大特征值所对应的特征向量作为所述第一主成分向量。
作为上述方案的改进，所述水平加速度序列计算单元具体用于：
将所述水平加速度矩阵与所述第一主成分向量相乘，从而得到降维后的所述水平加速度序列。
作为上述方案的改进，所述运动方向校正单元具体用于：
将所述水平加速度序列和滤波后的Z轴加速度序列进行对比，若滤波后的Z轴加速度序列的相位领先于所述水平加速度序列，则智能设备的运动方向与所述第一主成分向量方向相同；若所述水平加速度序列的相位领先于滤波后的Z轴加速度序列，则智能设备的运动方向与所述第一主成分向量方向相反。
与现有技术相比，本发明公开的一种智能设备的实时定位方法及系统根据预设截取周期(时间划窗)依次对由加速度传感器、重力加速度传感器和磁场传感器实时采集并计算得到世界坐标系下的X、Y、Z轴加速度进行截取，从而得到当前截取周期下的X、Y、Z轴加速度序列，然后基于当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波，将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量，利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能 设备的运动方向，然后结合用户步伐数量和智能设备的运动方向计算得到当前截取周期内智能设备的位移和当前位置。在本发明中，预设截取周期(时间划窗)包括6～8个步伐周期，所设定的截取周期使得在后续步骤中利用快速傅里叶变换判断Z轴加速度序列的主频率时较为准确，且保证了计步的实时性。而且本发明并不是采用固定通带的带通滤波器，而是根据对当前截取周期所截取的Z轴加速度序列(因为加速度数据在Z轴上振幅最大)进行快速傅里叶变换后得到主频率fstep，并根据这个主频率fstep设计通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波，因此本发明采用的带通滤波器的通带会根据所采集的Z轴加速度数据做相应自动调整，相比现有的固定通带的带通滤波器能够更有效的对干扰数据进行过滤。另外，为了避免在非步行状态下，由于智能设备的轻微晃动所产生加速度数据的波动以及相应的波峰波谷而带来这部分数据导致的误判，本发明还加入了步伐判断阈值这一条件，即将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而提高用户步伐计算的准确度，进而提高了智能设备的实时定位的精确度。
为了解决上述的技术问题，本发明的目的是提供一种手机运动姿态的判定方法。本发明解决其技术问题所采用的技术方案是：
一种手机运动姿态的判定方法，包括：
S1、根据加速度计、重力加速度计和磁场传感器采集的实时数据，计算获得手机运动时在世界坐标系下的实时三轴加速度；
S2、使用带通滤波器对实时三轴加速度进行滤波，并分别对水平方向的加速度和垂直方向的加速度进行处理后，获得手机的步行方向向量并进行步伐计数；
S3、计算当前步伐的前进方向向量，并将水平方向的加速度在该向量方向上进行降维处理；
S4、对当前步伐的前进方向向量进行校正后，结合步伐计数计算获得手机当前的位置。
进一步，所述步骤S1，包括：
S11、根据重力加速度计和磁场传感器采集的手机坐标系下的实时数据，计算获得手机坐标系与世界坐标系之间的转换矩阵；
S12、将转换矩阵与加速度计采集的手机坐标系下的实时三轴加速度相乘，获得手机运动时在世界坐标系下的实时三轴加速度。
进一步，所述步骤S2，包括：
S21、使用通带为1～3Hz的带通滤波器对实时三轴加速度进行滤波；
S22、对预设时间滑窗内的水平方向的加速度做主成分分析后，获得手机的步行方向向量；
S23、根据垂直方向的加速度数据的周期性进行步伐计数；
所述预设时间滑窗是指时间长度为一个步伐周期的时间窗口。
进一步，所述步骤S22，其具体为：
将预设时间滑窗内的水平方向的X轴和Y轴的加速度数据构成水平加速度矩阵后，计算该矩阵的协方差矩阵，进而计算获得协方差矩阵的特征向量及其相应的特征值，最后获得协方差矩阵的最大特征值所对应的特征向量作为手机在当前时刻的步行方向向量。
进一步，所述步骤S3，包括：
S31、计算一个步伐周期内的步行方向向量的平均值作为当前步伐的前进方向向量；
S32、将水平方向的加速度在前进方向向量的方向上进行降维处理。
进一步，所述步骤S3中所述将水平方向的加速度在该向量方向上进行降维处理的步骤，其具体为：
将一个步伐周期内的水平方向的X轴和Y轴的加速度数据构成采样矩阵，进而将采样矩阵与前进方向向量相乘后，将获得的结果作为降维后的水平方向的加速度。
进一步，所述步骤S4，包括：
S41、计算垂直方向的加速度和降维后的水平方向的加速度之间的相位差，并根据该相位 差对当前步伐的前进方向向量进行校正，获得手机的步行前进方向向量；
S42、结合步行前进方向向量和步伐计数，获得手机当前的位置。
进一步，所述步骤S41，其具体为：
计算垂直方向的加速度和降维后的水平方向的加速度之间的相位差，并根据该相位差对当前步伐的前进方向向量进行校正，若判断垂直方向的加速度的相位领先于降维后的水平方向的加速度的相位，则将前进方向向量作为手机的步行前进方向向量，反之，将前进方向向量反相后作为手机的步行前进方向向量。
所述步骤S42，其具体为：
结合步行前进方向向量和步伐计数，根据下式，计算获得手机当前的位置：
Figure PCTCN2016073895-appb-000007
上式中，scurrent表示手机当前的位置，sinitial表示手机的初始位置，N表示当前的步伐计数总值，i为自然数，si表示手机在第i个步伐的位移，Li表示第i个步伐的长度，vgi表示第i个步伐的步行前进方向向量。
进一步，所述步骤S11中所述转换矩阵为：
Figure PCTCN2016073895-appb-000008
其中，R表示转换矩阵，
Figure PCTCN2016073895-appb-000009
分别表示方位角、俯仰角、滚动角，且满足下式：
Figure PCTCN2016073895-appb-000010
上式中，g表示重力加速度，Gx、Gy、Gz分别表示重力加速度计采集的手机坐标系下的X轴、Y轴和Z轴的实时数据，Bx、By、Bz分别表示磁场传感器采集的手机坐标系下的X轴、Y轴和Z轴的实时数据。
本发明的有益效果是：本发明的一种手机运动姿态的判定方法，包括：S1、根据加速度计、重力加速度计和磁场传感器采集的实时数据，计算获得手机运动时在世界坐标系下的实时三轴加速度；S2、使用带通滤波器对实时三轴加速度进行滤波，并分别对水平方向的加速度和垂直方向的加速度进行处理后，获得手机的步行方向向量并进行步伐计数；S3、计算当前步伐的前进方向向量，并将水平方向的加速度在该向量方向上进行降维处理；S4、对当前步伐的前进方向向量进行校正后，结合步伐计数计算获得手机当前的位置。本方法综合利用智能手机内置的加速度计、重力加速度计和磁场传感器采集的数据来进行使用者步行状态下手机运动姿态的判定，通过对手机在水平以及垂直方向上的加速度数据进行处理和比较，计算出在使用者步行状态下手机的运动姿态，有效的避免了累计误差，精确度高，测量稳定性强，不需要用户以固定姿态持握手机，并且对于不同的用户具有良好的适用性。
附图说明
图1是本发明实施例1中一种智能设备的实时定位方法的流程示意图。
图2是图1中的步骤S1的流程示意图。
图3是本发明实施例1中从智能设备坐标系到世界坐标系的三维旋转示意图。
图4是本发明实施例1中当前截取周期所截取的Z轴加速度序列的频域示意图。
图5是本发明实施例1中滤波后的Z轴加速度序列的时域波形图，并显示了设定的步伐判定阈值。
图6是图1中的步骤S5的流程示意图。
图7是本发明实施例1中滤波后的X轴加速度序列和Y轴加速度序列在水平面上的散点分布图。
图8是本发明实施例1中水平加速度序列和滤波后的Z轴加速度序列的时域波形对比图。
图9是本发明实施例2中一种智能设备的实时定位方法的流程示意图。
图10是本发明实施例3中一种智能设备的实时定位方法的流程示意图。
图11是本发明实施例4中一种智能设备的实时定位系统的结构框图。
图12是图11中的三轴加速度计算模块的结构框图。
图13是图11中的运动方向计算模块的结构框图。
图14是本发明实施例5中一种智能设备的实时定位系统的结构框图。
图15是本发明实施例6中一种智能设备的实时定位系统的结构框图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
参见图1，是本发明实施例提供的一种智能设备的实时定位方法的流程示意图。该方法包括步骤S1～S6：
S1、根据智能设备中的加速度传感器、重力加速度传感器和磁场传感器实时采集的数据，计算获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度。
其中，参考图2，所述步骤S1具体包括步骤S11～S12：
S11、根据重力加速度传感器和磁场传感器实时采集的智能设备坐标系下的数据，计算得到智能设备坐标系与世界坐标系之间的旋转矩阵。
其中，图3为从智能设备坐标系到世界坐标系的三维旋转示意图。智能设备内置加速度传感器中参照的智能设备坐标系的坐标轴定义如图3所示101：X轴平行于屏幕并指向屏幕右侧，Y轴平行于屏幕并指向屏幕上方，Z轴垂直于屏幕并指向屏幕外。为了将加速度数据转换到世界坐标系下，首先利用重力加速度传感器确定垂直于水平面向上的方向，再利用磁场传感器确定磁北方向。然后再计算出世界坐标系下三轴加速度的值。世界坐标系三轴定义如图1所示102：X、Y、Z三轴分别指向东、北、天三个方向。
具体的，由重力加速度确定垂直水平面向上的方向，由磁通量确定指向地磁北极的方向。根据这两个方向，可计算出方位角
Figure PCTCN2016073895-appb-000011
俯仰角θ、滚动角γ，并根据以下公式得到三维旋转矩阵R：
Figure PCTCN2016073895-appb-000012
其中，R为转换矩阵，
Figure PCTCN2016073895-appb-000013
分别表示方位角、俯仰角、滚动角，且满足公式(2)：
Figure PCTCN2016073895-appb-000014
  公式(2)
其中，GX、GY、GZ分别表示重力加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴重力加速度，BX、BY、BZ分别表示磁场传感器实时采集的智能设备坐标系下的X、Y、Z轴磁通量。将公式(2)计算得到的
Figure PCTCN2016073895-appb-000015
分别代入转换矩阵的公式中即可计算得到R。
S12、将加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴加速度与所述旋转矩阵相乘，从而获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度。
其中，设加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴加速度为:
Figure PCTCN2016073895-appb-000016
将智能设备坐标系下的X、Y、Z轴加速度与所述旋转矩阵R相乘后，得到的智能设备运动时在世界坐标系下的X、Y、Z轴加速度为：
Figure PCTCN2016073895-appb-000017
其中，ACCx表示世界坐标系下的X轴加速度，ACCy表示世界坐标系下的Y轴加速度，即ACCx和ACCy为世界坐标系下的水平方向的加速度；而ACCx表示世界坐标系下的Z轴加速度，即为世界坐标系下的垂直方向的加速度。
S2、根据预设截取周期依次截取并获得所述世界坐标系下的X、Y、Z轴加速度序列；其中，所述截取周期包括6～8个步伐周期，每一所述步伐周期是指携带所述智能设备的用户的平均步长。
经大量测试，人步行的频率最主要集中在1.5-2Hz。本发明中截取周期(即，时间滑窗)设定为4秒。若本发明的加速度传感器的采集频率为50Hz，那么，该时间滑窗的长度L为200。这个长度的时间窗会包括大约6～8个步伐周期，每一所述步伐周期是指携带所述智能设备的用户的平均步长，即每个平均步伐时间。该时间滑窗的长度L的获取，为后续步骤中利用快速傅里叶变换判断主频率时能提高准确度，且保证了计步的实时性。具体实施时，可设计一个长度为600的空数组(3轴加速度、200时间滑窗长度)，加速度传感器每次采集数据并旋转坐标系后按xyz轴顺序填充数组，该空数组填充满后进行后续处理并同时重启一个空数组进行下一段的数据记录。
S3、对当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波。
在利用时间长度为L的时间滑窗分别对所述世界坐标系下的X、Y、Z轴加速度进行截取，以获得当前截取周期的X、Y、Z轴加速度序列后，需要对该当前截取周期的X、Y、Z轴加速度序列进行后续处理。首先，就是对该当前截取周期的X、Y、Z轴加速度序列进行滤波处理。
由于加速度数据在Z轴上振幅最大，从上述空数组中提取出当前截取周期的Z轴加速度序列并对其做快速傅里叶变换，得到的频率响应图中在步伐频率上的幅值会明显高于其他频段，如图4所示。在1.855Hz处的频率幅值明显高于其他频率，即可确定该组数据采集时的 步伐频率(即，主频率fstep)为1.855Hz。获得的步伐频率幅值将用于滤波器设计。检测出其主频率fstep后，设计一个通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对该当前截取周期的X、Y、Z轴加速度序列进行滤波。本发明中选用了延迟较低的IIR滤波器，稳定性更高的FIR亦可。由于滤波会产生相位延迟，采用同一滤波器滤波分别对该当前截取周期的X、Y、Z轴加速度序列进行滤波可保证滤波后数据时间上的一致性。
通过本发明设计的带通滤波器对当前截取周期的X、Y、Z轴加速度序列进行滤波后，得到当前截取周期滤波后的X、Y、Z轴加速度序列：
Figure PCTCN2016073895-appb-000018
其中，fACCx表示当前截取周期滤波后的X轴加速度序列，fACCy表示当前截取周期滤波后的Y轴加速度序列，即fACCx和fACCy为当前截取周期滤波后的水平方向的加速度序列；而fACCx表示当前截取周期滤波后的Z轴加速度序列，即为当前截取周期滤波后的垂直方向的加速度序列(垂直加速度序列)。
S4、将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量。
滤波后的Z轴加速度序列fACCz，即为垂直加速度序列ACCV，根据设定的步伐判断阈值，进行步伐检测并记录用户步数。本实施例根据步伐频率设计一个带通滤波器用以过滤原始加速度数据中的环境噪声。滤波后的垂直加速度序列ACCV的波形呈现类似正弦波的波形，判断其周期性并参考步伐阈值来检测用户步伐数量。如图5所示，滤波后的垂直加速度序列的波形301相当平滑，很大程度上降低了步伐误判的几率。但是，在非步行状态下，智能设备的轻微晃动也会产生加速度数据的波动以及相应的波峰波谷，所以为了避免这部分数据导致的误判，本实施例加入了步伐判断阈值这一条件。在滤波后的Z轴加速度序列的时域波形图中，每个在步伐判定阈值302之上的波峰被才判定为一个步伐，从而得到当前截取周期的用户步伐数量N。
在本实施例中，所述预设步伐阈值为0.5m/s2。
S5、利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向。
具体的，参考图6，所述步骤S5具体包括步骤S51～S54：
S51、将滤波后的X轴加速度序列和Y轴加速度序列构成一个水平加速度矩阵；
其中，由于时间滑窗长度为L，被截取且滤波后的X轴加速度序列fACCx和Y轴加速度序列fACCy均为长度为L的序列。将这两个序列组成一个L*2的矩阵，记作为水平加速度矩阵X。
S52、利用主成分分析算法得到所述水平加速度矩阵的第一主成分向量；
其中，所述步骤S52具体包括步骤S521～S522：
S521、根据公式(3)计算出所述水平加速度矩阵的协方差矩阵:
C＝E{(X-E[X])(X-E[X])T}  公式(3)
其中，X为水平加速度矩阵，C为协方差矩阵；
S522、根据公式(4)计算得到协方差矩阵C的特征向量v及其相应的特征值λ：
Cv＝λv  公式(4)
其中，将计算得到的最大特征值λmax所对应的特征向量vmax作为所述第一主成分向量。
参考图7，图7是本发明实施例1中滤波后的X轴加速度序列和Y轴加速度序列在水平面上的散点分布图。在本实施例中，加速度传感器采集原始数据采集于一次向西步行过程中。而滤波后的水平加速度序列构成一组二维散点值401分布在X轴附近。对水平加速度序列做主成分分析，得到的第一主成分向量402即指向原数组方差最大方向。
S53、将所述水平加速度矩阵在所述第一主成分向量上进行降维，得到水平加速度序列；
具体的，将所述水平加速度矩阵X与所述第一主成分向量vmax相乘，从得到降维后的所述水平加速度序列ACCH。
S54、根据所述水平加速度序列和滤波后的Z轴加速度序列之间的相位差对所述第一主成分向量进行校正，从而得到智能设备的运动方向。
图8是本实施例中的水平加速度序列和滤波后的Z轴加速度序列(即垂直加速度序列)的时域波形对比图。在图8中，水平加速度序列501和垂直加速度序列502具有相同的频率和不同的相位。根据两组加速度序列波形之间的相位差校正第一主成分的向量方向，得到智能设备的运动方向。具体的，将所述水平加速度序列ACCH和滤波后的Z轴加速度序列fACCz(即垂直加速度序列ACCV)进行对比，若垂直加速度序列ACCV的相位领先于所述水平加速度序列，则智能设备的运动方向与所述第一主成分向量方向相同，即，vmove=vmax；若所述水平加速度序列的相位领先于垂直加速度序列ACCV，则智能设备的运动方向与所述第一主成分向量方向相反，即，vmove=-vmax。
S6、利用公式(1)计算当前截取周期内智能设备的位移，并结合上个截取周期所记录的位置得到并记录当前位置：
S＝N*d*Vmove；       公式(1)
其中，S为当前截取周期内智能设备的位移，N为当前截取周期的用户步伐数量，Vmove为当前截取周期内智能设备的运动方向，d为用户的平均步长。
结合步骤S4中获得的当前截取周期内的用户步伐数量N和步骤S5中获得的当前截取周期内的智能设备的运动方向Vmove，即可计算出当前截取周期内的位移，并通过上个截取周期所记录的智能设备位置S上得到并记录当前智能设备位置S当。假设用户的平均步长为d，则当前截取周期内的位移为S＝N·d·Vmove，当前智能设备位置S当=S上+S。
图9是本发明实施例2中一种智能设备的实时定位方法的流程示意图。在该实施例中，该方法在实施例1的基础上，还包括步骤：
S41、在滤波后的Z轴加速度序列所构成的时域波形图中，若存在M对大于预设步伐阈值的相邻两个波峰之间的时间间隔小于预设的最短步伐间隔，则将所述用户步伐数量减M，其中，M≥1。
即，在本实施例公开的智能设备的实时定位方法中，在实施例1的基础上还增加了对计算得到的用户步伐数量进行纠正的步骤。通过带通滤波器对加速度数据进行了滤后，但数据仍然会有不规则的波形，例如一个周期里有2个小波峰。这种情况需要设定一个最短步伐间隔来降低误判，如果两个波峰出现的时间间隔小于最短步伐间隔，则只记为一个步伐。因此，本实施例为了有效避免不规则波形对计步结果的影响，在滤波后的Z轴加速度序列所构成的时域波形图中，每对相邻波峰之间的时间间隔将被检测。如果该时间间隔比预设的最短步伐间隔(如0.5秒)更短，则这对波峰只记作一个步伐。在滤波后的Z轴加速度序列所构成的时域波形图中，若存在M对大于预设步伐阈值的相邻两个波峰之间的时间间隔小于预设的最短步伐间隔，则将计算得到的所述用户步伐数量N进行减M处理，其中，M≥1，即纠正后的用户步伐数量N纠=N-M。
在得到纠正后的用户步伐数量N纠后，则基于纠正后的用户步伐数量N纠来计算当前截取周期内智能设备的位移，即将N纠替换到实施例1中公式(1)中的N以计算当前截取周期内智能设备的位移S，从而提高计算结果的精确度。
图10是本发明实施例3中一种智能设备的实时定位方法的流程示意图。在该实施例中，该方法在实施例1的基础上，还包括步骤：
S7、将步骤S6得到的当前位置与上个截取周期所记录的位置连接，从而得到智能设备的运动轨迹并显示。
在本实施例中，除了能够实时定位智能设备的当前位置外，还能将智能设备的运动轨迹 实时显示给用户。具体的，通过将步骤S6可以依次得到每个截取周期的当前位置，将每个截取周期的当前位置与上个截取周期所记录的位置连接，即可得到智能设备的运动轨迹，并通过显示屏进行显示给用户。
图11是本发明实施例4中一种智能设备的实时定位系统的结构框图。该系统包括模块41～46：
三轴加速度计算模块41，用于根据智能设备中的加速度传感器、重力加速度传感器和磁场传感器实时采集的数据，计算获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度；
三轴加速度序列截取模块42，用于根据预设截取周期依次截取并获得所述世界坐标系下的X、Y、Z轴加速度序列；其中，所述截取周期包括6～8个步伐周期，每一所述步伐周期是指携带所述智能设备的用户的平均步长；
具体的，本发明中截取周期(即，时间滑窗)设定为4秒。若本发明的加速度传感器的采集频率为50Hz，那么，该时间滑窗的长度L为200。
滤波模块43，用于对当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波；
步伐计数模块44，用于将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量；
运动方向计算模块45，利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向；
位置计算及记录模块46，利用公式(1)计算当前截取周期内智能设备的位移，并结合上个截取周期所记录的位置得到并记录当前位置：
S＝N*d*Vmove；  公式(1)
其中，S为当前截取周期内智能设备的位移，N为当前截取周期的用户步伐数量，Vmove为当前截取周期内智能设备的运动方向，d为用户的平均步长。
具体的，图12是图11中的三轴加速度计算模块的结构框图。该三轴加速度计算模块41具体包括单元411～412：
旋转矩阵计算单元411，用于根据重力加速度传感器和磁场传感器实时采集的智能设备坐标系下的数据，计算得到智能设备坐标系与世界坐标系之间的旋转矩阵；
三轴加速度计算单元412，用于将加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴加速度与所述旋转矩阵相乘，从而获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度。
具体的，图13是图11中的运动方向计算模块的结构框图。该运动方向计算模块45具体包括单元451～454：
水平加速度矩阵构成单元451，用于将滤波后的X轴加速度序列和Y轴加速度序列构成一个水平加速度矩阵；
第一主成分向量计算单元452，利用主成分分析算法得到所述水平加速度矩阵的第一主成分向量；
水平加速度序列计算单元453，用于将所述水平加速度矩阵在所述第一主成分向量上进行降维，得到水平加速度序列；
具体的，将所述水平加速度矩阵与所述第一主成分向量相乘，从而得到降维后的所述水平加速度序列。
运动方向校正单元454，用于根据所述水平加速度序列和滤波后的Z轴加速度序列之间的相位差对所述第一主成分向量进行校正，从而得到智能设备的运动方向。
具体的，将所述水平加速度序列和滤波后的Z轴加速度序列进行对比，若滤波后的Z轴加速度序列的相位领先于所述水平加速度序列，则智能设备的运动方向与所述第一主成分向量方向相同；若所述水平加速度序列的相位领先于滤波后的Z轴加速度序列，则 智能设备的运动方向与所述第一主成分向量方向相反。
具体的，所述第一主成分向量计算单元452具体包括：
协方差矩阵计算单元，用于根据公式(3)计算出所述水平加速度矩阵的协方差矩阵:
C＝E{(X-E[X])(X-E[X])T}  公式(3)
其中，X为水平加速度矩阵，C为协方差矩阵；
最大特征向量计算单元，用于根据公式(4)计算得到协方差矩阵C的特征向量v及其相应的特征值λ：
Cv＝λv  公式(4)
其中，将计算得到的最大特征值所对应的特征向量作为所述第一主成分向量。
本实施例的智能设备的实时定位系统的工作原理及过程请参考实施例1，在此不再赘述。
图14是本发明实施例5中一种智能设备的实时定位系统的结构框图。该系统在实施例4的基础上，还包括：
步伐数量纠正模块441，用于在滤波后的Z轴加速度序列所构成的时域波形图中，若存在M对大于预设步伐阈值的相邻两个波峰之间的时间间隔小于预设的最短步伐间隔，则将所述用户步伐数量减M，其中，M≥1。
本实施例的智能设备的实时定位系统的工作原理及过程请参考实施例2，在此不再赘述。
图15是本发明实施例6中一种智能设备的实时定位系统的结构框图。该系统在实施例4的基础上，还包括：
运动轨迹显示模块47，用于将位置计算及记录模块得到的当前位置与上个截取周期所记录的位置连接，从而得到智能设备的运动轨迹并显示。
本实施例的智能设备的实时定位系统的工作原理及过程请参考实施例3，在此不再赘述。
综上所述，本发明公开的一种智能设备的实时定位方法及系统根据预设截取周期(时间划窗)依次对由加速度传感器、重力加速度传感器和磁场传感器实时采集并计算得到世界坐标系下的X、Y、Z轴加速度进行截取，从而得到当前截取周期下的X、Y、Z轴加速度序列，然后基于当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波，将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量，利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向，然后结合用户步伐数量和智能设备的运动方向计算得到当前截取周期内智能设备的位移和当前位置。在本发明中，预设截取周期(时间划窗)包括6～8个步伐周期，所设定的截取周期使得在后续步骤中利用快速傅里叶变换判断Z轴加速度序列的主频率时较为准确，且保证了计步的实时性。而且本发明并不是采用固定通带的带通滤波器，而是根据对当前截取周期所截取的Z轴加速度序列(因为加速度数据在Z轴上振幅最大)进行快速傅里叶变换后得到主频率fstep，并根据这个主频率fstep设计通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波，因此本发明采用的带通滤波器的通带会根据所采集的Z轴加速度数据做相应自动调整，相比现有的固定通带的带通滤波器能够更有效的对干扰数据进行过滤。另外，为了避免在非步行状态下，由于智能设备的轻微晃动所产生加速度数据的波动以及相应的波峰波谷而带来这部分数据导致的误判，本发明还加入了步伐判断阈值这一条件，即将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而提高用户步伐计算的准确度，进而提高了智能设备的实时定位的精确度。
本发明提供了一种手机运动姿态的判定方法，包括：
S1、根据加速度计、重力加速度计和磁场传感器采集的实时数据，计算获得手机运动时在世界坐标系下的实时三轴加速度；
S2、使用带通滤波器对实时三轴加速度进行滤波，并分别对水平方向的加速度和垂直方向的加速度进行处理后，获得手机的步行方向向量并进行步伐计数；
S3、计算当前步伐的前进方向向量，并将水平方向的加速度在该向量方向上进行降维处理；
S4、对当前步伐的前进方向向量进行校正后，结合步伐计数计算获得手机当前的位置。
进一步作为优选的实施方式，所述步骤S1，包括：
S11、根据重力加速度计和磁场传感器采集的手机坐标系下的实时数据，计算获得手机坐标系与世界坐标系之间的转换矩阵；
S12、将转换矩阵与加速度计采集的手机坐标系下的实时三轴加速度相乘，获得手机运动时在世界坐标系下的实时三轴加速度。
进一步作为优选的实施方式，所述步骤S2，包括：
S21、使用通带为1～3Hz的带通滤波器对实时三轴加速度进行滤波；
S22、对预设时间滑窗内的水平方向的加速度做主成分分析后，获得手机的步行方向向量；
S23、根据垂直方向的加速度数据的周期性进行步伐计数；
所述预设时间滑窗是指时间长度为一个步伐周期的时间窗口。
进一步作为优选的实施方式，所述步骤S22，其具体为：
将预设时间滑窗内的水平方向的X轴和Y轴的加速度数据构成水平加速度矩阵后，计算该矩阵的协方差矩阵，进而计算获得协方差矩阵的特征向量及其相应的特征值，最后获得协方差矩阵的最大特征值所对应的特征向量作为手机在当前时刻的步行方向向量。
进一步作为优选的实施方式，所述步骤S3，包括：
S31、计算一个步伐周期内的步行方向向量的平均值作为当前步伐的前进方向向量；
S32、将水平方向的加速度在前进方向向量的方向上进行降维处理。
进一步作为优选的实施方式，所述步骤S3中所述将水平方向的加速度在该向量方向上进行降维处理的步骤，其具体为：
将一个步伐周期内的水平方向的X轴和Y轴的加速度数据构成采样矩阵，进而将采样矩阵与前进方向向量相乘后，将获得的结果作为降维后的水平方向的加速度。
进一步作为优选的实施方式，所述步骤S4，包括：
S41、计算垂直方向的加速度和降维后的水平方向的加速度之间的相位差，并根据该相位差对当前步伐的前进方向向量进行校正，获得手机的步行前进方向向量；
S42、结合步行前进方向向量和步伐计数，获得手机当前的位置。
进一步作为优选的实施方式，所述步骤S41，其具体为：
计算垂直方向的加速度和降维后的水平方向的加速度之间的相位差，并根据该相位差对当前步伐的前进方向向量进行校正，若判断垂直方向的加速度的相位领先于降维后的水平方向的加速度的相位，则将前进方向向量作为手机的步行前进方向向量，反之，将前进方向向量反相后作为手机的步行前进方向向量。
所述步骤S42，其具体为：
结合步行前进方向向量和步伐计数，根据下式，计算获得手机当前的位置：
Figure PCTCN2016073895-appb-000019
上式中，scurrent表示手机当前的位置，sinitial表示手机的初始位置，N表示当前的步伐计数总值，i为自然数，si表示手机在第i个步伐的位移，Li表示第i个步伐的长度，vgi表示第i个步伐的步行前进方向向量。
进一步作为优选的实施方式，所述步骤S11中所述转换矩阵为：
Figure PCTCN2016073895-appb-000020
其中，R表示转换矩阵，
Figure PCTCN2016073895-appb-000021
分别表示方位角、俯仰角、滚动角，且满足下式：
Figure PCTCN2016073895-appb-000022
上式中，g表示重力加速度，Gx、Gy、Gz分别表示重力加速度计采集的手机坐标系下的X轴、Y轴和Z轴的实时数据，Bx、By、Bz分别表示磁场传感器采集的手机坐标系下的X轴、Y轴和Z轴的实时数据。
以下结合一具体实施例对本发明做详细说明。
一种手机运动姿态的判定方法，包括：
S1、根据加速度计、重力加速度计和磁场传感器采集的实时数据，计算获得手机运动时在世界坐标系下的实时三轴加速度，具体包括步骤S11和S12：
S11、根据重力加速度计和磁场传感器采集的手机坐标系下的实时数据，计算获得手机坐标系与世界坐标系之间的转换矩阵；
本步骤中，定义世界坐标系三轴XW、YW、ZW分别为东、北、天方向，手机坐标系三轴XP、YP、ZP分别为沿屏幕向右、沿屏幕向上、垂直屏幕向外方向，
Figure PCTCN2016073895-appb-000023
分别表示方位角、俯仰角、滚动角，则世界坐标系到手机坐标系的转换矩阵为：
Figure PCTCN2016073895-appb-000024
因为世界坐标系下，重力加速度的理论值应为[0 0 g]T，重力加速度计采集的手机坐标系下的X轴、Y轴和Z轴的实时数据为[Gx Gy Gz]T，则可得
Figure PCTCN2016073895-appb-000025
磁场传感器采集的手机坐标系下的X轴、Y轴和Z轴的实时数据为[Bx By Bz]T，将其投射到水平面后，XW、YW轴的磁通量分别为
Figure PCTCN2016073895-appb-000026
根据上述几个式子可获得
Figure PCTCN2016073895-appb-000027
满足下式：
Figure PCTCN2016073895-appb-000028
将
Figure PCTCN2016073895-appb-000029
的表达式带入转换矩阵R的表达式即可求解获得转换矩阵R。
S12、将转换矩阵R与加速度计采集的手机坐标系下的实时三轴加速度acc相乘，获得手机运动时在世界坐标系下的实时三轴加速度ACC＝R·acc。其中，手机坐标系下的实时三轴加 速度为：
Figure PCTCN2016073895-appb-000030
世界坐标系下的实时三轴加速度为：
Figure PCTCN2016073895-appb-000031
其中ACCx、ACCy表示水平方向的加速度，ACCz表示垂直方向的加速度。
S2、使用带通滤波器对实时三轴加速度进行滤波，并分别对水平方向的加速度和垂直方向的加速度进行处理后，获得手机的步行方向向量并进行步伐计数，具体包括步骤S21～S23：
S21、使用通带为1～3Hz的带通滤波器对实时三轴加速度ACCx、ACCy、ACCz进行滤波后，得到滤波后的世界坐标系下的实时三轴加速度，分别记为fACCx、fACCy、fACCz；
S22、对预设时间滑窗内的水平方向的加速度做主成分分析后，获得手机的步行方向向量，具体为：将预设时间滑窗内的水平方向的X轴和Y轴的加速度数据fACCx、fACCy构成一个N*2的水平加速度矩阵X，其中N为预设时间滑窗内的样本数量，N的值由步伐周期和加速度采样频率决定。然后，根据公式C＝E{(X-E[X])(X-E[X])T}计算该矩阵X的协方差矩阵C，进而根据公式Cv＝λv计算获得协方差矩阵C的特征向量v及其相应的特征值λ，最后获得协方差矩阵C的最大特征值λmax所对应的特征向量vmax作为手机在当前时刻的步行方向向量；
S23、根据垂直方向的加速度数据fACCz的周期性进行步伐计数，具体的，滤波后的垂直加速度fACCz呈现出类似正弦波的周期性波形，根据其波形的斜率变化趋势可得到其波峰，每一个波峰记为一个新的步伐，相邻两个波峰之间的时间间隔记为一个步伐周期；
所述预设时间滑窗是指时间长度为一个步伐周期的时间窗口。
S3、计算当前步伐的前进方向向量，并将水平方向的加速度在该向量方向上进行降维处理，具体包括S31和S32：
S31、计算一个步伐周期内的步行方向向量的平均值作为当前步伐的前进方向向量；
S32、将水平方向的加速度在前进方向向量的方向上进行降维处理。
本步骤中，将水平方向的加速度在该向量方向上进行降维处理的步骤，其具体为：
将一个步伐周期内的水平方向的X轴和Y轴的加速度数据fACCx和fACCy构成一个M*2的采样矩阵P，M为该周期内获得的样本数，进而将采样矩阵P与前进方向向量相乘后，将获得的结果作为降维后的水平方向的加速度ACCh；
S4、对当前步伐的前进方向向量进行校正后，结合步伐计数计算获得手机当前的位置，包括步骤S41和S42：
S41、计算垂直方向的加速度fACCz和降维后的水平方向的加速度ACCh之间的相位差，并根据该相位差对当前步伐的前进方向向量进行校正，若判断垂直方向的加速度fACCz的相位领先于降维后的水平方向的加速度ACCh的相位，则将前进方向向量作为手机的步行前进方向向量，反之，将前进方向向量反相后作为手机的步行前进方向向量；
S42、结合步行前进方向向量和步伐计数，根据下式，计算获得手机当前的位置：
Figure PCTCN2016073895-appb-000032
上式中，scurrent表示手机当前的位置，sinitial表示手机的初始位置，N表示当前的步伐计数总值，i为自然数，si表示手机在第i个步伐的位移，Li表示第i个步伐的长度，vgi表示第i个步伐的步行前进方向向量。
本方法综合利用智能手机内置的加速度计、重力加速度计和磁场传感器采集的数据来进行使用者步行状态下手机运动姿态的判定。本方法通过对手机在水平以及垂直方向上的加速度数据进行处理和比较，计算出在使用者步行状态下手机的运动姿态，有效的避免了累计误差，精确度高，测量稳定性强，不需要用户以固定姿态持握手机，并且对于不同的用户具有良好的适用性，可用于室内定位、步行导航等各类应用中。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计 算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。
以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (34)

  1. 一种智能设备的实时定位方法，其特征在于,包括步骤：
    S1、根据智能设备中的加速度传感器、重力加速度传感器和磁场传感器实时采集的数据，计算获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度；
    S2、根据预设截取周期依次截取并获得世界坐标系下的X、Y、Z轴加速度序列；其中，所述截取周期包括6～8个步伐周期，每一所述步伐周期是指携带所述智能设备的用户的平均步长；
    S3、对当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波；
    S4、将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量；
    S5、利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向；
    S6、利用公式(1)计算当前截取周期内智能设备的位移，并结合上个截取周期所记录的位置得到并记录当前位置：
    S＝N*d*Vmove；          公式(1)
    其中，S为当前截取周期内智能设备的位移，N为当前截取周期的用户步伐数量，Vmove为当前截取周期内智能设备的运动方向，d为用户的平均步长。
  2. 如权利要求1所述的智能设备的实时定位方法，其特征在于，所述方法还包括步骤：
    S7、将步骤S6得到的当前位置与上个截取周期所记录的位置连接，从而得到智能设备的运动轨迹并显示。
  3. 如权利要求1所述的智能设备的实时定位方法，其特征在于，所述方法还包括步骤：
    S41、在滤波后的Z轴加速度序列所构成的时域波形图中，若存在M对大于预设步伐阈值的相邻两个波峰之间的时间间隔小于预设的最短步伐间隔，则将所述用户步伐数量减M，其中，M≥1。
  4. 如权利要求3所述的智能设备的实时定位方法，其特征在于，所述预设的最短步伐间隔为0.5秒。
  5. 如权利要求1所述的智能设备的实时定位方法，其特征在于，所述预设步伐阈值为0.5m/s2。
  6. 如权利要求1所述的智能设备的实时定位方法，其特征在于，所述预设截取周期为4秒；和/或，所述加速度传感器的采集频率为50Hz。
  7. 如权利要求1所述的智能设备的实时定位方法，其特征在于，所述步骤S1具体包括：
    S11、根据重力加速度传感器和磁场传感器实时采集的智能设备坐标系下的数据，计算得到智能设备坐标系与世界坐标系之间的旋转矩阵；
    S12、将加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴加速度与所述旋转矩阵相乘，从而获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度。
  8. 如权利要求7所述的智能设备的实时定位方法，其特征在于，所述旋转矩阵为:
    Figure PCTCN2016073895-appb-100001
    其中，R为转换矩阵，
    Figure PCTCN2016073895-appb-100002
    θ、γ分别表示方位角、俯仰角、滚动角，且满足公式(2)：
    Figure PCTCN2016073895-appb-100003
            公式(2)
    其中，GX、GY、GZ分别表示重力加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴重力加速度，BX、BY、BZ分别表示磁场传感器实时采集的智能设备坐标系下的X、Y、Z轴磁通量。
  9. 如权利要求1所述的智能设备的实时定位方法，其特征在于，所述步骤S5具体包括：
    S51、将滤波后的X轴加速度序列和Y轴加速度序列构成一个水平加速度矩阵；
    S52、利用主成分分析算法得到所述水平加速度矩阵的第一主成分向量；
    S53、将所述水平加速度矩阵在所述第一主成分向量上进行降维，得到水平加速度序列；
    S54、根据所述水平加速度序列和滤波后的Z轴加速度序列之间的相位差对所述第一主成分向量进行校正，从而得到智能设备的运动方向。
  10. 如权利要求9所述的智能设备的实时定位方法，其特征在于，所述步骤S52具体包括：
    S521、根据公式(3)计算出所述水平加速度矩阵的协方差矩阵:
    C＝E{(X-E[x])(X-E[X])T}         公式(3)
    其中，X为水平加速度矩阵，C为协方差矩阵；
    S522、根据公式(4)计算得到协方差矩阵C的特征向量v及其相应的特征值λ：
    Cv＝λv         公式(4)
    其中，将计算得到的最大特征值所对应的特征向量作为所述第一主成分向量。
  11. 如权利要求9所述的智能设备的实时定位方法，其特征在于，所述步骤S53具体包括：
    将所述水平加速度矩阵与所述第一主成分向量相乘，从而得到降维后的所述水平加速度序列。
  12. 如权利要求9所述的智能设备的实时定位方法，其特征在于，所述步骤S54具体包括：
    将所述水平加速度序列和滤波后的Z轴加速度序列进行对比，若滤波后的Z轴加速度序列的相位领先于所述水平加速度序列，则智能设备的运动方向与所述第一主成分向量方向相同；若所述水平加速度序列的相位领先于滤波后的Z轴加速度序列，则智能设备的运动方向与所述第一主成分向量方向相反。
  13. 一种智能设备的实时定位系统，其特征在于,包括：
    三轴加速度计算模块，用于根据智能设备中的加速度传感器、重力加速度传感器和磁场传感器实时采集的数据，计算获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度；
    三轴加速度序列截取模块，用于根据预设截取周期依次截取并获得世界坐标系下的X、Y、Z轴加速度序列；其中，所述截取周期包括6～8个步伐周期，每一所述步伐周期是指携带所述智能设备的用户的平均步长；
    滤波模块，用于对当前截取周期所截取的Z轴加速度序列进行快速傅里叶变换后得到主频率fstep，并利用通带为[fstep-0.5Hz,fstep+0.5Hz]的带通滤波器分别对当前截取周期内所截取的X、Y、Z轴加速度序列进行滤波；
    步伐计数模块，用于将滤波后的Z轴加速度序列所构成的时域波形图中大于预设步伐阈值的每个波峰判定为一个步伐，从而得到当前截取周期的用户步伐数量；
    运动方向计算模块，利用主成分分析算法对滤波后的X轴加速度序列和Y轴加速度序列进行计算，得到当前截取周期内智能设备的运动方向；
    位置计算及记录模块，利用公式(1)计算当前截取周期内智能设备的位移，并结合上个截取周期所记录的位置得到并记录当前位置：
    S＝N*d*Vmove；        公式(1)
    其中，S为当前截取周期内智能设备的位移，N为当前截取周期的用户步伐数量，Vmove为当前截取周期内智能设备的运动方向，d为用户的平均步长。
  14. 如权利要求13所述的智能设备的实时定位系统，其特征在于，所述系统还包括：
    运动轨迹显示模块，用于将位置计算及记录模块得到的当前位置与上个截取周期所记录的位置连接，从而得到智能设备的运动轨迹并显示。
  15. 如权利要求13所述的智能设备的实时定位系统，其特征在于，所述系统还包括：
    步伐数量纠正模块，用于在滤波后的Z轴加速度序列所构成的时域波形图中，若存在M对大于预设步伐阈值的相邻两个波峰之间的时间间隔小于预设的最短步伐间隔，则将所述用户步伐数量减M，其中，M≥1。
  16. 如权利要求15所述的智能设备的实时定位系统，其特征在于，所述预设的最短步伐间隔为0.5秒。
  17. 如权利要求13所述的智能设备的实时定位系统，其特征在于，所述预设步伐阈值为0.5m/s2。
  18. 如权利要求13所述的智能设备的实时定位系统，其特征在于，所述预设截取周期为4秒；和/或，所述加速度传感器的采集频率为50Hz。
  19. 如权利要求13所述的智能设备的实时定位系统，其特征在于，所述三轴加速度计算模块具体包括：
    旋转矩阵计算单元，用于根据重力加速度传感器和磁场传感器实时采集的智能设备坐标系下的数据，计算得到智能设备坐标系与世界坐标系之间的旋转矩阵；
    三轴加速度计算单元，用于将加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴加速度与所述旋转矩阵相乘，从而获得智能设备运动时在世界坐标系下的X、Y、Z轴加速度。
  20. 如权利要求19所述的智能设备的实时定位系统，其特征在于，所述旋转矩阵为:
    Figure PCTCN2016073895-appb-100004
    其中，R为转换矩阵，
    Figure PCTCN2016073895-appb-100005
    θ、γ分别表示方位角、俯仰角、滚动角，且满足公式(2)：
    Figure PCTCN2016073895-appb-100006
           公式(2)
    其中，GX、GY、GZ分别表示重力加速度传感器实时采集的智能设备坐标系下的X、Y、Z轴重力加速度，BX、BY、BZ分别表示磁场传感器实时采集的智能设备坐标系下的X、Y、Z轴磁通量。
  21. 如权利要求13所述的智能设备的实时定位系统，其特征在于，所述运动方向计算模块具体包括：
    水平加速度矩阵构成单元，用于将滤波后的X轴加速度序列和Y轴加速度序列构成一个水平加速度矩阵；
    第一主成分向量计算单元，利用主成分分析算法得到所述水平加速度矩阵的第一主成分向量；
    水平加速度序列计算单元，用于将所述水平加速度矩阵在所述第一主成分向量上进行降维，得到水平加速度序列；
    运动方向校正单元，用于根据所述水平加速度序列和滤波后的Z轴加速度序列之间的相位差对所述第一主成分向量进行校正，从而得到智能设备的运动方向。
  22. 如权利要求21所述的智能设备的实时定位系统，其特征在于，所述第一主成分向量计算单元具体包括：
    协方差矩阵计算单元，用于根据公式(3)计算出所述水平加速度矩阵的协方差矩阵:
    C＝E{(X-E[X])(X-E[X])T}        公式(3)
    其中，X为水平加速度矩阵，C为协方差矩阵；
    最大特征向量计算单元，用于根据公式(4)计算得到协方差矩阵C的特征向量v及其相应的特征值λ：
    Cv＝λv        公式(4)
    其中，将计算得到的最大特征值所对应的特征向量作为所述第一主成分向量。
  23. 如权利要求21所述的智能设备的实时定位系统，其特征在于，所述水平加速度序列计算单元具体用于：
    将所述水平加速度矩阵与所述第一主成分向量相乘，从而得到降维后的所述水平加速度序列。
  24. 如权利要求21所述的智能设备的实时定位系统，其特征在于，所述运动方向校正单元具体用于：
    将所述水平加速度序列和滤波后的Z轴加速度序列进行对比，若滤波后的Z轴加速度序列的相位领先于所述水平加速度序列，则智能设备的运动方向与所述第一主成分向量方向相同；若所述水平加速度序列的相位领先于滤波后的Z轴加速度序列，则智能设备的运动方向与所述第一主成分向量方向相反。
  25. 一种手机运动姿态的判定方法，其特征在于，包括：
    S1、根据加速度计、重力加速度计和磁场传感器采集的实时数据，计算获得手机运动时在世界坐标系下的实时三轴加速度；
    S2、使用带通滤波器对实时三轴加速度进行滤波，并分别对水平方向的加速度和垂直方向的加速度进行处理后，获得手机的步行方向向量并进行步伐计数；
    S3、计算当前步伐的前进方向向量，并将水平方向的加速度在该向量方向上进行降维处理；
    S4、对当前步伐的前进方向向量进行校正后，结合步伐计数计算获得手机当前的位置。
  26. 根据权利要求25所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S1，包括：
    S11、根据重力加速度计和磁场传感器采集的手机坐标系下的实时数据，计算获得手机坐标系与世界坐标系之间的转换矩阵；
    S12、将转换矩阵与加速度计采集的手机坐标系下的实时三轴加速度相乘，获得手机运动时在世界坐标系下的实时三轴加速度。
  27. 根据权利要求25所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S2，包括：
    S21、使用通带为1～3Hz的带通滤波器对实时三轴加速度进行滤波；
    S22、对预设时间滑窗内的水平方向的加速度做主成分分析后，获得手机的步行方向向量；
    S23、根据垂直方向的加速度数据的周期性进行步伐计数；
    所述预设时间滑窗是指时间长度为一个步伐周期的时间窗口。
  28. 根据权利要求27所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S22，其具体为：
    将预设时间滑窗内的水平方向的X轴和Y轴的加速度数据构成水平加速度矩阵后，计算该矩阵的协方差矩阵，进而计算获得协方差矩阵的特征向量及其相应的特征值，最后获得协方差矩阵的最大特征值所对应的特征向量作为手机在当前时刻的步行方向向量。
  29. 根据权利要求25所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S3，包括：
    S31、计算一个步伐周期内的步行方向向量的平均值作为当前步伐的前进方向向量；
    S32、将水平方向的加速度在前进方向向量的方向上进行降维处理。
  30. 根据权利要求25所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S3中所述将水平方向的加速度在该向量方向上进行降维处理的步骤，其具体为：
    将一个步伐周期内的水平方向的X轴和Y轴的加速度数据构成采样矩阵，进而将采样矩阵与前进方向向量相乘后，将获得的结果作为降维后的水平方向的加速度。
  31. 根据权利要求25所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S4，包括：
    S41、计算垂直方向的加速度和降维后的水平方向的加速度之间的相位差，并根据该相位差对当前步伐的前进方向向量进行校正，获得手机的步行前进方向向量；
    S42、结合步行前进方向向量和步伐计数，获得手机当前的位置。
  32. 根据权利要求25所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S41， 其具体为：
    计算垂直方向的加速度和降维后的水平方向的加速度之间的相位差，并根据该相位差对当前步伐的前进方向向量进行校正，若判断垂直方向的加速度的相位领先于降维后的水平方向的加速度的相位，则将前进方向向量作为手机的步行前进方向向量，反之，将前进方向向量反相后作为手机的步行前进方向向量。
  33. 根据权利要求31所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S42，其具体为：
    结合步行前进方向向量和步伐计数，根据下式，计算获得手机当前的位置：
    Figure PCTCN2016073895-appb-100007
    上式中，scurrent表示手机当前的位置，sinitial表示手机的初始位置，N表示当前的步伐计数总值，i为自然数，si表示手机在第i个步伐的位移，Li表示第i个步伐的长度，vgi表示第i个步伐的步行前进方向向量。
  34. 根据权利要求26所述的一种手机运动姿态的判定方法，其特征在于，所述步骤S11中所述转换矩阵为：
    Figure PCTCN2016073895-appb-100008
    其中，R表示转换矩阵，
    Figure PCTCN2016073895-appb-100009
    θ、γ分别表示方位角、俯仰角、滚动角，且满足下式：
    Figure PCTCN2016073895-appb-100010
    上式中，g表示重力加速度，Gx、Gy、Gz分别表示重力加速度计采集的手机坐标系下的X轴、Y轴和Z轴的实时数据，Bx、By、Bz分别表示磁场传感器采集的手机坐标系下的X轴、Y轴和Z轴的实时数据。
PCT/CN2016/073895 2015-06-30 2016-02-16 智能设备的实时定位方法及系统、手机运动姿态的判定方法 WO2017000563A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US15/504,025 US10641608B2 (en) 2015-06-30 2016-02-16 Method and system for real-time positioning of smart device, and method for determining the motion gesture of mobile phone

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
CN201510376966.XA CN104964685B (zh) 2015-06-30 2015-06-30 一种手机运动姿态的判定方法
CN201510376966.X 2015-06-30
CN201610051156.1 2016-01-26
CN201610051156.1A CN105758404B (zh) 2016-01-26 2016-01-26 智能设备的实时定位方法及系统

Publications (1)

Publication Number Publication Date
WO2017000563A1 true WO2017000563A1 (zh) 2017-01-05

Family

ID=57607718

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2016/073895 WO2017000563A1 (zh) 2015-06-30 2016-02-16 智能设备的实时定位方法及系统、手机运动姿态的判定方法

Country Status (2)

Country Link
US (1) US10641608B2 (zh)
WO (1) WO2017000563A1 (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110595501A (zh) * 2019-10-09 2019-12-20 成都乐动信息技术有限公司 一种基于三轴传感器的跑步距离校正方法
CN111780779A (zh) * 2020-06-15 2020-10-16 Oppo（重庆）智能科技有限公司 数据处理方法、装置、设备及存储介质
CN114739412A (zh) * 2022-04-11 2022-07-12 河海大学 一种基于智能手机的行人步态实时检测方法及设备

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10342462B2 (en) * 2014-10-26 2019-07-09 David Martin Application of gait characteristics for mobile
US10564177B2 (en) * 2016-08-01 2020-02-18 Samsung Electronics Co., Ltd. System and method for user activity recognition using accelerometer
CN107966161B (zh) * 2017-11-09 2020-04-28 内蒙古大学 基于fft的步行检测方法
CN111902695B (zh) * 2018-08-08 2022-07-12 华为技术有限公司 一种获取用户运动轨迹的方法及终端
CN109103953B (zh) * 2018-08-23 2021-07-20 广州市香港科大霍英东研究院 一种电池组主动均衡控制方法、系统及装置
CN111435083A (zh) * 2019-01-11 2020-07-21 阿里巴巴集团控股有限公司 行人航迹推算方法、导航方法及装置、手持终端及介质
CN112161630B (zh) * 2020-10-12 2022-07-15 北京化工大学 适用于大型仓储系统的agv在线无碰撞路径规划方法
CN113938853B (zh) * 2021-09-26 2023-05-12 众虎物联网(广州)有限公司 一种室内导航的定位方法、装置及系统

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101900563A (zh) * 2009-05-28 2010-12-01 三星电子株式会社 使用终端的航向信息的步行导航设备和方法
US20110054787A1 (en) * 2009-08-27 2011-03-03 Apple Inc. Context Determination to Assist Location Determination Accuracy
CN102278988A (zh) * 2011-07-01 2011-12-14 微迈森惯性技术开发（北京）有限公司 一种步行定位方法与设备
US8380433B2 (en) * 2008-09-24 2013-02-19 Texas Instruments Incorporated Low-complexity tightly-coupled integration filter for sensor-assisted GNSS receiver
CN102944240A (zh) * 2012-11-02 2013-02-27 清华大学 一种基于智能手机的惯性导航系统及方法
CN103940442A (zh) * 2014-04-03 2014-07-23 深圳市宇恒互动科技开发有限公司 一种采用加速收敛算法的定位方法及装置
CN104964685A (zh) * 2015-06-30 2015-10-07 广州市香港科大霍英东研究院 一种手机运动姿态的判定方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110307213A1 (en) * 2006-07-10 2011-12-15 Yang Zhao System and method of sensing attitude and angular rate using a magnetic field sensor and accelerometer for portable electronic devices
US8462109B2 (en) * 2007-01-05 2013-06-11 Invensense, Inc. Controlling and accessing content using motion processing on mobile devices
AU2011374152B2 (en) * 2011-08-04 2016-07-21 Google Llc Moving direction determination with noisy signals from inertial navigation systems on mobile devices

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8380433B2 (en) * 2008-09-24 2013-02-19 Texas Instruments Incorporated Low-complexity tightly-coupled integration filter for sensor-assisted GNSS receiver
CN101900563A (zh) * 2009-05-28 2010-12-01 三星电子株式会社 使用终端的航向信息的步行导航设备和方法
US20110054787A1 (en) * 2009-08-27 2011-03-03 Apple Inc. Context Determination to Assist Location Determination Accuracy
CN102278988A (zh) * 2011-07-01 2011-12-14 微迈森惯性技术开发（北京）有限公司 一种步行定位方法与设备
CN102944240A (zh) * 2012-11-02 2013-02-27 清华大学 一种基于智能手机的惯性导航系统及方法
CN103940442A (zh) * 2014-04-03 2014-07-23 深圳市宇恒互动科技开发有限公司 一种采用加速收敛算法的定位方法及装置
CN104964685A (zh) * 2015-06-30 2015-10-07 广州市香港科大霍英东研究院 一种手机运动姿态的判定方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110595501A (zh) * 2019-10-09 2019-12-20 成都乐动信息技术有限公司 一种基于三轴传感器的跑步距离校正方法
CN110595501B (zh) * 2019-10-09 2022-10-04 成都乐动信息技术有限公司 一种基于三轴传感器的跑步距离校正方法
CN111780779A (zh) * 2020-06-15 2020-10-16 Oppo（重庆）智能科技有限公司 数据处理方法、装置、设备及存储介质
CN114739412A (zh) * 2022-04-11 2022-07-12 河海大学 一种基于智能手机的行人步态实时检测方法及设备
CN114739412B (zh) * 2022-04-11 2024-05-03 河海大学 一种基于智能手机的行人步态实时检测方法及设备

Also Published As

Publication number Publication date
US10641608B2 (en) 2020-05-05
US20170234686A1 (en) 2017-08-17

Similar Documents

Publication Publication Date Title
WO2017000563A1 (zh) 智能设备的实时定位方法及系统、手机运动姿态的判定方法
CN109579853B (zh) 基于bp神经网络的惯性导航室内定位方法
CN105758404B (zh) 智能设备的实时定位方法及系统
US10508920B2 (en) Pedestrian dead reckoning position tracker
CN104964685A (zh) 一种手机运动姿态的判定方法
Yu et al. Comparison of pedestrian tracking methods based on foot-and waist-mounted inertial sensors and handheld smartphones
CN109115207A (zh) 行人步行轨迹检测方法、装置及系统
CN107917707B (zh) 一种任意姿态下行人方向的确定方法、装置及电子设备
Liu et al. A novel heading estimation algorithm for pedestrian using a smartphone without attitude constraints
US10533874B2 (en) Inertial positioning and navigation device featuring a novel walk detection method
US20210093917A1 (en) Detecting outdoor walking workouts on a wearable device
WO2015051656A1 (zh) 一种移动智能终端及其计步方法、系统
US11833392B2 (en) Methods and systems for swimming performance analysis
US20180344217A1 (en) Fitness tracking for constrained-arm usage
CN111435083A (zh) 行人航迹推算方法、导航方法及装置、手持终端及介质
CN106767770A (zh) 一种基于便携智能设备的用户行走方向检测与追踪方法
CN110657802B (zh) 一种gps失效情况下的智能手环导航方法
CN114739412B (zh) 一种基于智能手机的行人步态实时检测方法及设备
CN116943130A (zh) 一种计数方法、装置、设备及存储介质
US20220218230A1 (en) System and method of detecting walking activity using waist-worn inertial sensors
Truong et al. Nonparametric regression-based step-length estimation for arm-swing walking using a smartphone
Mikov The multi-mode inertial tracking system for unconstrained indoor positioning
Wang et al. An improved PDR localization algorithm based on particle filter
Wang et al. Posture recognition and adaptive step detection based on hand-held terminal
Suksuganjana et al. Improved step detection with smartphone handheld mode recognition

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 16816928

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 16816928

Country of ref document: EP

Kind code of ref document: A1