ES2348778T3

ES2348778T3 - Procedimiento de determinacion autonoma de niveles de proteccion para posicionamiento gnss basado en residuos de navegacion y en una relacion de confianza isotropica.

Info

Publication number: ES2348778T3
Application number: ES08380133T
Authority: ES
Inventors: Miguel Azaola Saenz
Original assignee: GMV Aerospace and Defence SA
Current assignee: GMV Aerospace and Defence SA
Priority date: 2008-04-30
Filing date: 2008-04-30
Publication date: 2010-12-14
Anticipated expiration: 2028-04-30
Also published as: DE602008002149D1; PT2113786E; SI2113786T1; EP2113786A1; PL2113786T3; EP2113786B1; US20090273515A1; ATE477508T1; US8203482B2

Abstract

Procedimiento de cálculo de un nivel de protección horizontal HPL y/o un nivel de protección vertical VPL, que acotan hasta un nivel de confianza dado 1 - α las componentes horizontales y verticales, respectivamente, del error de estimación de la posición δ de una solución de navegación GNSS basada en mínimos cuadrados con un número m de parámetros estimados, solución de navegación que se basa en n observaciones y en una matriz de observación H de n ×m , y cuya solución navegación también proporciona un vector de residuos r ; en el que el nivel de protección horizontal HPL y el nivel de protección vertical VPL se calculan según las siguientes expresiones: **Fórmula** en las que: r es la norma euclídea del vector de residuos de la estimación por mínimos cuadrados r ; **Fórmula** - EE, NN, EN, NE, EU, UE, UU, VU, UV son las componentes espaciales de la ) 1 (HT 20 matriz DoP ( H t.H) -1 de la estimación por mínimos cuadrados expresada en el sistema de coordenadas local horizontal de la posición estimada, donde los subíndices E , N y U representan las componentes este, norte y vertical, respectivamente, y sus combinaciones corresponden a sus correlaciones y correlaciones cruzadas: **Fórmula** y en las que - k es una relación de confianza isotrópica calculada resolviendo numéricamente la siguiente expresión Ec. 3 que vincula k con n , m y α : en la que: **Fórmula** - Γ es la función gamma de Euler''; y, Κm - indica un espacio vectorial real estándar m -dimensional.

Description

Campo de la invención

La presente invención está planteada en el campo de la navegación por satélite (o GNSS), pero puede aplicarse a cualquier campo en el que se utilice estimación por mínimos cuadrados. La invención responde a la necesidad de límites fiables para el error de posición en una diversidad de aplicaciones de navegación GNSS, que van desde la aviación civil al cobro electrónico (en carreteras), entre otras.

Antecedentes de la invención

El problema de la navegación GNSS es el problema de estimar la posición de un usuario GNSS mediante la información proporcionada por la señal GNSS cuando se recibe por un receptor de usuario GNSS. Hay diversas técnicas navegación GNSS estándar, la más común de las cuales es la navegación absoluta. En la navegación absoluta, el sistema de navegación calcula su posición absoluta sin más información que la contenida en las señales de satélite GNSS, mediante las denominadas medidas de seudodistancia (medidas con ruido de la distancia entre el receptor y los satélites GNSS basadas en la determinación del tiempo de recorrido de las señales GNSS desde los satélites al receptor). Con este fin, es necesario sincronizar el reloj del receptor con el tiempo del sistema GNSS (hacia el que se supone que los relojes del satélite GNSS deben estar orientados con una gran precisión). En otras palabras, en modo de navegación absoluta, el receptor debe estimar la desviación de su reloj para poder estimar su posición. Tanto la posición como la desviación del reloj se estiman normalmente simultáneamente en un proceso de estimación lineal por mínimos cuadrados común. El número de parámetros que deben estimarse en este modo de navegación son, por tanto, cuatro: tres coordenadas de posición más la desviación del reloj. Otras técnicas de navegación GNSS estándar incluyen navegación diferencial y navegación cinemática. En ambos casos, el receptor calcula su posición con respecto a una estación GNSS, aprovechando las observaciones de señal GNSS adquiridas por la estación. Esto elimina la necesidad de una sincronización de reloj de receptor precisa, ya que el receptor puede combinar las observaciones de la estación con sus propias observaciones de tal manera que la contribución de la desviación del reloj del receptor se anula, permitiendo así una estimación de estado de tres parámetros (en lugar de cuatro). Por lo que respecta a la presente invención, no hay diferencia entre las técnicas de navegación GNSS relativa, cinemática o incluso absoluta, en la medida en que están basadas en estimación por mínimos cuadrados. De modo que introduciremos brevemente la técnica de

5

10

20

25

30

2

navegación GNSS por mínimos cuadrados de forma genérica en la que la presencia o ausencia de la desviación del reloj como componente del vector de estimación vector es transparente. Esto es teoría estándar y puede encontrarse en la literatura (véase por ejemplo: “Understanding GPS: Principles y Applications”, Elliot D. Kaplan & Christopher J. Hegarty (editores), 2006), por lo que no profundizaremos en los detalles. En cualquier modo de navegación GNSS el problema de la estimación que ha de resolverse es no lineal, de modo que, puesto que van a aplicarse métodos de estimación lineal por mínimos cuadrados (que es algo estándar en la navegación GNSS), el problema de navegación debe

η

ante todo linealizarse. Póngase que es el vector de estado real de usuario (con tres o cuatro parámetros que dependen del modo de navegación). La linealización requiere un cálculo inicial

ηη0

de que se indicará como , alrededor del cual se diferenciarán las ecuaciones de observación GNSS no lineales. El problema de estimación lineal resultante puede escribirse como:

y = H ⋅ x +ε

[Eq. 1] 15

donde:

x =η −η

El vector x es la innovación de estado 0 , es decir, la diferencia entre el cálculo del η0 η

estado y el estado real , y por tanto es lo que debe estimarse para resolver el problema de navegación.

y

� El vector observación se forma con la diferencia entre las medidas reales (por ejemplo, medidas de seudodistancia en el caso de la navegación absoluta), que se obtienen a

η

partir de la posición definida por el estado real y las medidas (ficticias) que se obtendrían si

η

el receptor estuviera en la posición definida por el estado calculado 0 .

�: El vector de error ε es el vector de los errores de medida (por ejemplo, errores de seudodistancia).

�: La matriz de observación H (en ocasiones también denominada matriz geométrica) es la matriz (jacobiana) de derivadas parciales de la ecuación de observación GNSS no lineal en

η0 η0

el estado cálculo , y por tanto, relaciona pequeñas innovaciones del estado aproximado

η0

con pequeñas innovaciones de las medidas esperadas en . En lo que respecta a la presente invención, la forma real de la ecuación de observación GNSS no lineal o como se deriva la matriz de observación H a partir de la misma, no son cuestiones relevantes, por lo que no profundizaremos en estos detalles (que, por otro lado, son de uso estándar en navegación GNSS por mínimos cuadrados y puede aprenderse en numerosas fuentes de la literatura de GNSS como, por ejemplo, “Global Positioning System: Theory & Applications”, Bradford W. Parkinson & James J. Spilker (editores), 1996). El hecho importante es que H relaciona los vectores de estado innovación, observación y error tal como se indica en la ecuación Ec. 1.

y

ε

Obsérvese que los vectores y tienen tantas coordenadas como observaciones disponibles, por ejemplo, tantas como medidas de seudodistancia simultáneas disponibles en el momento en el que ha de calcularse la posición, en el caso de la navegación absoluta, o tantas como el número de medidas de fase doblemente diferenciadas en el caso de la navegación cinemática. Se supone que hay n observaciones disponibles. Entonces, los

y

ε n xmH

vectores y tienen coordenadas, mientras que tiene coordenadas y tiene el

m

tamaño n × (siendo m tres o cuatro dependiendo del tipo de navegación). Ha de observarse también que el vector y y la matriz H son conocidos, el vector de innovación de estado es el que debe estimarse y el vector de error ε quedará siempre sin conocerse (o de lo contrario se podrían corregir los errores y dejaría de haber errores, lo que es imposible). La estimación por mínimos cuadrados de x viene dada por la fórmula ampliamente conocida:

TT

xˆ =(H ⋅ H )−1 ⋅ H ⋅ y x

Proporciona una estimación ˆ de la innovación de estado x , que a su vez proporciona una ηˆ =η + xˆ η=η + x

estimación 0 del estado 0 . El error de estimación δ es la diferencia entre el estado estimado ηˆ y el estado real η que a su vez es igual a la diferencia entre la innovación de estado estimada ˆ y la innovación de estado real

xx :

δ =ηˆ −η = xˆ − x

Hasta ahora se han descrito los principios básicos de una técnica de navegación GNSS estándar (navegación GNSS por mínimos cuadrados). Ahora se definirá con precisión la noción que constituye el principal problema de la presente invención: el Nivel de protección.

1−α

El nivel de protección (PL) es un límite, hasta un nivel de confianza dado , para el error de la estimación de la solución de posición GNSS, es decir, un límite sobre la relación de tamaño entre el error de posicionamiento o una de sus componentes (por ejemplo, la componente vertical u horizontal) y el vector residual que resulta del propio cálculo de la posición.

1−α

Por tanto, un nivel de protección con confianza para la solución de navegación por mínimos cuadrados descrita anteriormente es u n número positivo PL tal que:

P(δ

≥ PL )≤α En la que P es el operador de probabilidad. Obsérvese que δ no se conoce nunca, y ese es el motivo por el que queremos acotarlo. El concepto de nivel de protección se particulariza normalmente a un subconjunto de las coordenadas del vector de estado; obsérvese que δ tiene tres componentes espaciales (más

la componente de desviación del reloj en el caso de la navegación absoluta). Si el sistema de coordenadas empleado es el sistema de coordenadas referidas al horizonte del observador local en ηˆ (y así se supondrá en lo sucesivo), entonces el vector δ también se expresa en coordenadas de horizonte local, y por tanto, sus tres componentes espaciales representan los

δδ

δ EN U

errores de estimación en las direcciones este, norte y vertical ( , y , respectivamente). No es inusual que una aplicación GNSS particular esté especialmente interesada en un límite

δU

para sólo una de las componentes (por ejemplo, la componente vertical , como en el caso de la aviación civil) o en un subconjunto de las mismas (por ejemplo las componentes

δδ N

horizontales E y , como en el caso del cobro electrónico en carretera).

1−α

Un nivel de protección vertical con confianza para la solución de navegación por mínimos cuadrados descrita anteriormente es un número positivo VPL tal que:

P(δU

≥ VPL )≤α 1−α

Un nivel de protección horizontal con confianza para la solución de navegación por mínimos cuadrados descrita anteriormente es un número positivo HPL tal que:

P(δ H ≥ HPL )≤α

⎡

δ E

⎤

δ=

H

δ

Este concepto de nivel de protección se planteó como el núcleo del concepto de integridad

⎢⎣

GNSS desarrollado por los sistemas de aumentación por satélite (SBAS), tal como el WAAS estadounidense o el EGNOS europeo entre otros, y se ha aplicado específicamente a estos sistemas así como a sistemas de aumentación basados en tierra (GBAS) tal como el LAAS. Tanto SBAS como GBAS son sistemas que proporcionan correcciones en tiempo real e información de integridad aplicable a observaciones GNSS. Aviones civiles equipados con unidades de navegación GNSS especiales pueden recibir estas correcciones y aplicarlas a sus observaciones GNSS, así como la información de integridad, en forma de límites estadísticos (es decir, límites hasta un cierto nivel de confianza estadística) para los errores de observación que quedan después de aplicar las correcciones. Por tanto, la unidad GNSS de abordo puede conseguir una estimación más precisa de su posición (gracias a las correcciones) y, además, puede calcular un nivel de protección, que es un límite estadístico para el error de posición que sigue habiendo (gracias a los límites estadísticos para los errores de observación emitidos por el sistema). Además, se han definido procedimientos autónomos para calcular un nivel de protección (autónomo significa que no dependen de correcciones o de información adicional procedente de un sistema de aumentación tal como un SBAS o GBAS) en el marco de los denominados

⎥⎦

es la componente horizontal del vector de error de estimación δ

N

en la que

.

procedimientos RAIM (Receiver Autonomous Integrity Monitoring, monitorización de integridad autónoma del receptor). El concepto RAIM pretende proporcionar una capa de integridad al proceso de navegación GNSS, implementando técnicas para detectar y aislar medidas incorrectas (es decir, medidas con un error excesivo) junto con los niveles de protección mencionados para limitar estadísticamente el error de estimación de la posición. Tales procedimientos de cálculo PL son, sin embargo, difíciles de justificar desde un punto de vista teórico, ya que se basan en hipótesis que raramente se cumplen en el mundo real. El documento US 7219013 B1 desvela un procedimiento de cálculo de los niveles de protección horizontal y vertical, basado en hipótesis relacionadas con la medición de fallos incorporados en la solución de navegación. La presente invención proporciona un procedimiento robusto y consistente, desde un punto de vista teórico, para calcular de manera autónoma niveles de protección basándose en hipótesis únicas y verosímiles, resolviendo así la debilidad de procedimientos previamente existentes.

Descripción de la invención

La presente invención se refiere a un procedimiento para el cálculo autónomo de niveles de protección para posicionamiento GNSS según la reivindicación 1. Realizaciones preferidas del procedimiento se definen en las reivindicaciones dependientes. Según el procedimiento de la presente invención, se calculan niveles de protección autónomos del receptor estadísticamente consistentes, tanto horizontal como vertical. Son aplicables a una solución de navegación GNSS basada en mínimos cuadrados, y como tal, son aplicables a diferentes técnicas de navegación que van desde el posicionamiento absoluto al diferencial y relativo, y también para múltiples soluciones de constelación. La presente invención proporciona un procedimiento para calcular niveles de protección que es autónomo (como el procedimiento para calcular PL previsto por los algoritmos RAIM), es decir, los niveles de protección se calculan sin información externa sobre la naturaleza estadística de los errores de medida (a diferencia de los PL calculados con GBAS o SBAS, que requieren una fuente externa de información y/o algunas hipótesis sobre el tamaño de la distribución estadística de los errores de medida individuales). No obstante, la presente invención difiere de los procedimientos RAIM en lo siguiente: para compensar la falta de información estadística de los errores de medida (que se proporciona en un sistema de aumentación externo tal como SBAS o GBAS), cuando se calculan los PL los procedimientos RAIM utilizan un modelo estadístico para tales errores de medida. Dicho modelo estadístico normalmente incluye, entre otras, las hipótesis de que cada error de medida individual es una variable aleatoria normalmente distribuida con varianza conocida o con un

límite conocido para esta varianza (y, en algunos casos, también se supone una media cero para tal distribución). El procedimiento de la presente invención no hace hipótesis alguna sobre el comportamiento estadístico de los errores de las medidas individuales, sino que se basa sólo en una hipótesis,

5 que se denominará isotropía, de que tales errores de medida se combinan en un vector de error que puede apuntar en cualquier dirección (del espacio de medida) con la misma probabilidad. En líneas generales, la isotropía puede entenderse como:

•: Los errores de medida son independientes estadísticamente por parejas.

•: Los errores de medida son similares estadísticamente por parejas.

10 Obsérvese que no estamos haciendo suposiciones sobre la distribución estadística de cada error de medida individual, sino sólo del modo en que se relacionan entre sí. Como tal, el procedimiento de la presente invención permite limitar el error de estimación de cualquier proceso de estimación por mínimos cuadrados cuyas observaciones cumplan de manera razonablemente la hipótesis de isotropía.

15 Un primer aspecto de la presente invención se refiere a un procedimiento para el cálculo de un nivel de protección horizontal HPL y de un nivel de protección vertical VPL, que acota hasta un

1 −α

nivel de confianza dado las componentes horizontales y verticales, respectivamente, del error de estimación de la posición δ de una solución de navegación GNSS basada en mínimos cuadrados, según las siguientes expresiones:

HPL = k ⋅

r

⋅ aH

20

,

VPL = k ⋅

r

⋅ aV

En las que:

r

-

es la norma euclídea del vector de residuos de la estimación por mínimos cuadrados r ; a

25 -H es el semieje mayor de la elipse bidimensional definida por las componentes HT )−1

