Ahmad Osman

Valérie Kaftandjian

Ulf Hassler

Improvement of X-ray castings inspection reliability by using Dempster-Shafer data fusion theory

Keywords: X-ray imaging, castings inspection, true defects, false alarms, features extraction, data fusion, confidence levels, mass value, Dempster-Shafer theory

The aim of this work is to improve the classification of defects in X-ray inspection by developing a new method based on Dempster-Shafer data fusion theory where measured features on the detected objects are considered as information sources. From the histogram of features values on a learning database of manually classified objects, an automatic procedure is proposed to define a set of mass functions for each feature. The spatial repartition of features is divided into regions of confidence with corresponding mass functions. A smooth transition between regions is ensured by using fuzzy membership functions. The whole process is carried out without any expert intervention. Validation takes place on a testing database. Data fusion leads to a significant improvement of classification performances with respect to the actual system.

Introduction:

X-ray inspection is increasingly used as a tool for non destructive testing of industrial parts, such as aluminium castings in the automotive sector. In this field, new materials or processes are developed for cost and weight optimization. As those materials are operated near their load limit, it is necessary to check each produced part thoroughly. One main difficulty in X-ray castings inspection is the detection of false alarms (or false defects), especially if very small and low contrasted defects have to be detected. Therefore, efforts have to be made to reduce the rejection rate of good parts without risking to miss true defects. Once the X-ray images are processed for defect detection, the problem can be reduced to a classification problem, in order to discriminate the true defects (TD) from the false defects (FD). * Corresponding author: Tel. +33 4 7243 8061. Email: valerie.kaftandjian@insa-lyon.fr (V. Kaftandjian). Several approaches have been proposed for the classification of true defects and false defects, using X-ray inspection.

Filtering techniques adapted to specific defect shapes such as bidirectional template matching techniques were used by [START_REF] Newman | Bidirectional template-matching for 3d cad-based inspection[END_REF] for finding gross casting defects. The method is attractive but limited to specific defect types and time consuming.

Neural networks were used by [START_REF] Lawson | Intelligent segmentation of industrial radiographic images using neural networks[END_REF] for the automated detection of discontinuities in X-ray images of welds. A big advantage of this approach is to learn the neural net on a limited number of defect samples. In another work [START_REF] Bonser | Defect detection in partially complete SAW and TIG welds using online radioscopy and image processing[END_REF] a defect enhancement technique is introduced using a laws filter matched for round defects and Kirsch oriented filter matched for longitudinal defects was reported. The authors concluded that the neural net is more sensitive to fine defects but gives a high false alarm rate, whereas the filtering approach is more reliable and faster but less sensitive to small defects. Additional investigations of neural networks are described in [START_REF] Dobrazanski | Computer aided classification of flaws occurred during casting of aluminum[END_REF] for flaws in aluminium castings. The optimum network type was selected among multilayer perceptron (MLP), linear network, radial base functions network (RBF), generalized regression neural network (GRNN) and probabilistic neural network (PNN).

The best results were obtained employing the PNN and MLP networks with performances of 90 % and 81 % respectively on the test vectors, and 85 % and 94.5 % for the validation vectors. Performance is here given by the classification correctness of different defect types, but the false alarms rate is not mentioned.

In [START_REF] Hernandez | Automated defect detection in aluminum castings and welds using neuro-fuzzy classifiers[END_REF] the extracted features of each potential flaw were analysed and assigned to one of the following classes: 'defect' or 'regular structure' by using a Neuro-Fuzzy method named ANFIS [START_REF] Jang | Anfis: Adaptive-network-based fuzzy inference system[END_REF] which is used for constructing and training a fuzzy inference model that best classifies new data. The advantage of Neuro-Fuzzy systems is the combination of both properties: non linear learning based on numerical data and handling uncertainties in data. However the fuzzy model identification is very sensitive to which attributes are presented and class distribution in the training data. The developed system was used in casting and welds defects detection. The best performance was obtained for aluminium castings with 95% of true defect classification (57 out of 60 instances of defect class) and less than 1% (199 non defects classified as defects out of 22876 available) of false alarms rate. However the number of true defects on which the method is validated is small. Statistical pattern recognition was used in [START_REF] Mery | Classification of potential defects in the automatic inspection of aluminum castings using statistical pattern recognition[END_REF] to classify castings defects. In statistical pattern recognition the classification is performed using the concept of similarity: similar patterns are assigned to the same class [START_REF] Newman | Bidirectional template-matching for 3d cad-based inspection[END_REF]. Five classifiers were used by Mery: Linear classifier, Threshold classifier, Nearest neighbour classifier, Mahalanobis classifier and Bayes classifier. For each detected region, 71 features were extracted to characterise the potential flaws. The selection of the relevant features is made using sequential forward selection. According to the authors, best performances were achieved with 100% of true defect classification and 10% of false alarms rate for the Bayes classifier, and 94.7% of true defects classification and 4.1% of false alarms respectively for the threshold classifier. However, the method is validated on only 76 defects. A flexible inspection system called ISAR (Intelligent System for Automated Radioscopy) is presented in [START_REF] Fuchs | A new system for fully automatic inspection of digital flat-panel detector radiographs of aluminium castings[END_REF]. Defect classification is here performed by individual evaluation of characteristic features of the segmented defect area, such as size, contrast, elongation. A special type of defect classification was implemented by [START_REF] Wenzel | Texture analysis for classification of cast defects in x-ray images". 2 nd workshop[END_REF]. Fourier coefficients are computed in order to get translation-and rotation-invariant defect features. Those features can be assigned to a certain group, each representing a certain proper structure of the image. By computing the minimum distance of an actual defect region to each of the groups, artefacts generated at regular object structures may be recognized. Our study aims to improve the ISAR system performance, especially concerning the false alarm rate, by introducing data fusion of objects features.

Data fusion theory was already applied by the authors in building a classification system, where the so-called confidence levels are assigned by a source to each detected object. The fusion between confidence levels (or mass functions) is based on Dempster-Shafer theory (DS) and the final decision on the sample acceptability was done by a threshold in the confidence level. This approach was used in improving the detection of weld defects in [START_REF] Kaftandjian | Uncertainty modeling using Dempster-Shafer theory for improving detection of weld defects[END_REF], and in castings inspection [START_REF] Lecomte | Combination of information from several x-ray images for improving defect detection performances-application to castings inspection[END_REF]. The obtained results showed a more precise and reliable decision by using the data fusion approach. However, the supervision of the expert was necessary to assign the confidence levels.

In the present work, a new mass value attribution procedure is presented, which is completely automatic, so that the expert supervision is no more necessary. Moreover, eleven features were computed for each segmented structure in the casting (instead of four features in the previous study).

Each feature is considered as a source of information. By using data fusion, the aim is to build an automatic classification system. To be able to fuse different sources, a common space where values are comparable is necessary. An automatic method is introduced to convert from the space of feature values into the mass value space. The rest of this paper is organized as follows. Section 2 introduces Dempster-Shafer data fusion theory. Sections 3 and 4 present the proposed classification method and the performance measures.

Afterwards experimental results are presented in section 5, discussed in section 6 before concluding.

Dempster-Shafer data fusion theory

Dempster-Shafer (DS) evidence theory was developed as an attempt to overcome the limitation of conventional probability theory by handling uncertain, imprecise and incomplete information [START_REF] Dempster | Upper and lower probabilities induced by multivalued mapping[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF]. It is also suited for combining information from different sources.

In DS theory, a fixed set of N mutually exclusive and exhaustive elements, called the frame of discernment, is defined that is symbolized by Θ = {H 1 , H 2 , …, H N }. The frame of discernment Θ defines the working space for the application being considered since it consists of all propositions for which the information sources can provide evidence through using the so-called mass values.

Information sources can also distribute mass values on subsets of the frame of discernment, A i ∈ 2 Θ .

Here, A i designates a single hypothesis H i or union of simple hypotheses (composite hypotheses).

An information source assigns mass values only to those hypotheses, for which it has direct evidence. That is, if an information source can not distinguish between two propositions A i and A j , it assigns a mass value to the set including both propositions (A i ∪ A j ). This point is precisely the reason for us to choose the DS theory because it reflects the hesitation between two hypotheses.

Obtaining the mass distribution or function m(A i ) (

( ) 1 0 ≤ ≤ i A m
) is the most important step since it represents the knowledge about the current application as well as the uncertainty and imprecision incorporated in the selected information source. The mass distribution for all the hypotheses has to fulfil the following conditions:

( ) 0 = φ m ( ) ∑ ∈ = Θ i A i A m 2 1 (1)
Mass distributions m 1 , m 2 from two different information sources are combined with Dempster's orthogonal rule. The result is a new distribution, m = m 1 ⊕ m 2 , which carries the joint information provided by the two sources:

( ) ( ) ( ) ( ) ∑ = ∩ × - = i A Aq Ap q 2 p 1 i A m A m K A m -1 1 (2) 4 where ( ) ( ) p q 1 2 Ap Aq K m A m A φ ∩ = = ∑
K is often interpreted as a measure of conflict between the two sources and introduced in Equation

(2) as a normalization factor. The larger the K, the more conflicting are the sources and the less sense makes their combination. As a consequence some authors, Smets in particular [START_REF] Smets | The combination of evidence in the Transferable Belief Model[END_REF], require the use of the Dempster combination rule without normalisation. However, in this case, the relation ( 1) is no more verified. In our case, the normalised rule is used, and the K value is taken into account in the decision process.

In case of N different information sources B 1 , B 2 ...B N , the DS rule is:

