Ultrasonic localization method based on receiver array optimization schemes

Abstract

Ultrasonic localization system using trilateration has been adopted in various types of full-enclosed and semi-enclosed space. For small-scale or enclosed scenarios, there is an important phenomenon, that is, the ultrasonic wave may rapidly reflect back and forth between the walls to form different normal waves during the process of propagation, which may break the inverse square law of sound intensity (the sound intensity in free space is inversely proportional to the square of the distance). The other problem is that the sound intensity of some local corner points are obviously stronger than those points near the acoustic source, which may confuse the ultrasonic receivers and bring significant differences in measurement errors. To avoid excessive errors influence for the positioning accuracy, in this article, a novel ultrasonic localization method based on receiver array optimization schemes is reported. The main purpose is to activate the ultrasonic receiver in the area as well as shield the ultrasonic receiver outside the area. The system error will be reduced based on this method, and most of the redundant data will be eliminated. The simulation results are visualized and compared with the existing methods. The positioning accuracy measured by our method is increased significantly.

Introduction

Location-based service is now ubiquitous with the high-speed Internet and has wide applications. ¹ In the outdoor environment, global satellite navigation, such as BeiDou Navigation Satellite System (BDS) or Global Positioning System (GPS), has been able to provide higher accurate positioning. However, in the complex indoor environment, satellite navigation cannot meet the requirements of localization service. Indoor positioning technology is used to solve the problems that the satellite signals are weak when they reach the ground and also cannot penetrate the buildings. With the fusion of indoor and outdoor positioning technologies, all azimuth should be taken into account.

The indoor localization technology is divided into two categories. One type of technology needs to sample the signal fingerprints in advance, by comparing the receiving fingerprints with the pre-generated signal fingerprint models to achieve the users’ location, such as Wi-Fi ² and geomagnetic location. However, this kind of method still has some drawback. For example, the acquisition of signal fingerprint characteristics in the communication system is high cost, and the storage of signal fingerprint features also needs a significant amount of disk space. Moreover, due to the small indoor size, sometimes the signal fingerprint characteristics of the communication system are influenced by multipath propagation. In addition, the penetration of electromagnetic waves is mighty, and adjacent rooms will interact with each other. Since the accuracy of commercial positioning is about meter or decimeter level, all the factors above will affect the positioning accuracy and produce significant errors. This technology may also receive other electromagnetic interferences and cannot be applied in the complex and sensitive electromagnetic environment.

The other location technology is to make use of pre-deployed additional base stations (BS) to achieve precise indoor positionings, such as radio frequency identification (RFID),^3,4 ultra-wideband (UWB),^5–9 Li-Fi, ¹⁰ ultrasound,^11–18 Bluetooth,^19,20 and ZigBee ²¹ technology. This type of technology can be used to increase accuracy by increasing the number of reference nodes, as shown in Figure 1; the reference nodes are arranged indoor and the positioning accuracy can be centimeter level. With the improvement of positioning accuracy and the reduction of location anchor cost, this kind of location is more and more popular. The RFID technology can get centimeter-level positioning information in a few milliseconds based on the proximity algorithm and the received signal strength indication (RSSI) algorithm. ⁴ Moreover, it can overcome the influence of line-of-sight (LOS) propagation, which means that the location label is tiny, as well as the cost is inexpensive. However, RFID technology has no communication ability, low anti-interference ability, and difficult to integrate into other systems. The safety of the user privacy protection and international standardization is imperfect. The UWB positioning method uses time difference of arrival (TDOA)^22,23 or time of arrival (TOA) algorithm to analyze and calculate the label position. ⁵ The multipath resolution ability is high, and the localization accuracy is at decimeter level. It does not need to use the carrier of traditional communication system but transmits and receives the very narrow pulse with nanosecond and below to transmit the data. Therefore, UWB has the bandwidth of the GHz magnitude. It is the mainstream technology in the indoor location application. Compared with RFID technology, UWB technology has the function of communication, anti-interference, and high security. ⁶ But, when the new blind spots are added, the active communication is needed, the power consumption will be high, and low-frequency band utilization is low.

Figure 1.

Ultrasonic radiation schematic.

Bluetooth for indoor location is also relatively mature, ²⁰ and positioning accuracy is higher than Wi-Fi, equivalent to UWB. The main advantage of Bluetooth is that it is easy to integrate into smartphones and other mobile devices. Also, it does not need additional location tag support and is very convenient to switch on or off. However, the anti-interference of Bluetooth is poor, and the distance of communication is short. ZigBee is a new wireless ad hoc network technology with the features of short range, low energy, low rate, and high efficiency. The data can be transmitted from one node to another by radio waves. ²¹ However, the signal transmission of ZigBee is affected by multipath effect and mobile. It is not suitable for the indoor mobile location. The machine vision positioning technology can be used in the indoor positioning. Laser radar or camera is used to modeling indoor scenes in simultaneous localization and mapping (SLAM)^24,25 and then the machine vision matching algorithm is used to estimate the position. This positioning method is only suitable for the use of specific scenes ²⁶ and has high time cost of positioning process.