( ⋅ H

horizontales de la matriz DoP como:

EE NN

a =

H

2

imagen1

⎛⎜ ⎝

imagen1

2 2

⎞⎟ ⎠

⎛⎜ ⎝

⎞⎟ ⎠

−

h h

EE NN

h

h h

EN NE

+

2

, y

a

-Ves el semieje mayor de la elipse unidimensional definida por las componentes 30 verticales de la matriz DoP como:

aV = imagen1 hUU

, donde

h , h , h , h , h , h , hh , h

EE NN EN NE EU UE UU NU UN

- son las componentes espaciales de matriz DoP

)−1

(HT ⋅ H

de dilución de precisión de la estimación por mínimos cuadrados expresada en el sistema de coordenadas horizontales local de la posición estimada, donde los subíndices E , NU

y simbolizan las componentes este, norte y vertical, respectivamente, y sus combinaciones corresponden a sus correlaciones y correlaciones cruzadas:

⎡

⎤

h hh

EE EN EU

⎢ ⎢ ⎢⎣

⎥ ⎥ ⎥⎦

T

(H ⋅ H )−1

Space

=

hhh

NE NN NU

h hh

UE UN UU

,

H n × m

-es la matriz de observación de la solución de navegación basada en mínimos cuadrados; y donde

k

10 -es una relación de confianza isotrópica calculada resolviendo numéricamente la

k nm α nm

siguiente expresión Ec. 3 que vincula con , y , siendo y el número de observaciones y el número de parámetros, estimados respectivamente, de la estimación por mínimos cuadrados:

mn−m−2

⎛⎜ ⎝

⎞⎟ ⎠

⎛⎜ ⎝

n

m

⎞⎟ ⎠

−

n

2

(1−)

2

π ⋅(1−α)

2

m ,2 ≤ k

∫∈ℜzz

Γ

⋅

⋅ =Γ

⋅

dz

2

z

2

1+k

[Ec. 3]

15 enlaque -Γ es la función Gamma de Euler; y, -ℜm indica un espacio vectorial real m -dimensional.

m = 4

En el caso en el que , es decir, en la solución de navegación GNSS absoluta, k puede calcularse resolviendo numéricamente la siguiente expresión Ec. 4, que también vincula k con nm

20 ,y α :

4− n 2−n

n − 22 n − 42

2

α= ⋅(1 + k )− ⋅(1 + k ) 2

2 2

[Ec. 4]. Tal como se indicó anteriormente, la presente invención es aplicable a técnicas de navegación GNSS relativa, cinemática o absoluta, en la medida en que estén basadas en estimación por mínimos cuadrados.

25 La estimación por mínimos cuadrados también puede ser una estimación por mínimos cuadrados ponderada, es decir, el vector r no es simplemente el vector de residuos de la estimación por mínimos cuadrados r , sino el resultado de normalizar el vector de residuos

r

0

(que se denominará en este caso para evitar la confusión) con una matriz no singular R ,y H

la matriz de observación no es simplemente la matriz de observación, sino el resultado de

5

10

15

20

25

30

8

H 0

normalizar la matriz de observación (que se denominará en el presente documento para evitar la confusión) con la misma matriz no singular R empleada para normalizar los residuos r

0

. De modo que:

r = R ⋅ r0−1 −

t tt

(H ⋅ H ) =(H ⋅ R ⋅ R ⋅ H ) 1

00

La matriz W = Rt ⋅ R se denomina normalmente matriz de ponderación en la teoría de la estimación por mínimos cuadrados ponderada estándar. La relación de confianza isotrópica k puede obtenerse consultando una tabla, tabla que incluye

k

valores de calculados previamente resolviendo numéricamente la ecuación Ec. 3 para un

nm k

conjunto dado de valores de , y α . Y por tanto también puede obtenerse mediante interpolación entre los valores calculados de k en dicha tabla. El procedimiento de la presente invención también puede comprender adicionalmente, para una generación de medida dada y para un nivel de confianza fijo (que se proporcionan como entrada): -calcular el nivel de protección horizontal HPL y/o el nivel de protección vertical VPL para

m + 1

todas las combinaciones posibles de o más de n observaciones disponibles; -buscar la solución cuyo nivel de protección horizontal o vertical sea el más pequeño de entre todas dichas combinaciones, y tomar tal solución y nivel de protección como salida. Es decir, es posible optimizar el nivel de protección obtenido: se comprueba si cualquier

m + 1

combinación de de las observaciones disponibles produce un nivel de protección, horizontal o vertical, que sea más pequeño que el nivel de protección producido por el conjunto completo de medidas disponibles. En tal caso, se toma como la solución el nivel de protección más pequeño de entre aquellos resultantes de todas estas combinaciones, junto con su posición asociada. Como se indicó anteriormente, el procedimiento de la presente invención está basado en la hipótesis de isotropía. Según la invención, la relación de confianza isotrópica k se define como un número positivo que cumple:

P( s

≥ k ⋅

r

)=α

o, de manera equivalente:

2 22

P(s ≥ k ⋅ r )=α

,

r = y − H ⋅ xˆ

donde r , , es el residuo de la estimación por mínimos cuadrados (tal como se indica en los Antecedentes de la invención, donde se describe la estimación por mínimos

s = H ⋅δ

cuadrados aplicada), y , la imagen a través de H del vector de error de la estimación

δ

.

(En la ecuación anterior y en lo sucesivo, la elevación de vectores al cuadrado, tales como r y s , deben entenderse como que toman el cuadrado de su norma euclídea). Se muestra ahora cómo esta relación de confianza isotrópica k , junto con el vector de HT )−1