( ) ( ) ( ) ... ... ( )
1 2 N i 1 1 2 2 N N B B B A i m B m B m B m A 1 -K ∩ ∩ ∩ = = ∑ (3)
where

( ) ( ) ... ... ( ) 1 1 2 N 1 1 2 2 N N B B B K m B m B m B φ ∩ ∩ ∩ = = < ∑
From a mass distribution, two functions can be evaluated that characterize the uncertainty about the hypothesis A i . The belief function Bel measures the minimum uncertainty value about A i whereas plausibility Pls reflects the maximum uncertainty value about this hypothesis. These two measures span an uncertainty interval [Bel(A i ), Pls(A i )], which is called "belief interval". The length of this interval gives a measurement of imprecision about the uncertainty value. Belief and plausibility functions are defined from the set of symbolic set of classes 2 Θ to the unit interval [0, 1]:

( ) ( ) ∑ = ⊆ i j A A j i A m A Bel (4) ( ) ( ) ∑ = φ ≠ ∩ i j A A j i A m A Pls (5)
These measures have been sometimes referred to as lower and upper probability functions.

If only two hypotheses are used, the credibility of each single hypothesis equals its mass value, and the credibility of the combined hypotheses (which represents the ignorance value) equals 1. This is the case for us when considering our frame of discernment: hypothesis H 1 ("this object is a true defect TD"), hypothesis H 2 ("this object is not a true defect", i.e. it is a false defect FD), and the ignorance is represented by the combined hypothesis (H 1 ∪H 2 =H 3 ).

Once data are combined, a decision rule must be applied. In our case, a threshold is applied on the mass value corresponding to the H 1 hypothesis. In the modelling step, specific care is taken in order to avoid conflicts between sources.

Proposed approach

The proposed new approach has the goal to automatically classify a detected object into either a TD or a FD, without expert supervision. This is done by analysing the features values extracted from each detected object, assigning to it a confidence level (or mass) and then combine mass values obtained from different features.

In the following, the term "object" refers to the result of the segmentation process, each object being potentially either a true or a false defect. Both the ISAR system and our new classifier are based on the same segmentation stage, so that the pre-processing does not influence the results (pre processing is explained in section 5.1).

The proposed approach consists of a learning stage, and a classification stage, which are detailed hereafter.

The learning is done on a population of known and manually classified objects. For those objects, a set of features is computed, and the histogram of each feature is built. From these histograms, regions of confidence are automatically defined with fuzzy transitions and masses are attributed, depending on the true defects proportion in the region. The complete procedure is detailed below.

After this learning stage, for each object, a set of mass values is computed for each of its features. It is then possible to consider each feature as a source of information about the set of hypotheses, and combine them. Final decision is obtained by a threshold applied on the combined mass associated to the hypothesis "defect" (H 1 ).

Learning

Estimation of true defect proportion

For a population of known detected objects (TD and FD), the set of values taken by a feature can be represented in the form of a spatial repartition. The histogram is used to divide this spatial repartition into several regions of confidence. To illustrate the method, we use an individual normal dis-tributed random variable x and two sets: set A representing the class of true defects TD, and set B representing the class of false defects FD (see figure 1). An instance of x can be classified into A or B. x can take any value of the codomain of a defect feature. This example will be used for illustration in each step.

The histogram of feature distribution is plotted in figure 2 where a number of 100 intervals was chosen in the feature values range in order to have a precise estimation.

Let I be the set of intervals forming the feature's histogram. For each interval i ∈ I, the percentage of instances of the class A present in this region is calculated using the following function:

, ( ) ( ) ( ) ( ) A A B A B h i P i h i h i = + ) (i h A
represents the number of instances of A inside i.

) (i h B
represents the number of instances of B inside i.

Elaboration of regions of confidence and associated mass functions

At first, each interval i is considered as a region, and a set of mass functions is attributed to this interval:

m(H 1 ) = , ( ) A B P i m(H 2 ) = 0 m(H 3 ) = , 1 ( ) A B P i - m(H 1 )
is the mass assigned to the hypothesis: instance of A or TD.

No mass is affected to the hypothesis H 2 (instance of B) to ensure that no conflict occurs between two sources of information in the fusion process.

m(H 3 ) is the mass assigned to the hypothesis ignorance (instance of A or B).

Secondly, subsequent intervals are congregated to form a region of confidence. To be able to merge two intervals, the variation of the function P A,B is computed. Let , ( )

A B P i ∆
be the obtained function (see figure 3).

If the difference between the P A,B values

, A B P ∆
for two adjacent regions is less than a fixed threshold T [first constraint], these two regions are merged and they will have in this case the same mass values:

If T ) i ( P ) 1 i ( P ) i ( P B , A B , A B , A < - + = ∆
: Region i is merged with region i+1.

Applying this method shows that the first constraint is not completely effective alone because some regions do not contain enough points to be considered as significant. For this reason a second constraint on the number of points existing in each region is imposed: a region should contain at least a certain percentage of points [second constraint] to be considered as having enough significance. Applying the merging procedure and ensuring enough points in each region, the obtained result is shown in Figure 4, where three regions of confidence are found. The m(H 1 ) is specified for each region.

Definition of fuzzy transitions between regions

It can be observed in figure 4 that the mass value from one region to another may change greatly.

Therefore it is important to be able to represent the transition between two regions of the graph in a continuous manner. The theory of fuzzy sets is used to solve this issue, by building membership functions which allow the continuous passage from one region to another [START_REF] Kaftandjian | Uncertainty modeling using Dempster-Shafer theory for improving detection of weld defects[END_REF][START_REF] Kaftandjian | The combined use of the Evidence Theory and Fuzzy Logic for improving multimodal Non Destructive Testing systems[END_REF].

Let fk be a source of information used to classify an object in class A or B according to its feature x, and R1, R2,. . . , Rs be the set of regions of confidence. k=1...N is the number of features considered as information sources.

Each object has a set of degrees of membership µ i (i=1...s) to each region i with:

Mass values attribution

The final mass attributed by a source f k to the object from its feature x is calculated by weighting the mass values of each region i by its degree of membership µ i :

1 1 ( ) ( / ) ( ) ( ) s k i i i m object H m x f x m R µ = ∈ = = × ∑ (6)
The final mass function obtained for the regions of figure 4 is represented in figure 5. The mass assigned to hypothesis H 3 is the complement to 1 of the mass assigned to H 1 .

3 1 ( ) 1 ( ) m object H m object H ∈ = - ∈ .
Afterwards the fusion process of different sources takes place.

Data fusion

The Dempster-Shafer rule is used to combine each two sources among N (see eq. 2), three most optimistic sources for each object and the N sources (see eq. 3).

Statistical fusion methods are also used by introducing: the mean mass value and the median mass value.

For an object, the mean mass is the mean value of all the mass values corresponding to this object.

It gives a sort of equal influence to the different opinions.

The median mass is the median value of all the mass values corresponding to an object. It allows to introduce a majority decision for the cases where one or few sources are not in accordance with the others.

Decision

To classify an object using the information source f k , a threshold S is applied on its mass value m(H 1 ) (which represents here the credibility of the hypothesis as said in section 2).

The object is classified as:

9

-a defect if m(H 1 ) ≥ S.

-unknown (defect or not) if m(H 1 ) < S.

Performance measurements and choice of optimal sources

Classification rates are introduced to measure the performance of the different combinations of information sources. These measures are calculated for different thresholds S applied on the mass values. For a certain threshold S, an object with mass m(H 1 ) is correctly classified if the source decision is in accordance with the human decision, this means:

-the object is a true defect (human decision) and its corresponding mass value is higher or equal to the threshold m(H 1 ) ≥ S.

-or it is a false defect (human decision) and its mass is lower than the threshold m(H 1 ) < S.

We introduce four measures:

-Correct decisions rate (PCD): The rate R was introduced in such a way to give the user the possibility to attribute more importance to the correct classification of either TD or FD. Usually more importance is given to TD detection.

Using these measures allows comparing the final performance of the method with the existing inspection system ISAR introduced in the next section.

Experimental results

This section is devoted to present the practical application of our new method. The pre-processing stage and features are described first, then the data base is presented and analysis of the results after data fusion is finally detailed.

Pre-processing stage and selected features

ISAR is a fully automatic system for radioscopic quality control currently in use in the production of castings. Within the database used in the present paper, evaluation is done by using a reference image, because thus a higher sensitivity to small defects is possible [START_REF] Fuchs | A new system for fully automatic inspection of digital flat-panel detector radiographs of aluminium castings[END_REF]. The data processing can be separated into five stages: registration, calibration, image processing, fault segmentation, and quality assessment.

• Registration Each inspected part that is processed during inline testing is not necessarily measured at exactly the same position. Due to mechanical instabilities, long-term changes of the moulds and slight differences between the moulds in use, an automatic registration step is needed, where the current test image is registrated to the appropriate reference image using an affine transformation.

• Calibration In order to allow for quantitative measurement of the depth of faults, i.e. to extract their length in direction of the penetrating radiation, a physical calibration object is evaluated. The necessary characteristic line is created using a step-wedge made out of the same material as the part and with typically 10 to 15 steps covering the maximum thickness of the objects to be tested.

• Image processing The reference image is subtracted from the current image. Thereby, all structures that do not appear within the fault-free reference part will become clearly visible, that is to say potential defects. In case of improper registration results however, edges due to the object structure may remain in the subtraction image, which is one cause of false defects.

• Segmentation The segmentation step is to assemble the single suspicious pixels into well defined fault regions appropriately. The analysis is based on the subtraction image. As result of this operation, we will have for each detected region several features that describe this object and which will be used for defect classification. Each detected object in this stage is classified by ISAR as true defect or false defect, and our new classification system starts with the same features.

• Quality assessment Quality criteria are applied by ISAR to the identified defects, which are defined by the end user for a final decision quality of the part (e.g. acceptance based on size or shape of defects). This postprocessing was not included in our system.

1

The features which are selected are described as following:

• Area [mm²]: is defined as the number of pixels of the detected object multiplied by the size of the pixel defined by the acquisition geometry. • Depth [mm]: the main objective of using the length image is to precise the size of the defect in X-ray penetration direction. Depth is defined as the mean value of the pixels of the detected object.

• Volume [mm³]: is calculated using the area occupied by this object multiplied by its depth.

• MaxDepth [mm]: is the highest value inside the defect region.

• DepthDev: is equal to standard deviation of the defect's depth in the defect region.

• InOutContrast [mm]: is also measured on the length image and is equal to InOutContrast = |mean(defect region) -mean(surrounding region)| • MaxElongation [mm]: is defined as the diameter of the minimum covering circle of the detected object. • Diff1stDev: This feature represents the difference, inside the defect region, of the mean gradient value within the reference image to the mean gradient value within the current image. • InOutContrastGV: is the same as InOutContrast but in grey values instead of mm.

• Thickness [mm]: is the thickness of the part inside the defect region. • Depth2Thickness: is the ratio between the depth of the object and the thickness of the part in the corresponding region.

Data base description

In this work a database is extracted from industrial images. The detected objects, also called potential defects, are classified manually into TD and FD.

The database is formed of 597 objects. It contains 382 true defects including oxides, gas voids and porosities and 215 false defects.

The database was divided into two parts:

-Learning database: formed of 115 false defects and of 243 true defects.

-Testing database: formed of 100 false defects and 139 true defects.

Application: Learning and classification

Spatial repartitions of four features (Area, Depth, InOutcontrast and Thickness) corresponding to the objects of the learning database are presented in figure 6. All the other features present their own spatial repartition with some overlapping between TD and FD classes. The Area graph exhibits the 1 better discrimination between classes with true defects mainly having small areas, although some false defects also.

First the estimations of the features histograms take place. Figure 7 illustrates the histograms of the four features represented in figure 6. As already mentioned, overlapping between TD and FD classes is important.

Using the feature's histogram the regions of confidence are built. The results for the four previous features are shown in figure 8 where different number of confidence regions are identified. The values on the right correspond to the mass values m(H 1 ) attributed for each point falling into the corresponding region. In addition, membership functions were built for each feature (see figure 9).

Using the estimated regions of confidence, their corresponding mass and membership functions, the single mass values of the learning database objects are computed (eq. 6). Mass functions obtained for the four previous features are also shown in figure 9.

Afterwards the fusion process takes place. As an example, the histogram of mass values obtained from the combination of Area and InOutContrast is presented in figure 10. It can be observed that the translation from feature values to mass values results into an easier separation between TD and FD. Moreover, the fact that the mass values are included between 0 and 1 facilitates the interpretation.

Subsequently, the classification of the learning database objects using different thresholds on the mass values given by each source (single and combined sources) is performed.

Performance results: Learning database

The performance measures for ISAR on the learning database are: PTD=0.978, PFD=0.723 and R=0.928 (with a=c=1, b=5). This choice of parameters a, b and c allows to give more importance to the classification of true defects, which is an industrial requirement.

Individual features performance after masses attribution is shown in Table 1. As it can be noticed, some of the features are better than the others (Area in particular, as was visible from its histogram), nevertheless, all of them will be kept for fusion due to the fact that they could bring some information.

1 ISAR has, in the tested configuration an excellent performance in the detection of true defects but it has a false detection rate of about 28 %. After the data fusion process, PTD, PFD and R are measured for all the single and combined sources (total of 70 sources). Several combined sources give better classification rates compared to the original ISAR decision (see figure 11). All these combinations present better PTD, PFD and R than ISAR. Grouped into a set L Σ , they are then used as sources of information to classify the objects of the testing database.

The best classification performance as measured by the overall rate R is obtained by the DS combination of the features Area and Depth with 96%, for a threshold on the mass values S=0.8, instead of 92% for ISAR. Table 2 gives the results obtained for PTD and PFD of this combination which are all higher than ISAR's measures.

When considering the performance of true defects classification only, the best source is the DS combination of Area and Thickness with PTD of 99%, for the threshold S=0.9, instead of 97% for ISAR. Here also PFD and R are both higher than ISAR's measures (see table 3).

The best performance for false defects classification only is obtained by the mean mass with a PFD of 98 %, for the threshold S=0.7, which is 28% better than ISAR's PFD. However, PTD is in this case lower than ISAR's PTD which is not acceptable for this application (see table 4), therefore the mean mass is not presented in figure 11. Nevertheless, this source will be considered for testing.

Validation: Testing database

The testing database is formed of 139 true defects and 100 false defects.

ISAR shows similar performance on the testing database with: PTD=0.980, PFD=0.725 and R=0.930.

After the learning step, the set of sources L Σ are used to automatically attribute mass values for the new objects of the testing database. Afterwards data fusion and classification take place.

All combined sources improve with an important percentage the classification of false defects and also are equal or even exceed ISAR performance on true defects detection (see figure 11 (b)). In the learning database, the optimal overall source was the fusion between Area and Depth. On the testing database, this combination presents a rate R of 96.37 % (for the threshold S=0.8), instead of 93.07% 1 for ISAR, with a very small difference from the optimal overall source on the testing database which is the fusion between Area and MaxElongation with R=96.86% (see table 5).

The optimal combination for the detection of true defects in the learning database (DS combination of Area and Thickness) is still the best combination for the detection of true defects in the testing database with a PTD of 0.99% for the threshold S=0.9 (see table 6). The improvement in the detection of true defects is 1 % relatively to ISAR, and 12 % in the detection of false defects. The optimal combination for the detection of false defects in case of testing database is still the mean mass value for the threshold S=0.7 (see table 7). The improvement in the detection of false defects is 25.8 % relatively to ISAR, but a loss of almost 4.5 % in the percentage of true defects detection occurs.

To illustrate two cases where ISAR fails, while data fusion gives the good result, two objects are chosen. Figure 12a shows a false defect classified as a true one, and figure 12b shows a true defect classified as a false (the object appears brighter than the background inside the rectangular zone).

Discussion

The first point to emphasize is that the results obtained during the testing phase confirm the robustness of the method. All sources selected during learning keep their high performance during testing.

In our previous work [START_REF] Kaftandjian | Uncertainty modeling using Dempster-Shafer theory for improving detection of weld defects[END_REF] only 4 features were selected as best candidates to discriminate between TD and FD. This selection of features (contrast to noise ratio, area, elongation and position) was done by the expert based on his knowledge of X-ray inspection. The drawback was to require the expert's supervision during learning. Here the approach is completely automated and all the features are considered without a priori information. The obtained performance cannot be compared with [START_REF] Kaftandjian | Uncertainty modeling using Dempster-Shafer theory for improving detection of weld defects[END_REF] due to the fact that the application was not the same.

When considering which features give the best classification rates, one can be surprised that the contrast feature is not among the best ones, although contrast is usually considered as having the main physical meaning. Indeed it is impossible to justify a posteriori why some features were better than the others when combined. This is due to the fact that our knowledge is only valid for TD, and the classification of FD cannot be based on a physical basis. This is precisely the difficulty of FD detection and the reason why even the expert find them difficult to classify. Therefore an important advantage of our method is that it adapts itself to all the objective data that are measured on the true 1 and false defects, without any a priori knowledge.

Another fact to mention is that the developed method is general in the sense that it can work for any inspection method. The translation from the feature space to the mass values space allows to combine any information sources independently of their origin (ultrasonic data, radiographic data...).

Conclusion

A new implemented method able to automatically detect regions of confidence from a feature's histogram and to affect degrees of confidence (mass values) to these regions was presented. The automatic learning process is based on a classified database. The method consists in translating a feature value into a confidence or mass value. These mass values are used within the fusion process.

The method was applied on a learning database and validated on a testing database for a casting application. As could be shown, several feature combinations are exceeding the current ISAR decision. Using the data fusion theory, while conserving at least the 98 % percentage of classification of TD given by ISAR, an improvement from 72 % to 92 % (see table 4) in the percentage of FD classification was achieved. An improvement of 3.8% of the overall rate R is obtained.

Some combinations can optimize even better either the TD or the FD classification. A PTD rate of 99% and PFD of 98.3% was obtained separately. Further work is needed to find an optimal decision rule which would associate the best TD and FD classification rates.

In a future work, a study of the robustness of this method for another application will be done.

Further investigation will also be conducted to measure the influence of the constraints on the derivative variation and on the percentage of points inside the regions of confidence on the results.

Table 1:

Performances of single features after mass attribution, in the learning database. Features are classified with respect to the overall rate R

Table 2:

Combination giving the highest overall rate R in the learning database.

Table 3:

Combination giving the highest true defects detection rate PTD in the learning database.

Table 4:

Combination giving the highest false defects detection rate PFD in the learning database.

Table 5:

Combination giving the highest overall rate R in the testing database.

Table 6:

Combination giving the highest true defects detection rate PTD in the testing database.

Table 7:

Combination giving the highest false defects detection rate PFD in the testing database.

2

Figure 5

 5 Figure 5 illustrates the fuzzy transitions obtained for the regions of figure 4. A fuzzy set is defined and a membership function (classical trapezoidal shape) is built for each region. The slope of the membership function is chosen to be proportional to the difference between the mass values of the adjacent regions.

Fig. 1 :

 1 Fig. 1: Two normal distributions where a star point represents a member of the class A and a cross point represents a member of class B.

Fig. 2 :

 2 Fig. 2: Feature's histogram of class A (full line) and class B (dotted line).

Fig

  Fig. 3:Illustration of

Fig. 4 :

 4 Fig. 4:Regions of confidence after merging and their corresponding mass values m(H 1 ).

Fig. 5 :

 5 Fig. 5: Membership functions corresponding to the regions of figure 4 (full line) and mass functions (dotted line).

Fig. 6 :Fig. 7 :

 67 Fig. 6: Spatial repartition of True and False defects for Area (high left), Depth (high right), InOutContrast (low left) and Thickness (low right) Fig. 7: Histograms for Area (high left), Depth (high right), InOutContrast (low left) and Thickness (low right).

Fig. 8 :

 8 Fig. 8: Regions of confidence and corresponding mass values m(H 1 ) for Area (high left), Depth (high right), InOutContrast (low left) and Thickness (low right).

Fig. 9 :Fig. 10 :

 910 Fig. 9: Membership functions (solid line) and mass functions (dotted line) corresponding to the confidence regions of Area (high left), Depth (high right), InOutContrast (low left) and Thickness (low right) Fig. 10: Histogram of mass values obtained from the combination of Area and InOutContrast. The majority of true defects have mass values close to one.

Fig. 11 : 1 False

 111 Fig. 11: Classification rates obtained on the learning database (a) and on the testing database (b). The selected sources are those giving a better overall classification rate than ISAR, and the same set of sources are still better on the testing database

Acknowledgements

This work was co-financed by the European Union and the Free State of Bavaria within the project ZeLuR.


Zhe He
email: zhehe@ucalgary.ca

Valé Rie Renaudin
email: valerie.renaudin@ifsttar.fr

Mark G Petovello
email: mark.petovello@ucalgary.ca

Gé Rard Lachapelle
email: gerard.lachapelle@ucalgary.ca

Use of High Sensitivity GNSS Receiver Doppler Measurements for Indoor Pedestrian Dead Reckoning

Keywords: high sensitivity GNSS, indoor multipath, pedestrian dead reckoning, tight integration, Doppler measurements, direct vector processing

Dead-reckoning (DR) algorithms, which use self-contained inertial sensors combined with gait analysis, have proven to be effective for pedestrian navigation purposes. In such DR systems, the primary error is often due to accumulated heading drifts. By tightly integrating global navigation satellite system (GNSS) Doppler measurements with DR, such accumulated heading errors can usually be accurately compensated. Under weak signal conditions, high sensitivity GNSS (HSGNSS) receivers with block processing techniques are often used, however, the Doppler quality of such receivers is relatively poor due to multipath, fading and signal attenuation. This often limits the benefits of integrating HSGNSS Doppler with DR. This paper investigates the benefits of using Doppler measurements from a novel direct vector HSGNSS receiver with pedestrian dead-reckoning (PDR) for indoor navigation. An indoor signal and multipath model is introduced which explains how conventional HSGNSS Doppler measurements are affected by indoor multipath. Velocity and Doppler estimated by using direct vector receivers are introduced and discussed. Real experimental data is processed and analyzed to assess the veracity of proposed method. It is shown when integrating HSGNSS Doppler with PDR algorithm, the proposed direct vector method are more helpful than conventional block processing method for the indoor environments considered herein.

Introduction

A wide range of commercial applications such as emergency services and cell phone location-based services (LBS) have driven the development of pedestrian navigation technology over the past several years. With a demand for low cost and high reliability, attention has been given to using additional sensors or devices such as Wi-Fi, inertial sensors and ZigBee radios to integrate with GNSS receivers. In particular, inertial sensors that can provide DR information have proven to be of great potential [START_REF] Foxlin | Inertial head-tracker sensor fusion by a complementary separate-bias Kalman filter[END_REF][START_REF] Babu | Ultra-Tight Integration of Pseudolites with INS[END_REF][START_REF] Pany | Performance of a Partially Coherent Ultra-Tightly Coupled GNSS/INS Pedestrian Navigation System Enabling Coherent Integration Times of Several Seconds to Track GNSS Signals Down to 1.5 dB-Hz[END_REF][START_REF] Li | Real-Time Ultra-Tight Integration of GPS L1/L2C and Vehicle Sensors[END_REF]. For pedestrian navigation, the PDR algorithm is often utilized because it makes best use of the fact that users are most likely to move on foot, and the costs of the required inertial sensors are relatively low. Much research has thus been directed towards improving PDR algorithms, either on reliable step length detection or improved heading estimation, such as [START_REF] Ladetto | Two Different Approaches for Augmented GPS Pedestrian Navigation[END_REF][START_REF] Randell | Personal Position Measurement Using Dead Reckoning[END_REF][START_REF] Afzal | Multi-magnetometer based perturbation mitigation for indoor orientation estimation[END_REF][START_REF] Renaudin | Step length estimation using handheld inertial sensors[END_REF].

On the other hand, GNSS receivers are often used together with PDR algorithms due to the fact that their errors are not accumulated. Many researchers have investigated integrating PDR algorithms with global positioning system (GPS) receivers. The feasibility and performance of using a low-cost motion sensor integrated with GPS and differential GPS (DGPS) was assessed in [START_REF] Salychev | Low Cost INS/GPS Integration: Concepts and Testing[END_REF]. The performance of using pedestrian dead-reckoning with a micro-electro-mechanical system (MEMS) inertial measurement unit (IMU) to aid high sensitivity GPS in harsh environment was reported in [START_REF] Lachapelle | Pedestrian navigation with high sensitivity GPS receivers and MEMS[END_REF]. Generally speaking, the absolute accuracy of the integrated system is governed by the accuracy of the GNSS receivers.

In order to enhance the performance of the GNSS receiver, a typical method is to increase the coherent integration, as shown in [START_REF] Pany | Performance of a Partially Coherent Ultra-Tightly Coupled GNSS/INS Pedestrian Navigation System Enabling Coherent Integration Times of Several Seconds to Track GNSS Signals Down to 1.5 dB-Hz[END_REF]. However, extremely long coherent integration requires that stringent requirements be met, such as to compensate for user motion, which is not the focus of this paper. Instead, the spatial gains obtained by non-coherent integration among satellites are explored. As signal degradation is inherent for indoor GNSS signals, even if high sensitivity receivers are still able to generate measurements in such challenged environments, their quality is usually poor. For example, pseudorange observations in shopping malls or tower block buildings are largely biased and can result in 20 to 60 m horizontal root mean squared errors, even with commercial HSGNSS [START_REF] Kjaergaard | Indoor positioning using GPS revisited[END_REF], and in these indoor environments, the benefits of integrating PDR with HSGPS are often limited due to multipath and fading. The performance of PDR integrated with conventional HSGPS using Doppler measurements in various indoor scenarios was assessed in [START_REF] Renaudin | Characterization of the Impact of Indoor Doppler Errors on Pedestrian Dead Reckoning[END_REF]. Results showed that the performance improvement of integrating conventional HSGPS's Doppler measurements with PDR was bottlenecked by the quality of Doppler measurement. It also indicated that HSGPS Doppler which uses block processing techniques [START_REF] Van Graas | Comparison of Two Approaches for GNSS Receiver Algorithms: Batch Processing and Sequential Processing Considerations[END_REF][START_REF] Borio | Doppler measurement accuracy in standard and high-sensitivity GNSS receivers[END_REF] in some indoor environments cannot provide beneficial Doppler measurements.

With this in mind, the major objective of this paper is to investigate a new method of generating HSGNSS Doppler measurements with the goal of improving PDR implementation in certain degraded signal scenarios. A direct vector processing method is thus proposed and developed. First, the velocity maximum likelihood estimate (MLE) is obtained. Then, Doppler measurements are generated based on such velocity MLE. The advantage of this approach is its reliability in harsh indoor environments where line of sight (LOS) and/or non-LOS (NLOS) signals are present. Subsequently, the benefit of these measurements for improving PDR algorithms indoors is investigated. The methodology proposed here is analyzed based on the indoor signal and multipath models, which are intrinsically related with the distribution of multipath statistics. Real experimental data is then presented to further verify the effectiveness of the proposed methodology.

The contributions of the paper are two-fold. First, a new direct vector processing receiver architecture is introduced and developed, which is shown to provide a more reliable velocity solution as well as Doppler measurements. Second, by using the new Doppler measurements integrated with PDR, the results are shown to improve the horizontal velocity accuracies by factors of more than 9% over the tradition implementation. Thus the effectiveness and benefits of the proposed Doppler estimation method are demonstrated and validated.

The paper is organized as follows: in Section 2, the signal and multipath models are introduced. After reviewing the architecture of conventional HSGNSS receivers, the proposed direct vector receiver is introduced. Then the velocity and Doppler estimation with direct vector processing in indoors are discussed in detail. In Section 3, the HSGNSS/PDR tight integration algorithm used in this paper is introduced. In Section 4, real indoor data is processed and analyzed. PDR-only solution, HSGPS/PDR tight integration with conventional Doppler and proposed Doppler measurements solutions are shown, compared, and discussed. Finally, conclusions are drawn in Section 5.

Direct Vector Processing in Indoor Multipath Environments

In this section, an indoor signal and multipath model is first introduced. The model is used to analyze how indoor multipath signals affect conventional HSGNSS Doppler estimation. After that, the proposed direct vector receiver architecture is introduced and discussed with comparison to the conventional HSGNSS receiver.

Signal and Multipath Model

The environment considered herein is indoors with dense multipath, where the multipath delay spread is usually smaller than one chip duration, or equivalently, the coherence bandwidth is much larger than the signal bandwidth (spreading code bandwidth in GNSS case). Under this scenario, a non-frequency selective channel or flat-fading channel is usually assumed which implies the multipath time-delay is non-resolvable [START_REF] Proakis | Digital Communications[END_REF].

Once the radio frequency signal is received by the antenna, the receiver down-converts it to near baseband. At this point, the general complex signal envelope can be expressed as:

1 ( ) ( ) ( ) ( ), 1 sat N l i i i i sat i r t a x t h t w t N        (1)
In Equation ( 1), N sat is the number of satellites in view, a i is the signal amplitude, τ i is the line-of-sight (LOS) signal time delay, h i (t) is the channel gain series. x i (t) = d i (t)p i (t) is defined as the product of the spreading code or pseudo random noise (PRN) code, p i (t) and the navigation data bits, d i (t). w(t) is input additive white Gaussian noise. Here the noise w(t) is assumed independent of the signal and has a flat power spectrum over the pre-correlation bandwidth. In the indoor case, the channel gain series can be further broken down as follows [START_REF] He | Modeling Hsgps Doppler Errors in Indoor Environments for Pedestrian Dead-Reckoning[END_REF]:

' , , 0 0 0 cos cos ( cos cos ) 1 ,, 11 ( ) ( ) 1 1 (
) ( ) i D Max m m m Di max i i i i i i M j t j jt i im m i ii LOS i NLOS i K h t e A t e K KM h t h t                  (2)
From Equation (2), The NLOS channel gain has M i multipath components, ,,