Compared with the technology scheme mentioned above, the location system based on ultrasonic and wireless sensor network (WSN) has obvious advantages. For example, the slow speed of voice transmission (about 340 m/s) makes the structure of ultrasonic receiving sensor simple and low cost. Also, when using the distance from transmitter to receiver signal (time-of-flight, ToF) to calculate distance, and the time delay $(\tau )$ data are transmitted back to the central calculating unit, ultrasonic technology can tolerate lower time synchronization accuracy than UWB. Ultrasonic signal has slight influence on environmental noise, strong diffraction, and strong directionality, which makes ultrasonic technology more practical and valuable. Furthermore, the miniaturized ultrasonic sensor can be easily integrated into smart phone, and it is easier to promote.

In this article, we propose a new method based on receiver array optimization schemes to improve the measurement accuracy. The main purpose is to solve the problems about the rapid multiple reflection and attenuation of the ultrasonic waves in the confined space, which results in the irregular distribution of the indoor wave intensity. The key technology is to divide the sound intensity points of indoor space into the active point spaces using graphics methods. Then, the ultrasonic receivers in the active space will be activated. The receiver data outside the active space will be shielded. The important advantages about this method are that the accuracy of measured data can be optimized, and at the same time, the large offset caused by too low sound intensity is excluded. The main theory and algorithm used in this method include Delaunay’s excellent triangulation feature,^27,28 TDOA ranging method based on CHAN algorithm,^29,30 and Kalman linear filter ³¹ to perform error calibration.

The rest of this article is organized as follows. In section “Related works,” the related works are introduced, including research status and system architecture of ultrasonic positioning. In section “Active area integration and optimization based on DT,” active area integration and optimization method is reported based on Delaunay triangulation (DT). In section “The CHAN algorithm and optimization of trilateral localization,” the CHAN algorithm and optimization of trilateral localization schemes are presented. And in section “Experiments and results,” the simulation experiment and the result analysis are given. The last section is “Conclusion.”

Related works

The technology of ultrasonic indoor positioning has been developed for a long period. In the same distance range, the ultrasonic frequency (20–40 kHz) is much lower than the radio wave (30 GHz), making the ultrasonic ranging more accurate. At present, ultrasonic positioning technology has been able to achieve the accuracy of 0.01 m. ³² However, most of the systems are one-dimensional (1D) or two-dimensional (2D), and few can achieve three-dimensional (3D) high-precision positioning experiments. In order to achieve high-precision 3D positioning system, M Cobos et al. ³³ used Beepbeep strategy to eliminate synchronization problems and provided high-accuracy location values. Loudspeaker and microphone are used as the system models for positioning calculation. But the problem of indoor not line-of-sight (NLOS) propagation cannot be considered. This problem is solved by AD Angelis et al. ³⁴ Besides using fixed anchor points, several mobile Beacons are set. These Beacons can be moved indoors and can be accurately positioned in complex environments. The centimeter-level positioning accuracy within the scope of 4 m can be obtained. Concerning joint parameter estimation and precise calibration for broadband acoustics, A Kumar and J Mcnames ³⁵ used 5 anchors and 13 reference points to achieve the average positioning accuracy of 0.86 mm in the test space of 0.5 m × 0.5 m × 0.25 m.

Based on the research of ultrasonic positioning system, we have established the indoor acoustic attenuation model according to the basic properties of acoustics. In this article, the ultrasonic receiver array and the handhold ultrasonic transmitting device are used to realize the high-precision 3D location of the full range multi-sound source in the general indoor.

The basic characteristic of ultrasonic

The ultrasonic transmitters are much smaller in size than the space-average free path and are considered as the point source. The point source is also in line with the strong directionality of the ultrasonic beam angle. ³⁶ It is also vital to choose the suitable frequency. It determines the detectable distance and the beam angle of the ultrasonic to a certain extent. The higher the frequency, the farther the measuring distance, and the smaller the beam angle. As shown in Figure 2, after emitting the ultrasonic wave, it is exposed outwardly in the form of a cone, θ is the beam angle of the ultrasonic, h is the detectable distance of the ultrasonic wave, and the shaded area is the area of ultrasonic radiofrequency. ³⁷ The radius r is the circumference at the furthest distance that travels. The relationship between r and h is as formula (1)

$$r=h\times \tan \theta$$

(1)

With the same amplitude, the vibrational energy of an object is proportional to the vibration frequency. The frequency of the ultrasonic wave is greater than 20 kHz. When it propagates in the air, the vibration frequency of the medium particle is very high, and the energy is much enormous. According to the formula (2), the wavelength $(\lambda )$ of ultrasound is very short and the diffraction power is very poor. It can be linear propagation in a homogeneous medium with highly directional

$$\lambda =\frac{u}{f}$$

(2)

Figure 2.

Ultrasonic radiation schematic: (a) beam angle and (b) beam surface.

In previous literature, ultrasonic transmitters were placed on the ground and receivers were placed on ceilings. However, such placement is not suitable in actual positioning applications, and transmitters are placed in an arbitrary position. The direct propagation direction may be arbitrarily changed. Regarding the stereo surround sound effect of a gorgeous movie theater, the receivers are arranged on all of the planes to obtain a stereoscopic receiving array.

Then, the rectangular reverberation room needs to be established. In a chamber, all surface reflections and absorbed sound fields form a reverberation room. Direct ultrasonic reaches all the walls in all directions in succession. After a part of the energy is lost through reflection, countless wavelets with different amplitudes and phases are formed in the room, forming numerous sound fields. These strong and weak sound fields make up a Normal Wave. The cost of direct sound to Normal Wave ³⁸ is so instantaneous that people can hardly feel it. For example, in a standard living room (e.g. $7\times 5\times 3{\mathrm{m}}^3$ ), the average free path is only about 3 m. Even after five times reflections, normal waves are established only spending 0.005 s. This process is a moment and can be considered as an initial reverberation.