( ⋅ H

residuos r y la matriz DoP son medios para obtener los niveles de protección PL

2 22

P(s ≥ k ⋅ r )=α

5 deseados. Se supone que ya se ha hallado un número k que cumple . Si interesa el nivel de protección horizontal (HPL), se define una región de protección horizontal (HPR) como:

⎧

⎫

⎡

⎤

x

E

⎪⎪

⎬

⎢⎢⎢⎢

⎥⎥⎥⎥

⎡

⎤

x

E

N U

[x x xUx ]⋅ HT

22

:

∃

≤

k

⋅

⎨

HPR

=

C

|

H

⎢⎣

⎥⎦

xx

U

r

,

E N C

x

⎪ ⎪⎩

⎪ ⎪⎭

N

⎣

⎦

x

C

Consideramos ahora la componente horizontal δ H del estado error de estimación δ ,

⎡

⎤

δ Eδ=

H

δ

probabilidad de que H esté dentro de la región de protección horizontal es:

⎢⎣

⎥⎦

(δ )

TT 22

δ

α

⋅

≤

⋅

−

P

H

k

= 1

r

N

. Se sabe por construcción que

. Por lo tanto, la

⎛

⎞⎟ ⎟ ⎟ ⎟⎟

⎡

⎤

δ E

⎢ ⎢ ⎢ ⎢⎣

⎜ ⎜ ⎜ ⎜⎜⎠

δ

⎥⎥⎥⎥

⎛⎜⎜

⎞⎟⎟

⎡

⎤

δ δ

⎢⎣

E

N U

[δ xU x ]⋅ HT

22

δ

∈

∃x

≤

⋅

P

HPR

P

|

H

k

⎥⎦

=

x

r

,

U C E N C

⎝

x

N

⎦

⎝

⎠

x

C

Por tanto,

⎛⎜ ⎜ ⎜ ⎜⎜

⎞⎟ ⎟ ⎟ ⎟⎟

⎡⎢⎢⎢⎢

⎤

δ E

⎥⎥⎥⎥

⎛⎜⎜

⎞⎟⎟

⎡

⎤

⎢⎣

δ δ δ δU δ C

E N

[δδ δδ]⋅ HT

22

α

∈

≥

≤

k

⋅

−

P

HPR

P

H

= 1

⎥⎦

r

E N U C

⎝

⎠

N

⎣

⎦

⎝

⎠

15 De manera análoga, una región de protección vertical (VPR) puede calcularse como:

⎧

⎫

⎡

⎤

x

E

x

N

⎪⎪

⎬

⎢⎢⎢⎢

⎥⎥⎥⎥

[x x xUx ]⋅ HT

22

x

U

:

∃

≤

k

⋅

⎨

VPR

=

C

|

H

xxx

r

,

E N E N C

x

U

x

C

⎪ ⎪⎩

⎪ ⎪⎭

⎣

⎦

La VPR cumple la inecuación de probabilidad análoga:

P(δU ∈VPR )≥ 1−α

En el caso horizontal, la región de protección es una elipse centrada, y la norma máxima es el 20 semieje mayor de la elipse. En el caso vertical, la región de protección es un segmento centrado, y la norma máxima es la mitad de la longitud del segmento. En ambos casos, la

HT )−1

