The paper describes the methodology and experimental results for revealing similarities in thermal dependencies of biases of accelerometers and gyroscopes from 250 inertial MEMS chips (MPU-9250). Temperature pro les were measured on an experimental setup with a Peltier element for temperature control. Classi cation of temperature curves was carried out with machine learning approach. A perfect sensor should not have thermal dependency at all. Thus, only sensors inside the clusters with smaller dependency (smaller total temperature slopes) might be pre-selected for production of high accuracy inertial navigation modules. It was found that no uni ed thermal pro le (\family" curve) exists for all sensors in a production batch. However, obviously, sensors might be grouped according to their parameters. Therefore, the temperature compensation pro les might be regressed for each group. 12 slope coe cients on 5 degrees temperature intervals from 0oC to +60oC were used as the features for the k-means++ clustering algorithm. The minimum number of clusters for all sensors to be well separated from each other by bias thermal pro les in our case is 6. It was found by applying the elbow method. For each cluster a regression curve can be obtained.
fabricated with high yield and no additional or subsequent assembly. The integration process may be implemented as hybrid integration using conventional wire bonding and ip-chips, or monolithic integration [11]. Monolithic integration o ers superior system integration performance to hybrid systems but at an overall higher e ort in upfront and non-recurring engineering (NRE) costs in terms of involved technology and processing [1].
Despite the many similarities between IC (integrated circuit) and MEMS fabrication, MEMS fabrication methods are signi cantly more complex, especially in case of production of multisensor modules [7,11]. Unlike an ordinary IC factory, which performs one or two standard processes, a MEMS factory normally performs a wide variety of processes.
The physical characteristics of the base material, accuracy of topology building over the wafer and the complexity of MEMS fabrication are the major factors, which have an impact on random and systematic errors of sensor readouts. The error mechanisms a ecting the accuracy of MEMS sensors originate from sensitivities and stability of dimensions of sensor parts and sensitivities and stability of electronics. Mechanical stress also could be a source of inaccuracies.
Under static conditions major part of errors can be considered, measured and excluded. However, in dynamics it is hard to separate uctuating zero o sets from the true physical signal.
As is typical of inertial sensors, thermal e ects are a primary driver of inaccuracies [13,14]. For example, the sti ness coe cient of beams, the damping ratio and other MEMS material parameters change with temperature and a ect the gyroscope. These phenomena are subject to analysis, modelling [8,21] and compensating [4,12,22].
Consequently, MEMS chips manufacturers denote large error intervals for produced sensors. For example, gyroscope zero-rate output (ZRO) variation over temperature range of 40oC + 85oC is declared as 30o=s [15].
Inertial sensors may be embedded in dead reckoning systems, which calculate the current position by using a previously determined position (also called x), and advancing that position based upon known or estimated linear or rotational speeds or accelerations over elapsed time and trajectory.
One possible way to improve the accuracy of the estimation positioning or orientation is to carry out periodical self-calibration of inertial modules. However, during long periods of continuous motion under conditions of varying temperature it is impossible to run autocalibration algorithms.
For this, developers use temperature compensation data obtained in preliminary tests in thermal chambers. The compensation functions might be presented in either \family" curves (tables) or unique pro les for the chips. Both dependencies are pre-measured for the manufactured inertial measurement unit on a factory side, not by a customer, which leads to both a limited visibility of the internal processes in the modules and a signi cant dependency from the module vendor. This is restricting the achievable performance as the chip vendor will typically choose a compromise between accuracy and computational complexity. If the algorithms were open, system developers could more exibly trade-o these two objectives and reach better accuracy with the same physical modules through more complex compensation algorithms.
So the goal is either to nd the \family" curve for all sensors in a batch, or to divide them into groups and to obtain separate regression curves for each group. In order to achieve this goal, we started this cluster analysis of thermal curves of MEMS sensors, which is - to our best knowledge - the rst public discussion of these aspects.
The remainder of the paper is organized as follows. In chapter 2 we give a short overview of the related works for the state of the art. In chapter 3, we describe our ongoing project, sensor chips used as devices under test (DUT) and our experimental setup and proposed processing algorithms. In Section 4 we present clustering groups of obtained bias thermal pro les and discuss possible correlation between them. Section 5 concludes the paper. 2
Machine learning is widely applied to separate signal spectra in groups [18]. However, the papers devoted to the statistical analysis of sensors compensation curves are rare, since researchers usually do not have statistically signi cant amount of sensors. One of the examples of applying polynomial regression for nding a \family" curve for MEMS barometric sensors is presented in [9].
One of the problems concerning clustering time series data or any other functional data is that each data point contains a noise component. Thus, additional data preprocessing is required. T. Tarpey investigated the curve clustering performance depending on the quality of how the curves are t to the data [20]. Another approach was demostrated by A. Antoniadis et al. [3]. They represent time series data on electricity consumption as a set of wavelets and then tried to group the processed data by means of K-centroid algorithm.
In practice, the clustering problem can be solved by various algorithms: Kmedoids, hierarchical clustering, density-based clustering, etc. [10].
It this paper the k-means++ algorithm was applied [2,5]. The algorithm is based on minimizing the within-cluster sum of squares. Euclidean distance is used as a metric in data space. We also implemented additional procedure of initialization of the cluster centers as presented in [6]. It is known that the kmeans++ algorithm converges to a local minimum. Therefore, if it is appropriate, the exact optimum might be found by the brute force. 3 3.1
Our ongoing project is devoted to the development of the high precision autonomous self-calibrating MEMS multisensor inertial module unit MIMU2.5, presented in Fig. 1 [17,19].
Each MIMU2.5 consists of ve 9DOF MPU-9250 chips (3D accelerometer, 3D wope, 3D magnetometer and temperature sensor) [15]. MPU-9250 chips are mounted on the inner surface of the top aluminium cover at di erent angles to each other.
The pilot production batch of MIMU2.5 modules consisted of 50 pieces. Thus, approximately 250 MPU-9250 chips can be used as devices under test (DUTs).
According to the datasheet of MPU-9250 [15], it contains two dice: a MPU6500 accelerometer and gyroscope produced by InvenSense and an AK8963 magnetometer of Asahi Kasei Microdevices Corporation. 3.2
To register temperature pro les for all axes of accelerometers and gyroscopes an experimental setup was constructed. It consists of a thermoelectric cooler (TEC, Peltier element) controlled by a microcontroller and a cooling system comprising a heat sink and a fan (Fig. 2a). The aluminium ange of MIMU module was attached to the cooling side of the TEC.
First, the inertial module was cooled to the minimum temperature of approximately 0oC, which depends on the ambient temperature. The achieved temperature was monitored by the sensors embedded in MPU9250 chips. Acceleration and rotation rate values along with the temperature data were registered by the external computer connected to the MIMU module via RS-232 interface.
Then, the polarity of the supplied voltage was changed to inverted and heating started. After reaching 60oC the experiment ended. 3.3
We propose the following method for clustering bias thermal pro les. First, thermal dependencies are registered as a set of measurement points. Then, the whole temperature range is divided into n intervals. In each interval the regressed slope coe cients for all sensors are determined. Thus each sensor curve is characterized by a set of n slope coe cients. Then, slope coe cients are standardized and used as features in machine learning. Finally, bias thermal pro les are distributed in a number of clusters according to their sets of n slope coe cients.
From the one hand, the number of intervals should not be large in order to improve classi cation. From the other, it should not be small, since the obtained piece-wise function should well describe the curved pro le.
In our investigation we chose n = 12, leaving 5 measurement points per interval for linear regression. The number of clusters was varied from 5 to 12.
Both accelerometer and gyroscope bias thermal pro les are handled similarly. Further in formal representation of the algorithm we consider measurements from a triaxial accelerometer.
Let I be the index set used for numbering accelerometers.
Consider a vector ai(t) = (aiX (t); aiY (t); aiZ (t)) consisting of temperature
pro les registered for X-, Y-, and Z-axes of i-th accelerometer, i 2 I. We assume
that
aiX (t) = aX +
aiY (t) = aY +
aiZ (t) = aZ +
iX (t) + iX
iY (t) + iY
iZ (t) + iZ ;
(
i(t) = ( iX (t); iY (t); iZ (t)) is the bias (zero o set) vector of the accelerometer i at a temperature t, i = ( iX ; iY ; iZ ) is the vector of normally distributed random variables with zero means, characterizing the noise of measurements.
The functions ai(t) are assumed di erentiable in the mean-square sense.
Let the temperature range T = [t1; t2] is divided into N = mn points with equal intervals: 0 = t1; 1 = 0 + h; : : : ; N = t2. The interval h is selected as follows: h = t2Nt1 , where m is the size of the temperature interval, n is the number of features for the subsequent clustering of sensors.
Thus we divide thermal pro les aiX (t); aiY (t); aiZ (t) into a certain number of intervals n, see Fig. 3 (n = 12).
F1
F 3 F2
F4
F11
F12
(
We will use a~iXj as the sample mean values of aiX (t) on every intervals [ j 1; j ) to estimate the expectation of the values of this function, j = 1; : : : ; N .
We will use the following method to estimate the derivatives of the expectation value of the function aiX (t).
For each k = 1; : : : ; n we construct linear regressions f ( ) = BiXk + CXk by sets of m points
The angular coe cients BiXk will serve as estimates for (E[aiX (t)])0 on the intervals [ (k 1)m; km) for k = 1; : : : ; n.
Before clustering, the resulting values of BiXk attributes for accelerometers need to be standardized.
Let BXk be the sample mean and iXk be the sample standard deviation for the numbers BiXk, i 2 I.
We introduce the standardized features:
DiXk =
BiXk
BiXk ;
iXk
i 2 I; k = 1; : : : ; n:
(
Standardization involves the preprocessing of data, after that every feature has an average value of 0 and a variance of 1.
For clustering we use the k-means++ algorithm [6]. Though it converges to a local minimum, visual control of plots with grouped pro les approved the possibility of this approach. Also we checked the results by applying hierarchical clustering method. It yielded similar results for clustering bias thermal pro les.
For the axes Y and Z, clustering is carried out in a similar way. 3.4
We estimate the similarity between clusters belonging to di erent sensors (either between di erent axes of the single sensor or between the accelerometer and the gyroscope in the certain chip) by means of the Jaccard index J (A; B), which measaures the similarity and diversity between two nite sample sets A and B.
The Jaccard index is de ned as the size of the intersection divided by the size of the union of the sets A and B:
J (A; B) = jA \ Bj :
jA [ Bj
(
The closer the value of the index to 1, the more similar the sets of data. If the Jaccard index is 1, then the sets are identical. If the Jaccard index is zero, then the sets do not contain common elements. This index was calculated for all pairs of sets of obtained clusters. 4
The bias thermal pro les for 250 MPU-9250 chips were obtained according the procedure described in Sec. 3.2. The overall per-axis data for accelerometers and gyroscopes are presented in Fig. 4. At rst, it seems that no uni ed \family" curve can be applied for all sensors in the batch to compensate the temperature dependencies of biases. From the other side, it is possible to split thermal pro les for each axis into groups with similar shape of the curves inside each group. Thus, the procedure of clustering was applied. As the temperature range for the obtained bias dependencies was 0 : : : 60 C, the size of the temperature interval 0.6 0.4 0.6 0.4 0.2 0.0 −0.2 −0.4 −0.6 3 2 1 0 −1 −2 −3 −4 −20 0 20
T, [oC] 40 60 was chosen of 5oC, that corresponded to 12 features { 12 slope coe cients for each bias curve.
Then the procedure of clustering was repeated with di erent number of clusters. The typical result of clustering into 5 groups for accelerometers and gyroscopes are shown in Fig. 5.
It can be mentioned that ve clusters are not enough for e ective clustering in our case, since a cluster might contain curves of di erent shapes. For example, top and bottom curves on AX4 inset (Fig. 5) are slightly di erent from other 4 pro les.
The example of clustering of accelerometers x-axis bias thermal pro le in 9 groups are shown in Fig. 6.
AX1
AX2
AX3
AX4
AX5
AX6
AX7
AX8
AX9
The distribution in 9 clusters seems to be more adequate than in 5 clusters as presented in Fig. 5. Though clusters AX1, AX4 and AX8 are very similar. GX1 )K1000 J( ion 800 t cnu 600 ft so 400 C 200 GX2 GX3 GX4 GX5
Also, clustering the data on Z-axis of gyroscope (bottom right inset in Fig. 4) produce good separation only in the case of at least 9 groups.
The number of clusters depends on the diversity of the pro les, which in turn depends on the production lot, complexity of the circuit, used materials and production technology, etc.
The elbow method [16] might be applied to determine the appropriate number of clusters. For the investigated set of MPU-9250 chips the results are presented in Fig. 7. Cost function is the inter cluster sum of squared distances from all cluster elements to the cluster centroid.
K = 9 5 10
15 20 25 Number of clusters K 30 35
According to the plot presented in Fig. 7 we recommend to use 6-10 clusters for each axis.
To check the correlation between clusters belonging to di erent measurement axes (X-, Y- and Z-axes of both accelerometers and gyroscopes), the Jaccard coe cient was calculated for all possible pairs of clusters (Sec. 3.4).
A typical result is shown in Table 1 demonstrating the data for accelerometer X-Y pair. In headers, the number in brackets represents elements in certain cluster. The number in brackets near a Jaccard index represents the number of common elements (intersection) between two clusters. s i x a X
AX1
(37)
AX2
(20)
AX3
(
AY1 (35)
The maximum value of Jaccard index is only 0.14 whereas the typical values are close to the zero. Consequently, there is no signi cant correlation between the clusters obtained. This indicates that the measurement axes of each sensor have their own independent temperature pro les of biases.
All other correlations between di erent axes of accelerometers and gyroscopes, for example, between AX and AZ, GX and GZ, AX and GY etc., were investigated and the corresponding tables of Jaccard indexes were created. However, no signi cant correlation observed.
It should be noted that the number of variations in the shape of the temperature pro le is nite, and each temperature pro le can be attributed to a particular cluster by measuring the bias only at several temperatures instead of measuring over the entire temperature range. Therefore, signi cant amount of time might be saved while obtaining thermal pro les.
As discussed above, MEMS chips manufacturers denote large error intervals for produced sensors. Even though it might seem that there is no \family" curve for a batch of sensors to compensate the thermal dependency, we observed certain similarity of these curves for a set of chips.
K-means++ clustering algorithm was used to reveal the groups of thermal curves. It was shown that sensors can be divided in several clusters by the features of slopes on 12 temperature intervals. The number of intervals was chosen experimentally. Further analysis will be done to tune the clustering parameters for optimum results.
The number of clusters depends on the complexity of the bias thermal pro le. Good group separation for each sensor axis for 250 MPU-9250 chips was achieved for a number of clusters from 6 to 10.
Also it was shown that there is no correlation between temperature pro les neither for the di erent axes (X, Y and Z) of one sensor, nor for di erent axes of an accelerometer and a gyroscope within a certain MPU-9250 chip. It means that temperature dependencies of biases of MEMS formed even on a single die are not similar to each other.
This research is supported by the grant 333GR/24464 (IRA-SME program).
8. Dai, G., Li, M., He, X., Du, L., Shao, B., Su, W.: Thermal drift
analysis using a multiphysics model of bulk silicon MEMS capacitive
accelerometer. Sensors and Actuators A: Physical 172(
In: Proceedings of SENSOR 2015. pp. 495{499. Association for Sensors and Measurement (2015) 10. Estivill-Castro, V.: Why so many clustering algorithms: a position paper. ACM
SIGKDD explorations newsletter 4(
New York (2011) 12. He, J., Xie, J., He, X., Du, L., Zhou, W.: Analytical study and compensation for temperature drifts of a bulk silicon mems capacitive accelerometer. Sensors and Actuators A: Physical 239, 174{184 (2016) 13. IEEE: Standard Speci cation Format Guide and Test Procedure for Linear, Single
Axis, Non-Gyroscopic Accelerometers. IEEE Std 1293 (1998)
14. IEEE: Standard Speci cation Format Guide and Test Procedure for Coriolis
Vibratory Gyros. IEEE Std 1431 (2004)
15. InvenSense: MPU-9250 Product Speci cation Revision 1.1. InvenSense Inc (2016)
16. Ketchen, D.J., Shook, C.L.: The application of cluster analysis in strategic
management research: an analysis and critique. Strategic management journal 17(