Calculating a Normal Wave should consider the wave equation initially. According to the continuity equation as formula (3), Euler momentum equation as formula (4), and state equation in thermodynamic as formula (5), the wave equation of acoustic is as formula (6)

$$\frac{\partial \rho }{\partial t}+\nabla (\rho ·\mathbf{U})=0$$

(3)

$$\frac{\partial \rho }{\partial t}+{\rho }_0\nabla ·\mathbf{U}=0$$

(4)

$$\frac{\partial p}{\partial t}=c^2\frac{\partial \rho }{\partial t}$$

(5)

$$\frac{{\partial }^{2}p}{\partial t^2}=c^2{\nabla }^{2}p$$

(6)

The wave equation is the linearized, lossless wave equation for the propagation of sound in fluids. All fluids in this article are air fluids. c is the speed for acoustic wave in air. The acoustic pressure $p(x,y,z,t)$ is a function of the label $x,y,z$ and time t. In rectangular coordinates, ∇ and ${\nabla }^2$ are the 3D Laplacian operator, where $\nabla =(\partial /\partial x)+(\partial /\partial y)+(\partial /\partial z)$ and ${\nabla }^2=({\partial }^2/\partial x^2)+({\partial }^2/\partial y^2)+({\partial }^2/\partial z^2)$ .

In this article, we use the single-chip microcomputer to produce more than 20 kHz square wave pulse signal as the ultrasonic sound source, ³⁹ as a set of discrete values. Due to the difference in sampling frequency, the sound source information capacity of unit time is different. The influence of sound source on the receiving point at a certain time is calculated as formula (7)