i i i m m m
A  are the weighting factor, azimuth and elevation angles for m i th multipath component. K i is the Ricean factor for i th satellites, which is the ratio between LOS and NLOS signal powers. The subscript 0 represents LOS component, while subscript m represents one NLOS path. It is noted that the channel gain series is decomposed into two components; one for the LOS signals and one for all NLOS signals. The first term on the right hand side in the above equation is the LOS channel gain series. The term involving the summation is the channel gain due to NLOS signals. For convenience, the total channel gain series is defined to have unity power.

Having presented the basic signal model with dense multipath, attention is now given to how this signal is handled within a GNSS receiver, and how it affects the conventional HSGNSS Doppler estimation. The conventional block processing technique for Doppler measurement is discussed in [START_REF] Borio | Doppler measurement accuracy in standard and high-sensitivity GNSS receivers[END_REF] and is based on the Doppler frequency MLE. With the notation introduced above, the correlator output for i th satellite with code delay and Doppler frequency   , , , ,

i j D i k f  can be expressed as: ,, , , , 2 ( , ) ,
[ ] [ ] ( )

D i k s i j D i k j f nT i f i s i j y n r n x nT e      (3)
i j D i k i j i D i k s ij D i k s N D i i f kK n N j f nT i s i j kK n N j f nT i i kK n i i D i k i i D i k kK f y n N r n x nT e N a h n e N a I f a P f                            1: ) kK  (4)
In Equation ( 4), it is apparent that the single Doppler frequency estimation will be directly affected by the channel gain series statistics. More precisely, this test statistic is I i (f), which denotes periodogram of the channel gain series. This statistic is an asymptotically unbiased spectral estimator [START_REF] Stocia | Introduction to Spectral Analysis[END_REF] of the true power spectral density (PSD) of channel gain series i.e., P i (f). Due to the fact that all estimators are based on time-limited sequences, the actual spectrum is the PSD of channel gain series convolved with spectrum of the time window. However, the periodogram is a reasonable approximation of actual PSD. So throughout the paper, the ideal PSD of channel gain series (P i (f) for i th satellite) is investigated instead of I i (f).

