BRPI0515452B1 - método de determinação de parâmetros de navegação de um objeto em movimento - Google Patents

método de determinação de parâmetros de navegação de um objeto em movimento Download PDF

Info

Publication number
BRPI0515452B1
BRPI0515452B1 BRPI0515452A BRPI0515452A BRPI0515452B1 BR PI0515452 B1 BRPI0515452 B1 BR PI0515452B1 BR PI0515452 A BRPI0515452 A BR PI0515452A BR PI0515452 A BRPI0515452 A BR PI0515452A BR PI0515452 B1 BRPI0515452 B1 BR PI0515452B1
Authority
BR
Brazil
Prior art keywords
adr
measurements
propagation
state vector
computation
Prior art date
Application number
BRPI0515452A
Other languages
English (en)
Inventor
Draganov Alexander
Original Assignee
Exelis Inc
Itt Mfg Enterprises Inc
Itt Mfg Enterprises Llc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Exelis Inc, Itt Mfg Enterprises Inc, Itt Mfg Enterprises Llc filed Critical Exelis Inc
Publication of BRPI0515452A publication Critical patent/BRPI0515452A/pt
Publication of BRPI0515452B1 publication Critical patent/BRPI0515452B1/pt

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42Determining position
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42Determining position
    • G01S19/43Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42Determining position
    • G01S19/45Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement
    • G01S19/47Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement the supplementary measurement being an inertial measurement, e.g. tightly coupled inertial
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/52Determining velocity

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)
  • Navigation (AREA)
  • Burglar Alarm Systems (AREA)

Abstract

processamento de delta de distância acumulado de gps melhorado para aplicações de navegação. técnicas para navegação por gps usadas para a determinação da posição e da velocidade de um objeto em movimento. medições de pseudodistância (pr) e medições de deita de distância acumulada (adr) são feitas no objeto a partir de sinais de gps recebidos. diferenças são computadas entre medições de adr que são separadas por um intervalo de tempo que é maior do que um intervalo de tempo entre medições de adr consecutivas. os parâmetros de navegação (por exemplo, posição, velocidade e relógio) são estimados a partir das medições de pr e das diferenças de adr. as equações de medição de adr estabelecidas aqui são formuladas de uma forma muito mais acurada, de modo que o intervalo de tempo entre as medições de adr usado para a computação de uma diferença de adr possa ser muito maior do que aquele usado para as técnicas atuais de diferenciação de adr em aplicações de navegação por gps. conseqüentemente, as diferenças de adr são mais acuradas, o que se traduz em uma solução de navegação muito mais acurada. além disso, a técnica de diferenciação de adr contribui para encurtar os tempos de convergência do processamento de filtro de kalman, e, desse modo, melhorar a acurácia da navegação de uma nave espacial. também são providas técnicas para extensão destes algoritmos de processamento de adr altamente acurados para integração de aplicações de navegação por gps/imu, onde os dados de imu são usados em um modelo de propagação acurado para propagação do vetor de estado.

Description