$$W(S,R,t)=\{\begin{array}{cc}A\lceil t_0\times f\rceil ,\hfill & t>t_0\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{otherwise}\hfill \end{array}$$

(7)

where W is the amplitude of the acoustic wave under the influence of the sound at the current time, and A is the array acoustic amplitudes at the point where the sound source propagates for a period of time. The acoustic wave is determined by the source position S, receiving position R, and the vibration time t. The propagation time between S and R can be represented by the formula $t_0=l/v$ , where l is the length of the distance. The numeric value of l can be achieved by the following formula $l=|S-R|$ .

To achieve simple localization, in this article, we propose a transmitter-based indoor positioning system. In this system, a transmitter-based method that achieves ToF measurement without clock synchronization is introduced for indoor acoustical positioning. Moreover, an adaptive technique of adjusting the signal level according to the measurement environment for accurate ToF measurement is proposed and combined with a transponder-based technique. In this method, a terminal is localized as follows:

• A terminal transmits a request signal to anchors.

• When the anchors received the request signal, they first measure the request signal level, and transmit a response signal to the terminal immediately. The response signal level from each anchor is adaptively normalized so that the terminal can receive multiple response signals with appropriate levels.

• The terminal processes the response signals and obtains multiple ToFs. By combining multiple ToFs and the position of anchors, the terminal calculates its position by a trilateration method.

Based on the above ultrasonic properties, a reverberation room model for indoor acoustic attenuation is constructed. The ultrasonic is a high-frequency $(f)$ wave and it attenuates very rapidly in the reverberation room. Due to the sensitivity of temperature and humidity, the experimental environment must also ensure constant pressure and heat insulation.

Similarly, the coordinates of the sound source that are relative to other boundaries can be obtained.

Active area integration and optimization based on DT

At first, the Delaunay subdivision is applied to ultrasonic location. The receiver array is optimized according to the nature of the Delaunay partition to improving the positioning accuracy.

DT construction and formation in the 3D space

The distribution of indoor active sound field can show the sound intensity characteristics of the local region in a sense. However, due to the presence of the line of sight in 3D space, and the point set cannot intuitively represent the overall sound intensity distribution in closed rooms. In this article, the method of graphics is used to make the discrete point set continuous regionalization, the active space and active surface are formed, and the indoor sound intensity distribution can be visualized.

It has always been one of the hot topics in graphics research to aggregate discrete points into surfaces or spaces in 3D space. The DT algorithm is superior to other algorithms based on its mature theory and criterion. On one hand, DT has a complete mathematical theoretical foundation, which can ensure DT convergence for any complex model of triangulation. On the other hand, the algorithm is efficient and easy to implement, and the quality of generated cells is good.

The DT network, which was first proposed by Russian mathematician B. Delaunay in 1934, is a geometric dual graph of the Voronoi graph, which has a strict mathematical definition and a complete theoretical basis. He proposed that there must be a subdivision algorithm, so that the minimum inner angle of all triangles is the largest, which is the DT.

For the active points that have been measured, we use the DT to generate the corresponding active space. In 2D, the DT has the features of the maximum circumscribed circle and minimum angle maximization, which can reduce the emergence of deformed triangles as much as possible. We can use these properties in the 3D space. To obtain the accuracy of ultrasonic positioning to the millimeter level, the size of the generated subdivision point set may reach tens of millions of levels, which poses a technical challenge to the existing algorithms. In this article, a rapid method to generate Delaunay tetrahedron is presented, which can be used to meet the needs of indoor instant location service. The relationship between DT and Voronoi diagram is shown in Figure 3. As shown in Figure 4, the DT will be established quickly.

Figure 3.

Dual relationship between Delaunay (solid lines) and Voronoi graphs (dashed lined).

Figure 4.

The incremental insertion algorithm: (a) check the three sides of the insertion point, (b) check the criterion of the circumcircle, and (c) connect the insertion point and the rest of the nodes.

In this article, we choose 3D space scattered point set. ⁴⁰ Similar to the point-by-point insertion algorithm on the plane, the algorithm used in the 3D space-scattered point set is: First of all, create a very large triangular network, such as a cube, which contains five tetrahedral. Next, select a point at the remaining point and find the tetrahedron that contains the point. Starting from the tetrahedron that contains the point, the related topological information of the body record utilizes the empty outer ball detection to find out the outer tetrahedron set where the outer ball contains the current insertion point. The outer boundary of the tetrahedron set is the nuclear cavity to be found, which is formed by some triangular patches. ⁴¹ Finally, reconstructing the tetrahedral mesh in the nucleus to create a 3D DT. The main steps of construction DT is described in Algorithm 1

$$3\{\begin{array}{c}\overleftarrow{n}·\overleftarrow{\mathrm{AP}}=u(x-x_1)+v(y-y_1)+w(z-z_1)\hfill \\ \multicolumn{1}{c}{\overleftarrow{n}·\overleftarrow{\mathrm{BP}}=u(x-x_2)+v(y-y_2)+w(z-z_2)}\\ \multicolumn{1}{c}{\overleftarrow{n}·\overleftarrow{\mathrm{CP}}=u(x-x_3)+v(y-y_3)+w(z-z_3)}\\ \multicolumn{1}{c}{\overleftarrow{n}·\overleftarrow{\mathrm{PO}}=u(x-x_4)+v(y-y_4)+w(z-z_4)}\end{array}$$

(8)

Algorithm 1 The Main Steps of Construction Delaunay Triangulation

Set: Build initial triangulation. Get $X_{\min },Y_{\min },Z_{\min },X_{\max }$ , $Y_{\max },Z_{\max }$ of a given set of points, and create a cube box. Proceed as follows iteration until all data points have been processed. Input: Pick a point $p(x,y,z)$ from a given set of points and search for the tetrahedron containing the point in the established tetrahedral mesh. 1: Find a normal vector $\overleftarrow{n}=(u,v,w)$ of one of the three sides ABC 2: Calculate formula (8) 3: If $\overleftarrow{n}·\overleftarrow{\mathrm{AP}}\le 0$ , $\overleftarrow{n}·\overleftarrow{\mathrm{BP}}\le 0$ , $\overleftarrow{n}·\overleftarrow{\mathrm{CP}}\le 0$ , and $\overleftarrow{n}·\overleftarrow{\mathrm{OP}}\le 0$ simultaneously, the Tetrahedron OABC contains a point P that can locate the Tetrahedron OABC 4: Otherwise, continue the above steps until a tetrahedron is located. 5: End for Output: Fully divided Delaunay Space.

In Algorithm 1, $A(x_1,y_1,z_1)$ , $B(x_2,y_2,z_2)$ , and $C(x_3,y_3,z_3)$ are the closest three points closest to P and $O(x_4,y_4,z_4)$ constitutes a tetrahedral cavity that does not include any known nodes.

The experimental environment was carried out in a closed room of $3\times 5\times 7{\mathrm{m}}^3$ , and the experiment was conducted in a room of the Guilin University of Electronic Science and Technology Library. The next experimental environment will be expanded to large shopping malls like Wanda Plaza. According to the ranging experiment, the 20 kHz ultrasonic wave used in this experiment can stably propagate indoors for 20 m. In areas over $20\times 20{\mathrm{m}}^2$ , they can be arranged in parallel, and the positioning result can be queried in a polling manner in a positive triangle array. It is undeniable that the human body will also block. This part of the NLOS situation can be eliminated by pedestrian dead reckoning (PDR). ⁴²

The Tetrahedron ABCD is the tetrahedron that needs to be detected. The center O of the ball is the center of the ball of the ABCD outer tetragon. The radius r is also the corresponding radius. The point P is the current insertion point. If $d^2(O,P)-r^2<0(>0\mathrm{or}=0)$ , then it can be judged that the point P is on the inner side (outer side or spherical surface) of the tetrahedron ABCD. Finally, reconstruction of the grid—deleting the tetrahedron in the core cavity and re-establishing the mesh in the core cavity according to the Delaunay criterion—is relatively simple. Only the triangular patch in the core cavity needs to be connected with the current insertion point in turn, and the resulting tetrahedral mesh is Delaunay tetrahedron grid.

The CHAN algorithm and optimization of trilateral localization

To solve the signal synchronous problem of the ultrasonic receiver, TDOA algorithm is selected. The Chan Algorithm (CA) is used to solve the positioning equations based on the weighted least square (WLS) method, and to obtain the noniterative closed solution. It has the best estimation performance for the TDOA measurement error. The measurement error mainly includes two aspects: one is the distance error, caused by the frequency deviation between the transmitters and the receivers in the wireless transmission process and the multipath effect and so on. Another is positioning error caused by the nonlinearity of the positioning equations. In order to get the accurate 3D position information, we use Kalman filter to correct the calculation error and to obtain the actual value in this article.

In the 3D indoor space, the receiver BS are activated according to the Delaunay subspace. The estimated location of the measured point P is $(x,y,z)$ , the location of the receiver base station BS _i is $(x_i,y_i,z_i)$ , and the distance between P and BS _i is expressed as formula (9)

$$R_i=\sqrt{{(x_i-x)}^2+{(y_i-y)}^2+{(z_i-z)}^2}$$

(9)

The basic principle of ultrasonic distance measurement is trilateration positioning method, that is, according to the three BS, we can obtain the distance $d_i$ between $B{S}_i$ and P. By sending an ultrasonic burst through the air from a transmitter to a receiver, the traveling time is $t_k$ , and the distance $(d)$ is $d=v{t}_k$ where v is the speed of sound. To accurately estimate the parameter $t_k$ , the transmitter and the receivers shared a common clock and the received time-stamp of the incoming signal that was transmitted at a known time instance was measured as $t_k$ . Depending on the system synchronicity, the measured ranges are combined with trilateration or multilateration algorithms to estimate the 3D position of the target.

For the activated base station $\mathrm{B}{\mathrm{S}}_i$ $(i=1,2,\dots ,n)$ , if the $\mathrm{B}{\mathrm{S}}_1$ is the reference anchor, we can define the $R_{i,1}(i=2,3,\dots ,n)$ in formula (10)

$$\begin{array}{cc}\multicolumn{1}{c}{R_{i,1}=R_i-R_1}& =\sqrt{{(x_i-x)}^2+{(y_i-y)}^2+{(z_i-z)}^2}\hfill \\ \multicolumn{1}{c}{}& -\sqrt{{(x_1-x)}^2+{(y_1-y)}^2+{(z_1-z)}^2}\hfill \end{array}$$

(10)

Let $X_{i,1}=x-x_i,Y_{i,1}=y-y_i,Z_{i,1}=z-z_i$ , the formula (11) can be defined as follows

$$R_i^2-R_1^2=2\left[\begin{array}{ccc}{X}_{i,1}\hfill & Y_{i,1}\hfill & Z_{i,1}\hfill \end{array}\right]\left[\begin{array}{c}x\hfill \\ \multicolumn{1}{c}{y}\\ \multicolumn{1}{c}{z}\end{array}\right]+K_i-K_1$$

(11)

where $K_i=x_i^2+y_i^2+z_i^2$ . The four BS are needed to be detected at least. When the number of BS is greater than five, the amount of measured values are more than unknown parameters. In this case, we could process redundant data by WLS. The method will obtain superior estimates. First of all, TDOA matrix is linearized, and WLS is used to get the initial value. Then, we use the known constraint condition such as the estimated coordinates and additional variables obtained in above to get the WLS estimation in second time. Accordingly, the improved estimated coordinates are obtained. Suppose that ${\mathbf{z}}_a=[{\mathbf{z}}_p^t{R}_1]$ is unknown parameter, when ${\mathbf{z}}_p=[xyz{]}^{\mathrm{T}}$ , then the linear equation with noise is established, as formula (12)

$$\psi =\mathbf{h}-{\mathbf{G}}_a{\mathbf{z}}_a$$

(12)

where

$$\mathbf{h}=\frac{1}{2}\left[\begin{array}{c}{R}_{2,1}^2-K_2+K_1\hfill \\ \multicolumn{1}{c}{R_{3,1}^2-K_3+K_1}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{R_{n,1}^2-K_n+K_1}\end{array}\right]$$

$${\mathbf{G}}_a=\left[\begin{array}{cccc}-x_{2,1}\hfill & -y_{2,1}\hfill & -z_{2,1}\hfill & -R_{2,1}\hfill \\ \multicolumn{1}{c}{-x_{3,1}}& -y_{3,1}\hfill & -z_{3,1}\hfill & -R_{3,1}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{-x_{n,1}}& -y_{n,1}\hfill & -z_{n,1}\hfill & -R_{n,1}\hfill \end{array}\right]$$

and ${\mathbf{z}}_a^0$ is actual location of unknown point.

Assume that each element in ${\mathbf{z}}_a$ is independent of each other, ${\mathbf{z}}_a$ is defined as formula (13)

$${\mathbf{z}}_a=({\mathbf{G}}_a^{\mathrm{T}}{\mathrm{\Psi }}^{-1}{\mathbf{G}}_a){\mathbf{G}}_a^{\mathrm{T}}{\mathrm{\Psi }}^{-1}\mathbf{h}$$

(13)

where $\mathrm{\Psi }$ is covariance matrix of error matrix $\psi$ , $\mathrm{\Psi }=E[\psi {\psi }^{\mathrm{T}}]=c^2\mathbf{BQB}$ , and $\mathbf{Q}=\mathrm{diag}\{{\sigma }_{2,1}^0{\sigma }_{3,1}^0\dots ,{\sigma }_{n,1}^0\}$ is covariance matrix of measured values, and the matrix $\mathbf{B}=\mathrm{diag}\{R_2^0,R_3^0,\dots ,R_n^0\}$ . Formula (13) can be calculated approximately as preliminary results, and the results are based on the condition that each element in ${\mathbf{z}}_a$ is independent of other. However, because the variable R is related to $(x,y,z)$ in ${\mathbf{z}}_a$ , some errors occurred. To obtain the more accuracy positioning results, further evaluation should be done. Next, we calculate the covariance matrix of the estimated positioning ${\mathbf{z}}_a$ , in noisy case as formula (14)

$$\{\begin{array}{c}{R}_{i,1}=R_{i,1}^0+c{n}_{i,1}\hfill \\ \multicolumn{1}{c}{{\mathbf{G}}_a={\mathbf{G}}_a^0+\mathrm{\Delta }\mathbf{G}}\\ \multicolumn{1}{c}{\mathbf{h}={\mathbf{h}}^0+\mathrm{\Delta }\mathbf{h}}\end{array}$$

(14)

As a result of ${\mathbf{h}}^0={\mathbf{G}}_a^0{\mathbf{z}}_a^0$ , replace ${\mathbf{z}}_a$ with ${\mathbf{z}}_a^0$ , according to formulas (14) and (12), we can obtain formula (15)

$$\begin{array}{cc}\multicolumn{1}{c}{\psi }& =\mathbf{h}-{\mathbf{G}}_a{\mathbf{z}}_a^0={\mathbf{h}}^0+\mathrm{\Delta }\mathbf{h}-{\mathbf{G}}_a{\mathbf{z}}_a^0\hfill \\ \multicolumn{1}{c}{}& ={\mathbf{G}}_a^0{\mathbf{z}}_a^0+\mathrm{\Delta }\mathbf{h}-({\mathbf{G}}_a^0+\mathrm{\Delta }\mathbf{G}){\mathbf{z}}_a^0\hfill \\ \multicolumn{1}{c}{}& =\mathrm{\Delta }\mathbf{h}-\mathrm{\Delta }{\mathbf{G}}_a{\mathbf{z}}_a^0\hfill \end{array}$$

(15)

When ${\mathbf{z}}_a={\mathbf{z}}_a^0+\mathrm{\Delta }{\mathbf{z}}_a$ , then $\mathrm{\Delta }{\mathbf{z}}_a$ and its covariance matrix can be expressed in formulas (16) and (17)

$$\mathrm{\Delta }{\mathbf{z}}_a=({\mathbf{G}}_a^{\mathrm{T}}{\mathrm{\Psi }}^{-1}{\mathbf{G}}_a){\mathbf{G}}_a^{\mathrm{T}}{\mathrm{\Psi }}^{-1}\mathbf{Bn}$$

(16)

$$\mathrm{cov}({\mathbf{z}}_a)=E[\mathrm{\Delta }{\mathbf{z}}_a\mathrm{\Delta }{\mathbf{z}}_a^{\mathrm{T}}]=({\mathbf{G}}_a^{0\mathrm{T}}{\mathrm{\Psi }}^{-1}{\mathbf{G}}_a{)}^{-1}$$

(17)

According to the error estimate

$${\mathbf{z}}_a=\left[\begin{array}{c}{x}^0+e_1\hfill \\ \multicolumn{1}{c}{y^0+e_2}\\ \multicolumn{1}{c}{z^0+e_3}\\ \multicolumn{1}{c}{R^0+e_4}\end{array}\right]$$

in which $e_1,e_2,e_3,e_4$ are assumed error of ${\mathbf{z}}_a$ . The system of linear equations is defined in formula (18)

$$\psi \prime =\mathbf{h}\prime -{\mathbf{G}}_a^{\prime }{\mathbf{z}}_a^{\prime }$$

(18)

In formula (18), $\psi \prime$ is error vector of ${\mathbf{z}}_a$

$$\begin{array}{c}\mathbf{h}\prime =\left[\begin{array}{c}{(x^0+e_1-x_1)}^2\hfill \\ \multicolumn{1}{c}{{(y^0+e_2-y_1)}^2}\\ \multicolumn{1}{c}{{(z^0+e_3-z_1)}^2}\\ \multicolumn{1}{c}{{(R^0+e_4)}^2}\end{array}\right],{\mathbf{G}}_a^{\prime }=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill \end{array}\right]\hfill \\ {\mathbf{z}}_a^{\prime }=\left[\begin{array}{c}{(x-x_1)}^2\hfill \\ \multicolumn{1}{c}{{(y-y_1)}^2}\\ \multicolumn{1}{c}{{(z-z_1)}^2}\end{array}\right]\hfill \end{array}$$

The same as formula (13), we could obtain formula (19) for ${\mathbf{z}}_a^{\prime }$

$${\mathbf{z}}_a^{\prime }=(\mathbf{G}{\prime }_a^{\mathrm{T}}\mathrm{\Psi }{\prime }^{-1}{\mathbf{G}}_a^{\prime })G{\prime }_a^{\mathrm{T}}\mathrm{\Psi }{\prime }^{-1}h\prime$$

(19)

where $\mathrm{\Psi }{\prime }^{-1}$ is covariance matrix of estimated error, $\mathrm{\Psi }\prime =E[\psi \prime \psi {\prime }^{\mathrm{T}}]=4\mathbf{B}\prime \mathrm{cov}({\mathbf{z}}_a)\mathbf{B}\prime$ , and $\mathbf{B}\prime =\mathrm{diag}\{x^0-x_1,y^0-y_1,z^0-z_1,R_1^0\}$ . Finally, we obtain the ultimate prediction position ${\mathbf{z}}_p$ of the unknown node as formula (20)

$${\mathbf{z}}_p=\pm \sqrt{{\mathbf{z}}_a^{\prime }}+\left[\begin{array}{c}{x}_i\hfill \\ \multicolumn{1}{c}{y_i}\\ \multicolumn{1}{c}{z_i}\end{array}\right]$$

(20)

Kalman filter in estimation and prediction in a series of positionings

Kalman filter is used to further alleviate the range of error between actual value and calculated value in a series of positionings. The great success of the Kalman filter is due to its small computational requirement, elegant recursive properties, and its status as the optimal estimator for 1D linear systems with Gaussian error statistics. The estimated value by TDOA at a time t evolved from the prior state at time $t-1$ according to formula (21)

$$\overline{{\mathbf{z}}_t}=\mathbf{F}{\mathbf{z}}_{t-1}+\mathbf{B}{\mathbf{u}}_t$$

(21)

where $\overline{{\mathbf{z}}_t}$ is predicted current position by Kalman filter, ${\mathbf{z}}_{t-1}$ is last position, and $\mathbf{F}$ is the state transition matrix between time $t-1$ and t. The covariance of predicted value can be calculated by formula (22)

$$\overline{{\mathbf{P}}_t}=\mathbf{F}{\mathbf{P}}_{t-1}{\mathbf{F}}^{\mathrm{T}}+\mathbf{E}$$

(22)

Experiments and results

In our experiments, we describe the space model based on the rigid cuboid and boundary reflection. First, we used a virtual cuboid room that is with $7\times 5\times 3{\mathrm{m}}^3$ to describe the indoor acoustic field border. Then, the image sound source method was used to simulate the indoor acoustic.

By the method mentioned above, we take a single source and dual source as the test sample separately. The location of the sound sources are $\mathrm{Scr}1=[1.5\mathrm{m},2.0\mathrm{m},1.5\mathrm{m}]$ , $\mathrm{Scr}2=[3.0\mathrm{m},2.0\mathrm{m},1.5\mathrm{m}]$ , and $\mathrm{Scr}3=[4.5\mathrm{m},2.0\mathrm{m},1.5\mathrm{m}]$ .

The active ultrasound points distributed in the room are generated by the Delaunay 3D space algorithm described above. We can set the point $\mathrm{Scr}2=[3.0\mathrm{m},2.0\mathrm{m},1.5\mathrm{m}]$ , shown in Figure 5(a)–(c), and set points $\mathrm{Scr}1=[1.5\mathrm{m},2.0\mathrm{m},1.5\mathrm{m}]$ and $\mathrm{Scr}3=[4.5\mathrm{m},2.0\mathrm{m},1.5\mathrm{m}]$ , shown in Figure 5(d)–(f). From Figure 5, we could find that there are few variances in the acoustic intensity distribution between single-source and dual-source after a period of time $(>0.3\mathrm{s})$ , which can be regarded as the same situation. In a very short time, ultrasound fills the whole chamber, when $t>0.3\mathrm{s}$ , due to the absorption and reflection of the wall, the ultrasonic decay rapidly, the active point dropped sharply, the distribution of active point of ultrasonic wave in the chamber are obviously.

Figure 5.

Distribution of active point of acoustic field: (a), (b), and (c) are single-source active point distributions (Scr2). (d), (e), and (f) are dual-source active point distributions (Scr1 in blue, Scr3 in red), in which $t=0.003\mathrm{s}$ in (a) and (d), $t=0.1\mathrm{s}$ in (b) and (e), and $t=0.3\mathrm{s}$ in (c) and (f).

Next, set $t=0.3\mathrm{s}$ as the sampling time, the active points near the left and right walls are generated by DT classification points, as shown in Figure 6(a). If DT is applied in all of these points to partition the active points, it is found that the active tetragonal inclusion region will be generated, which cannot reflect the real active region of ultrasound due to the excessive distance between the left and right walls, shown in Figure 6(b). To solve this problem, we apply K-means clustering to optimize the entire active point through the Euclidean distance between two points in space, find out the position of each centroid iteratively, and mark the data with similar distances as the data set tuple. Then, DT algorithm is used to partition the active points. The result is shown in Figure 6(c). In Figure 6(c), we can find that the deeper the color value, the stronger the ultrasound intensity.

Figure 6.

The active point distribution of dual sound source at $t=0.3\mathrm{s}$ and Delaunay results: (a) active point distribution, (b) Delaunay triangulation of active points, and (c) Delaunay triangulation after clustering optimization.

Secondly, the active subspace is formed by DT subdivision. Figure 7 shows the different subspaces of two different ultrasonic receiving arrays arranged and formed. The receivers in Delaunay district are activated, which can reduce the calculation and improved the real-time of positioning.

Figure 7.

Ultrasonic receiving array schematic layout in the active space: (a) traditional receiver array and (b) optimized receiver array, in which the solid points are in the Delaunay District and activated, and the hollow points are not in the area and shielded.

The indoor positioning test in a real environment is designed and finished. Three handhold ultrasonic transmitters (Scr1, Scr2, and Scr3) are used and continuously emit ultrasonic signals. The indoor temperature $(\phi )$ was ${20}^{°}\mathrm{C}$ , and ultrasonic velocity $(v)$ was calculated as 343.30 m/s using the formula $v=(331.3+0.6\phi )\mathrm{m}/\mathrm{s}$ . Since the effect of humidity for the speed of sound is much smaller than for temperature, the effect of humidity on sound velocity can be ignored. The positioning coordinates are measured 10 times for each positioning point, and the average result is shown in Figure 8. The estimation error $\mathrm{\Delta }=\sqrt{{(\overline{x}-x)}^2+{(\overline{y}-y)}^2+{(\overline{z}-z)}^2}$ , $\overline{x},\overline{y},\overline{z}$ is the calculation positioning coordinate, and the $x,y,z$ is the transmitters positioning coordinate. In Figure 8(a), $\mathrm{\Delta }=7.23\mathrm{cm}$ , and in Figure 8(b), ${\mathrm{\Delta }}_1=4.84\mathrm{cm},{\mathrm{\Delta }}_2=3.99\mathrm{cm}$ .

Figure 8.

Illustration of the positioning of real in red and calculation in blue: (a) single point Scr2 and (b) double points Scr1 and Scr3.

Figure 9 shows the comparison of the results between method of activating all receivers, method of integral translation, and our method. In Figure 9, the blue cross line is the positioning coordinate obtained by normal receivers which are all of activating, resulting in excessive error mixing in the calculation. The green trig line is the positioning coordinate obtained by method of integral translation, and the red square line is our method used in this article. The black points are the actual position of ultrasonic transmitters. In Figure 9(a) and (b), we fix $z=1.5\mathrm{m}$ to measure the X and Y coordinates. Figure 9(a) show the LOS condition. We simulate an obstacle in $x=2.5\mathrm{m}$ and measure its impact on localization. It can be found that an enormous error fluctuation appears at the obstacle in method of activating all receivers and method of integral translation. However, for NLOS condition, measurement discarding method is utilized in our method. Combined with Kalman filter, it can remove NLOS error primely. As you can see from Figure 9(b), obstacle only has little impact on our localization algorithm. In Figure 9(c) and (d), we fix $y=2.0\mathrm{m}$ to measure the X and Z coordinates. It can be seen from Figure 9(c) and (d) that the error generated by the Z-axis is obviously greater than that of the Y-axis. Besides this, the rest of the case is similar to that of Figure 9(a) and (b).

Figure 9.

Comparison of measured values. Fixed $z=1.5\mathrm{m}$ in (a) and (b) and fixed $y=2.0\mathrm{m}$ in (c) and (d). (a) and (c) are LOS (line of sight) conditions and (b) and (d) are NLOS (not line of sight) conditions.

Figure 10 shows the comparison of positioning error among the method of activating all receivers, method of integral translation, and our method. We can see clearly in Figure 10(a) and (b) that the method reported in our article is obviously better than the other two methods in existence. In Figure 10(a), the accuracy of activating all receivers fluctuates greatly and the overall accuracy is unsatisfactory and the same as method of integral translation. As for Figure 10(b), although the method of integral translation has a certain correction for activating all receivers, the positioning accuracy is still far less than our method.

Figure 10.

Comparison of positioning error between other methods and our method: LOS for (a) and NLOS for (b).

Table 1 shows the comparison of absolute location errors (NLOS and LOS). It can be seen that in both NLOS and LOS cases, our method has higher accuracy and more stable error. Using DT to activate the space and activate the receiver’s positioning, the algorithm in this article will directly improve the indoor positioning accuracy to centimeter level.

Table 1.

Comparison of absolute location errors (NLOS and LOS) (cm).

Technique	Method of activate all receivers	Method of integral translation	Our method
Average error	15.2281 and 7.8873	12.0290 and 6.2244	1.8212 and 1.6575
Standard deviation	12.3258 and 5.8547	11.5319 and 3.9541	1.6047 and 1.3937
RMSE	18.2460 and 8.6819	15.3534 and 6.5412	2.2269 and 1.9753

NLOS: non-line of sight; LOS: line of sight; RMSE: root mean square error.

Due to the inevitable noise and environmental factors, there are still some errors. It can be seen that although the overall error has become smaller, the error between the initial point $(0\mathrm{cm},200\mathrm{cm})$ and $(300\mathrm{cm},200\mathrm{cm})$ still cannot be ignored. To solve this problem, consider replacing the walls of the experimental environment so that the wall reflections are as close as possible to the mirror effect, or ignore the calculated values of receivers that are too close to each other, and only use signals from the receivers on other walls for calculation.

Conclusion

In this article, an ultrasonic localization method based on receiver array optimization schemes is presented. We use the excellent properties of DT to partition the active indoor sound field space and combine theory of room acoustics and computer graphics to filter out large-error data. The main highlights with the help of graphic method can improve the indoor positioning accuracy efficiently. Specially, this method is suitable for general homes and indoor areas such as shopping malls and factories, which integrates a small ultrasonic sensor network in each direction into a large-scale sensor network. The other advantage in this method is that the receiver array can be relatively loosely disposed, which means that this arrangement will not only form a modular, distributed topology but also improve system scalability and security.

In the future work, we will consider the problem about absorption and damping effects of furniture and different materials of floors and ceilings in the ultrasonic localization. At the same time, multimodal data integration with Wi-Fi and iBeacon positioning is needed to research. The positioning algorithm and Delaunay generation algorithm can also be further simplified, achieving cost-effectiveness and commercial results.