In the vector form, the Doppler MLEs for all available satellites will be:

sat sat NN D i i D i k i i D i k k K k K ii a I f a P f     f (5)
From Equation ( 4) and ( 5), it is also observed that for individual Doppler MLE, each channel or satellite is processed independently.

When a channel is time varying, the so called -Doppler spread‖ is sometimes introduced [START_REF] Fontan | Modelling the Wireless Propagation Channel[END_REF]. This term is defined as frequency spread in the signal spectrum due to different signal propagation paths contributing to a single fading channel. Consequently, in some non-symmetric ring scattering environments (i.e., where multipath from certain directions is much stronger than other direction), the NLOS components dominate the signal. Then the PSD of the channel gain series will be significantly distorted by certain channel parameters, and Doppler MLE can be biased. The most important channel parameters in this regard are the multipath angle-of-arrivals (AOAs). As such, once the distribution of the multipath AOA is specified, the ideal PSD of the channel gain can be derived [START_REF] He | Modeling Hsgps Doppler Errors in Indoor Environments for Pedestrian Dead-Reckoning[END_REF]. In a more compact form, if the NLOS carrier power spectrum is expressed as a function of multipath AOAs, then the overall carrier power spectrum has the following form: ,, ,

i i NLOS i LOS i ii i NLOS i avg i i avg i i D i ii K P f P f P f KK K g f f f KK             , , , , , 1 ( ) ( ) ( ) 11 1 ( , , , , ) ( ) 11
In Equation ( 6), K i is the Ricean factor, f D,i is the LOS Doppler, ,

i i i    , , , , , , avg i avg
  respectively represent averaged multipath azimuth, elevation angle and beamwidth for azimuth and elevation angles. Function g NLOS,i describes the major power distribution of all multipath components for the i th satellite. Its shape and centre are determined by the above mentioned multipath AOA parameters. In Equation ( 6), when the Ricean factor K is extremely small, multipath will dominate the signal and the Doppler MLE for this satellite may be biased. It is noted that as coherent integration time is extended, the LOS signal power might be accumulated to the point where it is usable, and the multipath power might decrease due to exceeding the coherence time of fading channel. However it is not practical to unlimitedly increase the coherent integration time, which usually calls for stringent requirements for both oscillator and inertial systems. In real scenarios, it is more common to have combined LOS and NLOS signals or NLOS-only signals present in the correlator outputs, which are the focus of the paper.

Conventional and Direct Vector GNSS Software Receiver Architectures

The architecture of a conventional high sensitivity GNSS software receiver that uses block processing technique is illustrated in Figure 1. The received intermediate frequency (IF) data is first fed to a Doppler removal and correlation (DRC) block. After the integration module, which might have coherent or non-coherent form, the despread and demodulated signals are then processed by a block processing technique. The pseudorange and Doppler measurements are estimated by selecting the maximum power in the correlator outputs (code phase and Doppler domain). With such measurements, the user position and velocity can be estimated by either least squares or Kalman filtering. When it comes to weak signal conditions, these measurements are often themselves biased and lead to very large residuals. An effective way to deal with such measurements is to weight them according to their accuracies. Due to the fact that the Doppler estimation accuracy is inversely proportional to the received signal power, once accurate signal power can be estimated, such weighting can be applied.

However, the C/N 0 estimation will get worse under weak signal conditions, which limits the benefits of using conventional HSGNSS Doppler measurements. Direct vector receiver architecture is thus proposed in Figure 2. With scrutiny of the figure, one can observe that in direct vector receiver, the signal power in the Doppler domain for all satellites is projected onto the velocity domains. In this way, the estimated velocity solutions are based on all the information contained in all correlator outputs across all satellites which is essentially the velocity MLE.

As for conventional high sensitivity GNSS receivers, the velocity estimation is only based on the Doppler MLEs, which discards some information before final solution is made. Another major benefit of the direct vector receiver is that weighting is automatically performed according to the received signal strength [START_REF] He | Improved Velocity Estimation with High-Sensitivity Receivers in Indoor Environments[END_REF]. In such direct vector processing, all the power in correlator domain is preserved Reference NavSolution

Nav-Domain MLE

i i I jQ   2 2
( )

i i I Q   Integration Correlation Map DRC NCO Channel #K i i I jQ   2 2
( )

i i I Q   DRC NCO Received IF
i i I jQ   2 2
( )

i i I Q   Integration Correlation Map DRC NCO Channel #K i i I jQ   2 2
( )

i i I Q   DRC NCO Received IF Data NavSolution Reference NavSolution
Bit Aiding

Receiver Tracking Channels Dopplers

Channel #1 until a navigation solution is obtained. For example in this paper, once a valid velocity solution is obtained, the Doppler measurement for each satellite can be generated accordingly. Details of how the velocity MLE in the presence of multipath signals is obtained are discussed in the following section.

Direct Vector Processing in Indoor Multipath

It is known that the measured Doppler frequency from GNSS receivers has the following relationship with the user velocity:

, , , , , ,, 1 ( ) ( ( ) ( ) 1 (
)) , 1, 2...,

D i i x u,x sat,i x i y u,y sat,i y i z u,z sa s t at ,i z f u v -v u v -v u v -v cdt i N       v ( 7
)
where λ is the wavelength of transmitted RF signal,

i i x i y i z u u u  u is
the direction unit vector from the satellite to the receiver, cdt is the receiver clock drift. As such, the user velocity vector is assumed to also include the clock drift term, i.e., , ,

.

u xyz u x u y u z v v v cdt  v Due
to the invariance property of maximum likelihood estimates [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF], the velocity MLE can then be easily obtained by: 22 11 ˆarg max ( ) arg max ( ) .