( ⋅ H

norma máxima de la región puede obtenerse a partir de la matriz DoP . La matriz

ηˆ

DoP puede expresarse en coordenadas de horizonte local en como:

imagen1

en la que los subíndices E , N y U representan las componentes este, norte y vertical,

5 respectivamente, y sus combinaciones corresponden a sus correlaciones y correlaciones cruzadas, y los puntos representan el hecho de que puede haber otras filas y columnas correspondientes a otros parámetros estimados (por ejemplo, la desviación de reloj en el caso de la navegación absoluta). A los efectos de la presente invención, sólo son relevantes las componentes espaciales de la

10 matriz DoP:

⎡

⎤

h hh

EE EN EU

⎢ ⎢ ⎢⎣

⎥ ⎥ ⎥⎦

T

(H ⋅ H )−1

Space

=

hhh

NE NN NU

h hh

UE UN UU

. En estas condiciones puede mostrarse (por ejemplo, mediante el teorema de multiplicadores de Lagrange) que la norma máxima en la región horizontal (por ejemplo, el semieje mayor de la elipse) puede obtenerse como:

a =

H

2

imagen1

imagen2

2 2

⎞⎟ ⎠

⎛⎜ ⎝

⎞⎟ ⎠

−

h h

EE NN

h

h h

EN NE

+

2

y en la región vertical como:

imagen3

a =

V

Por lo tanto el HPL y VPL deseado puede obtenerse como:

HPL = k ⋅

r

⋅ aH VPL = k ⋅

r

⋅ aV

,

, 20 enlaque k es la relación de confianza isotrópica para un nivel de confianza dado 1−α , y los números n y m de observaciones y parámetros que deben estimarse, respectivamente.

2 22

P(s ≥ k ⋅ r )=α

La condición que ha de mantenerse es o, de manera equivalente:

⎛

⎞

22 2

sr

=

α

≥

⋅

P

k

⎜⎜⎝

⎟⎟⎠

ε 2 ε 2

Para que se cumpla la ecuación anterior (que se refiere a una medida de probabilidad) es 25 necesaria cierta información sobre la distribución de probabilidad de los elementos implicados.

Para ello, se supone una única hipótesis sobre la que se basa todo el procedimiento, que es la

ε ε~ = ε

k

hipótesis de isotropía (que da su nombre al valor que debe obtenerse): el vector

(n −1)

define un punto aleatorio uniformemente distribuido en la esfera unitaria -dimensional (es decir, en la esfera unitaria del espacio n -dimensional).

5 En otras palabras, la hipótesis de isotropía afirma que, dentro del espacio n -dimensional en el que se encuentra, el vector ε puede apuntar en cualquier dirección con la misma probabilidad. Puede observarse que:

r = y − H ⋅ xˆ = H ⋅ x +ε − H ⋅ xˆ =ε − H ⋅(xˆ − x)=ε − H ⋅δ =ε− s ε= r + s

, es decir, . Por otro lado, r y s son ortogonales:

−1

T T

T T TTT

s ⋅ r =(xˆ − x)⋅ H ⋅(y − H ⋅ xˆ)= (xˆ − x)⋅ H ⋅(I − H ⋅(H ⋅ H )⋅ H )⋅ y =

10

=(xˆ − x)T ⋅ HT ⋅ y −(xˆ − x)T ⋅ HT ⋅ y = 0 rs ε

Por tanto, y definen una descomposición ortogonal de , y por tanto puede

2 2 2

r + s =ε

invocarse el teorema de Pitágoras: . Definiendo:

ε rs

ε~ = ~ = ~ =

r

s

ε ε ε

;

2 2 2

~ ~~~

15 De ello se deriva que r + s = 1 y que ε= 1. Además, por la hipótesis de isotropía, εes un vector unitario aleatorio que se distribuye de manera informen la esfera unitaria del espacio de medida n -dimensional. Ahora, la condición que queremos que se cumpla puede escribirse como:

~2 2~2

P(s ≥ k ⋅ r )=α

~2 ~ ~

r + s 2 = 1 r

20 Puesto que , puede eliminarse de la condición:

2 ⎛ k 2 ⎞

~2 ~2~2

α= P(s ≥ k ⋅(1− s )) = P⎜ s ≥ ⎟

⎜2 ⎟

1+ k

⎝⎠

k 2

2

ρ= imagen1 1+ k ρ ∈[0,1)

Digamos que . Obsérvese que . Nuestra condición se traduce en:

~

P( s

≤ ρ)= 1−α

[Ec. 2]

25

ρ α

Nuestro objetivo ahora es obtener una fórmula (más explícita) que relacione y .

ρ ∈[0,1)

Calculemos el lado izquierdo de la ecuación anterior para un genérico.

Sd (τ )

Supongamos que denota la esfera de radio τ y dimensión d (genérica) e supongamos V

30 que denota el operador volumen que corresponde a la dimensión de su operando (por

V (S (τ ))= 2π ⋅τ V (S (τ )) = 4π ⋅τ 2

ejemplo, 1 es el perímetro del círculo y 2 es la superficie de la 2-esfera, que es la esfera en el 3-espacio).

~

Puesto que εse distribuye uniformemente en la esfera unitaria del espacio n –dimensional (es decir, en la esfera Sn−1 (1) de dimensión n −1 y radio unitario), la probabilidad de que ε~ se

D

n−1 n −1

encuentre en cualquier región de la esfera puede calcularse como el ( )-volumen de

n −1

esa región normalizado con (es decir, dividido por) el ( )-volumen total de la esfera:

~

∈D

∫ε

~

dε

~

(ε )

−1

∈

n

P

D

=

−1

V (S () 1 )

n n−1

~~

s

El vector es simplemente la proyección ortogonal de εen un subespacio m -dimensional, que es la imagen del subespacio a través de H del espacio de estado m -dimensional (donde

m

toma el valor 3 o 4 dependiendo del tipo de técnica de navegación que se esté utilizando) de modo que:

~~

P( s

≤ ρ)= P(ε ∈ D )

n−1

donde

D

n

−1

=

⎧⎨⎩

⎫⎬⎭

=

⎧⎨⎩

⎫⎬⎭

⎡

⎤

⎡

⎤

~

r

~

r

()1 :0 ≤ 2 ≤ 2

222 2

~

ρ

∈

S

:0

≤

−

1

⎢⎣

⎥⎦

⎢⎣

⎥⎦

,

r

=

s

n

−1

~

s

~

s

⎫⎬⎭

Consideremos la siguiente región de la imagen del espacio de estado m -dimensional:

2

D ={~ s :0 ≤ ~ s ≤ρ 2 }

m

D

La región n−1 puede describirse como:

)