PEDIDO RELACIONADO [001] Este pedido reivindica prioridade para o Pedido Provisório U.S. N° 60/610.609, depositado em 17 de setembro de 2004, intitulado GPS ADR Processing for Navigation of a Spacecraft, cuja totalidade é incorporada aqui como referência.
CAMPO DA INVENÇÃO [002] A presente invenção refere-se a tecnologias de navegação de sistema de posicionamento global (GPS) e, mais particularmente, a algoritmos de navegação que usam medições de delta de distância acumulado (ADR) para resultados de navegação mais acurados.
ANTECEDENTES DA INVENÇÃO [003] O GPS é um meio bem estabelecido de navegação para espaçonaves, aeronaves, bem como veículos ou pessoas baseados em terra. Tipicamente, um usuário recebe sinais de satélites de GPS e usa medições de pseudodistância (PR) para determinar sua posição, velocidade e outros parâmetros. O termo usuário refere-se ao objeto cuja posição é computada, com base em GPS ou outros sinais de medição de distância que ele recebe. Se uma aplicação requerer uma navegação de alta acurácia, as medições de PR serão processadas por um filtro de Kalman, ao invés de por um algoritmo de solução pontual. Um filtro de Kalman usa medições de PR e modelos de propagação sofisticados para a estimativa da posição de usuário, da velocidade e de outros parâmetros de vetor de estado. Os modelos de propagação são projetados para computarem o estado do usuário no tempo N, se o estado for conhecido no tempo N-1. Eles permitem que o filtro de Kalman use o conhecimento prévio do estado do usuário para a computação da estimativa atual, desse modo melhorando a acurácia.
Petição 870180135007, de 27/09/2018, pág. 5/85
2/70 [004] Infelizmente, embora os modelos de propagação possam ser bastante acurados, isto não garante a acurácia do estado propagado. O modelo de propagação usa a estimativa de vetor de estado na época prévia como uma entrada, e a acurácia do último pode ser um fator limitante. A acurácia do estado propagado é particularmente vulnerável a erros na estimativa da velocidade do objeto. A velocidade é estimada indiretamente pelo filtro de Kalman, devido a sua correlação com a posição.
[005] A obtenção de uma navegação de alta acurácia está ligada proximamente à estimativa de velocidade de forma acurada. Isto cria uma dependência circular entre acurácia de posicionamento e acurácia de velocidade. Na prática, o processo de estimativa é gradual, onde melhoramentos na acurácia de velocidade e de posição facilitam uns aos outros ao longo de algum tempo. O processo inteiro é referido como a convergência do filtro de Kalman. Em um caso típico, a convergência de filtro para aplicações de espaço leva de várias horas a um dia ou mesmo mais. Isto pode ser um problema, especialmente após uma manobra de espaçonave, onde o filtro tem de reconvergir.
[006] Há espaço para melhoria da acurácia de um processamento de navegação por GPS, em termos das soluções de navegação que são produzidas e do tempo de convergência das computações de filtro de Kalman.
SUMÁRIO DA INVENÇÃO [007] Brevemente, são providas técnicas para navegação por GPS usadas para a determinação da posição e da velocidade de um objeto em movimento. Medições de pseudodistância (PR) e de delta de distância acumulado (ADR) são feitas no objeto a partir de sinais de GPS recebidos. Diferenças são computadas entre medições de ADR que são separadas por um intervalo de tempo que é maior do que um intervalo de tempo entre medições de ADR consecutivas. Parâmetros
Petição 870180135007, de 27/09/2018, pág. 6/85
3/70 de navegação (por exemplo, posição, velocidade e relógio) são estimados a partir de medições de PR e das diferenças de ADR.
[008] Uma analogia para as técnicas descritas aqui é conforme se segue. Uma pessoa viajando em um carro começa um cronômetro conforme o carro começa a viajar. Após algum intervalo de tempo, a pessoa pára o cronômetro e mede a distância que o carro viajou. A cada vez em que o cronômetro é começado e parado um erro é introduzido. Se um período de tempo curto for usado para a medição das distâncias, a estimativa da velocidade do carro provavelmente não será muito acurada. Inversamente, se um carro se mover a uma velocidade constante, e se a pessoa realizar uma medição por um período de tempo longo, a estimativa da velocidade se torna muito mais acurada. Embora a segunda opção seja preferível, nem sempre é praticável, uma vez que a hipótese da velocidade constante pode ser inválida. Mais ainda, a medição da velocidade ao longo de um intervalo de tempo relativamente longo produz a velocidade média por este intervalo de tempo. A velocidade média é de pouco uso prático, uma vez que a maioria das aplicações requer conhecimento do estado atual.
[009] Estas técnicas descritas aqui combinam a vantagem de medição de uma velocidade de objeto por intervalos de tempo relativamente longos, com a capacidade de se mapearem de forma acurada tais medições para uma velocidade instantânea no tempo de processamento. Isto é realizado pelo processamento de medições de ADR diferenciadas e pela aplicação de uma formulação matemática que mapeia de forma acurada essas diferenças para o vetor de estado atual do usuário.
[0010] As equações de medição de ADR estabelecidas aqui são formuladas de uma forma muito mais acurada, de modo que o intervalo de tempo entre as medições de ADR usado para a computação de uma diferença de ADR possa ser muito maior do que aquele usado
Petição 870180135007, de 27/09/2018, pág. 7/85
4/70 para técnicas de diferenciação de ADR atuais em aplicações de navegação por GPS. Conseqüentemente, as diferenças de ADR são mais acuradas, o que se traduz em uma solução de navegação muito mais acurada. Além disso, a técnica de diferenciação de ADR contribui para encurtar os tempos de convergência do processamento de filtro de Kalman e, desse modo, melhorar a acurácia de navegação.
[0011] A dificuldade no uso de medições de ADR é que elas refere-sem a quantidades as quais são acumuladas/integradas por algum intervalo de tempo, ao passo que o filtro de Kalman requer medições instantâneas. Esta dificuldade é eliminada com as técnicas descritas aqui, por uma formulação matemática que mapeia de forma acurada medições de ADR em medições instantâneas. Mais ainda, essas medições instantâneas são primariamente para velocidade e, portanto, provêem um meio direto para estimativa da velocidade do objeto.
[0012] Há muitos benefícios do uso de medições de ADR com as formulações matemáticas descritas aqui. As estimativas de velocidade combinadas com um modelo de propagação imediatamente resultam em um estado propagado mais acurado, o que permite uma redução dramática no tempo de convergência do filtro de Kalman. O uso de medições de ADR além de PR aproximadamente dobra o número total de medições usadas para o processamento. Além disso, as medições de ADR têm ruído muito mais baixa do que medições de PR. Mais medições e ruído mais baixo resultam em acurácia geral melhor.
[0013] Também são providas técnicas para extensão destes algoritmos de processamento de ADR altamente acurados para aplicações de navegação em que um modelo de propagação confiável pode não estar disponível para um movimento de objeto portador. Por exemplo, modelos de programação precisa não são disponíveis para certa aeronave, veículos em terra, pessoas, etc. Contudo, se estes objetos tiverem a capacidade de fazer medições inerciais, tais como aceleração e
Petição 870180135007, de 27/09/2018, pág. 8/85
5/70 taxa de rotação, as medições inerciais podem ser usadas para ajudarem no processamento de diferença de ADR.
[0014] Outros objetivos e vantagens tornar-se-ão prontamente evidentes quando uma referência for feita à descrição a seguir tomada em conjunto com os desenhos associados.
BREVE DESCRIÇÃO DOS DESENHOS [0015] A FIGURA 1 é um diagrama de blocos de um sistema de navegação que inclui satélites de GPS e uma unidade de navegação que reside em um objeto em movimento, e o qual usa medições de PR e de ADR para a computação de parâmetros de navegação, com base unicamente em sinais de GPS ou com base em sinais de GPS e dados de IMU.
[0016] A FIGURA 2 é um digrama que mostra as medições e os parâmetros que são introduzidos em um filtro de Kalman para as computações de navegação, com base em medições de PR e de ADR.
[0017] A FIGURA 3 é um fluxograma que mostra o processo de computação de navegação geral com base em medições de PR e de ADR.
[0018] A FIGURA 4 é um diagrama que representa de forma pictórica como as diferenças de ADR são computadas a partir de duas medições de ADR separadas no tempo por múltiplas épocas.
[0019] A FIGURA 5 é um fluxograma que mostra as etapas de processamento para as computações de navegação apenas de GPS, com base em medições de PR e de ADR.
[0020] A FIGURA 6 é um fluxograma que mostra as etapas de processamento para o uso de dados de IMU para propagação do vetor de estado para computações de navegação de GPS - IMU integradas.
[0021] A FIGURA 7 é um diagrama que contém gráficos de erro de navegação versus tempo para as computações de navegação apenas de GPS representadas pelo fluxograma da FIGURA comparado com o
Petição 870180135007, de 27/09/2018, pág. 9/85
6/70 processamento de GPS convencional.
[0022] A FIGURA 8 é um diagrama que contém os gráficos de erro de navegação versus tempo para os dois primeiros minutos de processamento de ADR, magnificados a partir dos gráficos mostrados na FIGURA 7.
[0023] A FIGURA 9 é um diagrama que contém gráficos de erros de navegação de raiz quadrada da soma dos quadrados dos erros individuais como função do tempo usando-se uma propagação de primeira ordem de dados de IMU para medições de PR apenas.
[0024] A FIGURA 10 é um diagrama que contém gráficos de erros de navegação de raiz quadrada da soma dos quadrados dos erros individuais como função do tempo usando-se uma propagação de terceira ordem de dados de IMU para medições de PR apenas.
[0025] A FIGURA 11 é um diagrama que contém gráficos de erros de navegação de raiz quadrada da soma dos quadrados dos erros individuais como função do tempo usando-se uma propagação de terceira ordem de dados de IMU para medições de PR e de ADR.
DESCRIÇÃO DETALHADA [0026] A FIGURA 1 ilustra um ambiente em que as técnicas de navegação descritas aqui são úteis. Um objeto 10 é mostrado, cuja posição é para ser computada a partir de sinais recebidos a partir de uma pluralidade de satélites de sistema de posicionamento global (GPS) 20(1) a 20(N). O objeto 10 pode ser uma espaçonave, tal como um satélite ou um veículo espacial, um veículo em terra ou no ar (por exemplo, um carro, um caminhão, um tanque, avião, helicóptero, navio, etc.). O objeto 10 tem um receptor 12 que é capaz de receber sinais a partir dos satélites de GPS 20(1) a 20(N), um processador 14 que faz medições a partir dos sinais recebidos e computa a posição, a velocidade, a direção, etc. do objeto, e uma interface de objeto 16 (tais como um teclado e um visor) para o transporte da informação de naPetição 870180135007, de 27/09/2018, pág. 10/85
7/70 vegação para um ser humano, se desejado. Além disso, o objeto pode incluir uma unidade de medição inercial (IMU) 18 que extrai dados de medição para aceleração e taxa de rotação do objeto. O processador 14 realiza computações em medições derivadas a partir dos sinais de GPS recebidos.
[0027] O receptor 14 faz medições de pseudodistância (PR) e medições de delta de distância acumulado (ADR) a partir dos sinais de GPS recebidos. O processador 14 então faz computações de filtro de Kalman sobre as medições de PR e de ADR, usando equações de medição, e atualiza o vetor de estado usando computações de modelagem de propagação ou usando os dados de IMU.
[0028] A FIGURA 1 também mostra que, embora os satélites de GPS 20(1) a 20(N) sejam comumente usados para a determinação da posição, há outras fontes de medição de distância que transmitem sinais de medição de distância. O número de referência 22 é pretendido para a representação destas outras fontes de cálculo de distância, tais como sistemas baseados em terra ou sistemas de acompanhamento de posição internos. Os sinais de GPS devem ser apenas um exemplo do tipo de sinais de cálculo de distância com os quais as técnicas descritas aqui podem ser empregadas.
[0029] Um filtro de Kalman é um estimador linear ótimo que estima um vetor de estado (por exemplo, posição, velocidade e relógio) a partir de múltiplas medições de vários tipos (por exemplo, medições de PR e de ADR). As medições de PR e de ADR são feitas a partir dos sinais de GPS recebidos em uma base repetida, onde cada evento de medição é denominado uma época. O filtro de Kalman produz uma estimativa de variância mínima no sentido de mínimos quadrados. Além do estado de navegação, ele geralmente estima erros no estado de navegação geral. O filtro extrai uma medida da acurácia de sua estimativa de vetor de estado de erro. Esta medida de acurácia é reprePetição 870180135007, de 27/09/2018, pág. 11/85
8/70 sentada por uma matriz de covariância, a qual é uma matriz de segundos momentos centrais dos erros da estimativa de estado.
[0030] Tipicamente, um filtro de Kalman é aplicado em um ciclo, onde o algoritmo alterna entre o processamento de medições disponíveis (com um modelo de medição representado por equações de medição) e propagação da solução do vetor de estado para a próxima época (com um modelo de propagação representado por equações de propagação). O modelo de propagação é um processo que descreve como o vetor de estado muda no tempo. O modelo de medição define a relação entre o vetor de estado e quaisquer medições processadas pelo filtro. A acurácia da solução depende da qualidade e da quantidade das medições, da acurácia de equações de medição e da acurácia de rotinas de propagação. Para problemas não lineares, tal como uma estimativa de um estado de satélite, as equações de medição tipicamente são linearizadas.
[0031] A linearização inclui uma computação de derivadas parciais da quantidade medida com respeito aos componentes do vetor de estado. Uma computação dessas derivadas parciais para medições de PR é direta, e é o fundamento para algoritmos de navegação convencionais. Contudo, uma computação de derivadas parciais para medições de ADR é mais complexa. A dificuldade primária é a não combinação de tempo: o vetor de estado representa o estado de satélite em um tempo em particular, ao passo que as medições de ADR representam uma acumulação de medições de fase por algum intervalo de tempo. Uma técnica é mapear medições de ADR a partir de um intervalo de tempo até um instante de tempo. Esse mapeamento necessariamente será uma aproximação e introduzirá um erro na solução de navegação, além de erros de ruído de medição.
[0032] Conceitualmente, o filtro de Kalman pesa as contribuições relativas das medições e do comportamento dinâmico do vetor de esPetição 870180135007, de 27/09/2018, pág. 12/85
9/70 tado. As medições e o vetor de estado são ponderados por suas respectivas covariâncias. Se as medições forem não acuradas (grandes variâncias), quando comparadas com a estimativa de vetor de estado, então, o filtro desenfatizará as medições. Quando as medições são muito acuradas (variâncias pequenas), quando comparadas com a estimativa de vetor de estado, o filtro ponderará as medições de forma pesada, de modo que sua estimativa de estado previamente computada contribua pouco para a última estimativa de estado.
[0033] Mesmo uma abordagem ingênua para o processamento de dados de ADR pode produzir equações de medição as quais são mais acuradas do que as medições de ADR subjacentes (isto é, um erro nas equações de medição é menor do que o ruído de medição). Ainda o efeito de um erro aparentemente insignificante em equações de medição na acurácia de navegação pode ser bastante grande. A razão para isso é que este erro está altamente correlacionado no tempo e, portanto, não é removido por filtração.
[0034] Um risco dos erros correlacionados de longa duração em equações de medição pode ser apreciado pelo exame de efeitos similares para a propagação. O processamento de medição e a propagação são dois estágios de um algoritmo de filtro de Kalman, e são executados em cada ciclo. É bem sabido que pequenos erros de aceleração não considerados em um modelo de propagação podem se acumular ao longo do tempo e afetar a acurácia de navegação. Este efeito forçou o desenvolvimento de modelos de propagação sofisticados e precisos.
[0035] É intuitivamente claro que efeitos similares podem ser verdadeiros para equações de medição. De fato, um erro nas equações de medição tenderá a puxar o objeto (por exemplo, uma espaçonave) em cada época no estágio de atualização (processamento de dados) de processamento, e teria um efeito similar àquele a partir de uma
Petição 870180135007, de 27/09/2018, pág. 13/85
10/70 aceleração não modelada no estágio de propagação. O efeito desse erro sobre a solução de navegação pode ser substancial. Neste caso, usar medições com erros de longa duração pode ser pior do que não usar medições de ADR de forma alguma. Assim, o uso de medições de ADR em navegação é altamente dependente da capacidade de mapeá-las para o vetor de estado sem ou com poucos erros correlacionados de longa duração. Esta exigência de acurácia pode ser quantificada conforme se segue: erros correlacionados nas equações de medição devem estar abaixo ou ser da ordem de erros nas equações de propagação. Neste caso, usar equações de medição para ADR mais provavelmente tornará melhor a acurácia.
[0036] Isto é possível definitivamente para navegação relativa, por exemplo, vôo de formação de satélite. Uma diferenciação e uma diferenciação dupla removem erros sistemáticos em relação a um valor de referência e permitem uma acurácia extraordinária em soluções de navegação relativa. Contudo, é menos evidente como desenvolver um algoritmo acurado para aplicação de medições de ADR para aplicações de navegação absoluta. Isto pode ser a razão para não se usarem medições de ADR para navegação absoluta de uma espaçonave em muitas aplicações prévias.
[0037] Um aspecto importante das técnicas descritas aqui é uma formulação matemática, a qual permite um processamento de medição de ADR acurado para navegação absoluta de uma espaçonave. Este processamento de ADR é complementar a e usado em conjunto com um processamento de PR convencional. O recurso primário de distinção do algoritmo a seguir a ser descrito e sua acurácia. As equações de medição contêm erros não contabilizados, mas a magnitude destes erros é reduzida à mesma ordem que aquela de erros de propagação não modelados.
[0038] Com referência à FIGURA 2, uma representação pictórica
Petição 870180135007, de 27/09/2018, pág. 14/85
11/70 de alto nível do processamento de navegação é mostrada. O objeto em movimento usa os sinais de GPS recebidos para a derivação de medições de PR e diferenças de ADR (a partir de duas medições de ADR separadas no tempo) usando técnicas descritas abaixo em conjunto com a FIGURA 3. As medições de PR e as diferenças de ADR são supridas para o filtro de Kalman representado pelo bloco 70. Também são supridos para o filtro de Kalman um modelo de propagação, se conhecido, para o objeto 10 e/ou dados de IMU, se o objeto 10 incluir uma IMU. Os dados de IMU são úteis em circunstâncias em que computações de navegação de GPS - IMU integradas devem ser realizadas, tal como quando um modelo de propagação confiável não é conhecido para o objeto 10.
[0039] Voltando-nos para a FIGURA 3, as etapas básicas para o algoritmo de navegação apenas de GPS ou de GPS - IMU serão descritas. Na etapa 80, o objeto 10 recebe sinais de cálculo de distância, tais como sinais de navegação por GPS a partir de satélites de GPS ou sinais de cálculo de distância a partir de outras fontes de cálculo de distância, conforme mostrado na FIGURA 1. Em seguida, na etapa 82, o objeto faz medições de PR e de ADR a partir dos sinais de cálculo de distância recebidos na etapa 80. Na etapa 84, é computada uma diferença entre duas medições de ADR em duas épocas (por exemplo, instantes de tempo) que são separadas por um intervalo de tempo que é maior do que um intervalo de tempo entre medições de ADR consecutivas. Esta técnica de diferenciação de ADR é descrita em maiores detalhes em conjunto com a FIGURA 4. Na etapa 86, as diferenças de ADR são processadas usando-se equações derivadas parciais modificadas no filtro de Kalman, e o vetor de estado é propagado usando-se o modelo de propagação (e as medições de PR) ou pela propagação do vetor de estado usando-se dados de IMU obtidos pelo objeto. As equações de derivada parcial modificadas são descritas em maiores
Petição 870180135007, de 27/09/2018, pág. 15/85
12/70 detalhes aqui adiante em relação com a FIGURA 5. As técnicas para propagação do vetor de estado usando-se os dados de IMU são descritos em conjunto com a FIGURA 6. Novamente, as técnicas da FIGURA 6 são úteis em aplicações de navegação em que o objeto cuja posição é para ser medida não tem uma aceleração previsível e, assim, um modelo de propagação confiável não é conhecido.
I. Equações de Processamento de Medição de ADR para Navegação por GPS Apenas e GPS - IMU Integrados
A. Cancelando Erros Sistemáticos em Relação a um Valor de Referência [0040] Os dados de ADR têm um erro sistemático em relação a um valor de referência constante desconhecido (se não houver lapsos de ciclo). Devido a este erro sistemático em relação a um valor de referência, as medições de ADR não podem ser aplicadas da mesma forma que as medições de PR. Um método bem conhecido de se lidar com o erro sistemático em relação a um valor de referência é estimálo, por exemplo, pelo uso de ADR para atenuação de medições de PR. Uma técnica que é superior à atenuação é descrita aqui para processamento de medições de ADR e eliminação do erro sistemático em relação a um valor de referência.
[0041] Uma forma geral de exclusão de erro sistemático em relação a um valor de referência é usar uma medição na forma a seguir:
P = J P(t\ W(t) dt (1.1) t1 onde P(t) é a medição de ADR e a função de peso W(t) é t2 selecionada de forma tal que j W (t )dt = 0. Pode-se verificar por substit1 tuição direta que a equação (1.1) cancela um erro sistemático em relação a um valor de referência constante contido na medição de ADR P(t). Os pesos W (t) podem ser selecionados para a minimização de
Petição 870180135007, de 27/09/2018, pág. 16/85
13/70 erros na medição resultante P e para melhoria dos resultados de seu processamento pelo filtro de Kalman.
[0042] Para ilustração dos benefícios de uso da equação (1.1), um caso especial de função de ponderação W(t) é considerado, o qual produz simplesmente uma diferença de valores de ADR P = ap = P(tn)-P(tn1) em duas épocas tn_1, tn. Isto provavelmente é a expressão mais simples (embora não a ótima) que cancela o erro sistemático em relação a um valor de referência, e produz a medição a seguir:
ΔΡ = [fe - rs,n )2 l2 - [(rn-1 - rs, n-1 + fτd tn-1 (1.2) onde rn é o vetor de posição do objeto, rs n é o vetor de posição de um satélite de GPS, τ é a diferença entre a leitura atual e uma de referência em um relógio de objeto, e os subscritos n, n -1 denotam duas épocas (não necessariamente consecutivas). A natureza da equação (1.2) é aquela de velocidade relativa entre o objeto e o satélite de GPS, integrada por algum intervalo de tempo. Em técnicas de processamento de ADR conhecidas, as duas instâncias de tempo tn-1, tn para as quais valores de ADR são diferenciados correspondem a duas épocas consecutivas. Em contraste, nas técnicas de processamento de ADR descritas aqui, a qualidade de medição de ADR pode ser grandemente melhorada se estas duas instâncias de tempo forem separadas por um intervalo de tempo mais longo. Em uma formulação mais geral, isto corresponde a uma faixa mais ampla de t onde a função de ponderação W (t) é não nula. O caso mais simples de uso da equação de medição (1.2) pode trazer benefícios substanciais. Um melhoramento adicional pode ser possível pelo uso de medição na forma da fórmula (1.1).
B. Escolha do Intervalo de Tempo entre Medições de
Petição 870180135007, de 27/09/2018, pág. 17/85
14/70
ADR para a Computação da Diferença de ADR [0043] A equação (1.2) não especifica a duração do intervalo de tempo entre as duas épocas para o qual a computação de diferença de ADR é feita. A escolha da duração deste intervalo de tempo é muito importante para se atingir o pleno potencial do método sugerido abaixo.
[0044] Com base nas considerações qualitativas apresentadas acima, a forma básica aproximada da estimativa de velocidade derivada a partir de duas medições de ADR é:
V = (1.3) ^n-1 n [0045] A descrição a seguir especifica a equação de medição a ser usada no filtro de Kalman. Para as finalidades desta descrição, é necessário apenas saber que o modelo de propagação será usado para formulá-lo, e erros no modelo de propagação contribuirão para os erros na estimativa de velocidade juntamente com o ruído de medição nas medições de ADR. Assim, a extensão do intervalo de tempo entre medições de ADR para computação de uma diferença de ADR é dependente da acurácia do modelo de propagação. Um modelo de propagação mais acurado permitirá um intervalo de tempo mais longo entre medições de ADR para computação da diferença de ADR. De modo similar, dados de IMU mais acurados e o uso de dados de IMU para fins de propagação (conforme descrito aqui adiante em conjunto com a FIGURA 6) permitirão um intervalo de tempo mais longo entre medições de ADR para a computação da diferença.
[0046] Os erros no modelo de propagação õP são assumidos como sendo devido a acelerações não modeladas, e podem ser estimados conforme se segue:
õ - 2Sa-(t. -1,-7 (1.4) onde õa é a aceleração não modelada. Os ruídos de mediPetição 870180135007, de 27/09/2018, pág. 18/85
15/70 ção nas duas medições de ADR são contabilizados e são assumidos como sendo não correlacionados. Suas variâncias são assumidas como sendo iguais e são denotadas por σ}.
[0047] A variância total do erro de estimativa de velocidade pode ser aproximada como se segue:
σ2 = Ρ +1 Sa2 •(í„_1 - tn )4 (t„-1 - C )2 (1.5) [0048] É direto encontrar o mínimo desta variância com respeito à duração do intervalo de tempo. Isto define a duração ótima a ser usada na aplicação desta técnica, conforme descrito aqui:
(t„-1 - )opt = 24 · [0049] Se a qualidade do modelo de propagação for tal que acelerações não modeladas sejam da ordem de Sal « 10-6 m/s2, e se o ruído de ADR for da ordem de 10-2 m, então, a duração ótima do intervalo entre as duas épocas para as quais a computação de diferença de ADR é feita é da ordem de 100 s, onde cada época é de aproximadamente 1 s.
[0050] Com referência à FIGURA 4, é mostrado um exemplo para ilustrar como as diferenças de ADR são computadas. Ao passo que as técnicas atuais de navegação por GPS computam uma diferença entre medições de ADR consecutivas, foi descoberto que uma computação de uma diferença entre medições de ADR separadas por um intervalo de tempo que é maior e, em particular, substancialmente maior, do que o intervalo de tempo entre medições de ADR, pode produzir soluções de navegação muito mais acuradas, e ajudar na convergência do filtro de Kalman mais rapidamente. Por exemplo, a FIGURA 4 mostra que uma diferença de ADR é computada entre uma medição de ADR na época E1 e uma medição de ADR na época E100 (onde cada época
Petição 870180135007, de 27/09/2018, pág. 19/85
16/70 é separada por aproximadamente 1 s), entre uma medição de ADR na época E2 e uma medição de ADR na época E101, entre uma medição de ADR na época E3 e uma medição de ADR na época E102, e assim por diante. Assim, com este exemplo, as diferenças de ADR são computadas para cada nova medição de ADR, conforme ela for recebida. [0051] Não é necessário que uma diferença de ADR seja computada para cada nova medição de ADR com respeito a uma medição de ADR anterior a precedendo em muitas épocas. As diferenças de ADR podem ser computadas entre medições de ADR, separadas por um período de tempo adequado (por exemplo, um número de épocas) que não necessariamente usam toda e qualquer medição de ADR que seja feita. Por exemplo, uma diferença de ADR pode ser computada entre uma medição de ADR na época E1 e uma medição de ADR na época E100, e uma medição de ADR na época E101 e uma medição de ADR na época E200, e assim por diante. Ainda uma outra possibilidade é computar uma diferença de ADR entre a época E1 e uma medição na época E100, entre a época E10 e uma medição na época E111, entre a época E2 e uma medição na época E121, etc. Estas diferenças de ADR são processadas usando-se equações de derivada parcial modificadas descritas abaixo.
[0052] Note que intervalos de tempo pequenos (por exemplo, de 1 s) ou muito grandes (por exemplo, de 104 s) resultariam em erros maiores na estimativa de velocidade. Para a duração ótima, os erros na estimativa de velocidade podem ser estimados pela substituição do valor ótimo para duração na equação (2), para produzir:
Ή * 10 ' m/s (1.6) [0053] Esta estimativa para erro é pequena, se comparado com outros meios de estimativa de velocidade. É particularmente impressionante como resultado do processamento de uma única medição. Obviamente, o potencial para obtenção dessas medições de alta qualidaPetição 870180135007, de 27/09/2018, pág. 20/85
17/70 de pode ser percebido apenas usando-se equações de medição com uma acurácia comparável. As seções a seguir apresentam uma derivação para as equações de medição que são úteis para processamento de diferenças de ADR computadas conforme descrito acima.
C. Mitigação de Erros de Relógio de Objeto [0054] A derivação a seguir assume que um bom modelo de propagação esteja disponível para a posição e a velocidade do objeto (isto é, a posição e a velocidade do objeto 10), o que tipicamente é o caso para objetos de satélite. Infelizmente, o modelo de propagação para o relógio pode ser menos acurado (se um relógio barato for usado) e, portanto, o último termo na equação (1.2) é difícil de se estimar de forma acurada. Se parciais forem computadas para o lado direito da equação (1.2), um dominante poderia ser devido à diferença entre a leitura atual e uma de referência em um relógio, e pode criar problemas para uma estimativa acurada das outras componentes do vetor de estado (principalmente, a velocidade).
[0055] Em um processamento de filtro de Kalman de medições de ADR, a estimativa da diferença entre a leitura atual e uma de referência em um relógio é separada da estimativa de velocidade. Uma separação totalmente limpa não é praticável nem necessária; a meta é se ser capaz de estimar a velocidade na presença de um relógio instável. Foi observado que o último termo na equação (1.2) tem o mesmo valor para todos os satélites em visualização e, assim, pode ser cancelado por medições de diferenciação. Assim, o residual de medição para um par de satélites (p, q) a ser usado no filtro de Kalman é conforme se segue:
(1.7) [0056] Este residual pode ser prontamente computado a partir das medições de ADR disponíveis e a partir das posições conhecidas ou
Petição 870180135007, de 27/09/2018, pág. 21/85
18/70 estimadas de satélites de GPS e do objeto.
[0057] O residual de medição na forma da expressão (1.7) pode ser usado para estimativa da velocidade. Para a diferença entre a leitura atual e uma de referência em um relógio, é usado qualquer uma das medições únicas na forma dada pela equação (1.2).
[0058] Para se usar uma medição no filtro de Kalman, a matriz de parciais deve ser computada. O desafio na computação da matriz de parciais é que um modelo de medição instantâneo é necessário, ao passo que a equação (1.2) é para velocidade integrada por um período de tempo. Esta dificuldade é eliminada pela metodologia de processamento descrita abaixo em conjunto com o fluxograma da FIGURA 5. A meta das etapas de processamento mostradas na FIGURA 5 e descritas abaixo é computar derivadas das diferenças de ADR com respeito a uma ou mais componentes no tempo atual do vetor de estado (por exemplo, as componentes de velocidade e de posição de vetor de estado), para a produção de estimativas de velocidade instantânea que são usadas para a atualização da componente de velocidade do vetor de estado. As medições de PR também são processadas no filtro de Kalman para atualização da componente de posição, e um modelo de propagação conhecido ou de dados de IMU (conforme descrito aqui adiante em conjunto com a FIGURA 6) é usado para a propagação do vetor de estado. Estas técnicas, conforme se tornará evidente aqui adiante, têm uma acurácia de ordem mais alta do que as técnicas de navegação por GPS da técnica anterior com respeito à separação de tempo entre medições de ADR, por exemplo, o erro é de ~x4, onde x é um parâmetro pequeno proporcional à diferença de tempo entre medições, se comparado com o erro ~x de técnicas de navegação da técnica anterior.
[0059] Este método compreende a computação de derivadas parciais da posição de objeto na época passada com respeito ao vetor de
Petição 870180135007, de 27/09/2018, pág. 22/85
19/70 estado na época atual. Esta computação é uma aproximação de ordem mais alta (por exemplo, de 4a ordem) com respeito à diferença de tempo entre épocas, usando-se o modelo de propagação como um agente subjacente.
D. Parciais com Respeito à Posição e Velocidade no Quadro de Referência Inercial [0060] Uma vez que a equação (1.2) é um escalar, ela pode ser formulada usando-se qualquer quadro de referência. Nesta subseção, um quadro de referência inercial é usado. Para a computação de parciais em um quadro fixo na Terra (ECEF), elas são primeiramente computadas no quadro inercial, o qual é definido como tendo a mesma orientação de eixo geométrico que o quadro de ECEF na época n. Os resultados serão usados mais tarde para a computação de parciais no quadro de ECEF.
[0061] Para deixar as equações menos trabalhosas, as parciais são obtidas usando-se a equação (1.2) que correspondem a um dos satélites no par, isto é, computando-se:
(1.8) onde X é o vetor de estado de objeto que compreende a posição, a velocidade, o relógio e a diferença entre a leitura atual e uma de referência em um relógio do objeto.
[0062] Assim, na etapa 100, o vetor de estado pode ser convertido a partir do quadro de referência de ECEF para um quadro de referência inercial. A computação 100 é opcional. Uma transformação para e a partir do quadro de ECEF é opcional.
[0063] As derivadas parciais na expressão (1.8) são computadas apenas com respeito às componentes de posição e velocidade do vetor de estado na época n. Esta divisão no vetor de estado é denotada
Petição 870180135007, de 27/09/2018, pág. 23/85
20/70 como xn = {rn, vn}. As parciais com respeito a relógio e diferença entre a leitura atual e uma de referência em um relógio são computadas separadamente.
[0064] Uma vez que a posição de objeto rn é explicitamente usada na expressão (1.8), ela cria a impressão de que é necessário computar apenas as parciais com respeito a ela. Contudo, a posição de objeto rn é o resultado da propriedade de vetor de estado a partir da época n -1. Assim, rn é uma função de xn-1; daí, xn-1 é uma função (inversa) de rn, e rn-1 (que é uma parte de xn-1) é uma função de rn. Isto resulta em uma computação mais complexa das derivadas:
X(xn - Xs,n ) [AT ( rn-1 ~ Xs, n-1 )] [(rn - rs,n )2 11 [( rn-1 - Xs, n-1
AJ [(r
SXn n - rs,n )2 -[(rn-1 - rs,n-1 )' ]* > = n-1 ' s, n-1 n-1 ' s, n-1 - r n-1 s, n-1 , (1.9) onde õr.
A=Ã.
n-1,i (1.10) e χ é uma matriz 6x6, a qual contém tudo zero, exceto pelos três primeiros elementos da diagonal principal, os quais são iguais a 1. O lado direito da equação (1.9) é uma divisão da matriz de parciais, a qual pode ser usada no filtro de Kalman para processamento de medição.
[0065] Na equação (1.9), a matriz A é usada. Conforme foi explicado acima, esta matriz é devido à dependência entre vetores de estado em duas épocas, a qual, por sua vez, é devido à propagação do vetor de estado.
II. Processamento de ADR Adicional
A. Navegação Apenas por GPS [0066] Para a propagação de um vetor de estado, todas as componentes do vetor de estado são necessárias. Se a componente de
Petição 870180135007, de 27/09/2018, pág. 24/85
21/70 velocidade estiver faltando, uma propagação não pode ser feita. Note que a matriz Λ usa uma posição apenas na época n -1 e ambas a posição e a velocidade na época n. Isto torna necessário propagar para trás, isto é, da época n para a época n-1, ou, em outras palavras, por um intervalo de tempo que separa as medições de ADR para as quais as diferenças são computadas.
[0067] Assim, na etapa 110, o vetor de estado é propagado para trás um incremento usando-se uma rotina de propagação, tal como o esquema de Runge-Kutta. A propagação é governada por equações diferenciais ordinárias. Muitas rotinas de propagação usam o assim denominado esquema de Runge-Kutta para garantir uma alta acurácia. Quando da computação da matriz Λ, uma acurácia de nível alto é mantida usando-se um método de Runge-Kutta de 4a ordem para o algoritmo de propagação subjacente. Este algoritmo de propagação não tem de combinar exatamente com um usado para a propagação do vetor de estado (há vários esquemas de Runge-Kutta igualmente bons). Um exemplo de um esquema de Runge-Kutta útil é representado pelo conjunto de equações a seguir:
dx r.
d = f (íX) Xn+1 - Xn + T(k1 + 2k2 + 2k3 + k4 ) k1 - h · f (t n , xn ) , _ 1 + 2hXn + 7k1 +2h (1.11)
- 1 ,η “ 2 1) xn +1 k2 | n 2 2) + k3 ) + h, Xn onde h é o incremento da propagação. O conjunto de equações em (1.11) é escrito em uma forma canônica, isto é, quando as equações são propagadas a partir da época n para a época n +1. Uma vez que é necessário propagar na direção oposta, h é substituído por (- h) no conjunto de equações (1.11).
Petição 870180135007, de 27/09/2018, pág. 25/85
22/70 [0068] Para a computação da matriz Λ, as derivadas parciais óx —n são computadas:
ÓX„ j ÓXn+1, ÓXn,j = δ· +1_í dki,' + 2 Ók+ 2 dk3,' +Ó6 [õXnp ÓXnj 'X ÓXn,j J
1,i (1.12) n, j n, j n, j onde δν é a função de δ de Kronecker, e óktl ÓXnj ók 2,/ dXnp = h Óf {tn , Xn ) ÓXnj Óf [ tn + 2 h Xn + 2 k1 J h Óxnj (1.13) 'n,j óa = h.
ÓX„ j Óxn j 'Pj xF-kPn-k} ÓXn,j [0069]
Por brevidade, as notações a seguir são introduzidas: Óf t {tn, Xn ) ψ1,j ψ {tn> Xn ) = ψ, 1tn +1 h, xn +1 k | t} [ n 2 n 2 1J = ψ„ | tn +1 h, xn +1 k2 | ,} [ n 2 n 2 2 J = ψ {t n + h Xn + k3 ) ψ {tn > Xn ) =
l.i/' ψ
T 2,j ψ
T3,y
4,í/' H (1.14) [0070] A substituição do conjunto de equações (1.14) no conjunto de equações (1.13) e a computação de derivadas parciais leva ao seguinte:
Petição 870180135007, de 27/09/2018, pág. 26/85
23/70
Sk,.
= h ·Ψ1 ü íx 1 j SXn,J ík íXnj í
· δ· + Sk3, íXnj
Sk íXnj __ 4, im m
Sku ' i 2 ) f +11 í
· Õmj + (1.15) [0071] 2 )
Sk3,m ' δΧη, j ;
Pela substituição de íki. íXnj a partir da primeira equação do conjunto de equações (1.15) na segunda, e pela continuação dessas substituições em cadeia até a última equação, as derivadas k com respeito à xn são obtidas, as quais são expressas através de matrizes
Ψ1_4. Então, todas essas equações são substituídas na equação (1.12) para a obtenção do seguinte:
ÍX 1 j 6V 1 2 34 h2 Ψ 2 Ψ _ + h2 Ψ 3 Ψ 2 + h2 Ψ 4 Ψ 3 + 3 1 3 (1-16
- h3 Ψ 3 Ψ 2 Ψ, + - h3 Ψ 4 Ψ 3 Ψ 2 + 3 2 1 2 4 3 2
- h1Ψ 4 Ψ s Ψ 2 Ψ 11 [0072] Esta é uma fórmula trabalhosa contendo múltiplos produtos de matrizes 6x6. Felizmente, uma simplificação considerável pode ser feita, devido à forma especial das matrizes Ψ.
[0073] Note que xn+1,xn na equação (1.16) corresponde a rn_1i, x na equação (1.9). O vetor rn_1 é expresso no sistema de referência inercial. Mais ainda, é assumido que o movimento do objeto em órbita seja primariamente devido a um campo gravitacional potencial U (t,x). Isto significa que a função f(t,X) na primeira equação do conjunto (1.11) tem a forma a seguir:
Petição 870180135007, de 27/09/2018, pág. 27/85
24/70 fr (t, X ) = V fv (t, X )==^ õr (1.17) isto é, as três primeiras componentes de f (t, X) (as quais são responsáveis pelas derivadas da posição) são simplesmente aquelas da velocidade, e as três últimas componentes (as quais são responsáveis pelas derivadas de velocidade) são funções da posição apenas. Este fato é usado para a simplificação do lado direito da equação (16).
[0074] As derivadas de f (t, X) (conforme necessário para o conjunto de equações (1.14)) são como se segue:
(1.18) onde 0 é uma matriz nula 3χ3, I é a matriz identidade 3x3 e φ é uma matriz 3x3 definida como se segue:
= õ 2U (t, x) y õxt õx.
1 J (1.19)
Os produtos de matrizes Ψ agora podem ser simplificados:
• ! 0 Ί a_l___ _0_! ! A !<Pa\ 0 r
(1.20) onde os subscritos a,b,c,d assumem valores de 1 a 4. A substituição da equação (1.20) na equação (1.16) produz a expressão a seguir:
ÕXn+1 = I + 1 rΦ rr _! Φ rv õXn 6 lΦ vr ! Φ vv J (1.21) onde as matrizes Φ têm as dimensões de 3x3 e são definiPetição 870180135007, de 27/09/2018, pág. 28/85
25/70 das como se segue:
Φ rr
Φ rv
Φ vr Φ vv vv = h 3φλ + h 3φ2 + h 3φ3 +1 h Λφ3φ1 3 1 3 = 6hI + — h φ2 + — h φ3 2 2 (1.22)
3 1 3 = hφ1 + 2hφ2 + 2hφ3 + hφ4 +— h φ3φ1 + — h φ4φ2 = h Ίφ2 + h Ίφ3 + h 2φ4 +1 h 4φ4φ3
Õx [0075] A matriz n na equaçao (1.21) contem produtos de nao ÕXn mais do que pares de matrizes 3x3, se comparado com a expressão original na equação (1.16), a qual continha produtos de ate quatro matrizes 6x6. Note que a dimensão da matriz Λ e de 6x3, e, portanto, apenas duas das quatro matrizes 3x3 são necessárias no lado direito da equação (1.21).
[0076] As estimativas aproximadas apresentadas em uma das seções previas mostram que o intervalo de tempo entre duas epocas associadas à computação de diferença de ADR pode ser bastante longo. Embora derivadas parciais (1.21) sejam computadas para a acurácia de quarta ordem com respeito a h, a acurácia do resultado pode não ser suficiente para intervalos de tempo longos. Uma solução natural para melhoria da acurácia seria empregar vários subincrementos e computar uma derivada parcial (1.21) como um produto de matrizes, as quais correspondem a estes subincrementos.
[0077] As matrizes cpk são a chave para a computação da equação (1.21). Elas podem ser computadas no caso de um potencial gravitacional J2:
= 3(rr μδ1} (J2 R2e )j‘ r5 r3 2
30(r · n)^(nirj + nj.ri
). 6nn r5
15riri 105 (r · nf rir.
15 (r ·n)2 δ (1.23)
Petição 870180135007, de 27/09/2018, pág. 29/85
26/70 onde n) é o vetor unitário na direção z, r = |7|, e re é o raio da Terra.
[0078] Um modelo de propagação freqüentemente inclui outros efeitos além do campo gravitacional da Terra. Alguns destes efeitos causam aceleração de espaçonave (objeto), o que é uma função da posição da espaçonave apenas, por exemplo, uma aceleração devido ao Sol e à Lua, e pressão de radiação solar. Neste caso, a derivação de parciais para processamento de ADR é análoga àquela acima. Em alguns outros efeitos (por exemplo, arrasto atmosférico), a aceleração é uma função de posição e velocidade. Neste caso, não se faz uma simplificação nas equações e se usa a forma completa (1.16), ou se tratam estes efeitos como uma perturbação com uma acurácia de ordem mais baixa.
[0079] Assim, na etapa 120, a matriz cpk é computada para cada um dos subincrementos de Runge-Kutta usando-se a equação (1.23), ou fórmulas para termos gravitacionais de ordem mais alta, além de J2, os quais são providos separadamente abaixo. Após as matrizes !k serem computadas, na etapa 130 elas são substituídas no conjunto de equações (1.22), e os resultados do conjunto de equações (1.22) são usados no conjunto de equações (1.21) para a computação da matriz Λ.
[0080] A matriz Λ é substituída na equação de parciais (1.9).
1. Parciais com Respeito à Posição e Velocidade no Quadro de Referência Fixo na Terra [0081] A subseção prévia apresenta uma computação de parciais de um quadro de referência inercial. Uma vez que as posições de satélite e os termos gravitacionais de ordem mais alta estão disponíveis no ECEF, é necessário computar parciais neste quadro. Uma vez que a equação (1.9) pode ser aplicada diretamente à computação de parciais com respeito à posição e à velocidade do objeto no quadro inerPetição 870180135007, de 27/09/2018, pág. 30/85
27/70 cial, na etapa 140, o resultado da computação na equação (1.9) é convertido para ECEF usando-se a equação de transformação (1.24) abaixo para a computação de parciais com respeito à posição e à velocidade (para pares de medições).
[0082] Para a conversão do resultado para ECEF, a fórmula a seguir é usada:
C\l>
Λ cx õAP
CX' CX _=1 Gxn,i õxn, _ (1.24) onde xEit é a iésima componente do vetor de posição e de velocidade no quadro de ECEF na época n, e xJ n, _ é a j-ésima componente do mesmo vetor no quadro inercial. O elemento J_ = ·· cX E ^Xn,i na matriz jacobiana J pode ser computado explicitamente, de modo que:
β:) (1.25) onde β = [[*]], * é o vetor da velocidade angular da Terra, e [[*]] é a matriz anti-simétrica 3x3 formada a partir de *.
[0083] É notado que as parciais são computadas primeiramente no quadro inercial. Isto é devido à forma especial da matriz Ψ no quadro inercial na equação (1.18), o que simplifica grandemente o resultado, e torna mais econômico computar parciais no quadro inercial primeiro, e, então, convertê-las para ECEF usando-se as equações (1.24). Contudo, isto não é uma exigência, mas apenas uma conveniência.
2. Parciais com Respeito ao Relógio e à Diferença entre a Leitura Atual e uma de Referência em um Relógio [0084] Algumas implementações de espaçonave usam um relógio com estabilidade relativamente baixa. Assim, uma estimativa do erro de relógio e, preferencialmente, de diferença entre a leitura atual e uma de referência em um relógio é necessária. Na etapa 150, as parPetição 870180135007, de 27/09/2018, pág. 31/85
28/70 ciais computadas acima com respeito à posição e à velocidade são aumentadas com parciais com respeito a relógio e diferença entre a leitura atual e uma de referência em um relógio, conforme descrito agora.
[0085] Isto começa pela estimativa da diferença entre a leitura atual e uma de referência em um relógio. Uma das medições de ADR na forma dada pela equação (1.2) é usada, e sua derivada é tomada com respeito ao erro de relógio e à diferença entre a leitura atual e uma de referência em um relógio. A diferença entre a leitura atual e uma de referência é assumida como sendo constante durante o intervalo de tempo entre as épocas n -1 e n. Isto produz as parciais a seguir:
(1.26) onde h é a diferença de tempo entre épocas e de longe é o maior termo. Assim, o processamento desta medição atualiza primariamente o valor da diferença entre a leitura atual e uma de referência em um relógio. Esta atualização pode não ser muito acurada, se um bom modelo de propagação não estiver disponível para o relógio (em um modelo, a diferença entre a leitura atual e uma de referência é assumida constante entre épocas). Inacurácias na estimativa da diferença entre a leitura atual e uma de referência em um relógio tornarão esta medição totalmente inútil para a estimativa da velocidade. Para estimativa da velocidade, precisa-se de uma medição diferenciada.
[0086] A computação de parciais para estimativa da velocidade do objeto com medições de diferença foi apresentada acima; contudo, nem todas as parciais foram computadas. Especificamente, as parciais das medições diferenciadas com respeito ao relógio do objeto e à diferença entre a leitura atual e uma de referência em um relógio permaPetição 870180135007, de 27/09/2018, pág. 32/85
29/70 necem por serem derivadas.
[0087] O numerador da equação (1.8) não contém o erro de relógio ou a diferença entre a leitura atual e uma de referência explicitamente, ainda há uma dependência implícita. O erro de relógio na época n é assumido como sendo τ, e a diferença entre a leitura atual e uma de referência em um relógio é τ. Isto significa que a n-ésima medição de época foi tomada em um tempo, o qual difere em -τ segundo inteiro, e a medição na (n -1)-ésima época foi tomada em um tempo, o qual difere em -τ + τ h do segundo inteiro. Uma vez que quaisquer mudanças em τ e τ afetam o tempo do instantâneo, elas também afetam a ADR medida. A computação das derivadas parciais correspondentes produz as fórmulas a seguir:
[(rn - Í n )2 11 - [(rn-1 - Ps, n-1 >
X.
n-1 ' s, n-1 (v - v )-(r - r ) (v 1 - v , )-(r . - r , ) \ n s, n J y n s, n / \ n-1 s, n-1 / \ n-1 s, n-1 / [(rn - rs,n )2 l2 [(rn-1 - P, n-1 ]
5τ n-1 ' s, n-1
Õ [(rn - rs,n )2 11 - [(rn-1 - P n-1 11 >
X.
n-1 ' s, n-1
(vn-1 Vs, n-1 )- (rn-1 P, n-1 ) , / >
---r-------------ü--h (1-27) [(rn-1 - rs, n-1 )2 f n-1 ys, n-1 ) χ n-1 ' s, n-1 * - r )2 12 n-1 f s, n-1) J [0088] Nos lados direitos da equação (1.27), os valores de posição e velocidade podem ser usados no quadro inercial ou de ECEF, uma vez que as fórmulas são escalares e válidas para ambos.
[0089] Usando-se o resultado da equação (1.27), uma matriz parcial é criada e é usada na etapa 160 para a formação das equações de medição de filtro de Kalman para a geração a partir das diferenças de ADR de uma estimativa de velocidade instantânea. As equações derivadas parciais descritas acima são formuladas de forma tal que o filtro de Kalman produza a velocidade instantânea a partir das diferenças de ADR. Assim, na etapa 170, esta velocidade instantânea é usada para a atualização da componente de velocidade do vetor de estado, durante a fase de medição do processamento de filtro de Kalman, como rePetição 870180135007, de 27/09/2018, pág. 33/85
30/70 presentado na etapa 170.
[0090] Assim, o filtro de Kalman processará um residual de medição de ADR para um par de fontes de sinal de cálculo de distância (por exemplo, satélites de GPS) na forma da equação (1.7), usando-se a matriz de parciais correspondentes. A matriz de parciais compreende duas divisões. A primeira divisão é para a posição e a velocidade (do objeto, por exemplo, uma espaçonave), e é definida por uma diferença dos dois termos, cada um correspondente a uma fonte de sinal de cálculo de distância (por exemplo, um satélite de GPS) e definido pelo lado direito da equação (1.9). A segunda divisão é para o relógio e a diferença entre a leitura atual e uma de referência em um relógio (de objeto, por exemplo, espaçonave), e é definida pelos lados direitos da equação (1.27). Se o vetor de estado da espaçonave incluir componentes adicionais (por exemplo, aceleração de relógio, coeficiente de reflexão de radiação solar, etc.) parciais correspondentes poderão ser derivadas e usadas de uma forma similar.
[0091] O precedente apresenta uma derivação de equações de medição para processamento de dados de ADR. As equações de medição para processamento de PR são bem conhecidas na técnica, e não são descritas aqui. O filtro de Kalman também processará as medições de PR para a atualização da componente de posição do vetor de estado na etapa 170.
[0092] As equações de medição de ADR estabelecidas acima são formuladas de uma forma muito mais acurada, de modo que o intervalo de tempo entre as medições de ADR usadas para a computação de uma diferença de ADR pode ser muito maior do que aquele usado para as técnicas atuais de diferenciação de ADR para aplicações de navegação por GPS. Uma analogia que ilumina as vantagens desta técnica é como se segue.
[0093] Considere a tarefa de medição de uma dimensão (compriPetição 870180135007, de 27/09/2018, pág. 34/85
31/70 mento, largura ou altura) de uma construção grande, tal como um prédio de escritórios. Uma forma de feitura destas medições é com uma régua relativamente curta. Devido ao fato de a régua ser curta, numerosas medições são requeridas com a régua para se completar uma medição de uma dimensão do prédio. Os erros de medição acumularse-ão para a pluralidade de medições feitas com a régua curta. Por outro lado, uma outra forma de feitura destas medições é usar uma trena muito longa. De fato, se a trena for longa o bastante, apenas uma medição pode ser necessária para uma dimensão do prédio, e uma medição muito mais acurada é obtida. As técnicas de navegação por GPS de ADR da técnica anterior são afins para se usar uma régua curta para a medição de uma distância significativamente grande. O erro de medição se acumula para cada computação de diferença de ADR porque o intervalo de tempo entre cada medição de ADR é curto (por exemplo, de 1 s), desse modo análogo ao problema da régua curta. Esses erros de medição de ADR diminuem a acurácia geral das computações de navegação, por causa do seu efeito cumulativo.
[0094] Embora os modelos de propagação usados nos filtros de Kalman tenham sido desenvolvidos e afiados para serem muito acurados, isto não foi o caso quanto às equações usadas para dados de ADR. As equações de medição de ADR estabelecidas acima atingem uma acurácia de medição que é comparável com a acurácia obtida por muitos modelos de propagação conhecidos. Assim, quando da atualização (por exemplo, propagação) de um vetor de estado com novos dados de medição, as equações de medição de ADR provêem dados de medição altamente acurados para as equações de propagação, o que melhora significativamente a acurácia das computações de navegação.
[0095] Uma outra analogia é começar um cronômetro e começar a viajar de carro e, após algum intervalo de tempo, parar o cronômetro e
Petição 870180135007, de 27/09/2018, pág. 35/85
32/70 medir a distância que o carro percorreu. A cada vez em que o cronômetro é começado e parado, um erro é introduzido. Se alguém usar um período de tempo curto, a estimativa da velocidade provavelmente não será muito acurada. Inversamente, se o carro se mover a uma velocidade constante, e se alguém realizar uma medição por um período de tempo longo, a estimativa de velocidade se tornará muito mais acurada. Embora a segunda opção seja preferível, nem sempre é praticável, uma vez que a hipótese de velocidade constante pode ser inválida. Neste caso, a medição de velocidade por um intervalo de tempo longo produzirá apenas a velocidade média por este intervalo. A velocidade média é de pouco uso prático, uma vez que a maioria das aplicações requer o conhecimento do estado atual.
[0096] Estas técnicas combinam a vantagem de medição da velocidade do objeto por intervalos de tempo relativamente longos com a capacidade de se mapearem de forma acurada tais medições para uma velocidade instantânea no momento do processamento. Isto é realizado pelo processamento de medições de ADR diferenciadas e pela aplicação de uma formulação matemática, a qual mapeia de forma acurada essas diferenças para o vetor de estado atual do usuário. [0097] Uma forma de explicação do método desta invenção é conforme se segue. Note que na analogia de velocidade do carregamento, a distância percorrida (a qual é similar à diferença de ADR no problema em questão) é diretamente proporcional à velocidade do carro. Portanto, a derivada da distância com respeito à velocidade do carro é igual ao intervalo de tempo da medição. O filtro de Kalman estima de forma canônica a velocidade pelo uso da medição em si (na analogia do carro, a distância percorrida), e sua derivada parcial com respeito ao desconhecido (na analogia do carro, esta derivada é o intervalo de tempo da medição). Nesta invenção, um método de alta acurácia é divisado para a computação das derivadas parciais requeridas pelo
Petição 870180135007, de 27/09/2018, pág. 36/85
33/70 filtro de Kalman, as quais, combinadas com o intervalo de tempo longo de cada uma das medições, produz uma estimativa de velocidade de alta acurácia para o objeto.
3. Propagação [0098] São estabelecidas abaixo as equações de propagação usadas pelo filtro de Kalman.
[0099] O vetor de estado em si pode ser propagado usando-se um modelo de força e uma rotina de integração (por exemplo, o esquema de Runge-Kutta). Novamente, este é bem conhecido na técnica e não é descrito aqui. A propagação da matriz de covariância é mais desafiadora. O ponto de começo é a derivada parcial do vetor de estado na época (n +1) com respeito ao vetor de estado na época n. Se estas derivadas parciais forem computadas, a matriz de covariância pode ser propagada usando-se fórmulas padronizadas. As derivadas requeridas são:
dxE 1.
= (1.28) BxE ϋ
[00100] O vetor de estado tem divisões para posição, velocidade e relógio/diferença entre a leitura atual e uma de referência. Assim, a matriz Ξ tem divisões para derivadas de posição com respeito à velocidade, posição com respeito ao relógio, velocidade com respeito à posição, etc. Estas divisões são consideradas separadamente abaixo.
a. Derivadas as Quais Não Envolvem Componentes de Relógio [00101] Uma comparação com a equação (1.10) mostra que esta matriz é similar à matriz Λ, e é assumido que a computação de Λ ajudará na computação de Ξ. Contudo, as componentes do vetor de estado na equação (1.28) são expressas no ECEF, e isto requer mais do que apenas a aplicação de Λ:
Petição 870180135007, de 27/09/2018, pág. 37/85
34/70 dx
E n+1,i õxE dXn, j =Σ k ,m õx
E n+1,i õx
I n+1,k õx^ õx^ uxn+1,k Uxn,m õx1 õxE n,m n, j (1.29)
[00102] õx1 Na equação (1.29), reconhece-se que —n = Jmj , conforme õxEJ
definido na equação (1.25), e õx
I n+1,k õx
I n,m +
km com as expressões para A+ km dadas pela equação (1.21). O sobrescrito + indica que, neste caso, a posição e a velocidade são propagadas para frente, em oposição à propagação para trás na definição para Akm na equação (1.10). (Estes dois casos correspondem ao uso de valores positivos e negativos para h na equação (1.21)).
[00103] õxE A quantidade remanescente a ser computada é —. Nesõxn+1,k
ta quantidade, o numerador e o denominador têm a mesma estampa de tempo, e nenhuma propagação está envolvida. Todas as derivadas são devido às conversões entre quadros de referência.
[00104] É assumido na equação (1.29) que o sistema inercial tem
seus eixos alinhados com os eixos de ECEF na época n. Ele mantém a mesma orientação para a época (n +1). Assim, a conversão de x[+1ki em xf+1>i. pode ser vista como compreendendo duas operações: conversão a partir do quadro inercial para o quadro de ECEF na época n e, então, conversão a partir do ECEF na época n para aquele na época (n +1). Esta conversão pode ser expressa conforme se segue:
õx
E n+1,i õx
I n+1,k = (ω · J 4.· (1.30) onde Ω = ’ e M é a matriz de rotação 3x3 a partir do ECEF na época n aquele na época (n +1). J ' pode ser formada a partir de J pela mudança de ) para -). A substituição de todas as
Petição 870180135007, de 27/09/2018, pág. 38/85
35/70 matrizes em (1.29) produz a expressão a seguir para as divisões Ξ, as quais não envolvem componentes de relógio:
Ξ r = M .(λ*„ +Λ*„ ·))
Ξ = M -Λ r t x Z x (1.31)
Ξ vr = -M )\Krr - V ))+M · (λ, - Λ-, ))
Ξ w = M-(Λ-/,-)·Λ+ν )
b. Derivadas de Posição/Velocidade com Respeito a Relógio/Diferença Entre a Leitura Atual e uma de Referência em um Relógio [00105] Esta subseção descreve as computações das derivadas
Grl. and õt õt õi õi
Quando da computação destas derivadas, é notado que a posição/velocidade na época n +1 depende da posição/velocidade na época n, e a última, por sua vez, depende do erro de relógio e da diferença entre a leitura atual e uma de referência em um relógio. Assim, há uma dependência implícita da posição velocidade na época n +1 do relógio/da diferença entre a leitura atual e uma de referência na época n, apesar da posição/velocidade na época n. Contudo, precisa-se computar derivadas parciais, e esta dependência implícita deve ser excluída. O resultado é como se segue.
[00106] Para equações conservativas (conforme é o caso com di-
nâmicas de satélite no campo gravitacional da Terra no quadro de ECEF), não há uma dependência os erro de relógio (isto é, as derivadas são nulas), mas há uma dependência ligeira da diferença entre a leitura atual e uma de referência em um relógio. É dada pela fórmula a
seguir: õfE η+λ = C · h õt (1.32) ÕVn-1 = aE · h õt
onde ãE+1 é a aceleração do objeto na época n +1 no quaPetição 870180135007, de 27/09/2018, pág. 39/85
36/70 dro de ECEF (a aceleração é computada usando-se modelos de força).
c. Derivadas do Relógio/da Diferença entre a Leitura Atual e uma de Referência com Respeito a Outras Compo- nentes
[00107] Estas são triviais e são dadas pelo seguinte: dt„+i δτη+ι 0τη+ι 0τη+ι 0 drn E õvE n õrE dvE Ô F = 1 Ôtn dTn+1 = 1 (1.33) = h Ôtn p = 0 4. Termos Gravitacionais de Ordem Mais Alta
[00108] Esta seção representa uma computação de matrizes φ pa-
ra termos gravitacionais de ordem mais alta. O potencial gravitacional de termos de ordem mais alta tipicamente é dado em coordenadas polares (veja M. Kaplan, Modern Spacecraft Dynamics & Control, John Wiley & Sons, 1976):
ii / k C R õk í \ u(τ,θ,φ}=ΣΣΙI P(sin.)-(Qcosj- + Sksin J-} (1-34) r k=2 j=0 \ r ) [00109] A fórmula (1.34) difere da fórmula (7.16) em Kaplan, conforme se segue: ele usa Θ para latitude (note que este não é o ângulo polar), - para longitude, e Ck 0 no lugar de - Jk.
[00110] As primeiras e segundas derivadas parciais de u com respeito às coordenadas polares são:
Petição 870180135007, de 27/09/2018, pág. 40/85
37/70 τ - t (. )*(k + n
IHt Σ ( r-E Pk(sin.) ( Σ Σ (—Ί p(sin .)·j (- cksin j-+Skcos j—— õ- rk=2j=0 i r ) õ 2 = (ΣΣ(—} (k + 1>(k + 2)· Pj (srn.)-(c/ õr r k=2 i r )
W-(l Σ (RE)‘Pk (sM-j >&
g. = (ΣΣfR' \P 1 (sin.) cos2 θ-Pk1 (sin.) sin.] (ckcos j-+Sksin j—) Õ. r k=2 j =0 i r ) (1.35) õ 2U õrõ.
õU cos j-+Sksin j-) cos. (c/ cos j- + Slj sin j—) õ 2U
COS j- + Sk sin j-) cos j- + Sksin j-) = -(Σ Σί ~ Í(k + 1)-P?'(sin.) r k=2 j=0 i r ) θ — = (ΣΣf—Ί Pkj(sinθ) cosθ j-(-cksin j-+Skcosj—— õ.õ- rk=2 j=0 i r )
Οο— = ~(Σ ΣfR' (k +1 )P (sin.j ( cksin j- + Sk cos j—— õrõ— r k =2 j =0 i r ) cos. (ck cos j— + S} k sin j—) [00111] As derivadas parciais com respeito às coordenadas cartesianas são computadas como se segue:
Õ U = Σ Õ Pm ÕU + Σ ÕPm ÕPn . Õ U (1 3θ) dxt dXj m ÕV ÕPm mn Õx j Ox 1 ÕPm ÕPn ' onde o vetor P é definido como P = {r,.,—}. Note que (A3) contém as derivadas de primeira e de segunda ordem de coordenadas polares com respeito às coordenadas cartesianas. Estas derivadas são dadas por:
õr x õr y õr z õx r’ õy r ’ õz r ’
õ. zx cos2
õ. = zy cos2 θ õ. = cos õx r2 3 õz r2 õ— = cos— õ— = y õy r2 õz rj ’
Petição 870180135007, de 27/09/2018, pág. 41/85
38/70 d2 r = 1 x2 dx2 r r3 ’ d2 r xy dxdy r3 ’ d 2Θ z cos2 dx2 d 2Θ dy y d 2Θ dz2 dxdz d2 r = 1 y2 dy2 ~ r r3’ d2 r xz dxdz r3( 3 zx cos r 3 + 3 r n 2 2 z cos Θ zy cos
--— + 3 r2 cosΘ· sinΘ d2 r = 1 _ z2 . dz2 r r3 ’ d2 r zy dzdy r3 ’ _ zx cos Θ · sin + 2-------3---r2 n zycos.· sin + 2-------5 3 r2 r2 dz ’
6Θ.
dx ’ (1.37)
d..
dy ’ cos^· sin δΘ x · cos2 dx r[ cos.· sin.
r2 zy cos Θ · sin d 2Θ = . dydz df== 2 dxdy d 2φ = 2yx d 2φ =sin dy2r± d2φ1 —— = y + 2^ dxdyr d 2φd dz2 r2
0Θ y · cos2 Θ dy rl d Θ 3 zxy · cos2 dx dzdx dzdy d- y·cos dy rl ~y2.
r24
77=o
Em resumo, as parciais usadas no filtro de Kalman para o [00112] processamento de dados de diferença de ADR são computadas usando-se a metodologia a seguir:
1. Converter o vetor de estado a partir do ECEF em um quadro inercial, o qual é momentaneamente alinhado com ECEF (mas não girando).
2. Propagar o vetor de estado atual de volta um incremento, usando um esquema de Runge-Kutta de 4a ordem.
3. Em cada subincremento, computar as matrizes de elemento φ usando a equação (1.23) ou as fórmulas para termos gravitaPetição 870180135007, de 27/09/2018, pág. 42/85
39/70 cionais de ordem mais alta.
4. Substituir estas matrizes no conjunto (1.22), e o resultado no conjunto (1.21). Isto computa Λ.
5. Substituir Λ na equação (1.9) e, então, aplicar uma equação de transformação (1.24) para computação das parciais com respeito à posição e à velocidade. Usá-las para pares de medições.
6. Usar a equação (1.26) para uma medição única e a equação (1.27) para pares para aumentar as parciais para as componentes de relógio e de diferença entre a leitura atual e uma de referência em um relógio.
B. Propagação de Vetor de Estado Usando-se Dados de IMU para uma Navegação de GPS - IMU Integrados [00113] O que vem a seguir descreve algoritmos de navegação para processamento de dados de ADR em sistemas de GPS - IMU integrados. Uma IMU não provê quaisquer medições absolutas de posição ou velocidade de objeto e, portanto, pode ser usada primariamente como um auxílio na navegação. Uma estimativa real de posição e velocidade ainda é o resultado do processamento de medições de GPS. Contudo, os dados de IMU são importantes como um meio de propagação acurada, o que é crucial para a obtenção de uma performance de navegação alta, particularmente quando um modelo de propagação acurado ainda não está disponível para o objeto de interesse.
[00114] Para apreciação da importância de medições de IMU, é instrutivo comparar seu papel com aquele de um modelo de prensa de punho de quadro duplo em aplicações espaciais. Para a obtenção de uma boa acurácia em aplicações de navegação de espaçonave, um modelo de propagação sofisticado é vital. Seu papel é assegurar uma consistência entre medições coletadas por um período de tempo estendido. Um bom modelo de propagação permite a regulagem do ruído de processo para um nível correspondentemente baixo, o qual faz com
Petição 870180135007, de 27/09/2018, pág. 43/85
40/70 que cada medição passada de um histórico recente contribua implicitamente para a estimativa atual quase a par com qualquer medição atual. O resultado é como se houvesse mais medições para a obtenção de uma solução em qualquer dado tempo. (Obviamente, as medições passadas não são aplicadas diretamente no tempo atual; ao invés disso, a estimativa prévia é propagada para o tempo atual e é ponderada pesadamente na solução.) Por exemplo, se um modelo de propagação for bom para 1000 segundos, isso facilita usar aproximadamente 1000 vezes mais medições do que em uma época, o que deve reduzir o erro de navegação por um fator de V1000 « 30 vezes.
[00115] Há uma outra consideração importante: um erro no modelo de propagação tipicamente está altamente correlacionado no tempo (isto é, sistemático), em oposição a um erro randômico nas medições. Os erros sistemáticos não são removidos por um cálculo de média de filtro, e um pequeno erro sintomático pode causar mais dano do que um erro randômico muito maior. Isto adicionalmente enfatiza uma necessidade de um bom modelo de propagação. As fontes de erro no modelo de propagação devem ser examinadas e mitigadas. Ao se fazer isso, deve ser notado que uma implementação em computador de um modelo matemático pode introduzir erros adicionais, freqüentemente de um tipo sistemático. Por exemplo, uma integração de equações diferenciais de movimento para propagação de espaçonave deve usar um esquema de ordem mais alta, tal como Runge-Kutta. Um esquema de ordem baixa (por exemplo, de Euler) é análogo a ter uma aceleração não modelada devido a erros na representação de equações diferenciais em uma forma discreta.
[00116] Uma aplicação atraente para um bom modelo de propagação derivado seria estender métodos para estimativa de velocidade através de medições de ADR para processamento de navegação de GPS/IMU integrados. Se melhores modelos de propagação forem dePetição 870180135007, de 27/09/2018, pág. 44/85
41/70 senvolvidos para sistemas de navegação de GPS/IMU integrados, uma extensão de processamento de ADR para esses sistemas de navegação de GPS/IMU integrados tornar-se-á possível.
[00117] Em um sistema de GPS/IMU, a parte de IMU do sistema serve até certo grau como um meio de propagação de uma solução de navegação. O que vem a seguir apresenta um modelo de alta acurácia para a propagação de uma solução usando-se dados de IMU.
1. Quantidades Observáveis e o Vetor de Estado [00118] Uma IMU mede acelerações e taxas de rotação. Estas quantidades são derivadas do tempo de parâmetros de velocidade e atitude, e podem ser incluídas no vetor de estado. Contudo, pode haver um custo de processamento significativo associado a isto: uma IMU tipicamente extrai dados a de 100 Hz a 1000 Hz, e, em geral, é excessivo rotar o processamento de filtro a uma taxa tão alta como essa. Ao invés disso, os dados de aceleração e taxa de rotação extraídos por uma IMU são usados estritamente como uma ferramenta de propagação.
[00119] Mesmo quando os dados de IMU são usados como uma ferramenta de propagação, há um acoplamento entre a atitude, as taxas de rotação, a aceleração, a velocidade e a posição. Este acoplamento é devido a elementos fora de diagonal da matriz de covariância e serão explorados em detalhes aqui adiante. É análogo à estimativa de velocidade a partir de medições de posição: embora a posição/distância seja a quantidade medida, a velocidade pode ser estimada, porque ela está acoplada à posição através de elementos fora de diagonal da matriz de covariância.
[00120] Esta abordagem permite evitar o uso de aceleração e taxas de rotação no vetor de estado. Os parâmetros candidatos para inclusão no vetor de estado são conforme se segue:
1. Três componentes de posição.
Petição 870180135007, de 27/09/2018, pág. 45/85
42/70
2. Três componentes de velocidade.
3. Relógio de objeto e taxa de relógio.
4. Atitude.
5. Erros sistemáticos em relação a um valor de referência, desalinhamento, fatores de escala e outros erros de IMU sistemáticos. [00121] A representação específica de informação de atitude no vetor de estado será descrita aqui adiante. Por agora, é notado que a atitude é plenamente descrita pelos três parâmetros.
2. Impacto de Erros de Medição de IMU na Propagação [00122] Esta subseção descreve os efeitos de fontes de erro de IMU randômicos.
[00123] Uma propagação é considerada por algum período de tempo T (por exemplo, um segundo) usando-se medições de aceleração ai, as quais estão disponíveis no tamanho de incremento h (por T exemplo, h = 0,01 segundo). Há um total de N = — medições de aceh leração por época.
[00124] Se uma medição de aceleração tiver um erro 3, o resultado será o erro na velocidade, o qual pode ser estimado como qh. Este erro de velocidade contribui para a propagação da posição do objeto durante o tempo remanescente até o final da época, isto é, durante um intervalo de tempo (t - ti), onde ti = ih. O erro resultante na posição a partir desta medição em particular então é estimado conforme se segue:
δ =3ih -(t - ti) (2.1) [00125] A variância deste erro pode ser estimada pela elevação ao quadrado da equação (2.1). As variâncias a partir de todas as subépocas se adicionam para a formação da variância de erro total no final da época:
Petição 870180135007, de 27/09/2018, pág. 46/85
43/70 δχ2 =Σ3 - h2 -(T-itif (2.2) ii [00126] A computação da soma leva a:
δχ2 « j^2 - h - T3(2.3) [00127] De modo similar, uma computação de variância de erro para velocidade leva ao seguinte:
δ2 = 32 - h - T(2.4) [00128] Erros de acelerômetro não contribuem para erros na atitude do objeto.
[00129] Erros de medição de giroscópio contribuem para um erro de posição de objeto, um erro de veículo e para um erro de atitude. Se as taxas de rotação p forem medidas em cada época com um erro õpt, então, os erros de atitude, posição e velocidade serão dados por:
ÕC2 = Õp2 - h - T (2.5)
Õx2 = — Õp2 - h - a2 - T5 (2.6) õv = ^Õp2 - h - T3 (2.7) onde a é a aceleração do objeto durante esta época (ela varia com o tempo, mas, para as finalidades de computação da variância de erro, a aceleração em qualquer ponto durante a época, por exemplo, no começo, pode ser usada), e quando o erro na posição e na velocidade é não nulo nas direções perpendiculares à aceleração. [00130] Conforme mencionado acima, isto é um erro randômico, e pode ser contabilizado pela aplicação de um ruído de processo para posição, velocidade e atitude propagadas. Erros sistemáticos em relação a um valor de referência, desalinhamentos, fatores de escala e outros erros sistemáticos são tratados separadamente.
3. Atualização de Medição [00131] As equações de atualização de medição processam mediPetição 870180135007, de 27/09/2018, pág. 47/85
44/70 ções de GPS de PR e, opcionalmente, medições de ADR. As equações de atualização de medição de PR são bem conhecidas na técnica e, portanto, não são descritas aqui. As equações de atualização de medição de ADR são as mesmas que aquelas apresentadas acima para uma navegação apenas por GPS. As únicas componentes não nulas da matriz de medição são aquelas correspondentes à posição e ao erro sistemático em relação a um valor de referência de relógio. [00132] Para medições de ADR, métodos de diferenciação ou adaptativos são usados para a estimativa da velocidade. No método descrito aqui, as medições de ADR são diferenciadas no tempo para cancelamento do erro sistemático em relação a um valor de referência. Para se atingir o potencial pleno, a diferenciação de ADR é realizada por intervalos de tempo os quais cobrem múltiplas épocas, conforme descrito acima para a navegação apenas por GPS.
[00133] Para se processar uma medição de ADR diferenciada no tempo, como por um filtro de Kalman, a matriz de observação (H) envolvendo derivadas parciais da medição com respeito às componentes de vetor de estado é computada:
Õ\ADR
ÕX
-M Xs, n-M õx onde X é o vetor de estado de objeto compreendendo posição, velocidade, atitude e estados de erro de IMU.
[00134] Conforme descrito acima em relação às equações (8) e (9), uma vez que a posição de objeto Xn é explicitamente usada na equação (2.8), ela cria uma impressão de que nós temos apenas que computar parciais com respeito a ela; contudo, o estado de objeto Xn é o resultado da propagação de vetor de estado a partir da época n - M. Assim, Xn é uma função deXn_M; daí, Xn_M é uma função (inversa) de
Petição 870180135007, de 27/09/2018, pág. 48/85
45/70
Xn, e Xn-M (que é uma parte de Xn_M) é uma função de Xn. Isto resulta em uma computação mais complexa das derivadas:
XX(X n Xs,n ) [Λ ( Xn-M Xs, n-M )] (2.9) [( Xn-M Xs, n-M onde
Figure BRPI0515452B1_D0001
(2.10) e χ contém tudo zero, exceto pelos três primeiros elementos da diagonal principal, os quais são iguais a 1.
4. Propagação
a. Propagação para Frente [00135] Os algoritmos de propagação para frente são necessários para o vetor de estado e a matriz de covariância. Com referência à FIGURA 6, um método para uso de dados de IMU para a propagação de um vetor de estado computado a partir de medições de PR e diferenças de ADR será descrito. Na etapa 200, a matriz de observação H é computada a partir do vetor de estado, conforme indicado na equação 2.10.
i. Propagação de Atitude [00136] Em seguida, na etapa 210, a componente de atitude é propagada pela propagação de uma matriz de rotação C(t) que transforma um vetor a partir dos eixos de corpo de IMU no sistema de coordenadas de ECI. A meta é derivar um algoritmo de propagação para esta matriz, o qual usa dados de sensor de giroscópio.
[00137] Em cada subépoca, a matriz C(t) sofre uma pequena mudança. Esta mudança é devido à rotação da IMU durante esta subépoca, cuja taxa é medida pelo giroscópio no sistema de referência do corpo. Assim, se uma matriz C (t + dt) for necessária para a transformaPetição 870180135007, de 27/09/2018, pág. 49/85
46/70 ção de um vetor a partir do sistema de corpo para ECI no tempo (t + dt), isto poderá ser feito pela aplicação de duas etapas consecuti vas:
1. Girar a partir do sistema de corpo no tempo (t + dt) para aquele no tempo t
2. Girar a partir do sistema de corpo no tempo t para ECI pela aplicação da matriz C (t) [00138] A primeira rotação é pequena e, portanto, todos os três ângulos de rotação são comutativos. É o inverso da rotação a partir do sistema de corpo no tempo t para aquele no tempo (t + dt), o que é da forma:
ÕC (t ) =
- dtp3 dtp3
- dtp2 ' dtp1 = I + p (t) · dt (2.11) γ dtp2 - dtp1
Figure BRPI0515452B1_D0002
onde a matriz p = [p ]] é a matriz anti-simétrica formada a partir do vetor de três componentes de taxas de rotação, conforme medido pelo giroscópio. A inversa desta matriz é simplesmente a mesma matriz, mas com um sinal menos após I. Assim, a propagação da matriz C (t) é descrita por:
C (t + dt )= C (t )·(ΐ -p · dt) (2.12) [00139] Esta fórmula seria perfeitamente boa, se uma propagação fosse necessária por um intervalo de tempo infinitesimalmente curto dt . Contudo, é modificada para a propagação dela por uma época (por exemplo, 1 segundo). Uma generalização simples leva a:
C(t + T) = C (t )· Π (I - p(l, )· dt) (2.13) onde a ordem de multiplicação de matriz é tal que as últimas subépocas correspondam aos multiplicadores mais à direita no produto. (Esta expressão está proximamente relacionada à integral multiplicativa de Volterra.) Este esquema de propagação tem a acuráPetição 870180135007, de 27/09/2018, pág. 50/85
47/70 cia de primeira ordem em cada incremento (isto é, o erro da ordem de o(dt2)) e o erro da ordem de o(dt) após uma propagação pela época inteira.
[00140] Em aplicações em que uma alta acurácia não é requerida, esta generalização simples pode ser suficiente; ainda, há espaço para um melhoramento adicional. Em primeiro lugar, a propagação de segunda ordem é buscada em cada incremento.
Esquema de propagação de segunda ordem para matriz de rotação [00141] A primeira etapa é centralizar o esquema de propagação. (De fato, se uma integral aditiva for propagada ao invés de uma integral multiplicativa, isto seria o que seria necessário.) Por brevidade, faça o subscrito 0 se referir ao tempo t, e o subscrito 1 se referir ao tempo (t + dt). Então, um esquema de propagação é buscado na forma a seguir:
C, = Co ·[/- 2(1 +1 )· d/j + o(di2) (2.14) onde C, é o valor exato de matriz de rotação no tempo (t + dt). Todas as diferenças entre a propagação exata e nossa fórmula são contabilizadas pelo termo o(dt2).
[00142] Em primeiro lugar, considere uma propagação por dois intervalos de tempo consecutivos. Há duas formas de se fazer isso: um incremento de 2dt de comprimento, ou dois incrementos de dt cada. Uma vez que a matriz de rotação exata é computada, os resultados devem combinar:
C2 = C0 [ I-2 (l +12 )· 2dt | + 4o(dt t ) k 2 7 (2.15) C2 = C0 - 2 (Ρθ +P1 )· dt^[l-2 (Â +P2 )· dt^j + 2o(dt 2 ) [00143] Note que, no primeiro caso, o erro o(dt2) aumentou por um
Petição 870180135007, de 27/09/2018, pág. 51/85
48/70 fator de 4, se comparado com a fórmula de dt único. Isto é porque o(dt1) tem uma dependência quadrática de dt, e dt agora é duas vezes maior. No segundo caso, o erro aumentou por um fator de 2, se comparado com a fórmula de dt único. Isto é porque o mesmo erro é obtido duas vezes pela aplicação da propagação duas vezes.
[00144] Após usar uma expansão de Taylor para pt, a segunda equação em (2.15) pode ser aproximada conforme se segue:
C2 = Co I - p0dt -1 p0 dt2I - p2dt +1 p2 dt22o(dt2) (2.16) [00145] Agora, o produto na equação (2.16) pode ser computado e o resultado igualado à primeira equação em (2.15). O resultado é conforme se segue (para a segunda ordem de acurácia com respeito a dt ):
0p) = d-2-C„-pt,-p, (2.17) [00146] Esta estimativa pode ser substituída na equação (2.14) para o resultado final, acurado para a segunda ordem com respeito a dt:
C1 = C0 ]j- 2 !Po +  )· dt + 2 Po  dt 2 )+ o(dt 3 ) (2-18)
Esquema de propagação de terceira ordem para matriz de rotação [00147] A mesma técnica pode ser repetida para a obtenção de uma fórmula válida para a terceira ordem com respeito a dt. As equações para a matriz de rotação após um incremento de 2dt de comprimento e após dois incrementos de dt são, cada uma, conforme se segue:
Petição 870180135007, de 27/09/2018, pág. 52/85
49/70
C2 = C0 · -1 (p0 + p2 )· 2dt + 2p0 · p2dt28o(dt3) C2 = C0 - 2 (p0 +p1 )· dt + 2 p0 -p x dt )· - 2 (p1 +p2 )· dt + 2 P1 ·ρ2 dt 2 ) + 2o(dt 3 ) [00148] Após uma álgebra similar, obtém-se:
3 (2.19) (dt3 ) = 4 (pQ . p. )-2 po 1 +p2 )p +(/?0 ·/?, -po -P2 ) (2.ii. * * * * * * * * 20) [00149] Assim, o esquema de terceira ordem para a propagação da matriz de rotação é conforme se segue:
(2.21)
C(tn )= C(tn-1 )'âC(tn-1 ) onde:
SC (tn-1 )= I-2 (pn-1 +pn )· dt + 2 Pn-1 -Pn ’ dt2 + ~ 4 (pn-1 + Pn )- 2 Pn-1 (pn-1 + Pn )pn + (pn-1 - Pn - Pn-1 - Pn ) n-1 n-1 (2.22) ii. Propagação de Posição e Velocidade [00150] Na etapa 220, os algoritmos de propagação descritos acima são estendidos para posição e velocidade. A propagação envolve uma integração ao longo do tempo, usando-se esquemas de integração de ordem mais alta, de uma expressão que contém a matriz de atitude.
[00151] A integração é representada como uma soma finita, e um esquema é derivado, que é válido para a segunda ordem de acurácia em cada incremento. Isto requer uma computação da matriz de atitude para a terceira ordem para cada incremento. De fato, o número de incrementos precedendo ao atual é inversamente proporcional ao tamanho de incremento. Assim, o erro na matriz de atitude pode se acumular para se tornar da segunda ordem com respeito ao tamanho de incremento, o que não seria consistente com o restante do esquema. Ainda, esta consideração não se aplica para a aproximação da matriz
Petição 870180135007, de 27/09/2018, pág. 53/85
50/70 de atitude no tamanho de incremento para as finalidades de derivação de uma fórmula para posição, onde um esquema de segunda ordem poderia ser aplicado.
[00152] O que vem a seguir é uma derivação para propagação de posição; uma fórmula para velocidade pode ser derivada de modo similar, e com um pouco menos de esforço.
[00153] A posição no fim de uma época (a qual contém N subépocas) é dada por:
!n t x = X0 + V0 T + JJC(τ)·a[τ)άτάί (2.23) t0 t0 onde X0, 1' são a posição e a velocidade no começo da época, e a (...) é a aceleração nas coordenadas do corpo, conforme medido pelo acelerômetro. Uma integração por parte simplifica a fórmula acima para uma integral única:
τΝ
X = X0 + V0 T + JC(t) a(t)-(T-t)dt (2.24) t0 [00154] A meta é computar esta integral. O método de trapezóide pode ser usado para uma aproximação da integral por uma soma finita. Cada incremento fará a contribuição a seguir:
tn
J C (t ) a (t )· (T -1 )dl = tn-1 d -[C(tn-1 )· a (tn-1 }(T - t (2.25) )+ C(t„ ) a (t„ )(t - t n )]+ o(dt3 ) [00155] Note que isto é uma fórmula exata; todos os erros no método de trapezóide são incorporados no termo o(dt2). Uma fórmula similar pode ser derivada para propagação de velocidade; ela difere apenas de (2.25) pela ausência dos multiplicadores (t - tn-1) e (t - tn).
[00156] C(tn-1) é assumido como já propagado para o tempo tn-1, usando-se o esquema de terceira ordem (o qual pode produzir um erro de segunda ordem por n -1 subépocas). Para as finalidades desta dePetição 870180135007, de 27/09/2018, pág. 54/85
51/70 rivação apenas, a atualização da matriz de atitude no tempo tn pode ser estimada a partir daquela no tempo tn-1, usando-se o esquema de segunda ordem:
Cn = Cn-1 '(j - 2 (fin-1 +Pn )' dt + j Pn-1 'Pn ’ dt t + (Á^i3 ) (2.26) [00157] O restante da derivação é conforme se segue: a equação (2.26) é substituída na equação (2.25), e a fórmula resultante é aplicada a dois casos: uma propagação única por um incremento, ou duas propagações de meio incremento subseqüentes. Os resultados devem ser idênticos (uma vez que a fórmula (2.25) é exata). Após alguma álgebra, isto produz uma equação para o(dt1), a qual leva ao seguinte:
j C (t)· ã (t )· (t - t)dt = r i (2.27)
- · C- (tn-1 )< - tn-1) · a (tn-1 )+(T - tn )' SC )· ã (t „ )] + · C(tn-1 )· onde os valores de todos os parâmetros no termo de ordem mais alta podem ser computados em qualquer momento no tempo no intervalo (tn-1, tn) (o esquema ainda reteria a ordem) e onde õC(tn_1) é derivado anteriormente. O incremento total na posição é a soma por n da parte direita da fórmula (2.27).
[00158] Uma equação para propagação de velocidade é conforme (T -1)· ^—p · a + P · ã - — p1 · ã -— ã^+ (ã - p · ã) + o(dt4) se segue:
tn j C (t )· ã (t )dt = tn-1 ~ · CC(tn-1 ) · [ã (tn-1 )+ Õ C(tn-1 )ã (tn ) ] dt3 d 1 A r . R 1 .
— · C(tn-1 )· ~P · ã + P · ã P
2 2 (2.28) + o(dt4) +
r 1 R > · ã--ã
Novas notações são introduzidas, para simplificação:
Petição 870180135007, de 27/09/2018, pág. 55/85
52/70
V = Σ:' n
Χ = Σβ(tn-1 )'5 n
(2.29) onde 6 = d ·[ (t«-1 )+ δ C(tn-1 )· a (tn )]+
6 5 = (T - tn )'6n +7
A r . R 1 .2 — p · a + p · a--p r 1 R · a--a r dt2 \ dt3 A „ λ 7n =— · a (tn-1 ) 6 ·(«-p · a) iii. Desalinhamentos [00159] Até agora, o tratamento de erros de IMU está limitado à consideração de desalinhamentos no sensor de giroscópio. Conforme notado nos colchetes, este é um erro difícil de se modelar. Na etapa 230, um vetor S, conforme descrito abaixo, é definido, que caracteriza o desalinhamento e os erros sistemáticos em relação a um valor de referência no sensor de IMU para o vetor de estado de filtro de Kalman, e derivadas correspondentes são computadas.
[00160] Se o sensor de giroscópio de IMU estiver desalinhado, as taxas de rotação medidas serão diferentes das taxas de rotação verdadeiras no quadro de referência de corpo. O modelo o qual descreve este efeito é conforme se segue:
p = M ·ρ (2.31) onde pt é o vetor verdadeiro de taxas de rotação no quadro de corpo, p é o vetor medido (de saída), e Ml é alguma matriz, a qual gira p para pt. Como antes, a matriz p é definida por p = [[r]]. A matriz Mí é rotativa e está próxima da matriz identidade e, portanto, pode ser representada na forma:
Mí = I + S +1S2 (2.32) onde S = [[S]], e o vetor S caracteriza três parâmetros de desalinhamento. Estes parâmetros são geralmente desconhecidos no
Petição 870180135007, de 27/09/2018, pág. 56/85
53/70 começo de uma aplicação de navegação e, assim, devem ser estimados.
[00161] Para a incorporação deste modelo no algoritmo de navegação, o que vem a seguir é realizado:
1. A inclusão do vetor S no vetor de estado de filtro de Kalman e a computação das derivadas parciais correspondentes.
2. A propagação do vetor de estado com as taxas de giroscópio corrigidas pela aplicação das fórmulas (2.31) e (2.32), antes do uso de p nas equações de propagação.
3. A estimativa de parâmetros de desalinhamento. Eles são assumidos como sendo estacionários.
[00162] Para se alcançar o item (3), a matriz de propagação é computada, a qual contém derivadas parciais da posição e da velocidade com respeito a parâmetros de desalinhamento. Também, uma computação similar é aplicada para a computação de derivadas parciais de uma medição com respeito a estes parâmetros.
[00163] A computação de derivadas parciais é conforme se segue. São definidas
Figure BRPI0515452B1_D0003
e as fórmulas gerais para G) são derivadas, onde pt, pt, etc., podem ser ligados. Então, estes blocos de construção são usados para se encontrarem as derivadas parciais desejadas.
[00164] Uma derivação para G(pt) é apresentada aqui adiante, o que resulta no seguinte:
G.,(p,) = |/ = p,-S,p, + .[.58p], + S,p,-[Sxpp -S„p,) (2·34)
Petição 870180135007, de 27/09/2018, pág. 57/85
54/70 [00165] Note que a parte direita da fórmula (2.34) contém as taxas de rotação p, as quais não são corrigidas quanto aos desalinhamentos (daí, não têm o subscrito t). Isto é porque a derivação da fórmula (2.34) considera correções de desalinhamento explicitamente.
[00166] Para derivadas da matriz de atitude, a fórmula recursiva a seguir é provida (onde todo p e suas derivadas abaixo representam valores, os quais já estão corrigidos quanto a desalinhamentos):
Figure BRPI0515452B1_D0004
ou, de forma equivalente:
Figure BRPI0515452B1_D0005
onde
Figure BRPI0515452B1_D0006
(2.37)
Para derivadas da posição e da velocidade:
ÕV
Figure BRPI0515452B1_D0007
: + 'C/' 7, - dt6- C Ç-Ç G (p ) a.
ÕS 6 (2.38) onde
Figure BRPI0515452B1_D0008
iv. Matriz de Propagação para Frente [00167] A próxima etapa, 240, é para a derivação de uma matriz de propagação. Derivadas parciais do vetor de estado são computadas no final de uma época com respeito ao vetor de estado no começo da
Petição 870180135007, de 27/09/2018, pág. 58/85
55/70 época.
[00168] A escolha quanto à representação da atitude no vetor de estado pode ser mais adiada. A matriz de atitude em si não é adequada, porque tem 9 elementos, ainda sendo definida por apenas 3 parâmetros independentes. Os ângulos de Euler têm uma desvantagem de exibirem degeneração em certos valores. Finalmente, uma representação de quaternião tem duas desvantagens: ela contém 4 elementos (apenas 3 sendo independentes), e requer a computação da matriz de atitude em cada incremento. Uma aplicação subseqüente de duas rotações de quaternião não é trivial. Isto nos leva à escolha a seguir quanto à representação de atitude: em cada época, um pequeno ângulo de rotação é definido com três ângulos de Euler pequenos com respeito à matriz de atitude atual. Assim, a matriz de atitude é representada na forma:
C,(( ) = (/ + (t) (2.40) onde C(t) é a estimativa atual da atitude, e (i + [[<?]]) é uma pequena rotação na forma
Figure BRPI0515452B1_D0009
(2.41)
Figure BRPI0515452B1_D0010
[00169] Os ângulos φ1, φ2, φ3 são pequenos, e as rotações corres pondentes são comutativas. Estes ângulos formam a divisão de atitude do vetor de estado.
[00170] Um exame próximo desta abordagem releva dois problemas em potencial:
1. Após algum tempo, uma acumulação de erro pode resultar em φ15 φ2, φ3 não serem mais pequenos. A solução para este problema é o reequilíbrio periódico (por exemplo, a cada época) da matriz de atitude e dos valores de φ15 φ2, φ3. Em cada época, uma rotação de
Petição 870180135007, de 27/09/2018, pág. 59/85
56/70
Δά(ζ) para C(t) é aplicada e, assim, atualiza a matriz de atitude atual e reinicializa φ,, φ2, φ3 para zero. A transformação correspondente também deve ser aplicada às covariâncias de φ,, φ2, φ3. Isto é similar a uma linearização de medições de pseudodistância em cada incremento com respeito a uma nova posição estimada atualizada. A única diferença é que esta linearização para matrizes é multiplicativa, ao invés de aditiva.
2. Uma aplicação sucessiva de pequenas rotações pode produzir uma matriz a qual não é mais puramente rotativa (por exemplo, a condição de ortonormalidade CT(t)·C(t) = I se torna não acurada). A solução para este problema é periodicamente refundir a matriz de atitude para uma forma ortonormal.
[00171] As parciais são derivadas para propagação por uma época no sistema de referência inercial:
V.,= V+(I -II'’V0 (242) XI,1 = xi.o + Vi.oT + (l + |Ιφ]])· C0 ·X0
Computação Auxiliar:
C~ = Co · c (2.43) [00172] Se houver um desalinhamento, a matriz de atitude final diferirá do valor nominal:
~ ~ ~ ( 1C j
C, + dC1 = Co · Co + . · dS (2.44) l dS ) [00173] Por outro lado, esta mudança na atitude é modelada pelos parâmetros !. :
C, + dCj = (I + »7, (2.45) [00174] Uma comparação de (2.44) e (2.45) leva a:
Petição 870180135007, de 27/09/2018, pág. 60/85
57/70 όφ
ÕS = c RC1 = C ·C ·C (2.46)
ÕS ÕS [00175] (O lado direito de (2.46) é uma matriz anti-simétrica; isso se segue a partir de X(c0 · C0 )= 0.) [00176] Para um vetor de estado contendo posição, velocidade, atitude e desalinhamentos, as derivadas parciais são dadas pelo seguinte:
ÕXI ,1 = ÕVI ,1 = Õ! ,1 = ÕSI,! = I ^I ,0 ÕVj ,0 ι ,0 ÕSIfi
ÕXI 1 . = I · T
ÕV
1,0 ÕxI1i . Õ!0,j ÕVi,i,i Õ!0,j õ[! ]] = ~ ÕC70 · CT · ~T
R C0 R C0 C0
ÕS0ÕS = (I + fe ]])· C„-V
ÕS 0 v 0ÕS 'x.=(I+! ]])· C0
ÕS 0 0ÕS
C · X ^0 — 0 J ij
Co V 0
-IJ ij (2.47) [00177] Todas as outras componentes são nulas. As derivadas com respeito a relógio e diferença entre a leitura atual e uma de referência não são mostradas, mas elas não diferem das fórmulas para uma navegação baseada puramente em GPS.
[00178] Estas fórmulas parecem trabalhosas, mas há muitas oportunidades para se tornar uma computação mais eficiente. Por exemplo, múltiplos produtos de matrizes 3x3 1 podem ser simplificados, algumas computações podem ser reusadas, etc.
[00179] Esta matriz de propagação pode ser usada para a propagação de covariâncias. O vetor de estado ainda tem de ser propagado
Petição 870180135007, de 27/09/2018, pág. 61/85
58/70 usando-se as equações (2.29) e (2.30) para posição, onde a atitude é propagada usando-se as equações (2.21-2.22).
b. Propagação para Trás [00180] Em seguida, na etapa 250, algoritmos de propagação para trás são requeridos para uma propagação para trás, com uma matriz de propagação para trás, do vetor de estado por múltiplas épocas e para a computação de derivadas parciais do estado passado com respeito às componentes do estado atual. Uma propagação para trás é usada no processamento de ADR.
[00181] Embora as matrizes de estado e de covariância sejam propagadas para trás, isto não pode ser computado por incrementos no tempo para trás. Isto requereria o armazenamento de dados de IMU pela duração da propagação para trás em uma pilha. Se o tempo de propagação for de dezenas de segundos, e as taxas de IMU estiverem nas centenas de Hertz, as exigências de armazenamento e de processamento poderiam se tornar excessivas. Assim, uma rotina de propagação para trás é divisada, que trabalharia como se fosse para frente. Em outras palavras, a propagação para trás deve ser baseada na (e derivada a partir da) propagação para frente. Esta derivação é o assunto desta subseção. A propagação para trás é baseada nas fórmulas de propagação para frente (2.42), onde os valores de C e Vo, X são computados como uma função do estado no final do intervalo de propagação de época múltipla.
[00182] Em cada época, os valores de C0 e V0, x0 em (2.42) são para serem computados usando-se valores de parâmetros de desalinhamento S, os quais são válidos no fim do intervalo de propagação para trás de múltiplas épocas. Infelizmente, durante as épocas precedentes, este valor futuro é desconhecido. Para consideração desse valor futuro de S, devem-se estimar os valores de C0 e V0, x0 para
Petição 870180135007, de 27/09/2018, pág. 62/85
59/70 algum S fixo (por exemplo, S = 0) e as derivadas dos mesmos com respeito a S. Quando o valor de S se torna disponível ao final do intervalo de propagação, os valores de C0 e V0, X0 podem ser aproxi mados por uma expansão de Taylor de dois termos usando-se valores em S = 0 e derivadas com respeito a S.
i. Matriz de Propagação para Trás [00183] Para se usarem de forma ótima medições de ADR, a matriz de propagação para trás é computada por múltiplas épocas. Assim, tais quantidades são computadas como:
ZXi ,0 Ζχ ,0 ,0 ÕVj ,0 ÕXj ,0,i ÕVj ,0, i ,M ZVI ,M ZXI ,M ZVI ,M ΖΨΐ ,M, j ΖΨΐ ,M, j
Figure BRPI0515452B1_D0011
onde M é o número de épocas indo para trás.
[00184] Para se ir de M -1 para M, as equações modificadas (2.42) são usadas como um ponto de começo:
VI ,M VI ,M-1 XI ,M XI ,M-1 + VI ,M-1 T + (I + [!ll)· CM-1 · XM-1 (2.48) e computam-se derivadas destas fórmulas com respeito a componentes do vetor de estado na época M. O vetor de estado compreende quatro divisões: X = {r, ,φ, S}, onde a identidade zX ~ . = I. A matriz de atitude CM-1 é estampada no tempo na época
M -1; daí, para as finalidades de computação de derivada, deve ser considerado um resultado da propagação para trás da matriz de atitude a partir da época M para M -1. Antes de se obterem as fórmulas para a matriz de propagação para trás, duas fórmulas auxiliares são derivadas abaixo.
ii. Computação Auxilar: Derivadas da Matriz de Atitude [00185] Para a matriz de atitude, o que vem a seguir é usado:
Cm = Cm-,-Cm-1 (2.49)
Petição 870180135007, de 27/09/2018, pág. 63/85
60/70 onde CM_j é dado por (2.21) (note que (2.21) descreve uma propagação por uma subépoca, e deve ser repetida múltiplas vezes para a obtenção de CM _1, a qual descreve a propagação por uma época). Isto leva a:
c = c · CT CM-1 CM CM-1 (2.50) õCM-1 = C õS CM
M-1
FCm l õS j
M-1 , γ
M-1
M-1 M-1 (2.51) iii. Computação Auxiliar:
δφ ~dS [00186]
Isto é similar à computação prévia de δφ —, mas agora as õS derivadas são computadas para uma propagação para trás, usandose:
C = c · CT
-1 cm CM-1
M-1 (2.52) [00187] Se houver um desalinhamento, a matriz de atitude final diferirá do valor nominal:
CM-1 + dCM-1 CM ~ í = c · CT CM CM-1 l
, õCM-1 õS λ
· dS (2.53) j
[00188] Por outro lado, esta mudança na atitude é modelada pelos parâmetros φ:
Cu-, + dCu_,=(i+Hd«ÍM-iID· Cu - (2.54)
Comparar (2.53) e (2.54) leva a:
Figure BRPI0515452B1_D0012
T (2.55) [00189] Agora, as derivadas para a matriz de propagação para trás podem ser escritas conforme se segue:
Petição 870180135007, de 27/09/2018, pág. 64/85
61/70 ÕXI ,M-1
I ,M-1 Õ!I,M-1 = ÕSI,M-1 = I Õ!I ,M ÕSI ,M
I ,M-1
I ,M-1
Ί ,M-1
I ,M-1
C
CM-1
Vm-1
I ,M-1
ÕS õViM sm
ÕXI M 1 sm
I ,M-1
Ί ,M-1 CM-1 ’ (x
M-1 (±M-1 = -(i+IW)· C,
- T Vm -1) õCm-1 C T
ÕS CM-1
ÕV 'V.\í-1 - C C R CM-1 Cl
ÕSM-1 r L_M-1 - C C
M-1 CM-1 CM-1 ÕSM-1 õC UCM-1 . V
ÕS V tT
M-1 V VM-1
T - (í + !]]) Cm-1 fX^ - Cm-1 Cm-1 ÕSM ÕSM-1 (2.56)
C C xCM -1 x
CM-1 CM-1 X
ÕS
M-1 χ XM-1 [00190]
Todas as outras componentes são nulas.
[00191]
Armados com estas fórmulas, as fórmulas recursivas para ÕX podem ser reescritas
ÕXm
ÕX 0 = ÕX 0 ÕXm ÕXm ÕXm-1 ÕXm
M-1 (2.57)
M-1
c. Propagação para Trás do Vetor de estado
Um algoritmo recursivo é definido:
6M 6M-1 + C0 ’ CM-1 ’ V
M-1
M-1 !—M-1 XM Xm-1 + T 6-1 + CT C X ~C0 CM-1 x
M-1 —M-1 (2.58) 'M-1 onde 60 = 0, 50 = 0 . O produto CT CM-1 nada mais é do
M-1 que o produto 4Cm de matrizes individuais Cm.
m=0 [00192] A resolução da primeira equação de (2.58) para CM-1 VM-1 e a substituição do resultado na primeira equação de (2.48) leva ao seguinte:
Petição 870180135007, de 27/09/2018, pág. 65/85
62/70 VM - VM-1 = (l + Ml· C0 -(6M ~6M-1 ) (2·59) [00193] O somatório de (2.59) por M produz a fórmula de propagação para trás desejada para V0:
Vo = VM -(I+ M1D-Co (2.60) [00194] De modo similar, a resolução da segunda equação de (2.58) para CM-1 · xM-1 e a substituição do resultado na segunda equação de (2.48):
XM - XM-1 = TVM-1 +(I + Ml· C0 ·(5Μ ~Xm-1 - T 6M-1 ) (2-61) [00195] Em seguida, a equação (2.60) é usada, mudando-se o subscrito ali de M para M -1, resolvendo-se a equação resultante para VM-1 e substituindo em (2.61) para a obtenção de:
XM - XM-1 = T · V0 +(I + Ml· C0 ·(Χμ -Xm-1 ) (2-62) [00196] O somatório de (2.62) por M produz uma fórmula de propagação para trás desejada para x0:
X0 = XM -M · T· V0 -(I + Ml· C0 ·ΧΜ (2-63) [00197] As fórmulas (2.60) e (2.63), onde ΰΜ, xM são computados pela aplicação de (2.58), servem como um fundamento do algoritmo de propagação para trás.
[00198] Assim, conforme indicado na etapa 260, a matriz de propagação para trás é usada para a formação das equações de medição de filtro de Kalman, e, na etapa 270, o filtro de Kalman processa as medições de PR e as diferenças de ADR e atualiza o vetor de estado.
C. Ilustração de Performance de Navegação [00199] As FIGURAS 7 e 8 ilustram resultados de simulação de erro de navegação (erro de posição em metros, m, e erro de velocidade em metros/segundo, m/s) como uma função para algoritmos convencionais (usando apenas medições de PR), se comparado com técnicas de processamento de PR + ADR descritas acima em conjunto com a
Petição 870180135007, de 27/09/2018, pág. 66/85
63/70
FIGURA 5. As circunstâncias desta simulação são conforme se segue: [00200] Um satélite em uma órbita geoestacionária.
[00201] Os sinais de GPS recebidos com S/N0>27 dB*Hz. Isto corresponde ao acompanhamento de sinais no lobo principal da antena de transmissão de GPS e no pico do primeiro lobo lateral.
[00202] Um erro de acompanhamento de código (PR) é gaussiano, não correlacionado e calibrado para ter um desvio padrão de 1 m para um sinal no eixo ótico de uma antena direcional da antena de GPS. Para sinais mais fracos, o erro de acompanhamento é escalonado por 1/raiz quadrada de S/N0.
[00203] O erro de efeméride é gaussiano, altamente correlacionado com o tempo, tem um desvio padrão de 1,4 m para medições de pseudodistância.
[00204] Um erro de acompanhamento de portadora (ADR) é gaussiano, não correlacionado e calibrado para ter um desvio padrão de 1/100 daquele para o código.
[00205] O relógio de usuário tem diferença entre a leitura atual e uma de referência e um erro de passada randômico com um desvio padrão de 15 cm/s.
[00206] O traço 1 é o erro de posição para um processamento de PR convencional. O traço 2 é o erro de posição usando-se o processamento de PR + ADR. O traço 3 é o erro de velocidade para um processamento de PR convencional. O traço 4 é o erro de velocidade para um processamento de PR + ADR. O melhoramento no erro de navegação das técnicas de processamento de PR + ADR descritas aqui em relação a uma navegação por GPS usando uma medição de PR convencional apenas é aproximadamente de uma ordem de magnitude.
[00207] A FIGURA 8 mostra uma porção aumentada para os dois primeiros minutos dos gráficos mostrados na FIGURA 7. Na FIGURA
Petição 870180135007, de 27/09/2018, pág. 67/85
64/70
8, o processamento de ADR começa em t = 5. A queda aguda no erro de navegação de velocidade imediatamente após o começo do processamento de ADR (traço 4) é muito evidente.
[00208] Com respeito aos algoritmos de navegação de GPS - IMU descritos aqui, foram realizadas simulações segundo as condições a seguir:
[00209] A IMU produz dados ideais (ruído e erros sistemáticos em relação a um valor de referência são zero).
[00210] As medições de pseudodistância têm erro gaussiano com desvio padrão de 1 m.
[00211] As medições de ADR têm erro gaussiano com desvio padrão de 1 cm.
[00212] As FIGURAS 9 a 10 mostram os resultados destas simulações, onde os erros de navegação de raiz quadrada da soma dos quadrados dos erros individuais (RSS) são mostrados como uma função do tempo para posição (metros), velocidade (m/s) e atitude (componentes da matriz de atitude, aproximadamente corresponde a radianos). A FIGURA 9 mostra os erros para estes parâmetros, quando apenas medições de PR foram processadas e o algoritmo de propagação foi limitado à primeira ordem. A FIGURA 10 mostra o erro para estes parâmetros usando-se um algoritmo de propagação de segunda ordem usando apenas medições de PR. A FIGURA 11 mostra o erro para estes parâmetros usando-se um algoritmo de propagação de terceira ordem para medições de PR e diferença de ADR. Estas figuras mostram uma ordem dois de magnitude de melhoramento na atitude e estimativa de velocidade e uma ordem de magnitude de melhoramento na estimativa de posição, a qual usa um processamento de ADR em conjunto com algoritmos de propagação acurados.
[00213] Um algoritmo de propagação acurado é um fator chave para se possibilitarem novos métodos de processamento de ADR para
Petição 870180135007, de 27/09/2018, pág. 68/85
65/70 aplicações espaciais. Os métodos descritos aqui obtêm uma precisão extraordinária na medição da velocidade de uma espaçonave. As estimativas para uma espaçonave GEO prometem a obtenção de acurácias de medição da ordem de 0,1 mm/s para cada medição individual. [00214] As equações de propagação para um sistema de GPS/IMU integrados podem ser usadas de uma forma similar. Uma qualidade melhor de propagação permite que se use uma diferença em medições de ADR tomadas muitas épocas separadas, para uma acurácia maior. Uma diferença entre uma propagação no espaço e uma propagação direcionada por IMU é que no último caso as medições de ADR se aplicariam a uma estimativa de velocidade e de atitude.
III. Miscelânea de Derivações
A. Derivação de Fórmula para G
Por toda esta derivação, as identidades a seguir são usadas:
[[u ]· v = U v = -\μ x v]
Figure BRPI0515452B1_D0013
(3.1) [a x [b x c ]| = b (a c)- c (a b )
1. Expressão 1 [00215] Em primeiro lugar, é computado:
A = Φ ·ρ]|· x as onde x é um vetor arbitrário. Então,
I = a[(^ pV x] = as a[x x (s p)] a[x x [s x p]| as as
{- S ·(ρ· X )+p-(s x)} [00216] O vetor X é regulado igual ao k-ésimo vetor unitário do
Petição 870180135007, de 27/09/2018, pág. 69/85
66/70 quadro de referência. Então, as derivadas podem ser computadas conforme se segue:
õfrp I = õs,
-{- Sp, + p,St }= (3.4) ÕS,
S, Pi - S p,
2. Expressão 2 [00217] A próxima etapa é computar a expressão a seguir:
B .-]v pj., r uψxp
ÕS ÕS
-F S-(x-[s xp|)-[s xp|-(x ·S)}
Novamente, x é regulado igual ao k-ésimo vetor unitário do quadro de referência:
- P * {sjs xp|.-[S xpj-s. }(3.6)
ÕS ÕS l· (3.5)
Então, õ|[S Õc p|11 = S, [S x pj. + Sp, -S, · [S x pj - S,p, (3.7)
3. Expressão Final [00218] Com estes blocos de construção, o que vem a seguir é computado:
õ(p t )* ÕS õ
ÕS r ' r l '· r p + S-p + - S2-p (3.8)
Sp.-S„p, + Ife. -[s ,j + S,p,-S -[s\ -sp,}
B. Derivação de x, fórmulas de propagação de V [00219] Começando a partir de:
Petição 870180135007, de 27/09/2018, pág. 70/85
67/70 t„
In_1 = J(T ~α^ )· C (t )· a (t )· dt (3.9) tn-1 [00220] Esta integral pode ser usada para a derivação de ambas as fórmulas de propagação de posição e de velocidade. Para a posição, a = 1 no resultado final. Para a velocidade, faça a = 0, T = 1. Usando-se o método do trapezóide:
!n-1 = d[_ · C(tn-1 )^[(T - atn-1 ) a (tn-1 )+ (T - atn )’ ά(t n-1 )' a (t n )] + o(dt' ) (3'10) [00221] Em seguida, a mesma quantidade é estimada, mas usandose dois meios incrementos:
1 n-1 = · C (tn-1 )·[(Τ -a tn-1 ) ' a (t n-1 )+ 2(T - a t n-1 /2 ) ' / 2 (t n-1 ) ' a (t n-1 /2 )+(t n a (t n ) (T ' tn )] + — o(dt3) v 7 (3.11) onde âC1/2 (tn_1) é a atualização de matriz de atitude de tn-1 para tn-1/2. O segundo termo J é denotado e dividido em duas partes:
n-1/2 J 2(T atn-1/2 )<C1/2 (tn-1 )’ a (T - atn )<C1/2 (tn-1 )’ a n-1/2 H (T - atn-1 )âC1/2 (n-1 )’ a n-1/2 (3.12) [00222] São substituídas em J duas expressões separadas para as duas ocorrências de a n_1/2 (a partir deste ponto, omitindo os subs critos n, n -1/2, n -1 nos termos, onde o subscrito exato não importa no limite de acurácia de meta):
n-1/2 n-1 dt
H--a n , n-1 dt2 ..
H--a dt n-1/2
--a , n 2 n-1
3dt2 ..
---a (3.13) [00223] São substituídas duas expressões separadas para âC1/2(tn-1). A primeira expressão é a partir da expressão geral para âC:
~ / x T dt ( dt . õ dt2 2 â^1/2 (tn-1 )= 1 - ~ 2Pn-1 H ~ P n-1 8 P (3.14)
Petição 870180135007, de 27/09/2018, pág. 71/85
68/70 [00224] A segunda expressão para óC1/2(tn_.) é obtida a partir da identidade:
«V )Cp ) = SC(i„-t) (3.15)
Então:
C 2 7 )-£(t,_, ) [óC,,, (l„_,, 2 )]-1 (3.16)
Para [óC,,, (t_,,2)]', o que vem a seguir é usado:
[óC1/2 (tn-1/2 )^ - I + ~ {Pn-112 + Pn-112 ) + 8 P o
T dt í 3dt . j dtt· 1+7 [ 2 P /“T p (3.17) [00225] Todas estas substituições levam a:
τ h, \e-Pi \( T dt 3dt2 . dt2 2 J -(T-atn (t n-1 )·| 1 + — Pn-1 +— P + ~ P
1' / dt dt2 . dt2 2 Y dt .
+ (T - ;tn-1 ) 1 - — Pn-1--— Pn-1 + ~ P II X + ~ X-1 dt an -an-1 dt2 ,.j +--a | 8 ) (3.18)
Retendo os termos com dt1:
. \ í dt . 3dtX dt dt2 . 3dt2 . dt2 2 J -(T -atn 0 (tn-1 ) · I an - ~ an-1--~ a + ~ Pn-1 an---- Pn-1C. - 8 P^ — P ü, / d^ dt X dt dt2 dt2 . dt2 2 j + T -at . ) a . +--a . +--a--p .a .--p .a .--pa +--Pa |
X n-1/1 n-1 2 n-1 o 2' n-1 n-1 4 > n-1 n-1 8 ' 8 ' I (3.19) n-1 [00226] A extração do termo principal e de todos os outros termos estimados para a segunda ordem com respeito dt. Para esta computação, SC(tn-1 )« I - dtP, etc. é usado:
J - atn j.C (tn-1 ) · a n + (T - )a„-. }+ , , v , / dt . 3dtX dt dt2 . 3dt2 . dt2 2 j (3.20) (T ~atn-1 ~adt )(1 - dtP)I - — a„-1--— a + — P„-1a„ —— P„-a„-1^^- pá 8 P an I / T dt . dt X dt dt2 dt2 . dt2 2 j + (T - atn-1 )a„-1 +a - — Pn-1an-1 PXX1 + V P a I n-1 [00227] Mais uma vez, apenas os termos os quais têm a segunda ordem (ou menor) com respeito a dt são retidos:
Petição 870180135007, de 27/09/2018, pág. 72/85
69/70
J =j 1 (T -at, (tn-1 ) dt2 ( \ . 1 .. 1 .
(T -atn-l) pa- —a +—pa an + 2(T -;tn-1 )a d 1.. 1 . 1 2 Ί ' , -p an v 2 2 2 n) n-1 + a-(à - pa) (3.21) [00228] A substituição desta expressão na segunda equação (aquela obtida a partir da fórmula de dois meios incrementos) para in-1 leva a:
I„-1 = C (tn-1 ){ 2 (T-atn )W (tn-1) C(t n-1 (T ~atn-1 ~ 2 a + 2 /&a ,h-3 r Cw/ 8 &
[00229] a:
o(dt3 )= d3 V ’ 6 an + — n 2 ’ 1 2 — p an
2 2
Uma comparação com uma fórmula de um incremento leva (T ~atn-1 )an-1 j + ^ + a^(a - pa) +1 o(dt3) (3.22) dt3 AZ \ / . 1 .. 1 . 1 ~6~ C (tn-1 )· (T -atn-1 - 2 a + 2 p&a - 2 o que produz a expressão final:
I.-1 = C (t.-1 )·{ 2 (T-at. )W (t„-1 ) (tn-1 (T-atn-1 - 2 a + 2 &a -
2^2 an ^ + a-(<a ~pa) (3.23) n-1 • 3 an + 2 (T -atn-1 )a 1 2 Ί — p a + a
2 2 ) n-1 (a - pa) (3.24) [00230] Para resumir, é provido um método para a determinação dos parâmetros de navegação de um objeto em movimento, compreendendo a feitura de medições de pseudodistância (PR) no objeto a partir dos sinais recebidos a partir de fontes de sinal medição de distância; a feitura de medições de delta de distância acumulada (ADR) no objeto a partir de sinais recebidos a partir de fontes de sinal medição de distância; a computação de diferenças de ADR entre medições de ADR separadas por um intervalo de tempo que é maior do que o intervalo de tempo entre medições de ADR consecutivas; e a estimativa de pelo menos um parâmetro de navegação do objeto em movimento a partir das medições de PR e das diferenças de ADR.
[00231] De modo similar, é provido um método para a determinação
Petição 870180135007, de 27/09/2018, pág. 73/85
70/70 da posição e da velocidade de um objeto em movimento, que compreende a feitura de medições de pseudodistância (PR) no objeto a partir dos sinais recebidos a partir de fontes de sinal medição de distância; a feitura de medições de delta de distância acumulada (ADR) no objeto a partir de sinais recebidos a partir de fontes de sinal medição de distância; a computação de diferenças de ADR entre medições de ADR separadas por um intervalo de tempo que é maior do que o intervalo de tempo entre medições de ADR consecutivas; e a estimativa de componentes de um vetor de estado compreendendo pelo menos a posição e a velocidade do objeto em um estado atual a partir das medições de PR e das diferenças de ADR e do vetor de estado em um estado anterior.
[00232] Mais ainda, a metodologia descrita aqui pode ser concretizada em um processador ou um meio legível por computador que armazena instruções (por exemplo, um software) que, quando executadas por um processador ou por um computador, fazem com que o processador ou computador: obtenha medições de pseudodistância (PR) no objeto a partir dos sinais recebidos a partir de fontes de sinal medição de distância; obtenha medições de delta de distância acumulada (ADR) no objeto a partir de sinais recebidos a partir de fontes de sinal medição de distância; compute diferenças de ADR entre medições de ADR separadas por um intervalo de tempo que é maior do que o intervalo de tempo entre medições de ADR consecutivas; e estime pelo menos um parâmetro de navegação do objeto em movimento a partir das medições de PR e das diferenças de ADR.
[00233] O sistema e os métodos descritos aqui podem ser concretizados em outras formas específicas, sem se desviar do espírito ou das características essenciais dos mesmos. As modalidades precedentes devem ser consideradas, portanto, em todos os aspectos como ilustrativas e não têm por significado serem limitantes.

Claims (24)

  1. reivindicações
    1. Método de determinação de parâmetros de navegação de um objeto em movimento, caracterizado pelo fato de que compreende:
    (a) a feitura de medições de pseudodistância (PR) no objeto a partir dos sinais recebidos a partir de fontes de sinal medição de distância;
    (b) a feitura de medições de delta de distância acumulada (ADR) no objeto a partir de sinais recebidos a partir de fontes de sinal medição de distância;
    (c) a computação de diferenças de ADR entre medições de ADR separadas por um intervalo de tempo que é maior do que o intervalo de tempo entre medições de ADR consecutivas; e (d) a estimativa de pelo menos um parâmetro de navegação do objeto em movimento a partir das medições de PR e das diferenças de ADR.
  2. 2. Método de acordo com a reivindicação 1, caracterizado pelo fato de que ainda compreende a determinação de um intervalo de tempo ótimo entre medições de ADR a ser usado para a computação das diferenças de ADR.
  3. 3. Método de acordo com a reivindicação 2, caracterizado pelo fato de que a determinação do intervalo de tempo ótimo entre medições de ADR compreende a computação de um mínimo de uma variância de um erro de estimativa de velocidade para duas medições de ADR.
  4. 4. Método de acordo com a reivindicação 3, caracterizado pelo fato de que a determinação do intervalo de tempo ótimo é dependente de uma acurácia do modelo de propagação, de modo que o intervalo de tempo ótimo possa ser mais longo quando um modelo de propagação mais acurado for usado.
    Petição 870180135007, de 27/09/2018, pág. 75/85
    2/6
  5. 5. Método de acordo com a reivindicação 1, caracterizado pelo fato de que a estimativa compreende a computação de derivadas parciais das diferenças de ADR com respeito a pelo menos uma componente de velocidade de um vetor de estado que representa posição e velocidade do objeto para a derivação de estimativas de velocidade instantânea do objeto.
  6. 6. Método de acordo com a reivindicação 1, caracterizado pelo fato de que a estimativa compreende a computação de derivadas parciais das diferenças de ADR com respeito a pelo menos uma componente de velocidade de um vetor de estado que representa posição e velocidade do objeto para a derivação de estimativas de velocidade instantânea do objeto e propagação do vetor de estado com equações de propagação que computam o vetor de estado para o objeto no tempo N, com base no vetor de estado no tempo N-1.
  7. 7. Método de acordo com a reivindicação 6, caracterizado pelo fato de que a computação de derivadas parciais compreende a computação de derivadas parciais das diferenças de ADR com respeito à posição e à velocidade do vetor de estado; a computação de derivadas parciais das diferenças de ADR com respeito a um relógio e um desvio de relógio (diferença entre a leitura atual e uma de referência em um relógio) do objeto; e a computação de estimativas de velocidade instantânea do objeto pelo uso destas derivadas parciais em um filtro de Kalman para processamento das diferenças de ADR.
  8. 8. Método de acordo com a reivindicação 7, caracterizado pelo fato de que a estimativa da posição do objeto em movimento compreende a computação das componentes de posição e velocidade do vetor de estado usando-se as medições de PR e as medições de velocidade instantânea estimada.
  9. 9. Método de acordo com a reivindicação 8, caracterizado pelo fato de que a computação de derivadas parciais das diferenças
    Petição 870180135007, de 27/09/2018, pág. 76/85
    3/6 de ADR compreende a computação de derivadas da posição de objeto em uma época passada com respeito a um vetor de estado em uma época atual, onde o vetor de estado compreende componentes para pelo menos uma posição e uma velocidade do objeto em uma época.
  10. 10. Método de acordo com a reivindicação 9, caracterizado pelo fato de que a computação de derivadas da posição de objeto em uma época passada com respeito ao vetor de estado em uma época atual compreende a computação de uma aproximação de ordem mais alta com respeito a uma diferença de tempo entre épocas consecutivas.
  11. 11. Método de acordo com a reivindicação 10, caracterizado pelo fato de que a computação de uma aproximação de ordem mais alta compreende a propagação do vetor de estado para trás por um intervalo de tempo separando medições de ADR usando-se uma rotina de propagação que tem múltiplas subetapas; para cada uma das subetapas da rotina de propagação, a computação de uma matriz φ que corresponde a um elemento de uma matriz Λ; a computação da matriz Λ para cada uma das matrizes φ; usando-se a matriz Λ, a computação de parciais da posição de objeto em uma época passada com respeito à posição e à velocidade do vetor de estado.
  12. 12. Método de acordo com a reivindicação 11, caracterizado pelo fato de que ainda compreende a propagação do vetor de estado usando-se um modelo de força e uma rotina de integração, e a propagação da matriz de covariância pela computação de uma derivada parcial do vetor de estado na próxima época com respeito ao vetor de estado na época atual.
  13. 13. Método de acordo com a reivindicação 12, caracterizado pelo fato de que o vetor de estado tem partições para posição, velocidade, relógio e desvio de relógio do objeto, onde a propagação do vetor de estado compreende a computação de uma matriz que tem
    Petição 870180135007, de 27/09/2018, pág. 77/85
    4/6 partições para derivadas de posições com respeito à velocidade, derivadas de posição com respeito ao relógio, derivadas de posição com respeito ao desvio de relógio, derivadas de velocidade com respeito à posição, derivadas de velocidade com respeito ao relógio e derivadas de velocidade com respeito ao desvio de relógio.
  14. 14. Método de acordo com a reivindicação 1, caracterizado pelo fato de que (c) a computação compreende a computação das diferenças de ADR entre medições de ADR em duas épocas separadas pelo referido período de tempo, o qual é maior do que o período de tempo entre épocas consecutivas nas quais medições de ADR são feitas.
  15. 15. Método de acordo com a reivindicação 14, caracterizado pelo fato de que (c) a computação compreende a computação de diferenças de ADR entre medições de ADR em duas épocas separadas pelo mesmo período de tempo, o qual é várias ordens de magnitude maior do que o período de tempo entre épocas consecutivas nas quais as medições de ADR são feitas.
  16. 16. Método de acordo com a reivindicação 1, caracterizado pelo fato de que a estimativa ainda compreende a propagação do vetor de estado com base em dados de medição inercial produzidos por uma unidade de medição inercial (IMU) associada ao referido objeto.
  17. 17. Método de acordo com a reivindicação 16, caracterizado pelo fato de que a propagação compreende a computação de uma matriz de observação que compreende derivadas parciais das diferenças de ADR com respeito a um vetor de estado para o objeto, onde o vetor de estado compreende uma estimativa da posição, da velocidade, da atitude e estados de erro de dados de medição inercial produzidos por uma unidade de medição inercial (IMU) associada ao referido objeto.
  18. 18. Método de acordo com a reivindicação 17, caracteriza
    Petição 870180135007, de 27/09/2018, pág. 78/85
    5/6 do pelo fato de que a computação da matriz de observação compreende a computação de uma matriz de rotação que transforma um vetor representando dados de sensor de giroscópio em eixos geométricos de corpo da IMU em um sistema de coordenadas inerciais centralizadas na Terra (ECI) e a propagação da matriz de rotação por uma época, de modo a se propagar a componente de atitude do vetor de estado.
  19. 19. Método de acordo com a reivindicação 18, caracterizado pelo fato de que a computação da matriz de observação compreende a propagação das componentes de posição e de velocidade do vetor de estado.
  20. 20. Método de acordo com a reivindicação 19, caracterizado pelo fato de que ainda compreende a computação de um vetor de desalinhamento que caracteriza parâmetros de desalinhamento do sensor de giroscópio na IMU, incluindo o vetor de desalinhamento no vetor de estado de filtro de Kalman.
  21. 21. Método de acordo com a reivindicação 20, caracterizado pelo fato de que ainda compreende a computação de derivadas parciais do vetor de estado no final de cada época com respeito ao vetor de estado no começo da época, para a produção de uma matriz de propagação para frente, onde em cada época a partição de atitude do vetor de estado é representada por ângulos de Euler com respeito a uma estimativa atual da atitude.
  22. 22. Método de acordo com a reivindicação 21, caracterizado pelo fato de que ainda compreende a computação de uma matriz de propagação para trás por múltiplas épocas.
  23. 23. Método de acordo com a reivindicação 22, caracterizado pelo fato de que ainda compreende uma propagação para trás do vetor de estado com a matriz de propagação para trás usando um algoritmo recursivo.
    Petição 870180135007, de 27/09/2018, pág. 79/85
    6/6
  24. 24. Método de acordo com a reivindicação 1, caracterizado pelo fato de que (a) a feitura compreende a feitura de medições de PR no objeto a partir de sinais recebidos a partir de satélites de posicionamento global (GPS).
BRPI0515452A 2004-09-17 2005-08-10 método de determinação de parâmetros de navegação de um objeto em movimento BRPI0515452B1 (pt)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US61060904P 2004-09-17 2004-09-17
US11/129,423 US7490008B2 (en) 2004-09-17 2005-05-16 GPS accumulated delta range processing for navigation applications
PCT/US2005/028537 WO2006036321A1 (en) 2004-09-17 2005-08-10 Improved gps accumulated delta range processing for navigation applications

Publications (2)

Publication Number Publication Date
BRPI0515452A BRPI0515452A (pt) 2008-07-22
BRPI0515452B1 true BRPI0515452B1 (pt) 2018-12-11

Family

ID=35515678

Family Applications (1)

Application Number Title Priority Date Filing Date
BRPI0515452A BRPI0515452B1 (pt) 2004-09-17 2005-08-10 método de determinação de parâmetros de navegação de um objeto em movimento

Country Status (10)

Country Link
US (1) US7490008B2 (pt)
EP (1) EP1792201B1 (pt)
JP (1) JP4789216B2 (pt)
KR (1) KR101209667B1 (pt)
AT (1) ATE414916T1 (pt)
AU (1) AU2005290192B2 (pt)
BR (1) BRPI0515452B1 (pt)
CA (1) CA2580539C (pt)
DE (1) DE602005011158D1 (pt)
WO (1) WO2006036321A1 (pt)

Families Citing this family (30)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2888643B1 (fr) * 2005-07-18 2009-09-25 Airbus France Sas Procede et dispositif pour determiner la position au sol d'un mobile, particulier d'un avion sur un aeroport
US7551138B2 (en) * 2005-12-22 2009-06-23 L3 Communications Integrated Systems, L.P. Method and apparatus for signal tracking utilizing universal algorithm
US7535420B2 (en) * 2005-12-22 2009-05-19 L-3 Communications Integrated Systems L.P. Method and apparatus for signal tracking utilizing universal algorithm
US7489271B2 (en) 2006-03-22 2009-02-10 Lockheed Martin Corporation Optimized receive antenna and system for precision GPS-at-GEO navigation
US7612712B2 (en) 2006-04-25 2009-11-03 Rx Networks Inc. Distributed orbit modeling and propagation method for a predicted and real-time assisted GPS system
US8125382B2 (en) * 2006-04-25 2012-02-28 Rx Networks Inc. Autonomous orbit propagation system and method
GB2440572A (en) * 2006-08-01 2008-02-06 Roke Manor Research Method and apparatus for controlling a clock and frequency source at a receiver
US7564406B2 (en) * 2006-11-10 2009-07-21 Sirf Technology, Inc. Method and apparatus in standalone positioning without broadcast ephemeris
US7855678B2 (en) * 2007-05-16 2010-12-21 Trimble Navigation Limited Post-mission high accuracy position and orientation system
US20090093959A1 (en) * 2007-10-04 2009-04-09 Trimble Navigation Limited Real-time high accuracy position and orientation system
JP2009229065A (ja) * 2008-03-19 2009-10-08 Toyota Motor Corp 移動体用測位装置
US20090299929A1 (en) * 2008-05-30 2009-12-03 Robert Kozma Methods of improved learning in simultaneous recurrent neural networks
US8165728B2 (en) * 2008-08-19 2012-04-24 The United States Of America As Represented By The Secretary Of The Navy Method and system for providing a GPS-based position
KR101384487B1 (ko) 2010-05-14 2014-04-10 한국전자통신연구원 위성항법 수신기의 의사거리 검증 방법 및 장치
WO2012021760A2 (en) * 2010-08-12 2012-02-16 The Government Of The United States Of America As Represented By The Secretary Of The Navy Improved orbit covariance estimation and analysis (ocean) system and method
RU2460970C1 (ru) * 2011-04-04 2012-09-10 Федеральное государственное унитарное предприятие "Центральный научно-исследовательский институт машиностроения" (ФГУП ЦНИИмаш) Способ определения эфемеридной информации в аппаратуре потребителя и устройство для его осуществления
EP2555017B1 (en) * 2011-08-03 2017-10-04 Harman Becker Automotive Systems GmbH Vehicle navigation on the basis of satellite positioning data and vehicle sensor data
RU2474838C1 (ru) * 2011-08-19 2013-02-10 Открытое акционерное общество "Российская корпорация ракетно-космического приборостроения и информационных систем" (ОАО "Российские космические системы") Электронное устройство оперативного восстановления измерений псевдодальности
JP5977040B2 (ja) * 2012-02-17 2016-08-24 Necスペーステクノロジー株式会社 軌道位置推定方法、軌道位置推定装置及び軌道位置推定プログラム
US8457891B1 (en) 2012-06-19 2013-06-04 Honeywell International Inc. Systems and methods for compensating nonlinearities in a navigational model
US9128183B2 (en) * 2012-11-29 2015-09-08 Caterpillar Inc. Machine navigation system utilizing scale factor adjustment
WO2015047417A1 (en) * 2013-09-30 2015-04-02 Intel Corporation Filtering for global positioning system (gps) receivers
US9182236B2 (en) * 2013-10-25 2015-11-10 Novatel Inc. System for post processing GNSS/INS measurement data and camera image data
US10488505B2 (en) * 2014-05-30 2019-11-26 The Boeing Company Positioning in indoor locations and other GPS-denied environments
CN104915966B (zh) * 2015-05-08 2018-02-09 上海交通大学 基于卡尔曼滤波的帧率上变换运动估计方法及系统
EP3328731B1 (en) * 2015-07-28 2020-05-06 Margolin, Joshua Multi-rotor uav flight control method
EP3293549B1 (en) * 2016-09-09 2020-03-11 Trimble Inc. Advanced navigation satellite system positioning method and system using delayed precise information
US11747487B2 (en) * 2018-03-26 2023-09-05 Texas Instruments Incorporated GNSS receiver clock frequency drift detection
KR101982181B1 (ko) * 2018-08-30 2019-05-24 국방과학연구소 관성항법데이터를 이용한 비행정보 보정 장치 및 방법
CN111288993B (zh) * 2018-12-10 2023-12-05 北京京东尚科信息技术有限公司 导航处理方法、装置、导航设备及存储介质

Family Cites Families (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4646096A (en) * 1984-10-05 1987-02-24 Litton Systems, Inc. Enhanced global positioning system Delta-Range processing
US5041833A (en) 1988-03-28 1991-08-20 Stanford Telecommunications, Inc. Precise satellite ranging and timing system using pseudo-noise bandwidth synthesis
US5109346A (en) 1990-02-01 1992-04-28 Microcosm, Inc. Autonomous spacecraft navigation system
US5430657A (en) 1992-10-20 1995-07-04 Caterpillar Inc. Method and apparatus for predicting the position of a satellite in a satellite based navigation system
US5430654A (en) 1992-12-01 1995-07-04 Caterpillar Inc. Method and apparatus for improving the accuracy of position estimates in a satellite based navigation system
US5606506A (en) 1993-04-05 1997-02-25 Caterpillar Inc. Method and apparatus for improving the accuracy of position estimates in a satellite based navigation system using velocity data from an inertial reference unit
US6175806B1 (en) * 1993-07-16 2001-01-16 Caterpillar Inc. Method and apparatus for detecting cycle slips in navigation signals received at a receiver from a satellite-based navigation system
US5477458A (en) * 1994-01-03 1995-12-19 Trimble Navigation Limited Network for carrier phase differential GPS corrections
US5748651A (en) * 1995-05-05 1998-05-05 Trumble Navigation Limited Optimum utilization of pseudorange and range rate corrections by SATPS receiver
US5726659A (en) 1995-09-21 1998-03-10 Stanford University Multipath calibration in GPS pseudorange measurements
US5787384A (en) * 1995-11-22 1998-07-28 E-Systems, Inc. Apparatus and method for determining velocity of a platform
US5935196A (en) 1997-06-11 1999-08-10 Itt Manufacturing Enterprises Technique for the use of GPS for high orbiting satellites
US6608589B1 (en) 1999-04-21 2003-08-19 The Johns Hopkins University Autonomous satellite navigation system
US6134484A (en) * 2000-01-28 2000-10-17 Motorola, Inc. Method and apparatus for maintaining the integrity of spacecraft based time and position using GPS
US6268823B1 (en) 2000-03-08 2001-07-31 Trimble Navigation Ltd Unconventional range navigation system with efficient update process
US20020008661A1 (en) * 2000-07-20 2002-01-24 Mccall Hiram Micro integrated global positioning system/inertial measurement unit system
US6420999B1 (en) 2000-10-26 2002-07-16 Qualcomm, Inc. Method and apparatus for determining an error estimate in a hybrid position determination system
US6622091B2 (en) 2001-05-11 2003-09-16 Fibersense Technology Corporation Method and system for calibrating an IG/GP navigational system
US6664923B1 (en) 2002-09-24 2003-12-16 Novatel, Inc. Position and velocity Kalman filter for use with global navigation satelite system receivers

Also Published As

Publication number Publication date
JP4789216B2 (ja) 2011-10-12
DE602005011158D1 (de) 2009-01-02
US7490008B2 (en) 2009-02-10
WO2006036321A1 (en) 2006-04-06
AU2005290192B2 (en) 2009-04-09
EP1792201A1 (en) 2007-06-06
CA2580539A1 (en) 2006-04-06
CA2580539C (en) 2011-08-02
KR101209667B1 (ko) 2012-12-10
KR20070059105A (ko) 2007-06-11
ATE414916T1 (de) 2008-12-15
JP2008513775A (ja) 2008-05-01
BRPI0515452A (pt) 2008-07-22
US20060195262A1 (en) 2006-08-31
EP1792201B1 (en) 2008-11-19
AU2005290192A1 (en) 2006-04-06

Similar Documents

Publication Publication Date Title
BRPI0515452B1 (pt) método de determinação de parâmetros de navegação de um objeto em movimento
US20220082386A1 (en) Vision-aided inertial navigation
Indelman et al. Factor graph based incremental smoothing in inertial navigation systems
Xiong et al. Robust Kalman filtering for discrete-time nonlinear systems with parameter uncertainties
Hasberg et al. Simultaneous localization and mapping for path-constrained motion
Vu et al. Centimeter-accuracy smoothed vehicle trajectory estimation
CN104181574A (zh) 一种捷联惯导系统/全球导航卫星系统组合导航滤波系统及方法
CN105180938A (zh) 一种基于粒子滤波的重力采样矢量匹配定位方法
CN104713559B (zh) 一种高精度sins模拟器的设计方法
Wendel et al. Time-differenced carrier phase measurements for tightly coupled GPS/INS integration
CN112146655A (zh) 一种BeiDou/SINS紧组合导航系统弹性模型设计方法
CN104359496A (zh) 基于垂线偏差补偿的高精度姿态修正方法
CN110926499B (zh) 基于李群最优估计的sins捷联惯性导航系统晃动基座自对准方法
EP4047311A1 (en) Range image aided ins with map based localization
CN102830415B (zh) 一种降维度的基于Carlson滤波算法的快速组合导航方法
Rouzaud et al. Rigorous integration of inertial navigation with optical sensors by dynamic networks
Zhou et al. Applying quaternion-based unscented particle filter on INS/GPS with field experiments
CN109741372A (zh) 一种基于双目视觉的里程计运动估计方法
Liu et al. Implementation and analysis of tightly integrated INS/stereo VO for land vehicle navigation
CN116380038A (zh) 一种基于在线增量尺度因子图的多源导航信息融合方法
Jiang et al. Statistical modeling of long-range drift in visual odometry
EP4047310A1 (en) Range image aided ins
Dupuis et al. GPS-based preliminary map estimation for autonomous vehicle mission preparation
Gul et al. GPS/SINS navigation data fusion using quaternion model and unscented Kalman filter
Boey Odometry Error Reduction in Wheelchair Using More Than One Sensor

Legal Events

Date Code Title Description
B25D Requested change of name of applicant approved

Owner name: ITT MANUFACTURING ENTERPRISES LLC (US)

B25A Requested transfer of rights approved

Owner name: EXELIS INC. (US)

B06A Patent application procedure suspended [chapter 6.1 patent gazette]
B09A Decision: intention to grant [chapter 9.1 patent gazette]
B16A Patent or certificate of addition of invention granted [chapter 16.1 patent gazette]

Free format text: PRAZO DE VALIDADE: 10 (DEZ) ANOS CONTADOS A PARTIR DE 11/12/2018, OBSERVADAS AS CONDICOES LEGAIS.

B21F Lapse acc. art. 78, item iv - on non-payment of the annual fees in time

Free format text: REFERENTE A 18A ANUIDADE.

B24J Lapse because of non-payment of annual fees (definitively: art 78 iv lpi, resolution 113/2013 art. 12)

Free format text: EM VIRTUDE DA EXTINCAO PUBLICADA NA RPI 2735 DE 06-06-2023 E CONSIDERANDO AUSENCIA DE MANIFESTACAO DENTRO DOS PRAZOS LEGAIS, INFORMO QUE CABE SER MANTIDA A EXTINCAO DA PATENTE E SEUS CERTIFICADOS, CONFORME O DISPOSTO NO ARTIGO 12, DA RESOLUCAO 113/2013.