sat sat NN ML i i i i ii a I a P                vv v v v (8)
From Equation [START_REF] Renaudin | Step length estimation using handheld inertial sensors[END_REF], one can see that the final velocity solution is obtained by weighting available carrier power spectrums for all satellites. The weighting factors are automatically chosen as the actual signal powers (a i 2 ). The physical interpretation of velocity MLE can be considered as follows. First, the LOS signal power which is an impulse in the Doppler domain can be projected to the velocity domain on a satellite by satellite basis. Analogously, the NLOS signal power can also be projected into the velocity domain. The point in the velocity domain that has the greatest power will represent the final estimate. In the following, how the velocity MLE is obtained in the presence of indoor multipath signals will be shown in detail.

In order to project the signal power from the Doppler domain to the velocity domain, the relationship between the Doppler to the velocity can be used. For example, a small offset in the Doppler will cause a small offset in the velocity, these two terms are linearly related as shown in Equation ( 9

): , , , , , , , , 1 1 1 1 ,
, ,

T T e T D i i x i y i z u xyz i u xyz i n u enu i u enu f u u u C                 v a v a v e v (9)
In Equation ( 9), e i T is the projection vector from velocity to the Doppler for i th satellite. C n e is the rotation matrix from navigation frame to the earth frame.

, u enu v is the user velocity offset vector with respect to the velocity searching centre in navigation frame which also includes the clock drift term.

, Di f  is the Doppler frequency offset for this satellite relative to Doppler searching centre. In order to visualize the effect of multipath on velocity estimation, it is convenient to first consider the two dimensional case. Assuming vertical velocity and clock drift are already known or constrained, the Doppler offset is then only related with two horizontal velocity offsets ( E v  and N v  ). The carrier PSD over the Doppler domain is shown in the upper part of Figure 3. This is based on Equation ( 6).

The solid green line is the LOS signal power given by 2 1

ii i aK K  and is located at the LOS Doppler offset, , Di f  . On the other hand, the major NLOS signal power is distributed over a small region, and has a maximum power of 2 , ( )

1 i i DMP i i a g f K  
with a frequency offset of

, DMP i f 
. In the velocity domain, the LOS signal is equally distributed over the line or plane ,

T i u D i f    ev
. Similarly, for the NLOS power, the peak power is along the line or plane ,

T i u DMP i f    ev
, and powers of its adjacent region can be evaluated with different offsets according to Equation [START_REF] Randell | Personal Position Measurement Using Dead Reckoning[END_REF]. As more satellites are considered, the situation tends to that shown in Figure 4. In this case, there are three satellites considered. All three LOS powers intersect at the point A, so the power of point A denotes the total power for LOS signals. It is also shown that there are other intersection points, such as point B (two NLOS power together with one LOS), point C (two NLOS power only), and point C (one NLOS power and one LOS power). The more intersections located at a particular point, the more power will be observed, and the more probable it will be the final velocity estimate. In the velocity domain, the point that contains the maximum power will be considered as the final estimate of the velocity.

The velocity estimation procedure shown in Figure 4 relies on many factors. Under different scenarios, it may have totally different performance. Firstly, if all LOS components are dominating the signal, or equivalently, the Ricean factors for all three satellites are very large, then point

A will PSD f 2 , ( )
1 i Di aK ff K    , Di f  2 () 1 i i a gf K  , DMP i f  N v  E v  , T i D i f    ev , T i DMP i f    ev Power  Doppler Domain

Velocity Domain

Peak power Adjacent region naturally be the final estimate and should be similar to the least squares velocity estimation that processes Doppler MLEs independently. Secondly, if only one LOS component has a large Ricean factor, say the yellow curve, and the other two satellite signals have dominant NLOS components, then the power of point A may have less power than that of point B thus potentially making point B the final estimate. In this case, the result is not based on LOS signals only, but with partial LOS information and partial NLOS information. In other scenarios, it is equally possible that point C or point D might have the maximum power, which should be avoided. Consequently an admissible region is defined in order to discard such outliers, and the size of the regains varies according to different user activities. For example, the horizontal acceleration during walking in this paper is assumed within the range ±2 m/s 2 .

Figure 4. Velocity MLE with channel distortion (2D case).

On the other hand, as the number of satellites tracked increases, the LOS signal power in the velocity domain is accumulated without loss, while the NLOS signal power in velocity domain is usually dispersed. The reason is that the projection matrix ) is calculated for LOS signal propagation. In other words, D f is in the column space of the projection matrix. For the NLOS Doppler offset array, this is generally not true. The other reason is that the multipath statistics for each satellite are generally different and time varying, whereas the LOS signal is predictable and deterministic. If by accumulating power between satellites, the total power of LOS components continues to increase and if it surpasses all the dispersed NLOS power, a relatively good velocity estimate will occur. Another acceptable scenario is similar to point B in the Figure 4, where the power is highest with a combination of LOS and NLOS components. If more than two LOS powers intersect at point B along with other NLOS power, the final estimate then takes advantages of both LOS and NLOS signal components.

E v  N v  1 , 1
T D f    ev 2 , 2
T D f    ev 3 , 3
T D f    ev 2 , 2
T DMP f    ev 3 , 3
T DMP f    ev 1 , 1
T DMP f    ev A B C Admissible Range D Power 
The above section discusses the power projection in two dimensional spaces. In the following, the same principles are extended to the real scenario where four velocity states need to be estimated. Assuming there are N sat available satellites, the four-dimensional hyper-planes of LOS powers for all satellites are given by:

1 , , 1 2 , ,2 ,
,,

sat sat T u enu D T u enu D u enu D T u enu NN D f f or f                    ev ev E v f ev (10)
Each hyper-plane is associated with the LOS signal power of

2 22 1 1 2 2 12 , ,..., 1 1 1 sat sat sat NN N aK a K a K K K K    .
Analogously, the NLOS signal power is distributed over the hyper-planes and has a certain amount of power according to , () 6) with different frequency offsets:

NLOS i g f in Equation (
1 , , 1 2 , ,2 , ,, sat sat T
u enu DMP T u enu DMP u enu DMP T u enu DM NN P f f or f                    ev ev E v f ev (11)
For example, the peak powers associated with each hyper-plane of NLOS signal are

2 22 12 1 , 1 2 ,2 , 12 (
), ( ),..., ( ) 1 1 1

sat sat sat N DMP DMP M DMP N N a aa g f g f g f K K K      
. All hyper-planes described in Equations ( 10) and ( 11) finally form the velocity domain. The point in velocity domain associated with the maximum power will be considered as the velocity MLE.

The velocity MLE can be obtained by using the Equations ( 8), ( 9), [START_REF] Lachapelle | Pedestrian navigation with high sensitivity GPS receivers and MEMS[END_REF], and [START_REF] Kjaergaard | Indoor positioning using GPS revisited[END_REF]. The offset associated with maximum power in the velocity domain is considered as final estimate of velocity. Once the velocity MLE is obtained, the Doppler measurements for each satellite can then be computed according to Equation [START_REF] Salychev | Low Cost INS/GPS Integration: Concepts and Testing[END_REF]. In this direct vector processing method, the Doppler measurements are based on maximum use of LOS/NLOS signal power and maximum use of mutual information between each satellite. In this way, a certain amount of information shared among satellites (see Figure 4) is fully utilized before arriving at a result. In contrast, conventional block processing methods only use Doppler MLE, whose residual can be quite large in indoor environments, potentially causing the solution not to converge. In the following sections, the Doppler estimated from block processing and direct vector processing methods will be evaluated and compared by integrating with PDR algorithms with real experimental data.

HSGPS/PDR Tight Integration

In this section, the system model for the HSGNSS/PDR tight integration is first introduced. The measurement or observation models are then presented. Following this, the integration performance of using conventional and the proposed Doppler measurements with a PDR algorithm will be assessed in the next section.

System Model

The system state vector for the PDR filter is shown in Equation [START_REF] Renaudin | Characterization of the Impact of Indoor Doppler Errors on Pedestrian Dead Reckoning[END_REF] and includes (in order) the user's 3D coordinates in the local level frame, the horizontal and vertical speeds, the walking heading direction and the GNSS receiver's clock drift in units of range rate (i.e., multiplied by the speed of light, c):

, , , , , , hv E N U v v cdt     x (12)
The GNSS receiver clock bias is not present, since in this paper only the Doppler measurements are used to update the integration filter. In turn, this is because the paper focuses on assessing the benefit of the proposed Doppler measurements for PDR.

The system dynamic model of the pedestrian's position follows the equations of a classical PDR mechanization and is given by: sin( ) cos( )

h h v Ev Nv Uv      (13)
The velocities and the heading are further modeled as random walk processes:

h v hv vv v v        (14)
Similarly, the clock drift state is modeled as:

cdt cdt   (15)
In Equations ( 14) and ( 15),

h v  , u v
 ,   and cdt  are the white Gaussian driving noise of the corresponding state elements. The spectral density of the clock drift error noise is computed using a standard clock stability model, as found in [START_REF] Brown | Introduction to Random Signals and Applied Kalman Filtering[END_REF].

It is noted that discrepancies between PDR-derived velocity and GNSS-derived velocity are expected. In particular, due to the repetitive nature of the human gait, Doppler measurements typically exhibit oscillations over the course of a full gait cycle. In contrast, by its very nature, PDR velocities are effectively averaged over the course of a step and do not contain these oscillations, thus introducing an oscillatory discrepancy between the PDR-and GNSS-derived values. Correspondingly, these oscillations should be modeled to properly integrate the velocity information based on Doppler measurements.

For coping with these oscillations over a step, the classical PDR has been modified in this work. Normally, the filter state is only propagated when a step is detected. In the proposed approach, however, Doppler measurements are used at each GNSS measurement epoch, typically at a higher rate than the step frequency. As such, the PDR filter is propagated at this higher rate as well. However, PDR observations are still only used when a step is detected. It is expected that this approach will average the oscillating effects sensed by Doppler measurements over one step. The advantage of this asynchronous measurement update is that the integrated system could track the rapid changes of the heading sensed by the IMU or other heading sensors.

On the other hand, there are also some disadvantages with this new integration method. First, the computational load is increased because the filter is propagated more frequently. In this paper, a coherent integration time of 500 ms is used, such that the measurement update time for Doppler is 2 Hz, so computational load is not a major concern. Second, as the coherent integration time increases, the Doppler measurement actually conveys information about the average velocity (and thus attitude) during the integration interval. However, this type of averaged Doppler measurement is still usable since it can help to alleviate the long term heading drift of the PDR system.

Measurement Model

Having introduced the system model of the integration filter, the following discusses the measurement models. There are two types of measurement updates for the proposed tight integration; from the PDR and from the GNSS Doppler.

The PDR sensor update is composed of step length updated and heading update. The measured step length and walking directions are related to the user's position and velocity through the following equations:

h step SL obs s v t          ( 16
)
where s is the user's step length, step t  is the step duration and step  is the averaged heading over the last step.

The measured GNSS Doppler is related to the pedestrian's velocity via Equations ( 7) and [START_REF] Salychev | Low Cost INS/GPS Integration: Concepts and Testing[END_REF]. In order to be consistent with other observation models, it is listed again here:

, , , , , ,, 1 ( ) ( ( ) ( ) 1 (
))

D i i x u,x sat,i x i y u,y sat,i y i z u,z sat,i z f f u v -v u v -v u v -v cdt        v ( 17
)
where η f is the noise induced by the receiver. After linearization of Equation ( 17), then the observation equation is:

,1 ,1 ,1 ,2 ,2 , 3 1 , , , , 3 3 3 1 ˆ11 ˆ11 D D sat D sat sat sat sat f D D e D f D D M n M u DN fN Rx DN Filter h N N v f f v f f v cdt v f f v T v cdt                                                                  f H 1 H1        f η (18)
In Equation ( 18), ,

ˆˆˆT

D D D M Filter f f f    
are the filter predictions and ,1 ,2 ,

T D D D M Rx f f f    
are the measurements from the receiver. The subscript N sat before N sat was used denotes the number of satellites.

 

,1 3 cos( )sin( ), cos( ) cos( ), sin( )

H           
is the design matrix from velocities to pseudorange rates, ε and α are the elevation angles and azimuth angles for each satellite. λ is the wavelength of the transmitted GNSS signal, and matrix T 3×3 is defined as:

33 sin( ) cos( ) 0 cos( ) sin( ) 0 0 0 1 h h v Tv          (19)
It is noted that PDR only provides information in the horizontal plane rendering the vertical velocity unobservable except with Doppler observations. To further constrain the vertical component of the pedestrian's position, barometer records could be used, although this was not done here. Similarly, indoors, pedestrians are mainly walking on flat surfaces and only change their elevation when climbing/descending stairs or when taking an elevator.

Data Processing and Analysis

This section deals with real experimental data processing and analysis in order to assess the benefits of integrating conventional and proposed HSGNSS Doppler measurements with PDR sensors. First, the data collection is described briefly. Then the analysis and results are described. For simplicity, only GPS satellites are used, but it is expected that results would likely improve if a multiple GNSS constellation were used.

Data Collection Description

The experimental data was collected on the campus of the University of Calgary. The primary pieces of equipment were an NI front-end, a Novatel SPAN receiver and a LCI IMU [START_REF]Novatel SPAN TM Brochure[END_REF], all of which are shown in Figure 5, mounted on a backpack. The raw IF data was collected with the NI front-end at a rate of 5 Msps (complex). The reference trajectory is shown in Figure 6. This trajectory is generated using Novatel Inertial Explorer TM with a Novatel SPAN TM system. The SPAN TM system used here includes a tactical grade IMU (LCI) and a high precision GNSS receiver, i.e., SPAN TM receiver. The Inertial Explorer TM uses differential GNSS measurements to tightly integrate with IMU. However, the high precision SPAN TM receiver cannot provide GNSS measurements indoors due to the fact that the signal C/N 0 is too weak for standard tracking loops. In such cases, the solution is computed only using IMU measurements. In order to get better results, both forward/backward smoothing and coordinate updates are used. By using this approach, the estimated position standard deviation was better than 3 m at all times, and the estimated velocity standard deviation was better than 0.02 m/s.

The trajectory contains several different indoor environments, and two spots are chosen to emphasize how Doppler from direct vector HSGPS outperforms that from conventional HSGPS (scenario A and scenario B). The pedestrian carries a backpack containing an antenna, the NovAtel SPAN system and LCI IMU. The cable from the backpack is connected to the NI front-end.

One of the benefits of the proposed direct vector processing is its autonomous weighting by power (equivalent to C/N 0 ). Ideally, if it is possible to get accurate enough C/N 0 estimates in the conventional HSGPS, the results can be very close to the proposed approach. However, there are some difficulties for C/N 0 estimation in weak signal and multipath conditions. As shown in Figure 6, C/N 0 values are fluctuating during the indoor periods. This is caused by the user-satellite dynamic and multipath (or fading) phenomenon. The other fact is that, as the signal is weak along with multipath, the C/N 0 estimator will itself exhibit a larger variance or even be biased. By using these erroneous C/N 0 weights, the conventional HSGPS may perform even worse.

The sky plot of all available GPS satellites is shown in Figure 7.

To begin the test, the pedestrian walks in a circular path outside the building in order to align the inertial system. Then, the user walks through the inside of the building with periodic returns outside in order to maintain an accurate reference solution (details below). Finally, circular motion is repeated again at the end of the test in open sky scenario in order to facilitate backward processing of the data.

The raw IF data was processed using the GNSRx-ss™ software receiver [START_REF] O'driscoll | Dynamic duo-combined GPS/GLONASS receivers in urban environments[END_REF] in order to get the raw correlator values over a pre-defined Doppler search range from which Doppler MLE can be easily obtained. In the GSNRx-ss™ receiver, assistance information such as raw data bits, nominal trajectory and broadcast ephemerides are also provided as input. At each measurement epoch, the reference trajectory is used to compute nominal pseudorange and Doppler values for each satellite in view, which are then passed to the signal processing channels. Each channel then computes a grid of correlators around the nominal code phase and Doppler values. Finally, the maximum correlator outputs are used to generate pseudorange and Doppler measurements. In practice, a reference solution would not be available to control the center of the search range. The major impact of this is that the search range in the Doppler domain would likely have to increase in order to accommodate the resulting velocity errors and resulting in an increased computation load. As long as the searching space increases, the same estimation performance has been demonstrated. GSNRx-ss™ with a coherent integration of 500 ms is used. With this coherent integration time, it is expected that a desirable pre-detection SNR can be obtained. After projecting the correlator outputs onto the velocity domain, the velocity powers are computed using Equations ( 8) and [START_REF] Salychev | Low Cost INS/GPS Integration: Concepts and Testing[END_REF]. In Figure 6, the C/N 0 profiles of abovementioned two indoor scenarios are also shown. In scenario A, one can see that the power of three satellites is relatively stronger than the others, being in the range of 20 to 37 dB-Hz. The other five satellites have very weak signals in the range of 10 to 15 dB-Hz. Given the sky plot shown in Figure 7 and the picture of scenario A in Figure 6, one can assume that the satellite in the south part of the sky can penetrate the south-facing windows, resulting in larger C/N 0 values, such as PRN 15,18 and 29. In scenario B, one can see that C/N 0 values of six satellites increase to between 15 and 35 dB-Hz during epoch 40 to 100. Referring to the sky plot in Figure 7 and the picture in Figure 6 one can see that the signals transmitted by the satellites from the south and east portions of the sky penetrate the windows more easily, such as PRN 29 PRN 15.

Velocity MLE

Once correlator outputs from the receivers are available, the MLE velocities are computed first. The resulting velocity domain power distributions in the two indoor environments are depicted in Figures 8 and9. Vertical axes are the velocity offsets with respect to the reference solution computed by using the NovAtel SPAN TM system. Horizontal axes are the time axes. The white circles represent the maximum powers at a given epoch. And the color of any one pixel in the image represents instantaneous received signal power across all satellites in that candidate velocity solution. Due to the fact that the centre of the bins is aligned with the user's true velocity, if the LOS signals are dominant, the largest power should be located at the centre. In scenario A, it can be observed that around epoch 15 s, the total received power is the greatest and the dominant power is located very near the reference solution; the LOS signal appears to dominate the signal. However in the succeeding epochs, NLOS signals seem to dominate the received signal. The error statistics of the proposed direct vector and conventional HSGPS velocity solutions for this scenario are summarized in Table 1. One can see that the proposed algorithm outperforms the conventional one by 75%, 75%, and 89% in the east, north and up velocity axis, respectively. Since the multipath statistics are environment-dependent, it is useful to show that the proposed algorithm works in various indoor environments. In Figure 9, the distribution of the total received signal power over velocity domains is plotted. It can be observed that the power is concentrated and consistent from time 20 to 40 s. Along with C/N0 plots, it is likely that the subject is approaching the windows during this period. The RMS velocity errors of the proposed direct vector and conventional HSGPS are summarized in Table 2. One can see the proposed algorithm outperforms the conventional one by 24%, 72%, and 83% in the east, north and up velocity axis, respectively. Among all three navigation solutions, it can be seen that the navigation solution with the proposed Doppler measurements integrated with PDR is the nearest to the reference trajectory, which validates the effectiveness of the proposed method. In order to further show the error characteristics, the position and velocity errors are plotted as a function of time in Figures 11 and12 user displacements along the east-west axis are much larger than the displacements along north-south direction for the whole trajectory. Thus the position errors in north direction are not as easily observed as east direction. Similarly, for the east velocity errors, it is observed that the conventional HSGPS/PDR integrated solution also degrades performance. This is because the Doppler MLEs in the indoors are not as reliable and may even bias the results. However, the direct vector HSGPS/PDR integrated navigation solution shows that both the position and the velocity RMS errors have noticeable improvements as compared to PDR-only solution. For example, the east, north and RMS position errors of the proposed algorithm have improvements of 24% and 60% as compared to PDR-only solution. Similarly, the north velocity accuracy has an improvement of 13% as compared to PDR-only solution. One can also observe that there is 0.03 m/s degradation in the east velocity when integrating direct vector HSGPS with PDR. However, the direct vector HSGPS/PDR velocity still outperforms conventional HSGPS/PDR velocity by 16%, and 9% in the east and north axes, respectively. These results suggest that direct vector HSGPS Doppler actually reduces such errors on average, which, as discussed below, improves heading determination. For this data set, the benefit from direct vector HSGPS Doppler is more significant in the north direction. To further illustrate the improvement, the mean velocity errors are shown in Table 4. It can be seen that the proposed method generally outperforms the PDR only solution in both the east and north axis, and the biasness of east velocity has a slight improvement with respect to conventional approach. In Figure 14, the estimated headings are plotted. It can be seen that the error of the proposed method is generally smaller than with conventional HSGPS/PDR and does not exhibit the long-term drift seen in the PDR-only heading.

Conclusions

In this paper, the combined indoor signal and multipath model is first introduced. How multipath affects conventional HSGNSS Doppler measurements is then discussed. After that direct vector receiver architecture is proposed and compared to conventional high sensitivity GNSS receivers. The conventional method performs well in most good scenarios, and the proposed method is the maximum likelihood extension to the conventional method, which will asymptotically approach the performance bound. How velocity and Doppler measurements are obtained with direct vector processing performing in indoor multipath environments is also investigated. In order to evaluate the benefits of Doppler estimated with the proposed direct vector processing method over conventional block processing method, Doppler estimated with both approaches are tightly integrated in a detailed navigation filter, which follows a PDR strategy. Comparisons are made between PDR only solution and integrated solutions. From the experimental results, the following conclusions can be drawn:  In indoor environments, the direct vector processing makes the best use of LOS and NLOS signal power for each satellite. Additionally, the mutual information between satellites is further used before arriving at a solution.  Velocity estimation by using direct vector HSPGS outperforms the conventional HSGPS for two indoor scenarios considered herein. The HSGPS Doppler measurements obtained by using direct vector processing are thus more reliable and helpful than conventional HSGPS Doppler measurements for heading estimation.  In weak signal and multipath conditions, the benefit of using conventional HSGPS Doppler measurements is very limited due to the fact that proper weighting cannot be correctly set. However, the proposed direct vector processing method autonomously weighs the navigation solution according to power (equivalent to C/N 0 ), which is actually the maximum likelihood estimate of the navigation solution, and performs no worse or better than the conventional approach.

Figure 1 .

 1 Figure 1. Conventional high sensitivity GNSS receiver architecture.

Figure 2 .

 2 Figure 2. Direct vector GNSS receiver architecture.

Figure 3 .

 3 Figure 3. Power projection from Doppler domain to velocity domain (2D case).

  from user velocity to LOS Doppler offset array ( D f

Figure 5 .

 5 Figure 5. Primary equipments used in data collection. SPAN TM LCI IMU, CPT, and SPAN TM receivers are shown on the left; NI front-end is shown on the right.

Figure 6 .

 6 Figure 6. Data collection actual trajectory, indoor scenarios, and corresponding C/N0 plots (one epoch is 0.5 s)-scenario A picture taken facing east, scenario B picture taken facing south.

Figure 7 .

 7 Figure 7. Sky plots of all available GPS satellites.

Figure 8 .

 8 Figure 8. Power distributions in east, north and up velocity domains (top left, top right and bottom)-Scenario A.

Figure 9 .

 9 Figure 9. Power distributions in east, north and up velocity domains (top left, top right and bottom)-Scenario B.

Figure 10 .

 10 Figure 10. Trajectories of various navigation solutions.

  .

Figure 11 .

 11 Figure 11. Position errors of various navigation solutions.

Figure 14 .

 14 Figure 14. Interpolated heading errors.

Table 1 .

 1 Velocity RMS errors-indoor scenario A.

	Velocity RMS Errors	Direct	Conventional	Improvement
	East (m/s)	0.14	0.56	75%
	North (m/s)	0.32	1.29	75%
	Up (m/s)	0.52	4.96	89%

Table 3 .

 3 Position and velocity RMS errors.

		Position RMS Errors (m)	Velocity RMS Errors (m/s)
		East	North	Up	East	North	Up
	PDR-only	33.65	93.18	2.23	0.38	0.47	0.16
	PDR+HSGPS (conv.)	75.71	35.41	7.55	0.49	0.45	0.20
	PDR+HSGPS (dir.)	25.52	37.45	5.53	0.41	0.41	0.20

Table 4 .

 4 Mean velocity errors.

	Velocity

Mean errors (m/s)

		East	North	Up
	PDR-only	0.02	-0.06	0.00
	PDR+HSGPS (conventional)	0.02	0.01	0.00
	PDR+HSGPS (direct)	-0.01	-0.01	0.00

Acknowledgments

Financial support from Research In Motion (RIM), the Natural Science and Engineering Research Council of Canada, Alberta Advanced Education and Technology and the Western Economic Diversification Canada is acknowledged.

Comparing results of the above two scenarios, two major phenomena are observed. First, the accumulated power fluctuates more from epoch to epoch while indoors. Second, the location of the peak power in the velocity domain largely depends on the multipath statistics or environments.

HSGPS/PDR Integration Performance

In direct vector processing, the velocity MLE is first obtained, and then the corresponding Doppler measurements are computed. It is then convenient to assess the performance of the HSGPS/PDR tight integration using various types of observation. In the following paragraphs, the performance of the PDR only navigation solution and HSGPS/PDR tight integration with conventional and proposed Doppler measurements is assessed.

Here tight integration uses only Doppler measurements to update the user velocity and heading. No pseudorange measurements are used. In this way, the accumulating errors caused by inaccurate Doppler measurements will become more apparent. The other reason to choose HSGPS/PDR integration with conventional and proposed Doppler measurements is to make the results comparable. If the velocity MLE is directly integrated with PDR, it is then a loose integration scheme. This might obscure the benefits of the measurements and integration schemes. By using the same integration scheme, the benefits of better Doppler measurements will become evident.

For the DR algorithms, the step event is first detected by using pattern recognition techniques with a MEMS accelerometer [START_REF] Godha | Integrated GPS/INS System for Pedestrian Navigation in A Signal Degraded Environment[END_REF] and the step length is calculated by using an averaged reference velocity. In this way, the step length estimation error will be very small. The methods of step length estimation have reached a good and consistent level of accuracy, especially with foot and belt mounted sensors, thus it is not a concern in this paper. Only the long term heading accuracy is the focus of the paper. Regarding the PDR heading information, only one tactical grade vertical gyroscope rigidly attached to the backpack is used to provide heading angles with accumulated errors. And such heading angles are fed directly to the PDR filter.

The trajectories of PDR only and HSGPS/PDR tight integration with conventional and proposed Doppler measurements are shown in Figure 10. It can be seen that the PDR only solution is much smoother, which shows better short-term accuracy. However, as time evolves, the PDR only solution is rotated due to the large heading drift and is biased by hundreds of metres. When integrating PDR with conventional HSGPS, the shape of the trajectory is severely distorted. On the other hand, integrating PDR with the proposed approach gives a solution which alleviates the heading drift and better preserves the shape of the trajectory as compared to the other two cases. To further analyze the results, the cumulative histograms of the horizontal and vertical position errors are given in Figure 13. Two key things are worth noting. First, the errors are clearly non-Gaussian, which is not surprising given the operating environment and the high probability of highly-varying NLOS errors. Second, it can be observed that the proposed method has considerably better horizontal positioning performance compared to the conventional HSGPS/PDR solution. For the vertical position error, both solutions are similar until about 7 m, at which point the proposed method begins to outperform the conventional HSGPS/PDR solution. The corresponding RMS error statistics are summarized in Table 3. As shown, the east position errors of conventional HSGPS/PDR solution are nearly twice that of the PDR-only solution. Regarding the north position RMS errors, it can be seen that the improvement of the proposed solution over conventional approach is not as apparent as east direction. This is very likely due to the fact that the  By integrating the Doppler measurements from the proposed direct vector processing with PDR, the solutions show improvements in east and north velocity estimation of 16% and 9% for the indoor scenarios considered herein as compared to the conventional approach.

Currently only the GPS system in L1 band is used. In an on-going project, the GLONASS signals in L1 and L2 bands will also be included. It can be expected that with increased number of satellites, the proposed Doppler will have a tendency to be closer to the true value under certain circumstances. Various indoor environments will be assessed in order to show the strength and weakness of the methods proposed in this paper. Regarding the integration, other heading sensors, such as magnetometers, and low-cost MEMS gyroscopes will also be included, and the benefits of using HSGNSS Doppler measurements in various indoor environments will be further assessed and analyzed.