⎧⎨⎩

⎡

⎤

~

(

r

2

~~

~

∈

−

D

S

1

:

⎢⎣

⎥⎦

=

s

r

s

,

−1 −1−

n m n m

~

s

Por tanto, puede calcularse:

1 imagen1 ~2

imagen1 2 ⋅V (S ( 1− s ))⋅

n−5

~1−

∫ε D

~

∈

∫Ds

~

∈

~

ε

~

ds

d

=

n−1 m

s

1 1− ~ s 2

donde el factor

imagen1 representa la deformación del volumen debido a la proyección desde

DD

n−1 mn > m

una región esférica a una plana . Obsérvese que se asume implícitamente que (es decir, que hay más observaciones que parámetros a estimar). Esto se corresponde con la idea intuitiva de que la redundancia de las observaciones es necesaria para proporcionar los niveles de protección. De la ecuación anterior se desprende que:

2

imagen1 1 V (S ( imagen1 1− ~ s ))⋅ ds

n−1−m

~ ~

P( s

≤ ρ)= ⋅

∫~∈ m~

sD 2

1− sV (S ()

1 )

n−1

(d − 1)

La fórmula general para el volumen de una -esfera con radio τ es:

d

2⋅π 2 d −1

=

V (S ()

τ)

⋅τ

d −1

⎛⎜ ⎝

d

⎞⎟ ⎠

Γ

2

donde Γ es la función gamma de Euler. Por lo tanto,

⎛⎜ ⎝

⎞⎟ ⎠

n

imagen1 2 m Γ

V (S ( 1 − ~ s )) −

n−1−m

= 2 ⋅

−−1

n m

2

(1 )

~

π

⋅

−

2

s

V (S ()

1 )

⎛⎜ ⎝

⎞⎟ ⎠

−

n

m

−1

Γ

n

2

y por tanto se obtiene:

⎛⎜ ⎝

⎞⎟ ⎠

n

m Γ n−m−2

−

2

(1 )

(

)

2

∫~∈D

s

~

ds

2

⋅

ρ

=

π

≤

⋅

−

⋅

P

2

s

⎛⎜ ⎝

n

⎞⎟ ⎠

−

m

Γ

2

Dm ρ

Recuperando las definiciones y y sustituyendo la ecuación Ec. 2, se obtiene la ecuación que relaciona k , α , n y m :

−−

m

2

⎛⎜ ⎝

n

⎞⎟ ⎠

⎛⎜ ⎝

n

m

⎞⎟ ⎠

n m

−

2

(1−)

2

π ⋅(1−α)

2

1+ k

ℜm

∫z

Γ

⋅

⋅ =Γ

⋅

dz

2

z

2

[Ec. 3]

∈ℜ m , z 2 ≤ k

2

donde

2

10

m -dimensional (es decir, el espacio de

indica el espacio vectorial real estándar

m

coordenadas reales). Obsérvese que la integral en el lado izquierdo de la

vectores con

ecuación Ec.3 es una integral multivariante (de hecho, una integral m -variante). 15 La relación de confianza isotrópica k no puede despejarse a partir de la ecuación Ec. 3 como nm

una función de , y α . Sin embargo, la ecuación puede resolverse numéricamente en k nm

dados los valores de , y α . A efectos de eficacia del cálculo, tales soluciones numéricas

m

pueden calcularse previamente y tabularse para un conjunto de valores conveniente de n ,y α , conveniente para el contexto al que va a aplicarse el presente procedimiento de cálculo PL,

20 evitando así la necesidad de resolver la ecuación Ec.3 en aplicaciones en tiempo real. La resolución numérica de la ecuación Ec.3 se sustituiría entonces en tiempo real por la sencilla consulta o interpolación en la tabla.

Descripción detallada de las realizaciones preferidas

25 El procedimiento de la presente invención calcula niveles de protección, horizontal y vertical,

1 −α

que acotan, hasta un nivel de confianza dado , las componentes correspondientes del error de estimación de la posición de una solución de navegación GNSS basada en mínimos cuadrados. Y para ello, según una realización preferida:

� Recibe como entrada (o por configuración) el valor o valores del parámetro α que 1 −α

definen el nivel o niveles de confianza para el que van a calcularse los niveles de protección.

m

5 � Recibe como entrada (o por configuración) el número de parámetros que han de estimarse (por ejemplo, cuatro en navegación absoluta, tres en navegación relativa o cinemática).

n

� Recibe como entrada el número de observaciones empleadas en la estimación por mínimos cuadrados.

10 � Recibe como entrada el vector r de residuos de la estimación por mínimos cuadrados (posiblemente normalizado con una matriz de ponderación).

Ht )−1

