Introduction

A Micro-Electro-Mechanical-Systems (MEMS) device is a combination of integrated mechanical and electrical systems on a single chip that should be batch fabricated with high yield and no additional or subsequent assembly. The integration process may be implemented as hybrid integration using conventional wire bonding and flip-chips, or monolithic integration [11]. Monolithic integration offers superior system integration performance to hybrid systems but at an overall higher effort 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 significantly 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 affecting 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 fluctuating zero offsets from the true physical signal.

As is typical of inertial sensors, thermal effects are a primary driver of inaccuracies [13,14]. For example, the stiffness coefficient of beams, the damping ratio and other MEMS material parameters change with temperature and affect 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 [15].

−40 o C • • • + 85 o C is declared as ±30 o /s

Inertial sensors may be embedded in dead reckoning systems, which calculate the current position by using a previously determined position (also called fix), 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 profiles 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 significant 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 flexibly trade-off these two objectives and reach better accuracy with the same physical modules through more complex compensation algorithms.

So the goal is either to find 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 first 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 profiles and discuss possible correlation between them. Section 5 concludes the paper.

Related Works

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 significant amount of sensors. One of the examples of applying polynomial regression for finding 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 fit 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 k-means++ algorithm converges to a local minimum. Therefore, if it is appropriate, the exact optimum might be found by the brute force.

Devices Under Test and Measurements

Devices Under Test

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]. 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 MPU-6500 accelerometer and gyroscope produced by InvenSense and an AK8963 magnetometer of Asahi Kasei Microdevices Corporation.

Temperature Dependence Registration

To register temperature profiles 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 flange of MIMU module was attached to the cooling side of the TEC.

First, the inertial module was cooled to the minimum temperature of approximately 0 o C, 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 60 o C the experiment ended.

Clustering Methods for the Bias Thermal Dependencies

We propose the following method for clustering bias thermal profiles. 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 coefficients for all sensors are determined. Thus each sensor curve is characterized by a set of n slope coefficients. Then, slope coefficients are standardized and used as features in machine learning. Finally, bias thermal profiles are distributed in a number of clusters according to their sets of n slope coefficients. From the one hand, the number of intervals should not be large in order to improve classification. From the other, it should not be small, since the obtained piece-wise function should well describe the curved profile.

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 profiles 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 a i (t) = (a iX (t), a iY (t), a iZ (t)) consisting of temperature profiles registered for X-, Y-, and Z-axes of i-th accelerometer, i ∈ I. We assume that

a iX (t) = a X + ∆ iX (t) + δ iX a iY (t) = a Y + ∆ iY (t) + δ iY(1)a iZ (t) = a Z + ∆ iZ (t) + δ iZ ,

where a = (a X , a Y , a Z ) is the vector of true accelerations,

∆ i (t) = (∆ iX (t), ∆ iY (t), ∆ iZ (t)

) is the bias (zero offset) 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 a i (t) are assumed differentiable in the mean-square sense.

Let the temperature range T = [t 1 , t 2 ] is divided into N = mn points with equal intervals: τ 0 = t 1 , τ 1 = τ 0 + h, . . . , τ N = t 2 . The interval h is selected as follows: h = t2−t1 N , 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 profiles a iX (t), a iY (t), a iZ (t) into a certain number of intervals n, see Fig. 3 (n = 12). Similar thermal profiles of different sensors can be grouped. Every group has a typical set of derivatives (F 1 ...F 12 , Fig. 3).

We define ∆ i (t), i ∈ I as the bias thermal dependance for i-th accelerometer and carry out clustering accelerometers on the basis of the similarity of the derivatives of functions ∆ iX (t), ∆ iY (t), ∆ iZ (t).

We describe the clustering method for the X-axis. For the Y -and Z-axes the method works similarly.

Since the expectation value E[δ iX ] = 0, then

∆ iX (t) = E[a iX (t)] − a X .(2)

Further,

∆ iX (t) = (E[a iX (t)]) .(3)

We will use ãiXj as the sample mean values of a iX (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 a iX (t).

For each k = 1, . . . , n we construct linear regressions

f (τ ) = B iXk • τ + C Xk (4)

by sets of m points

τ j + h 2 , ãiXj | j = (k − 1)m + 1, . . . , km .(5)

The angular coefficients B iXk will serve as estimates for (E[a iX (t)]) on the intervals [τ (k−1)m , τ km ) for k = 1, . . . , n.

Before clustering, the resulting values of B iXk attributes for accelerometers need to be standardized.

Let B Xk be the sample mean and σ iXk be the sample standard deviation for the numbers B iXk , i ∈ I.

We introduce the standardized features:

D iXk = B iXk − B iXk σ iXk , i ∈ I, k = 1, . . . , n.(6)

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 profiles approved the possibility of this approach. Also we checked the results by applying hierarchical clustering method. It yielded similar results for clustering bias thermal profiles.

For the axes Y and Z, clustering is carried out in a similar way.

Similarity of Thermal Profiles in Different Clusters

We estimate the similarity between clusters belonging to different sensors (either between different 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 finite sample sets A and B.

The Jaccard index is defined as the size of the intersection divided by the size of the union of the sets A and B:

J(A, B) = |A ∩ B| |A ∪ B| .(7)

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.

Results and Discussion

The bias thermal profiles 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 first, it seems that no unified "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 profiles 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 was chosen of 5 o C, that corresponded to 12 features -12 slope coefficients for each bias curve.

Then the procedure of clustering was repeated with different 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 five clusters are not enough for effective clustering in our case, since a cluster might contain curves of different shapes. For example, top and bottom curves on AX4 inset (Fig. 5) are slightly different from other 4 profiles.

The example of clustering of accelerometers x-axis bias thermal profile in 9 groups are shown in Fig. 6.

AX1 AX2 AX3 AX4 AX5 AX6 AX7 AX8 AX9

Fig. 6: Clustering results for accelerometer X axis (9 clusters).

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.

AX1 AX2 AX3 AX4 AX5 GX1 GX2 GX3 GX4GX5

Fig. 5: Clustering thermal bias profiles into 5 groups for accelerometer X (top) axis and for for gyroscope X axis (bottom).

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 profiles, 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. 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 different measurement axes (X-, Y-and Z-axes of both accelerometers and gyroscopes), the Jaccard coefficient was calculated for all possible pairs of clusters (Sec. 3.4).

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 profile. 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 profiles neither for the different axes (X, Y and Z) of one sensor, nor for different 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.