CN104812063A - TOA (time of arrival) positioning method based on virtual sensors in indoor environment - Google Patents

Abstract

The invention belongs to the field of indoor wireless positioning, and relates to a TOA (time of arrival) positioning method based on virtual sensors in an indoor environment. The method includes the steps: deploying a plurality of solid sensors in the indoor environment; building one virtual sensor for each solid sensor; analyzing all barriers one by one, building a virtual sensor group corresponding to the i<th> solid sensor, and acquiring the length of a virtual path from each virtual sensor to an unknown object; acquiring virtual sensor groups corresponding to all M solid sensors; randomly selecting one virtual sensor from the virtual sensor groups corresponding to the solid sensors, and estimating all possible coordinates; screening obtained coordinate points and selecting final position estimation. Compared with a traditional indoor TOA positioning algorithm, the method has higher positioning precision in the indoor non-line-of-sight and multi-path environment.

Description

Time TOA localization method is reached based on the ripple of virtual-sensor under indoor environment

Technical field

The invention belongs to indoor wireless positioning field, under indoor environment, reach time (TOA) localization method based on the ripple of virtual-sensor.

Background technology

In recent years, due to the development of computer technology and wireless communication technology, context-aware services is come true by dream.The key realizing this service allows the environment of intelligent system understanding residing for service object, and then provide corresponding Intelligent Service to it, and the positional information of service object is one of most important environmental parameter, and this positional information not only comprises the position of outdoor but also has indoor.And under indoor environment for different buildings, indoor layout, material structure, the path loss that the difference of building yardstick result in signal is very large, meanwhile, the immanent structure of building can cause the reflection of signal, diffraction, refraction and transmission, form multipath (Multipath) and non line of sight (Non-Line-of-Sight NLOS) phenomenon, make the amplitude of Received signal strength, phase place and the time of advent change, cause the loss of signal, cause traditional outdoor wireless location algorithm precision when indoor application sharply to decline.

For above problem, the method for solution roughly has two classes: a class is statistical method.The error caused for non line of sight and multipath carries out modeling, describes error with certain statistical nature, but the shortcoming of this method needs model of error distribution known; A class is also had to be method of geometry.Mainly utilize the propagation characteristic of wireless signal, regarded as ray, suppress non line of sight and multipath error by the propagation path of analytic signal under specific environment.Such as utilize the TOA of redundancy to estimate, alleviated the impact of single non line of sight by the hybrid estimation of a large amount of non line of sight and sighting distance; Or utilize two-way TOA and weighting vector (AOA) to estimate to judge sighting distance and non line of sight situation, cast out multiple scattering, but receiving-transmitting sides all needs aerial array, hardware is comparatively complicated.Although virtual-sensor method is previously used, method before only considered the reflection of metope, and have ignored reflection and the diffraction situation of other barriers.

Summary of the invention

For the above-mentioned deficiency of prior art, the object of the invention is to propose a kind of ripple based on virtual-sensor and reach time (Time of Arrival, TOA) indoor orientation method, the error that minimizing non line of sight and multipath cause, improves the precision of indoor wireless location.

Technical scheme of the present invention reaches time TOA localization method based on the ripple of virtual-sensor under a kind of indoor environment, it is characterized in that comprising the following steps:

Step 1, is located in two-dimentional indoor environment and has N number of barrier, wherein jth barrier straight line A _jb _jrepresent, be expressed as with the linear equation of slope-intercept form:

$\left\{\begin{array}{cc}y=m_{j}x+a_j& m_j\&NotEqual;\&infin;\\ x=a_j& m_j=\&infin;\end{array}\right.$

Wherein, m _jrepresent slope, a _jrepresent y-axis intercept;

Step 2, disposes M entity transducer in indoor environment, and wherein i-th entity sensor coordinates is X _i=[x _i, y _i] ^t, the setting coordinate of unknown object to be positioned is X _t=[x _t, y _t] ^t, the time of advent of Received signal strength first observable crest is the TOA value of Article 1 detectable signal path FDP, is defined as τ, therefore, and the length of the signal path FDP measured be expressed as:

${\hat{l}}_i={\mathrm{\&tau;}}_i\&times;c=l_{\mathrm{FDPi}}+n_i$

Wherein c is velocity of wave, l _fDPifor the length of free from error actual FDP, n _irepresent measure error, obeying average is 0, and variance is σ _i ²normal distribution;

Step 3, sets up K to i-th entity transducer _iindividual virtual-sensor, has K _i>=N, wherein, kth _iindividual virtual-sensor coordinate definition is under indoor environment, virtual-sensor range-to-go is wherein, p _i=1,2,3 represent direct projection, reflection and diffraction respectively, and regulation Received signal strength experiences at most primary event or diffraction, and following steps 4 to step 6 is respectively the establishing method of the virtual base station in these three kinds of situations of direct projection, reflection and diffraction;

Step 4, when FDP be direct projection or transmission path time the path that measures be exactly the length of direct path:

l _FDPi＝‖X _i-X _t‖ ₂

Virtual-sensor is set up in the position of entity transducer there is virtual route length equal || || ₂represent 2 norms;

Step 5, when FDP is reflected signal, is arranged on entity sensors X by virtual-sensor _iabout straight line A _jb _jsymmetric points place, its positional representation is:

$X_i^{\mathrm{ki}}=\left\{\begin{array}{cc}{[m_j{x}_i,y_i/m_j-2a_j]}^T& m_j\&NotEqual;\&infin;\\ {[2a_j-x_i,y_i]}^T& m_j=\&infin;\end{array}\right.$

The path that reflected signal measures is equivalent to the distance of the direct path of signal between virtual-sensor and target, l _fDPibe expressed as:

$l_{\mathrm{FDPi}}={\left|\right|X_i^{\mathrm{ki}}-X_t\left|\right|}_2$

Virtual route length equal

Step 6, when FDP is diffraction path, Diffraction Point is the summit A of a jth barrier _j, then whole path is divided into two parts: a part be entity transducer to Diffraction Point, another part is that Diffraction Point arrives target, that is:

$l_{\mathrm{FDPi}}={\left|\right|X_i^{\mathrm{ki}}-X_t\left|\right|}_2+{\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2$

Wherein obtain from known indoor arrangement figure,

Virtual base station is arranged on an A _jplace, virtual route length be expressed as:

${\hat{l}}_{\mathrm{vi}}^3={\hat{l}}_i-{\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2;$

Step 7, according to step 4 ~ 6, until j=N from j=1, analyzes one by one to all barriers, sets up the virtual-sensor group corresponding to i-th entity transducer, and obtains the virtual route length of each virtual-sensor to unknown object;

Step 8, until i=M from i=1, analyzes one by one to all entity transducers, for a virtual-sensor group corresponding to each entity transducer, like this, obtains whole M the virtual-sensor group corresponding to entity transducer;

Step 9, selects all arbitrarily a virtual-sensor from the virtual-sensor group corresponding to each real sensor, the target location coordinate that first estimation one is possible the rest may be inferred, draws all possible K ₁× K ₂× ... × K _mindividual coordinate;

Step 10, screens the coordinate points obtained, if the virtual-sensor chosen is reflection foundation, and the condition that demand fulfillment reflection occurs, namely path needs to be greater than the distance between entity transducer to pip P, is expressed as:

${\left|\right|{\hat{X}}^{k_1,k_2,...,k_M}-X_i^{\mathrm{ki}}\left|\right|}_2{>\left|\right|P-X_i^{\mathrm{ki}}\left|\right|}_2$

Wherein P point is A _jb _jwith intersection point, its position coordinates is by calculating the method for above-mentioned two straight line intersection find intersections;

If the virtual-sensor chosen is diffraction foundation, the condition that demand fulfillment diffraction occurs, namely path needs to be greater than the distance between entity transducer to Diffraction Point, is expressed as:

${\left|\right|{\hat{X}}^{k_1,k_2,...,k_M}-X_i\left|\right|}_2{>\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2$

Step 11, choose meet following formula point as final location estimation result, that is:

$X_t={\hat{X}}^{a_1,a_2,...,a_M}$

$[a_1,a_2,...,a_M]=\underset{k_1\&Element;[1,K_1],k_2\&Element;[1,K_2],...,k_M\&Element;[1,K_M]}{\arg \min }\left\{D_{k_1,k_2,...k_M}\right\}$

$D_{k_1,k_2,...k_M}=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{({\left|\right|{\hat{X}}_t^{k_1,k_2,...,k_M}-X_i^{\mathrm{ki}}\left|\right|}_2-{{\hat{l}}_i}^{\mathrm{pi}})}^2.$

The present invention compared with prior art has the following advantages: (1) existing statistical property utilizing the method for statistical modeling to need to know in advance environment to the algorithm reducing non line of sight and multipath error, this is determined by a large amount of measurements with regard to needing, once environment change, this statistical property also can change needs and remeasure (2) existing method of geometry, or need the measurement by bulk redundancy, slow down the impact of non line of sight or judged the non line of sight degree of signal path by two-way TOA and AOA measurement, carrying out giving up or utilizing.First method needs a large amount of laying entity transducer, and its positioning precision depends on the entity number of sensors being in view distance environment; Second method signal transmitting and receiving both sides need aerial array, complex structure, and cost is higher.(3) existing virtual-sensor method only considered the reflection of metope, and have ignored reflection and the diffraction situation of other barriers.

Accompanying drawing explanation

Fig. 1 flow chart of the present invention.

The schematic diagram of virtual-sensor is reflected in Fig. 2 the present invention.

The schematic diagram of diffraction virtual-sensor in Fig. 3 the present invention.

Simulated environment in Fig. 4 the present invention.

The cumulative distribution function curve chart of Fig. 5 error, wherein, a) unknown object is positioned at A point; B) unknown object is positioned at B point; C) unknown object is positioned at C point.

Embodiment

With reference to Fig. 1, the concrete implementation step of TOA localization method based on virtual-sensor under a kind of indoor environment of invention is as follows:

Step 1, is defined in two-dimentional indoor environment and has N number of barrier, wherein jth barrier straight line A _jb _jrepresent, be expressed as with the linear equation of slope-intercept form:

$\left\{\begin{array}{cc}y=m_{j}x+a_j& m_j\&NotEqual;\&infin;\\ x=a_j& m_j=\&infin;\end{array}\right.$

Wherein, m _jrepresent slope, a _jrepresent y-axis intercept.

Step 2, disposes M entity transducer in indoor environment, and wherein i-th entity sensor coordinates is X _i=[x _i, y _i] ^t.The setting coordinate of unknown object to be positioned is X _t=[x _t, y _t] ^t.When measuring the TOA between entity transducer and unknown object, usually with Received signal strength first observable crest for benchmark, the time of advent of this crest is exactly the TOA value of Article 1 detectable signal path (first detectable path, FDP), is defined as τ.

Therefore, the length of the signal path FDP measured can be expressed as:

${\hat{l}}_i={\mathrm{\&tau;}}_i\&times;c=l_{\mathrm{FDPi}}+n_i$

Wherein c is velocity of wave, l _fDPifor the length of free from error actual FDP, n _irepresent measure error, obeying average is 0, and variance is σ _i ²normal distribution.

When there is barrier in indoor environment, signal path mainly contains following four kinds: direct projection, transmission, reflection and diffraction.Wherein, when FDP be direct path or transmission path when, signal be sighting distance (Light of Sight, LOS) propagate, l _fDPibe the Euclidean distance between i-th entity transducer and target.But when FDP be reflection or diffraction, signal be non line of sight (Non Light of Sight, NLOS) propagate, then l _fDPipath will be greater than Euclidean distance.

Step 3, in order to non line of sight problem is converted into line-of-sight problem, introduces virtual-sensor.To i-th entity transducer, set up K _iindividual virtual-sensor, has K _i>=N, namely for i-th entity transducer, needs to set up one or more virtual-sensor for each barrier.Wherein, kth _iindividual virtual-sensor coordinate definition is under indoor environment, virtual-sensor range-to-go is wherein, p _i=1,2,3 represent direct projection, reflection and diffraction respectively.Following steps 4 to step 6 discusses the establishing method that these three kinds of situations of direct projection, reflection and diffraction discuss virtual base station respectively.

Step 4, the path measured when FDP is direct path is exactly the length of direct path:

l _FDPi＝‖X _i-X _t‖ ₂

In order to the unification of whole model, need to set up virtual-sensor in the position of entity transducer obviously, virtual route length equal for transmission case, when signal can slow down through during barrier, energy attenuating and direction change, but compare with other non-market value, negligible.Therefore, transmission path is considered as direct path by us.

Step 5, after multiple reflections or diffraction, signal energy decay is very large, and not easily detect, therefore regulation Received signal strength of the present invention experienced by most primary event or diffraction.As shown in Figure 2, when FDP is reflected signal, virtual-sensor is arranged on entity sensors X _iabout straight line A _jb _jsymmetric points place, its position can be expressed as:

$X_i^{\mathrm{ki}}=\left\{\begin{array}{cc}{[m_j{x}_i,y_i/m_j-2a_j]}^T& m_j\&NotEqual;\&infin;\\ {[2a_j-x_i,y_i]}^T& m_j=\&infin;\end{array}\right.$

The path that reflected signal measures just is equivalent to the distance of the direct path of signal between virtual-sensor and target, at this moment l _fDPican be expressed as:

$l_{\mathrm{FDPi}}={\left|\right|X_i^{\mathrm{ki}}-X_t\left|\right|}_2$

Virtual route length equal

Step 6, as shown in Figure 3, when FDP is diffraction path, Diffraction Point is the summit A of a jth barrier _j, then whole path can be divided into two parts: a part be entity transducer to Diffraction Point, another part is that Diffraction Point arrives target, that is:

$l_{\mathrm{FDPi}}={\left|\right|X_i^{\mathrm{ki}}-X_t\left|\right|}_2+{\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2$

Wherein can obtain from known indoor arrangement figure.

Virtual base station is arranged on an A _jplace, virtual route length can be expressed as:

${\hat{l}}_{\mathrm{vi}}^3={\hat{l}}_i-{\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2$

Step 7, according to step 4 ~ 6, until j=N from j=1, analyzes one by one to all barriers, sets up the virtual-sensor group corresponding to i-th entity transducer, and this virtual-sensor group comprises K _iindividual virtual-sensor, and the virtual route length of each virtual-sensor to unknown object can be obtained.

Step 8, until i=M from i=1, analyzes one by one to all entity transducers, for a virtual-sensor group corresponding to each entity transducer.Like this, obtain whole M the virtual-sensor group corresponding to entity transducer, the virtual-sensor group wherein corresponding to i-th entity transducer comprises K _iindividual virtual-sensor.

Step 9, selects all arbitrarily a virtual-sensor from the virtual-sensor group corresponding to each entity transducer, corresponding virtual route is existing TOA location algorithm (as least square, maximal possibility estimation etc.) is utilized to draw a possible target location coordinate the rest may be inferred, can draw all possible coordinate, and total quantity of coordinate is K ₁× K ₂× ... × K _m.

Step 10, to the K obtained ₁× K ₂× ... × K _mindividual coordinate points is screened.

If the virtual-sensor chosen is reflection foundation, the condition that demand fulfillment reflection occurs, namely path needs to be greater than the distance between entity transducer to pip P, can be expressed as:

${\left|\right|{\hat{X}}^{k_1,k_2,...,k_M}-X_i^{\mathrm{ki}}\left|\right|}_2{>\left|\right|P-X_i^{\mathrm{ki}}\left|\right|}_2$

Wherein P point is A _jb _jwith intersection point, its position coordinates can by calculating the method for above-mentioned two straight line intersection find intersections.

If the virtual-sensor chosen is diffraction foundation, the condition that demand fulfillment diffraction occurs, namely path needs to be greater than the distance between entity transducer to Diffraction Point, can be expressed as:

${\left|\right|{\hat{X}}^{k_1,k_2,...,k_M}-X_i\left|\right|}_2{>\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2$

, for the coordinate point set after screening, if TOA estimates do not have error, so necessarily there is the coordinate points of a target in step 11 make:

$\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{({\left|\right|{\hat{X}}_t^{a_1,a_2,...,a_M}-X_i^{\mathrm{ai}}\left|\right|}_2-{{\hat{l}}_i}^{\mathrm{pi}})}^2=0$

Wherein a ₁∈ [1, K ₁], a ₂∈ [1, K ₂] ..., a _m∈ [1, K _m]

But measure error necessarily exists in practice, therefore choose make above formula value minimum point as final location estimation result, that is:

$X_t={\hat{X}}^{a_1,a_2,...,a_M}$

$[a_1,a_2,...,a_M]=\underset{k_1\&Element;[1,K_1],k_2\&Element;[1,K_2],...,k_M\&Element;[1,K_M]}{\arg \min }\left\{D_{k_1,k_2,...k_M}\right\}$

$D_{k_1,k_2,...k_M}=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{({\left|\right|{\hat{X}}_t^{k_1,k_2,...,k_M}-X_i^{\mathrm{ki}}\left|\right|}_2-{{\hat{l}}_i}^{\mathrm{pi}})}^2$

Effect of the present invention is described by following emulation:

(1) simulated conditions:

Set up simulated environment as shown in Figure 4.3 entity transducers are arranged in RS1 (6,8), RS2 (14,16), RS3 (24,5).Unknown object is arranged in 3 test points, A (13,16), B (6,12), C (16,1).FDP between each entity transducer and each test point as shown in the figure.P1, P2 are two people, like this, in some cases, even if direct path exists, and neither FDP.The selection of each test point, represents different non line of sight degree, and when target is positioned at A point, FDP is two direct path and a diffraction path; When being positioned at B point, FDP is a direct path and two reflection paths; When being positioned at C point, FDP is a diffraction path and two reflection paths, is all obstructed path.

(2) content is emulated:

The present invention is applied in the middle of two step weighted least-squares (TSWLS) and maximal possibility estimation (ML) two kinds of TOA location algorithms by we below, and compared with the positioning precision used before the present invention, in addition, that together compares also has Semidefinite Programming algorithm (SDP).

Be under the condition of 0dB in TOA noise power, we compare cumulative distribution function (CDF) curve of error, and result as shown in Figure 5.When unknown object is positioned at A point, relative to three entity transducers, be all substantially in view distance environment, therefore, five kinds of algorithm cumulative distribution are more or less the same, and CDF curve overlaps substantially.When unknown object is positioned at B point, non line of sight situation is aggravated, and utilize two kinds of algorithm advantages of the present invention to start to embody, cumulative distribution is starkly lower than its excess-three kind algorithm.Particularly when unknown object is positioned at C point, non line of sight situation is the most serious, and accuracy benefits of the present invention is also obvious.

In sum, the present invention can realize the TOA location under indoor non line of sight and multi-path environment, and precision is higher.

Obviously, those skilled in the art can carry out various change and modification to the present invention and not depart from the spirit and scope of the present invention.Like this, if within scope amendment of the present invention and modification being belonged to the claims in the present invention and equivalent technologies thereof, then the present invention is also intended to comprise these and changes and modification.

Claims (1)

1. reach a time TOA localization method based on the ripple of virtual-sensor under indoor environment, comprise the following steps:

Step 1, is located in two-dimentional indoor environment and has N number of barrier, wherein jth barrier straight line A _jb _jrepresent, be expressed as with the linear equation of slope-intercept form:

$\left\{\begin{array}{cc}y=m_{j}x+a_j& m_j\&NotEqual;\&infin;\\ x=a_j& m_j=\&infin;\end{array}\right.$

Wherein, m _jrepresent slope, a _jrepresent y-axis intercept;

Step 2, disposes M entity transducer in indoor environment, and wherein i-th entity sensor coordinates is X _i=[x _i, y _i] ^t, the setting coordinate of unknown object to be positioned is X _t=[x _t, y _t] ^t, the time of advent of Received signal strength first observable crest is the TOA value of Article 1 detectable signal path FDP, is defined as τ, therefore, and the length of the signal path FDP measured be expressed as:

${\hat{l}}_i={\mathrm{\&tau;}}_i\&times;c=l_{\mathrm{FDPi}}+n_i$

Wherein c is velocity of wave, l _fDPifor the length of free from error actual FDP, n _irepresent measure error, obeying average is 0, and variance is σ _i ²normal distribution;

Step 3, sets up K to i-th entity transducer _iindividual virtual-sensor, has K _i>=N, wherein, kth _iindividual virtual-sensor coordinate definition is under indoor environment, virtual-sensor range-to-go is wherein, p _i=1,2,3 represent direct projection, reflection and diffraction respectively, and regulation Received signal strength experiences at most primary event or diffraction, and following steps 4 to step 6 is respectively the establishing method of the virtual base station in these three kinds of situations of direct projection, reflection and diffraction;

Step 4, when FDP be direct projection or transmission path time the path that measures be exactly the length of direct path:

l _FDPi＝||X _i-X _t|| ₂

Virtual-sensor is set up in the position of entity transducer there is virtual route length equal || || ₂represent 2 norms;

Step 5, when FDP is reflected signal, is arranged on entity sensors X by virtual-sensor _iabout straight line A _jb _jsymmetric points place, its positional representation is:

$X_i^{\mathrm{ki}}=\left\{\begin{array}{cc}{[m_j{x}_i,y_i/m_j-2a_j]}^T& m_j\&NotEqual;\&infin;\\ {[2a_j-x_i,y_i]}^T& m_j=\&infin;\end{array}\right.$

The path that reflected signal measures is equivalent to the distance of the direct path of signal between virtual-sensor and target, l _fDPibe expressed as:

$l_{\mathrm{FDPi}}={\left|\right|X_i^{\mathrm{ki}}-X_t\left|\right|}_2$

Virtual route length equal

Step 6, when FDP is diffraction path, Diffraction Point is the summit A of a jth barrier _j, then whole path is divided into two parts: a part be entity transducer to Diffraction Point, another part is that Diffraction Point arrives target, that is:

$l_{\mathrm{FDPi}}={\left|\right|X_i^{\mathrm{ki}}-X_t\left|\right|}_2+{\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2$

Wherein obtain from known indoor arrangement figure,

Virtual base station is arranged on an A _jplace, virtual route length be expressed as:

${\hat{l}}_{\mathrm{vi}}^3={\hat{l}}_i-{\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2;$

Step 7, according to step 4 ~ 6, until j=N from j=1, analyzes one by one to all barriers, sets up the virtual-sensor group corresponding to i-th entity transducer, and obtains the virtual route length of each virtual-sensor to unknown object;

Step 8, until i=M from i=1, analyzes one by one to all entity transducers, for a virtual-sensor group corresponding to each entity transducer, like this, obtains whole M the virtual-sensor group corresponding to entity transducer;

Step 9, selects all arbitrarily a virtual-sensor from the virtual-sensor group corresponding to each real sensor, the target location coordinate that first estimation one is possible the rest may be inferred, draws all possible K ₁× K ₂× ... × K _mindividual coordinate;

Step 10, screens the coordinate points obtained, if the virtual-sensor chosen is reflection foundation, and the condition that demand fulfillment reflection occurs, namely path needs to be greater than the distance between entity transducer to pip P, is expressed as:

${\left|\right|{\hat{X}}^{k_1,k_2,...,k_M}-X_i^{\mathrm{ki}}\left|\right|}_2>{\left|\right|P-X_i^{\mathrm{ki}}\left|\right|}_2$

Wherein P point is A _jb _jwith intersection point, its position coordinates is by calculating the method for above-mentioned two straight line intersection find intersections;

If the virtual-sensor chosen is diffraction foundation, the condition that demand fulfillment diffraction occurs, namely path needs to be greater than the distance between entity transducer to Diffraction Point, is expressed as:

${\left|\right|{\hat{X}}^{k_1,k_2,...,k_M}-X_i\left|\right|}_2>{\left|\right|X_i^{\mathrm{ki}}-X_i\left|\right|}_2$

Step 11, choose meet following formula point as final location estimation result, that is:

$X_t={\hat{X}}^{a_1,a_2,...,a_M}$

$[a_1,a_2,...,a_M]=\underset{k_1\&Element;[1,K_1],k_2\&Element;[1,K_2],...,k_M\&Element;[1,K_M]}{\arg \min }\left\{D_{k_1,k_2,...k_M}\right\}$

$D_{k_1,k_2,...k_M}=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{({\left|\right|{\hat{X}}_t^{k_1,k_2,...,k_M}-X_i^{\mathrm{ki}}\left|\right|}_2-{\hat{l}}_i^{\mathrm{pi}})}^2.$