( ⋅ H

� Recibir como entrada la matriz DoP de la estimación por mínimos cuadrados (normalizada con la misma matriz de ponderación que el vector r si lo está r ) expresada en las coordenadas de horizonte local de la posición estimada.

15 � Obtiene una relación de confianza isotrópica k resolviendo numéricamente la siguiente m

ecuación, que relaciona k con n ,y α :

−−

m

2

⎛⎜ ⎝

⎞⎟ ⎠

⎛⎜ ⎝

⎞⎟ ⎠

n m

−

n

m

2

(1−)

2

π ⋅(1−α)

20

ℜm

∫zmdonde indicaelespaciovectorialrealestándar -dimensional(esdecir,elespaciode

m

vectores con coordenadas reales). Obsérvese que la integral en la lado izquierdo de la ecuación Ec. 3 es una integral multivariante (de hecho, una integral m -variante).

a

� Calcula el semieje mayor H de la elipse bidimensional definida por las componentes 25 horizontales de la matriz DoP como:

Γ

⋅

⋅ =Γ

⋅

dz

2

z

2

[Ec. 3]

∈ℜ m , z 2 ≤ k

2

1+ k

a =

H

2

imagen1

imagen2

2 2

⎞⎟ ⎠

⎛⎜ ⎝

⎞⎟ ⎠

−

h h

EE NN

h

h h

EN NE

+

2

.

a

� Calcula el semieje mayor Vde la elipse unidimensional (es decir, el intervalo) definido por las componentes horizontales de la matriz DoP como:

a =

V

r

30 � Calcula la norma euclídea del vector de errores residuos r .

� Y, finalmente, obtiene el nivel de protección horizontal HPL y el nivel de protección vertical VPL como:

imagen3

HPL = k ⋅ r

⋅ aH VPL = k ⋅ r

⋅ aV

,

m

Por desgracia no puede despejarse k a partir de la ecuación Ec. 3 como una función de n , y α . Sin embargo, la ecuación puede resolverse numéricamente en k dados los valores de n , m α

y . A efectos de eficacia de cálculo, tales soluciones numéricas pueden calcularse

nm

5 previamente y tabularse para un conjunto de valores conveniente de , y α , conveniente para el contexto al que va a aplicarse el presente procedimiento de cálculo de PL, evitando así la necesidad de resolver la ecuación Ec. 3 en aplicaciones en tiempo real. La resolución numérica de la ecuación Ec. 3 se sustituiría entonces en tiempo real simplemente consultando

o interpolando en la tabla. 10 Por ejemplo, si el procedimiento de la invención ha de aplicarse a la navegación GNSS

m = 4

absoluta ( ), suponiendo una constelación de satélites GNSS de hasta 60 satélites (por ejemplo GPS + Galileo) desde la que un usuario típico nunca verá más de 24 al mismo tiempo, la ecuación Eq. 3 puede resolverse para todos los n entre 5 y 24 y para un conjunto de valores

α α= 10 −ii = 1,2,K,7

de de interés (digamos , ). Si, por ejemplo, nos limitamos a la

α

15 navegación GNSS para aviación civil, pueden considerarse aquellos valores de que corresponden a los requisitos de aviación civil para las diferentes fases de vuelo (por ejemplo

−9 −7

α= 5⋅10 α= 10

para HPL en aproximación sin precisión, para VPL en aproximación con precisión).

m = 4 m = 3

Como ejemplo, considerando los casos y , que corresponden a problemas de

20 navegación GNSS absoluta y relativa/cinemática, respectivamente, la ecuación Eq. 3 se ha resuelto numéricamente para obtener los valores de k para todas las combinaciones de los

−1 −2 −7

n = 5,6,K,15 α= 10 ,10 ,K,10

valores y .

k

Los valores de resultantes se muestran en las siguientes tablas Tabla 0-1 y Tabla 0-2, respectivamente:

25

n \ α: 10-1 10-2 10-3 10-4 10-5 10-6 10-7

5: 14.94429 149.99444 1500.0 15000.0 150000.0 1500000.0 15000000.0

6: 4.299632 14.088958 44.704577 141.416177 447.211073 1414.21277 4472.13179

7: 2.668993 6.187072 13.520371 29.216247 62.984863 135.715594 292.399371

8: 2.026635 3.997127 7.309982 13.109950 23.374947 41.601852 73.999232

9: 1.678141 3.018864 4.986783 8.026302 12.797176 20.330060 32.251152

10: 1.456197 2.469589 3.823047 5.739841 8.510512 12.549471 18.459272

11: 1.300667 2.117498 3.134872 4.486298 6.324878 8.853090 12.347554

12: 1.184573 1.871641 2.682497 3.707618 5.038610 6.788587 9.104306

13: 1.093959 1.689456 2.362703 3.180923 4.204431 5.502947 7.162882

n \ α: 10-1 10-2 10-3 10-4 10-5 10-6 10-7

14: 1.020850 1.548462 2.124406 2.802147 3.624479 4.637763 5.897057

15: 0.960339 1.435688 1.939645 2.516931 3.199782 4.020948 5.017760

Tabla 0-1: valores de k como una función de n y α para

m = 4

n \

α: 10-1 10-2 10-3 10-4 10-5 10-6 10-7

5: 12.6798949 127.318718 1273.23902 12732.3954 127323.955 1273239.55 12732395.5

6: 3.70711279 12.1962823 38.7136918 122.469384 387.296721 1224.74436 3872.98319

7: 2.32180389 5.42740225 11.8789083 25.6776176 55.3602373 119.288334 257.007427

8: 1.7728918 3.53847099 6.49098539 11.6517568 20.7809015 36.9883721 65.7948352

9: 1.47366429 2.68997625 4.46334829 7.1954831 11.4796911 18.2415346 28.9407583

10: 1.28233412 2.21128223 3.44262318 5.18074906 7.6894481 11.3440256 16.6897284

11: 1.14780635 1.90314984 2.83641649 4.07131698 5.74811729 8.05159284 11.2338222

12: 1.04710248 1.68719352 2.4364028 3.37950465 4.60116217 6.20532948 8.32659438

13: 0.96830832 1.52664307 2.15265752 2.90996207 3.85476998 5.05158679 6.58012529

14: 0.90460043 1.40203201 1.94057301 2.57124221 3.33424766 4.27275924 5.43784769

15: 0.85177281 1.30210283 1.77567855 2.31547192 2.95200181 3.71595894 4.64215202

Tabla 0-2: valores de k como una función de n y α para

m = 3

5 La presente invención se refiere a un procedimiento para calcular niveles de protección horizontales y verticales autónomos para posicionamiento GNSS basado en mínimos cuadrados, basándose en residuos de navegación y en una relación de confianza isotrópica. Ha de entenderse que la descripción anterior es una ejemplificación de los principios de la invención y no limita la invención a dichas realizaciones.

10

Claims

REIVINDICACIONES

1. Procedimiento de cálculo de un nivel de protección horizontal HPL y/o un nivel de

1 −α

protección vertical VPL, que acotan hasta un nivel de confianza dado las componentes 5 horizontales y verticales, respectivamente, del error de estimación de la posición δ de una

m

solución de navegación GNSS basada en mínimos cuadrados con un número de

n

parámetros estimados, solución de navegación que se basa en observaciones y en una

n × m

matriz de observación H de , y cuya solución navegación también proporciona un vector de residuos r ;

10 en el que el nivel de protección horizontal HPL y el nivel de protección vertical VPL se calculan según las siguientes expresiones:

HPL = k ⋅

r

⋅ aH

,

VPL = k ⋅

r

⋅ aV

en las que:

r

15 es la norma euclídea del vector de residuos de la estimación por mínimos cuadrados r ;

EE NN

a =

H

2

imagen1

⎛⎜ ⎝

imagen1

2 2

⎞⎟ ⎠

⎛⎜ ⎝

⎞⎟ ⎠

−

h h

EE NN

h

h

h h

EN NE

+

+

+

+

2

2

a

-H se define como:

,

a = imagen1 h

-

a

V

se define como: V UU , en la que

hhhhhhhhh

-EE, NN, EN, NE, EU, UE, UU, VU, UV son las componentes espaciales de la

)−1

(HT ⋅ H

20 matriz DoP de la estimación por mínimos cuadrados expresada en el sistema de E , NU

coordenadas local horizontal de la posición estimada, donde los subíndices y representan las componentes este, norte y vertical, respectivamente, y sus combinaciones corresponden a sus correlaciones y correlaciones cruzadas:

⎡

⎤

h hh

EE EN EU

⎢ ⎢ ⎢⎣

⎥ ⎥ ⎥⎦

T

(H ⋅ H )−1

Space

=

hhh

NE NN NU

h hh

UE UN UU

;

25 yenlasque

k

-es una relación de confianza isotrópica calculada resolviendo numéricamente la m α :

siguiente expresión Ec. 3 que vincula k con n ,y

−−

m

2

⎛⎜ ⎝

⎞⎟ ⎠

⎛⎜ ⎝

n

m

⎞⎟ ⎠

n m

−

n

(1 )

(1 α)

2

∫∈ℜzz

2

m ,2 ≤ k

π 2

Γ

⋅

−

⋅ =Γ

⋅

⋅

−

dz

2

z

2

2

2

1+k

[Eq. 3]

en la que:

-Γ es la función gamma de Euler’; y,

ℜm -indica un espacio vectorial real estándar m -dimensional.

m = 4
2. Procedimiento según la reivindicación 1, en el que cuando en la solución de

k

navegación GNSS absoluta, se calcula resolviendo numéricamente la siguiente expresión

m

Ec. 4 que vincula k con n ,y α :

4− n 2−n

n − 22 n − 42

2

α= ⋅(1 + k )− ⋅(1 + k ) 2

22

[Ec. 4].
3. Procedimiento según cualquier reivindicación anterior, en el que la estimación por mínimos cuadrados es una estimación por mínimos cuadrados ponderada, lo que significa que

H

tanto la matriz de observación por mínimos cuadrados, que se usa para calcular la matriz DoP de dilución de la precisión, y el vector de residuos r por mínimos cuadrados se normalizan con una matriz no singular R .

k
4. Procedimiento según cualquier reivindicación anterior, en el que se obtiene

k

consultando una tabla que incluye valores de previamente calculados resolviendo

m α

numéricamente la expresión Ec. 3 para un conjunto dado de valores de n , y.

k
5.

Procedimiento según la reivindicación 4, en el que se obtiene interpolando entre valores calculados de k en dicha tabla.
6.

Procedimiento según cualquier reivindicación anterior, que, para una generación de medida y para un nivel de confianza dados comprende adicionalmente: -calcular el nivel de protección horizontal HPL y/o el nivel de protección vertical VPL para

m + 1

todas las posibles combinaciones de o más de las n observaciones disponibles; -buscar la solución cuyo nivel de protección horizontal o vertical sea el más pequeño de entre tales combinaciones, y tomar tal solución y nivel de protección como salida.