Detecção e seguimento de animais marinhos

marcados

Martin Jens Apel

Dissertação para obtenção do grau de Mestre em

Engenharia Electrotécnica e de Computadores

Orientadores: Doutor Paulo Jorge Coelho Ramalho Oliveira

Doutor Carlos Jorge Ferreira Silvestre

Júri

Presidente: Doutor João Fernando Cardoso Silva SequeiraOrientador: Doutor Paulo Jorge Coelho Ramalho OliveiraVogais: Doutor João Pedro Castilho Pereira Santos Gomes

Fevereiro 2015

Resumo

Em tempos de crescente preocupação com a preservação do meio ambiente é cada vez mais importante oestudo do meio ambiente e dos seres que nele habitam.Este trabalho encontrou motivação no âmbito do projecto MAST/AM que pretende obter o desenvolvi-mento de uma ferramenta para detecção e seguimento de animais marinhos previamente marcados comtags acústicos que satisfaça as necessidades dos biólogos durante alguns cenários típicos de tracking activo.Desta forma, este trabalho pretende descrever uma estratégia de descodificação de tags comerciais, usadaspara marcar animais marinhos, com uma ferramenta portátil (Portable Tool) a ser usada por um mergu-lhador. A detecção baseia-se num filtro passa-banda IIR e permite identificar as diferenças dos temposde chegada do sinal aos hidrofones configurados numa estrutura USBL. Assim, é possível encontrar umaestimativa da direcção de chegada do sinal ao Portable Tool.O trabalho pretende descrever também uma possível estratégia de processamento posterior, offline, dosdados recebidos pelo sistema para obter uma estimativa do percurso percorrido pelo animal marcado.Apresenta-se assim uma estratégia com Filtro de Kalman que combina as medições obtidas por umSurface Robotic Tool e por um Portable Tool com um modelo linear do movimento do animal marinhopara obter uma estimativa óptima da posição do animal.Por fim, uma simulação mostra o principio de funcionamento das estratégias de detecção e de seguimentosimulando também a atenuação imposta pela propagação no meio aquático. Consegue-se assim avaliaro desempenho do sistema proposto em relação à precisão na localização do animal bem como no quediz respeito ao alcance do sistema durante a detecção de animais marcados e compara-lo assim comum sistema de localização comercial para tracking passivo. Mantendo as vantagens de tracking activo osistema proposto apresenta um desempenho comparável com o do sistema comercial, sendo ao mesmotempo mais ágil e fácil na sua operação.

Palavras chave: tags acústicas; filtro passa-banda IIR; estimação de direcção de chegada; Filtro deKalman

II

Abstract

Nowadays as the protection of the environment is getting more and more important, studies on individualspecies and their natural habitat are also getting an increasing attention.This work was motivated by the MAST/AM project whose objective is to develop a tool to detect andtrack the movement of previously tagged marine animals to accomplish with the necessities of biologistsin some frequent sceneries.In this way, this work intends to describe an identifying and decoding strategy of commercial acoustictags by using a Portable Tool to be used by divers. Detection is based on the use of a band-pass filter(IIR) and allows the identification of the differences of time of arrival of the signal to the hydrophonesarranged in a USBL structure. Thus, it is possible to obtain an estimate on the arrival direction of thesignal to the Portable Tool.The objective of this work is also to describe a possible strategy of offline processing to get an estimateof the movement of the tagged animal. Therefore, a strategy is presented, in which a Kalman Filtercombines the measurements obtained by a Surface Robotic Tool and by a Portable Tool with a linearmodel of the animal’s movements in order to get an optimal estimate on the animal’s position.Finally, a simulation tries to demonstrate the previously described strategies of detection and trackingtaking as well into account the attenuation of the signal due to its propagation in the aquatic medium.Therefore it is possible to compare its performance in what concerns the localization precision as well asthe detection range with a commercial passive tracking positioning system. In conclusion the proposedsystem shows a similar performance as the commercial system while being more flexible and easy in itsoperation due to its active tracking carcteristics.

Key words: acoustic tags; IIR band-pass filters; Direction-of-arrival estimation; Kalman filters

IV

Conteúdo

Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IIAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IVLista de Figuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VILista de Tabelas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VIIIAcrónimos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . X1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Impacto da acústica subaquática e de tagging em animais marinhos . . . . . . . . . . . . . 11.2 Marcação de animais e telemetria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2.1 Tags . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.2 Receptores e estruturas de sensores . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.3 Acústica subaquática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Detecção e Descodificação das tags comerciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.0.1 Sinal emitido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.0.2 Estratégias de detecção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.0.3 Determinação do Time Difference of Arrival (TDOA) . . . . . . . . . . . . . . . . . . 282.0.4 Solução adoptada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.0.5 Descodificação do ID e Checksum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.0.6 Estimação da direcção de chegada do sinal . . . . . . . . . . . . . . . . . . . . . . . . 35

3 Filtro de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.1 Exemplo de funcionamento do Filtro de Kalman . . . . . . . . . . . . . . . . . . . . . . . . 383.2 O Filtro de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3 Aplicação do Filtro de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4 Simulação do Portable Tool em funcionamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474.1 Cenários de simulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.1.1 Cenário apenas com o Portable Tool . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.1.2 Cenário com Portable Tool e Surface Robotic Tool . . . . . . . . . . . . . . . . . . . . 48

4.2 Estrutura da simulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.3 Resultados de Simulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

5 Conclusões e Trabalho Futuro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63A Modelo de absorção acústica do meio aquático de Francois e Garrison . . . . . . . . . . . . . . . 66B Tags comerciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

B.1 Tags comerciais tipo I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67B.2 Tags comerciais tipo II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67B.3 Tags comerciais tipo III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69B.4 Sinais gravados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

C Código fonte da simulação em MatLab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73C.1 start_simulation.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73C.2 sistem.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75C.3 usbl.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92C.4 channel.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108C.5 eventList.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112C.6 event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114C.7 fish . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117C.8 fish_model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121C.9 detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

VI

Lista de Figuras

Fig. 1.1: Audiograma para alguns cetáceos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2Fig. 1.2: Audiograma para o american shad, o peixinho dourado, o salmão atlântico, o bacalhau

atlântico e de uma espécie de tubarão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3Fig. 1.3: Correspondência de espécies exemplares com o tipo de tag apropriado . . . . . . . . 5Fig. 1.4: Distância máxima de detecção de uma tag dependendo da frequência e da potência

de emissão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8Fig. 1.5: Desempenho do sistema comercial de localização . . . . . . . . . . . . . . . . . . . . 10Fig. 1.6: Perfil da probabilidade de detecção com a distância de uma tag comercial típica . . . 10Fig. 1.7: Propagação de som como onda esférica próximo da fonte e propagação como onda

plana longe da fonte de origem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11Fig. 1.8: Propagação de som de forma cilíndrica em água com pouca profundidade . . . . . . 13Fig. 1.9: Propagação multi-trajecto em águas de pouca profundidade e sinal detectado no

receptor com réplicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15Fig. 1.10: Refracção da onda acústica na interface de dois meios com características acústicas

diferentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15Fig. 1.11: Refracção devido à profundidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16Fig. 1.12: Perfil típico da velocidade de som dependendo da profundidade . . . . . . . . . . . . 17Fig. 1.13: Perfis genéricos da velocidade do som . . . . . . . . . . . . . . . . . . . . . . . . . . . 18Fig. 1.14: Criação de zonas de sombra pela variabilidade da velocidade do som na água . . . . 18

Fig. 2.1: Sinal emitido pelas tags comerciais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19Fig. 2.2: Amostragem e interpretação no domínio da frequência . . . . . . . . . . . . . . . . . 21Fig. 2.3: Amostragem de uma sinusóide a diferentes frequências de amostragem e suas implicações 22Fig. 2.4: Upsampling do sinal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Fig. 2.5: Diagrama de blocos do filtro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Fig. 2.6: A transformada de Fourier de um impulso unitário no tempo corresponde à função

Sinc no domínio da frequência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Fig. 2.7: Comparação entre o sinal à entrada do filtro e às saídas de um filtro adaptado e de

um filtro passa-banda a diferentes valores de SNR . . . . . . . . . . . . . . . . . . . . . . 25Fig. 2.8: Diagramas de Bode do filtro Infinite Impulse Response Filter (IIR) . . . . . . . . . . 27Fig. 2.9: Diagrama de blocos da implementação de um IIR de segunda ordem em estrutura

biquad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Fig. 2.10: TDOA entre pings correspondentes em canais diferentes e entre pings de um mesmo

canal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Fig. 2.11: Convolução com estratégia Overlap-add . . . . . . . . . . . . . . . . . . . . . . . . . 30Fig. 2.12: Saída do filtro adaptado para sinal de entrada sinusoidal na presença de ruído e

trajecto secundário com atraso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Fig. 2.13: Arquitectura da estratégia de detecção de sinal . . . . . . . . . . . . . . . . . . . . . 32Fig. 2.14: Diagrama de fluxo da estratégia de detecção para um canal . . . . . . . . . . . . . . 33Fig. 2.15: Incidência da onda acústica em dois receptores projectados no plano XY . . . . . . . 36

Fig. 3.1: Função de densidade de probabilidade para uma estimativa da posição x, dada amedida z1 e conhecendo o desvio padrão da medida �z1 . . . . . . . . . . . . . . . . . . . 39

Fig. 3.2: Função de densidade de probabilidade para uma estimativa x, dada a medida z2 emcomparação à função de densidade de probabilidade no instante t1 representada a tracejado 40

VII

Fig. 3.3: Função de densidade de probabilidade que descreve a estimativa x com valor esperadoµ e desvio padrão �, resultante da combinação das medidas z1 e z2 representadas a tracejado 40

Fig. 3.4: Devido ao modelo dinámico do sistema, o pico da função de densidade de probabili-dade desloca-se em direcção do eixo x, mas com aumento da incerteza o desvio padrão dadistribuição também aumenta [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Fig. 3.5: Esquema do funcionamento do Filtro de Kalman e da aplicação das suas equações . 44

Fig. 4.1: Cenário típico de missão para o projecto Advanced Tracking and Telemetry Metho-dologies to Study Marine Animals (MAST/AM) com Surface Robotic Tool e Portable Tool 49

Fig. 4.2: Visualização dos pontos de maior aproximação para duas rectas sem intersecção . . . 50Fig. 4.3: Visualização do referencial geodésico, o referencial Earth Centered, Earth Fixed Co-

ordinate System (ECEF) e do referencial North-East-Down Coordinate System (NED)local . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Fig. 4.4: O ponto p em relação aos referenciais A e B . . . . . . . . . . . . . . . . . . . . . . . 53Fig. 4.5: Ângulos de Euler e sequência de rotação yaw-pitch-roll . . . . . . . . . . . . . . . . . 54Fig. 4.6: Funcionamento esquematizado do programa de simulação . . . . . . . . . . . . . . . 56Fig. 4.7: Comparação entre o sinal emitido pela tag D (esquerda) e o sinal emitido pela tag B

(direita) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Fig. 4.8: Detecção de um sinal da tag G nos quatro canais do USBL . . . . . . . . . . . . . . 58Fig. 4.9: Detecção de um peixe marcado. (Posições:Surface Robotic Tool - [0;0;0], Portable

Tool - [800;100;150], Peixe inicialmente - [-350;100;150]) . . . . . . . . . . . . . . . . . . . 60Fig. 4.10: Detecção de um peixe marcado. (Posições:Surface Robotic Tool - [0;0;0], Portable

Tool - [800;100;150], Peixe inicialmente - [400;400;400]) . . . . . . . . . . . . . . . . . . . 61Fig. 4.11: Os pontos da trajectória do peixe que deram origem a uma collisão são marcados com

um circulo vermelho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

Fig. B.1: Sinal gravado da tag ID:A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Fig. B.2: Sinal gravado da tag ID:B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Fig. B.3: Sinal gravado da tag ID:C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Fig. B.4: Sinal gravado da tag ID:D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Fig. B.5: Sinal gravado da tag ID:E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Fig. B.6: Sinal gravado da tag ID:F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Fig. B.7: Sinal gravado da tag ID:G . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

VIII

Lista de Tabelas

Tab. 1.1: Mapa de códigos das Tags comerciais . . . . . . . . . . . . . . . . . . . . . . . . . . . 6Tab. 1.2: Taxa de falsas detecções para um conjunto de tags . . . . . . . . . . . . . . . . . . . 7Tab. 1.3: Distâncias máximas a que as tags comerciais podem ser detectadas . . . . . . . . . . 7

Tab. 2.1: Filtro Finite Impulse Response Filter (FIR) vs. filtro IIR . . . . . . . . . . . . . . . . 26

Tab. 4.1: Comparação do alcance e dos Signal to Noise Ratio (SNR)s do sinal no emissore no receptor para algumas tags simuladas considerando modelo de propagação esférica.(Esta simulação assume uma profundidade de 50m, uma salinidade de 35 Practical SalinityUnits (p.s.u.) e uma temperatura de 15 �C) . . . . . . . . . . . . . . . . . . . . . . . . . . 58

IX

Acrónimos

IIR Infinite Impulse Response FilterFIR Finite Impulse Response FilterDSP Digital Signal ProcessorTDOA Time Difference Of ArrivalAWGN Additive White Gaussian NoiseMAST/AM Advanced Tracking and Telemetry Methodologies to Study Marine AnimalsADC Analog Digital ConverterPAT Pop-up Satellite Archival TagSOFAR Sound Fixing And RangingRMS Root Mean SquarePOST Pacific Ocean Shelf TrackingOTN Ocean Tracking NetworkAATAMS Australian Acoustic Tracking And Monitoring SystemFADIO Fish Aggregating Devices as Instrumented Observatories of pelagic ecosystemsEKF Extended Kalman FilterUKF Unscented Kalman FilterECEF Earth Centered, Earth Fixed Coordinate SystemNED North-East-Down Coordinate SystemTOA Time Of Arrivalp.s.u. Practical Salinity UnitsUSBL Ultra Short BaselineUW-ASN Underwater Acoustic Sensor NetworkSNR Signal to Noise RatioDOA Direction-of-arrival

X

Capítulo 1

Introdução

A marcação de animais é uma prática comum entre cientistas e investigadores que possibilita a identi-ficação e assim um estudo intensivo durante um período de tempo alongado. A marca pode ser apenasuma identificação com tinta ou com um anel (afixada ao maxilar do peixe), ou pode ser uma marcaelectrónica que permite a utilização de telemetria ou de gravação de dados. Assim se pode recolher infor-mação acerca do animal ou do seu ambiente. Uma marca electrónica ajuda à detecção e identificação deanimais à distância podendo dar informação sobre a sua localização ou dados relevantes para o estudo daespécie no seu habitat natural. O conhecimento de padrões de migração podem servir para a definiçãode áreas de protecção e planos de pesca. Sensores de temperatura do corpo, do ritmo cardíaco ou donível PH no estômago contribuem para monitorizar períodos de actividade e de alimentação ou permitemtirar conclusões sobre o estado de saúde do animal marcado mantendo-o no seu habitat natural. Algunsanimais marinhos como as baleias percorrem longas distâncias nos oceanos podendo assim transportarsensores que recolhem dados do meio ambiente como a salinidade, a pressão e a temperatura da águaou os níveis da luminosidade. Assim podem ser criados perfis de temperatura da água ou ou podem serobtidas informações sobre uma possível existência de algas ou partículas nos oceanos.

1.1 Impacto da acústica subaquática e de tagging em animais

marinhos

Estudos recentes e dados analisados retrospectivamente mostram uma associação entre o uso de sonaresde grande potência e a desorientação de baleias que vão dar à costa. Embora não tenha sido possívelprovar uma ligação directa, a proximidade temporal e geográfica dos acontecimentos indicam pelo menosalguma relação causal entre o uso de sonares de alta potência e meio alcance durante exercícios militarese várias baleias bicudas que deram à costa . Alguns casos dramáticos de baleias que deram à costa comoem 1996 na Grécia e no norte das Bahamas em Março de 2000 fomentaram a discussão pública sobre osefeitos do ruído de origem humana em animais marinhos. [2]As fontes de ruído podem provir de sonares militares e comerciais, de plataformas de petróleo, de torrespara a produção de energia eólica ou podem ser causados pela locomoção de navios de todo o tipo. Oimpacto dos ruídos sonoros no comportamento dos animais marinhos pode variar bastante. Existemrelatos de que baleias da Gronelândia mudaram os seus padrões de respiração e de emersão quandosentiram ruído, o que biológicamente não apresenta uma influência negativa significativa. Noutros casosa reacção é mais relevante, como a de baleias francas em resposta à presença de barcos , baleias daGronelândia em resposta a ruídos industriais ou cachalotes em resposta a sequências curtas de pulsosde pingers acústicos , que diminuem ou cessam a comunicação, a qual faz parte da interacção social e éusada para a localização. [3]Apesar de baleias e golfinhos terem ouvidos resistentes e se ter observado alguma forma de habituaçãoem alguns casos, a variabilidade das reacções dos animais marinhos não depende apenas da espéciemas também de condições ambientais, de factores como a idade e o sexo do animal ou da situaçãocomportamental na qual se encontram. Golfinhos que estejam a descansar evitam barcos, enquanto

1

golfinhos que se estejam a alimentar os ignoram. Quando se encontram numa situação de interacçãosocial dentro do grupo, golfinhos poderão até procurar a proximidade de um barco.De qualquer modo, um animal apenas responderá directamente a um som se o conseguir detectar. O graude distúrbio que um ruído pode causar num animal depende do volume do mesmo que é função não só dasua intensidade no local em que ele se encontra, mas também da sensibilidade do animal às frequênciascorrespondentes.Os audiogramas desenvolvidos e os modelos dos sistemas auditivos de animais marinhos apresentam umagrande variação na sua capacidade auditiva e de sensibilidade às frequências. Mamíferos marinhos no seuconjunto têm um raio auditivo que vai de uma frequência de 10 Hz até 200 kHz. Este grupo pode-sesubdividir em cetáceos infrasónicos (sensibilidade geral de 15 Hz até 20 kHz, boa sensibilidade de 20Hz a 2 kHz), com sensibilidade geral a frequências sónicas até altas-frequências (100 Hz até 100 kHz,com grande variabilidade dos picos de sensibilidade), e espécies que dominam o espectro ultrasónico (200Hz até 200 kHz de sensibilidade geral, pico de sensibilidade de 16 a 120 kHz). [3] Alguns audiogramasestão apresentados na figura 1.1.Embora não haja muitos dados acerca de outras espécies marinhas para além dos mamíferos, relativamenteaos peixes e aos seláquios (tubarões e raias) existem para estes indicadores que apontam para o facto deque uma grande maioria dessas espécies poderá ter capacidade auditiva para detectar sons de frequênciasabaixo de 50 Hz (alguns mesmo a partir de 10 ou de 15 Hz) até frequências de 500 a 1,000 Hz. Algumasespécies conseguiram adaptar-se para poderem detectar frequências até 3 kHz ou mesmo mais comuma sensibilidade maior que outras espécies não especializadas em frequências mais baixas. O peixinhodourado e o American shad são exemplos de espécies especializadas, enquanto que o salmão atlântico e obacalhau atlântico são exemplos de espécies sem especialização. Alguns audiogramas estão apresentadosna figura 1.2. [3]

Figura 1.1: Audiograma para alguns cetáceos [3]

2

Figura 1.2: Audiograma para o american shad, o peixinho dourado, o salmão atlântico, o bacalhauatlântico e de uma espécie de tubarão [3]

Embora se conheça bastante sobre o repertório acústico de mamíferos marinhos, pouco se sabe sobre o seusignificado e sobre o modo como os seus padrões variam durante o dia, sazonalmente e dependendo da suasituação geográfica. Os estudos existentes abrangem apenas áreas confinadas e períodos de tempo muitolimitados pelo que instituições como o Committee on Potential Impacts of Ambient Noise in the Oceanon Marine Mammals, National Research Council em Wahington, D.C. e outras exigem que seja feitoum maior esforço no desenvolvimento de equipamento de tagging e em estudos que possam trazer maisinformações sobre o comportamento dos animais. [3] Com esse conhecimento podem ser desenvolvidasáreas de protecção e estratégias para melhorar a convivência entre humanos e animais marinhos.Instituições como a Society for Marine Mammalogy publicam recomendações [4] segundo as quais oscientistas se devem reger quando efectuam estudos científicos em animais marinhos como complementoàs directivas e à legislação nacional e internacional.Um estudo científico deve ser devidamente preparado e eventualmente precedido de estudos preliminarespara que possa ser garantida uma resposta cientificamente significativa, minimizando ao mesmo tempopossíveis impactos negativos sobre o animal individual ou sobre a sua população e eco-sistema.Para tal, é necessário tomar decisões relativamente à selecção do animal, às técnicas de estudos a utilizar,duração do estudo, ao número de amostras a tomar bem como ao grau de importância do ponto de vistaestatístico dos resultados obtidos.Para realizar os estudos há que seleccionar espécies que não estejam em perigo de extinção e cujo com-portamento se coadune melhor com os objectivos de investigação, especialmente quando se utiliza equi-pamento e instrumentos que estejam ainda a ser testados.A vulnerabilidade de um animal relativamente a possíveis distúrbios varia com a sua idade, sexo e estadode reprodução, com a sua posição social e com a sua localização ou com a profundidade da água. Deve-seter em atenção os possíveis distúrbios especialmente no que diz respeito a animais em estado de reproduçãoou com crias como no caso dos mamíferos.O grau de amostragem tem que ter suficiente significado estatístico para levar a conclusões válidas.Uma amostragem superior àquela que seja necessária produz distúrbios desnecessários num número maiselevado de animais e por outro lado possíveis perdas de tags não deverão pôr em causa o significadoestatístico do estudo.A escolha do grau de amostragem deve obedecer a um equilíbrio entre o aspecto científico, ético e logístico.

3

Qualquer tipo de intervenção física no animal em estudo, seja do ponto de vista de anestésicos ou daintrodução de ferramentas técnicas para a realização do estudo, terão naturalmente de ser realizadospor pessoal técnico especializado. O conhecimento actualizado da literatura científica específica, dosrequisitos e das normas legais são obviamente imprescindíveis.

1.2 Marcação de animais e telemetria

Nos diversos estudos sobre animais marinhos a grande variabilidade dos ambientes bem como as par-ticularidades de cada espécie exigem que a tecnologia envolvida se adapte ao máximo à situação emestudo. Consequentemente existem diferentes equipamentos que têm vindo a ser desenvolvidos para asmais variadas aplicações. Alguns desses equipamentos vão ser apresentados em seguida.

1.2.1 Tags

Os tags electrónicos podem ser classificados em duas categorias. Os transmissores enviam remotamentedados a partir do animal marcado para um receptor por via sonora ou por ondas de rádio. Dataloggersarmazenam os dados numa memória local e têm que ser recuperados do animal para obter os dados.Geralmente Dataloggers proporcionam uma maior quantidade de dados do que os transmissores mas arecuperação dos dataloggers por vezes não está garantida.Alguns tags electrónicos combinam num só aparelho as vantagens dos dois. Por exemplo no caso do Pop-up Satellite Archival Tag (PAT) um grande volume de dados pode ser gravado em memória e transmitidoa baixo volume para satélites quando o animal portador subir à superfície. No caso de que os PATspossam ser recuperados, os dados poderão ser obtidos na sua totalidade. Muitos dos tags electrónicosobtêm informação relativamente ao ambiente por meio de um conjunto de sensores. Os sensores maisusados medem a pressão (profundidade), a conductividade (salinidade) e os níveis de luminosidade. Estainformação ajuda os biólogos a melhor entender o ambiente em que vivem os animais. [5]

Tags comerciais

4

5

6

7

1.2.2 Receptores e estruturas de sensores

Depois de terem sido apresentadas as tags mais comuns na secção anterior, resta perceber quais osreceptores e as estruturas de sensores disponíveis para as detectar. Os sensores usados no meio aquáticopara detectar sinais acústicos são os hidrofones, que, dependendo da aplicação, se podem agregar paraestruturas de sensores maiores. Assim sendo, distingue-se aqui as seguintes estruturas que se diferenciamuma da outra principalmente pelo número e o distanciamento entre os sensores.

Ultra Short Baseline (USBL) apresenta um número reduzido de sensores com um posicionamentopróximo dos sensores. Uma estrutura USBL é constituída por pelo menos três sensores (normal-mente quatro por razões de redundância) montados de forma compacta com distâncias entre elesna ordem de centímetros ou de poucos metros. Assim é possível considerar que todos os receptoressejam afectados pelo meio aquático de forma semelhante. Para distâncias grandes em compara-ção ao distanciamento entre os sensores, é possível utilizar a simplificação da onda plana. Esta éusada para calcular uma estimativa de direcção de chegada do sinal a partir das Time DifferenceOf Arrival (TDOA), das diferenças nos tempos de chegada do sinal acústico aos vários receptores.No caso de se querer localizar um emissor, existe a possibilidade de que este envie uma estima-tiva de profundidade. Assim a intersecção entre a direcção de chegada do sinal do emissor como plano horizontal determinado pela estimativa da profundidade permite a localização do emissor.No entanto, uma estimativa de profundidade é geralmente baseada na medida da pressão e é poucoprecisa. Outro método mais robusto consiste em utilizar duas estruturas USBL e efectuar a inter-secção das rectas formadas pelas duas estimativas de direcção de chegada do sinal de cada uma dasestruturas USBL.A montagem e agregação de poucos sensores de forma compacta torna o sistema versátil e facilitaa sua utilização e instalação. Resulta daí uma complexidade reduzida na operação de um sistemacom configurações USBL, pois não existe a necessidade de efectuar instalações no fundo do mar.Desta forma uma estrutura USBL adequa-se a missões breves em que se pretenda também acedera dados recolhidos de forma online.Por outro lado uma estrutura USBL exige uma precisão adicional na estimação dos tempos dechegada do sinal aos receptores e é necessário conhecer-se a posição de cada estrutura USBL, bemcomo a sua atitude no espaço, o que normalmente exige a utilização de um sistema de navegaçãoinercial em conjunto com um sistema de navegação terrestre.

Underwater Acoustic Sensor Network (UW-ASN) consiste numa rede de sensores que se podeestender sobre áreas extensas fazendo uso de um número ilimitado de sensores individuais. Apesarde haver esforços de criar protocolos de comunicação entre os diferentes nós de uma rede de sensores[6], a abordagem tradicional para o funcionamento de uma rede de sensores consiste na instalaçãode receptores isolados em localizações fixas, onde recolhem dados que apenas se tornam disponíveisapós a recuperação do equipamento.Uma UW-ASN tradicional adequa-se a estudos prolongados com uma duração na ordem de váriosmeses em que os sensores recolhem dados de forma autónoma sendo que normalmente não é neces-sária qualquer manutenção. Isto significa também que a rede está sujeita às condições ambientais

8

sem que haja a possibilidade de reajustamento ou intervenção. Correntes na água ou outras in-fluencias externas poderão causar uma desconfiguração da rede que por vezes apenas poderá serdetectada durante a recuperação do equipamento ao final de uma missão. Em caso extremo istopoderá significar que a missão não concluirá os seus objectivos.Apesar da simplicidade do equipamento individual que constitui os nós da rede, o processo de ins-talação dos mesmos, bem como a análise dos dados recolhidos implicam uma elevada complexidadena operação de uma UW-ASN.

Sistema de localização comercial

9

10

1.3 Acústica subaquática

A acústica subaquática desempenha um papel fundamental no funcionamento de sistemas de telemetriano meio aquático. As ondas acústicas são originadas através da propagação de uma perturbação mecânicalocal que se afasta da sua fonte.A onda acústica é caracterizada pela amplitude de vibrações locais de cada partícula no meio de propa-gação à volta da sua posição de equilíbrio e pela pressão resultante dessas compressões e dilatações locaisdo meio. A variação da pressão local em relação à pressão média hidrostática, a pressão acústica, é aquantidade mais frequentemente usada na acústica subaquática e é expressa em micro Pascal [mPa]. Ascaracterísticas do meio de propagação vão definir a velocidade local de propagação. O módulo de elasti-cidade E do meio quantifica a variação relativa do volume ou da densidade do meio devido às variaçõesde pressão. A densidade da água ⇢ tem em média um valor aproximado de 1,030 kgm�3, dependendodos parâmetros da pressão, da salinidade e da temperatura. Assim a velocidade de propagação c de umaonda acústica atinge tipicamente um valor entre 1450 m/s e 1550 m/s e relaciona-se com a densidade eo módulo de elasticidade do seguinte modo

c =

sE

⇢(1.1)

As ondas de som podem ser caracterizadas pela sua frequência f (o número elementar de vibrações porsegundo [Hz]) ou pelo seu período T (duração de um ciclo de vibração elementar) sendo T = 1/f . Ocomprimento de onda � é a correspondente espacial à periodicidade temporal, ou seja, é a distânciapercorrida pela onda durante um período do sinal com a velocidade c. Logo

� = cT =

c

f(1.2)

Para descrever a propagação de uma onda acústica utilizam-se normalmente ondas planas e esféricas.As ondas planas são usadas quando o valor da amplitude se aproxima de uma constante e as frentes daonda apresentam uma curvatura negligível. Estas condições verificam-se a grandes distâncias da fonte.Para modelar processos locais próximos de uma fonte como um ponto e quando a diminuição do valor daamplitude durante a propagação não pode ser ignorada, usam-se as ondas esféricas.

Figura 1.3: Propagação de som como onda esférica próximo da fonte e propagação como onda plana longeda fonte de origem [7]

A propagação de uma onda sonora está associada à energia acústica. A intensidade acústica I é o valor

11

médio do fluxo de energia por unidade de superfície e de tempo. Para uma onda plana de amplitude p0e valor de Root Mean Square (RMS) prms = p0/

p2, a intensidade é dada por:

I =

p202⇢c

=

p2rms

⇢c(W/m2) (1.3)

A potência acústica P é a intensidade acústica recebida por uma superfície ⌃. No caso de uma ondaplana seria:

P = I ⇥ ⌃ =

p20⌃

2⇢c=

p2rms⌃

⇢c(W) (1.4)

Geralmente a pressão ou a intensidade de uma onda acústica é descrita em decibeis. Na acústica su-baquática a referência de pressão é o micro Pascal. Assim o nível de pressão absoluto é

pdB = 20 log

✓p

pref

◆(1.5)

e é expresso em dB relativamente a 1 mPa (dB re 1 mPa). É claro que p e pref devem estar definidosda mesma forma, ou seja ambos devem estar definidos em valores RMS ou ambos em valores de pico apico.

Perdas de propagação

A perda de intensidade é o efeito mais evidente na propagação de ondas acústicas devida à divergênciageométrica na propagação e à absorção da energia acústica pelo próprio meio de propagação.A propagação de uma onda acústica a partir de uma fonte faz divergir a energia acústica transmitida.Dado que a energia é conservada, a sua intensidade diminuirá inversamente proporcional ao aumento dasuperfície.No caso de uma fonte com radiação de forma esférica e envolvida por um meio homogéneo e infinito,a relação entre as intensidades em dois pontos diferentes é inversamente proporcional à relação dassuperfícies das esferas em que se encontram os pontos

I2I1

=

⌃1

⌃2=

4⇡r214⇡r22

=

✓r1r2

◆2

(1.6)

em que r1 e r2 são as distâncias radiais a partir da fonte. Logo, a intensidade diminui em 1/r2 e a pressãoem 1/r. As perdas devidas à divergência durante a transmissão TL consideradas a partir da distância dereferência (r1m = 1 m) podem ser expressas em dB:

TL = 20 log

✓r

r1m

◆(1.7)

A perda por divergência esférica corrigida para considerar perdas por atenuação do meio é geralmenteusada como primeira aproximação quando se avalia a perda de propagação e o desempenho de sistemassubaquáticos em condições em que não se preveja propagação por múltiplos trajectos ou refracção devidoa grandes variações no meio aquático.A perda de transmissão em dB é

TL = 20 log

✓r

r1m

◆+ ↵r (1.8)

Tem de se tomar atenção às unidades utilizadas na equação 1.8 em que r está expresso em metros enquantoque ↵ normalmente é dado em dB/km. Por esta razão as unidades deverão ser devidamente convertidas.Para frequências superiores a 1kHz este modelo geralmente não será adequado e ter-se-á de recorrer amodelos mais complexos [8].

12

No caso de águas pouco profundas não se verifica uma completa divergência esférica devido aos limitesconstituídos pela superfície da água e pelo fundo do mar que reflectem as ondas acústicas. A energiaacústica permanece entre esses limites propagando-se em múltiplos trajectos. Se a frequência da onda forsuficientemente alta, as oscilações no campo resultante podem ser desprezadas, considerando-se apenaso fluxo médio de energia que pode ser modelado como uma função da distância cilíndrica r.Como está representado na figura 1.4, a propagação é inicialmente esférica. Quando os limites de superfíciee do fundo do mar forem atingidos pela onda acústica a uma distância r0, a sua propagação pode sermodelada de forma cilíndrica. A sua intensidade diminui então em 10 log(r/r0) (em vez de 20 log(r/r0)com propagação esférica).

⇢TL = 20 log r + ↵r r < r0TL = 20 log r0 + 10 log

rr0

+ ↵r = 10 log(rr0) + ↵r r > r0(1.9)

Figura 1.4: Propagação de som de forma cilíndrica em água com pouca profundidade para distânciassuperiores a r0 [7]

A água do mar é um meio de propagação que absorve uma parte da energia de onda transmitida aqual é dissipada através de reacções químicas ou da sua viscosidade. A redução da amplitude local éproporcional à sua própria amplitude. Deste modo a pressão acústica diminui de forma exponencial coma distância. Estas perdas são adicionais às perdas por divergência.A diminuição exponencial da pressão dá uma perda em dB proporcional à distância de propagação o quepode ser expresso através de um coeficiente de atenuação em dB/m (frequentemente em dB/km)O grau de absorção depende fortemente do meio de propagação e da frequência. Na água do mar aabsorção deriva de

- viscosidade da água pura, cujo efeito aumenta com o quadrado da frequência;

- relaxação de moléculas de sulfato de magnésio (MgSO4) para frequências inferiores a 100 kHz;

- relaxação de moléculas de ácido bórico (B(OH)3) para frequências inferiores a 1 kHz.

As alterações de pressão local desfazem temporariamente as moléculas de sulfato de magnésio e de ácidobórico. Se o período da onda acústica for superior ao tempo necessário para a molécula se recompor, oprocesso repete-se em cada ciclo e está constantemente a dissipar energia. Assim, a atenuação devida aeste processo aparece em frequências mais baixas que a frequência de relaxação característica relevante.Tendo em consideração estes três factores pode-se modelar o factor de atenuação devido à absorção com

↵ = C1f1f2

f21 + f2

+ C2f2f2

f22 + f2

+ C3f2 (1.10)

13

em que os primeiros dois termos apresentam as contribuições resultantes dos processos de relaxação. Oúltimo termo corresponde à viscosidade apenas da água. Os coeficientes Ci dependem da temperatura,da pressão hidrostática e da salinidade. O modelo mais usado para o cálculo desses coeficientes é omodelo de Francois e Garrison. Para um estudo mais detalhado sobre este tema, consulte-se a literaturaidentificada na bibliografia no ponto [7] ou o anexo A.

Impacto do fundo e da superfície do mar na propagação de som

As características do fundo do mar condicionam a propagação de som. Assim uma onda acústica poderáser reflectida ou a sua energia dissipada conforme o contraste entre a impedância acústica da água e a dofundo do mar. Um fundo rochoso terá um coeficiente de reflexão mais elevado do que um fundo porosoou arenoso. Dependendo da frequência o perfil do fundo do mar poderá ser mais ou menos irregular emrelação à onda acústica o que tem impacto na criação de mais ou menos ecos.A superfície do mar apresenta uma interface reflectora devido ao grande contraste entre as impedânciasacústicas de água e de ar. A irregularidade da superfície relativa à onda acústica depende, como o perfildo fundo do mar, da frequência do sinal. No entanto, o seu perfil não é estático mas varia com o tempodependendo da ondulação existente.A ondulação cria também bolhas de ar próximo à superfície (geralmente o efeito de bolhas de ar poderáser desprezado em profundidades superiores de 10 a 20 m [7]) que criam um meio não homogéneo comgrande variabilidade local das características acústicas no meio de propagação. Possíveis efeitos são:

- atenuação na amplitude do sinal acústico devido a uma camada de absorção adicional;

- diferenças locais na velocidade de propagação de som dão origem a uma refracção na camadapróxima à superfície;

- criação de pequenos ecos parasitas que se sobrepõem ao sinal acústico e são criados devido à grandediferença de impedância acústica entre as bolhas de ar e o meio aquático envolvente;

Dependendo do tamanho e da densidade de bolhas de ar estes efeitos podem influenciar o desempenhode um sistema acústico de forma drástica ou mesmo inviabilizar o seu funcionamento por completo.

Propagação multi-trajecto

Dado que o meio de propagação está limitado pela superfície e pelo fundo do mar, os sinais, enquantose propagam, sofrem sucessivas reflexões nestes interfaces. Devido a estas reflexões um dado sinal podepropagar-se desde a fonte até ao receptor por trajectos distintos que correspondem a diferentes direcçõese durações de propagação. Para além do sinal principal que se propaga em linha recta, geram-se ecoscujas amplitudes diminuem com o número de reflexões sofridas.A formação de ecos pode ter um impacto no desempenho de um sistema acústico e o número de ecos quepodem interferir no funcionamento do sistema acústico é muito variável dependendo do atraso temporale da amplitude dos ecos. Na figura 1.5 está esquematizada uma propagação por múltiplos trajectos eapresenta-se o sinal recebido no receptor no domínio temporal. Cada eco incide sobre o receptor a partirde uma direcção diferente. No receptor o eco aparenta vir de uma imagem da fonte que se encontra noprolongamento do trajecto da última reflexão.

14

Figura 1.5: Encima: propagação por múltiplos trajectos em águas pouco profundas. (A) trajecto directo;(B) reflexão na superfície; (C ) reflexão no fundo do mar; (D) reflexão na superfície e no fundo do mar;(E) reflexão no fundo do mar e na superfície.Embaixo: Envolvente do sinal resultante do trajecto directo e de ecos causados por reflexões em domíniotemporal. As letras indicam a correspondência entre os picos de sinal visíveis e o trajecto correspondentepercorrido [7]

Refracção devido ao perfil de velocidade de som com a profundidade

O meio aquático não é homogéneo, mas encontra-se geralmente subdividido em várias camadas. Paraalem disso o aumento continuo da pressão com a profundidade, tem como consequência um aumentotambém da velocidade de som.Olhando para a interface entre duas camadas com velocidades de som c1 e c2, representada na figura 1.6surge o fenómeno de refracção. A onda acústica ao passar a fronteira entre as duas camadas altera asua direcção de propagação conforme descrito pela equação de Snell-Decartes, em que �1 é o angulo deincidência com a interface e �2 o ângulo de saída:

cos�1c1

=

cos�2c2

. (1.11)

Figura 1.6: Refracção da onda acústica na interface de dois meios com densidades ⇢i e velocidades desom ci diferentes (c2 > c1) [7]

15

A equação 1.11 apenas se aplica quando cos�2 6 1. Para ângulos rasos com valor menor que o angulocrítico dado por

�c = arccos

✓c1c2

◆, (1.12)

a onda acústica não passa da camada 1 para a camada 2, mas é reflectida na sua totalidade.Para um problema de N camadas de velocidade de som constante, identificadas pelos índices i = 1, 2, ...Ncomo representado na figura 1.7, pode-se escrever a lei de Snell-Decartes da seguinte forma:

cos�1c1

= · · · = cos�ici

=

cos�i+1

ci+1= · · · = cos�N

cN(1.13)

Generalizando para camadas infinitesimais num meio aquático em que a velocidade de som varia com aprofundidade, a expressão pode ser escrita como:

cos�ici

= constante (1.14)

Figura 1.7: Refracção devido à profundidade com N camadas (esquerda) e com aumento contínuo davelocidade de som (direita) [7]

Num meio em que a velocidade de som aumenta continuamente com a profundidade, o angulo que umaonda acústica, vinda da direcção da superfície, forma com a horizontal, vai diminuindo. No momento emque o ângulo crítico �1 seja atingido, a onda será reflectida para cima. Passando as camadas de águaagora na direcção oposta, o ângulo com a horizontal vai aumentando, conforme visível na figura 1.7.

Canais acústicos

A velocidade de propagação do som depende da temperatura, da pressão e da salinidade da água.A temperatura da água normalmente diminui da superfície para o fundo do mar, mas em camadas maispróximas da superfície factores como a mistura das camadas devido à ondulação, o aquecimento solar eoutros factores externos podem causar grandes variações no perfil da temperatura. Em profundidadesmais elevadas (em oceanos cerca de 1000 m e no caso de mares fechados como no mar mediterrânico cercade 100 a 200 m) a temperatura média permanece estável e diminui lentamente com a profundidade. [7]A pressão hidrostática dependendo da profundidade causa um aumento da velocidade do som devido avariações do módulo de elasticidade. Numa primeira aproximação o aumento da velocidade de som élinear e corresponde aproximadamente a 0.017 m/s [7].A percentagem da massa de sal dissolvida na água do mar é geralmente expressa em Practical SalinityUnits (p.s.u.). A salinidade nos grandes oceanos tem em média um valor de 35 p.s.u. mas pode localmenteter fortes variações. Em mares mais pequenos o valor médio de salinidade depende das condições locais

16

como a evaporação ou da entrada de água doce.A salinidade normalmente tem uma variação muitopequena em relação à profundidade (1 a 2 p.s.u.), excepto em águas doces menos profundas. [7]

Figura 1.8: Perfil típico da velocidade de som dependendo da profundidade [7]

A figura 1.8 mostra um perfil vertical idealizado da temperatura e da velocidade de som nos oceanos,típico para o canal Sound Fixing And Ranging (SOFAR). Este diagrama simplificado pode ser divididoem quatro camadas distintas.Na camada homogénea de superfície (Surface Mixed Layer), existe uma camada isotérmica e homogéneacom velocidade de som constante devido à mistura contínua causada pela ondulação.A camada de superfície (Surface Layer), que se pode estender até 200 m de profundidade, também elaisotérmica, a velocidade do som aumenta com a profundidade devido ao aumento da pressão hidrostá-tica. Esta camada é formada em condições de inverno, por introdução de água doce por rios ou peloderretimento de gelo à superfície.Na camada de termoclínio permanente a temperatura diminui com a profundidade e consequentementetambém diminui a velocidade de propagação de som. A base do termoclínio muda muito dependendo dalatitude mas é tipicamente encontrada a uma profundidade de 1000m. Esta camada pode ser permanenteou variar sasonalmente.Dentro da camada mais profunda, debaixo do termoclínio permanente, a mudança da temperatura émenor. Aqui, com o aumento da profundidade, verifica-se também um aumento contínuo da velocidadede som devido à subida de pressão. (tal como no Surface Layer).Na figura 1.9 encontram-se vários perfis típicos da velocidade de som.

17

Figura 1.9: Perfis genéricos da velocidade do som. A velocidade de som está representada no eixohorizontal, a profundidade no eixo vertical. Da esquerda para a direita estão representados os canais: (A)SOFAR de Verão; (B) SOFAR de Inverno; (C) mediterrânico de Inverno (isotermal); (D) mediterrânicode Verão; (E) atlântico Nordeste (com intrusão de águas do mar mediterrânico); (F) polar; (G) águapouco profunda no Inverno; (H) água pouco profunda no Verão; (I) água pouco profunda no Outono; (J)água pouco profunda com água doce à superfície. [7]

Devido à grande variabilidade da velocidade de propagação do som dependendo da temperatura, dapressão e da salinidade da água, existem diferenças espaciais e temporais consideráveis na velocidade depropagação do som. Estas poderão causar a formação de canais acústicos, como o Surface Channel e oDeep Channel assinalados na figura 1.8, que correspondem a um mínimo da velocidade de propagaçãode som. Dependendo da espessura do canal estes canais podem ser selectivos em relação à frequência. Aenergia acústica é focada e pode ser transmitida a longas distâncias com atenuação inferior à propagaçãode forma esférica. Por outro lado, tal como é visível na figura 1.10 , a focalização da energia acústicapode também dar origem a zonas de sombra que não sejam atingidas pela onda acústica.

Figura 1.10: Criação de zonas de sombra pela variabilidade da velocidade do som na água [8]

18

Capítulo 2

Detecção e Descodificação das tags

comerciais

Um dos objectivos deste trabalho é a detecção dos sinais emitidos pelas tags comerciais e a determinaçãoda direcção de origem do sinal. Como tal, é importante primeiramente analisar ao pormenor os sinaispor elas emitidos, para se poder elaborar uma estratégia de detecção.

2.0.1 Sinal emitido

19

20

21

22

2.0.2 Estratégias de detecção

A detecção do sinal emitido pelas tags comerciais requer, como previamente identificado, a verificação daexistência de um ping e a determinação do instante temporal em que surgiu relativamente aos restantes.Como o sinal no receptor foi corrompido pelos efeitos do canal acústico, sofrendo distorção e uma sobre-posição de ruído, convém a utilização de um filtro na sua recepção por forma a melhorar a relação sinalruído e assim melhorar as condições para o algoritmo de detecção.

Figura 2.1: Diagrama de blocos do filtro

Na figura 2.1, g(t) representa o sinal emitido após ter sido corrompido pelos efeitos do canal. Numambiente ruidoso assume-se n(t) como ruído gaussiano branco aditivo. À saída do filtro h(t) o sinal y(t)deverá mostrar melhores características ou seja estar mais próximo do sinal original s(t) que x(t) parapossibilitar uma boa detecção.Um impulso emitido como s(t) pode ser modelado de uma forma ideal por uma sinusóide infinita mul-tiplicada por uma função rectangular, ou impulso unitário, neste caso com uma duração de 10ms. Nodomínio da frequência esta operação reflecte-se numa convolução entre uma risca espectral centrada nafrequência da sinusóide e uma função Sinc que tem os seus zeros em múltiplos de 1

T , sendo T a largura doimpulso unitário, ver figura 2.2. A largura espectral do sinal esperado é de 2

T igual a 200Hz, equivalenteao lobo principal da função Sinc.

23

Figura 2.2: A transformada de Fourier de um impulso unitário no tempo corresponde à função Sinc nodomínio da frequência [9]

Como filtro foram considerados dois tipos de filtros diferentes:

- Utilização de um filtro adaptado

- Utilização de um filtro passa-banda

Em seguida estão descritas ambas as abordagens mais em detalhe e mencionadas algumas característicasque justificam a escolha tomada.

Utilização de um filtro adaptado

Um filtro adaptado é aquele em que a sua resposta a um sinal de entrada tem uma relação sinal-ruídomáxima [10]. Sendo o sinal de entrada x(t), t 2 [0;T ], corrompido por ruído branco gaussiano n(t), podeser aplicado um filtro adaptado, se a forma do sinal desejado s(t) for conhecida. O filtro adoptado édefinido por hadp na equação 2.1 e é uma versão do sinal desejado invertida no tempo, atrasada por T eescalada pela constante c.

hadp(t) = c · s(T � t) (2.1)

A saída do filtro adaptado está relacionada com a correlação entre o sinal esperado e o sinal de entrada.Desta forma a saída vai acusar um máximo no caso em que o sinal de entrada seja parecido ao sinaldesejado.O sinal esperado encontra-se guardado em memória num bloco de N elementos. Para poder utilizar umfiltro adaptado a implementação tem que seguir um processamento discreto por blocos. A saída do filtroadaptado y(m) em tempo discreto encontra-se descrito na seguinte equação [9].

hadp(m) = s(N � 1�m) , 0 m N � 1

y(m) =

N�1X

k=0

hadp(m� k)x(m)

(2.2)

24

Utilização de um filtro passa-banda

Dimensionamento de um filtro passa banda Um filtro digital pode ser implementado como umfiltro de resposta impulsional finita (FIR) ou infinita (IIR). Cada tipo de implementação tem as suas van-tagens e desvantagens, comparar com a tabela 2.1. Como o objectivo é apenas detectar a presença de umsinal, as implicações do filtro na fase não têm grande relevância face às vantagens de uma implementaçãomais eficiente no que diz respeito às operações de adição e multiplicação. A forma de implementaçãomais vantajosa do filtro é portanto um IIR cujo dimensionamento parte de uma especificação analógica.

Tabela 2.1: Filtro FIR vs. filtro IIR [11]

FIR IIR

não recursivo recursivosempre estável pode ser instável

fácil obtenção de fase linear fase difícil de controlarmais operações necessárias menos operações necessárias

25

O filtro passa-banda deve ser adaptado ao sinal desejado ajustando as frequências de corte inferior esuperior à largura de banda do sinal. Neste caso o espectro do sinal é apenas um pico estreito (200Hz).Um filtro adequado é um, em que ambas as frequências de corte coincidam, para que as frequências decorte inferior e superior sejam ambas iguais à frequência central do sinal de 69kHz. O filtro pode serdescrito pela função de transferência

H(s) =s

s2 + !0Q s+ !2

0

. (2.3)

No entanto pretende-se o processamento de forma digital. Aplicando a transformação bilinear 2.4consegue-se uma representação da função de transferência no domínio digital.

s =2

Ts

z � 1

z + 1

(2.4)

A transformação bilinear contém aproximações e por isso é necessário fazer um ajuste à frequência. Afrequência analógica !a ajustada deve ser usada para a parametrização do filtro no lugar de !0 para queo filtro se comporte da forma pretendida no domínio da frequência digital !d. A relação entre frequênciaanalógica e digital [11] é dada por:

j!a=2(ej!dTs � 1)

Ts(ej!dTs+ 1)

(2.5)

!a=2

Tstan(

!dTs

2

) (2.6)

(2.7)

Depois de aplicada a transformada bilinear, a função de transferência encontra-se no domínio da trans-formada Z.

H(z) =2Ts

z�1z+1

(

2Ts

z�1z+1 )

2+

!0Q

2Ts

z�1z+1 + !2

0

(2.8)

Y (z)

X(z)=

2(z � 1)Ts(z + 1)

4(z � 1)

2+

!0Q 2(z � 1)Ts(z + 1) + !2

0T2s (z + 1)

2(2.9)

Y (z){4z2 � 8z + 4 + 2

!0

QTs(z

2 � 1) + !20T

2s (z

2+ 2z + 1)}

= X(z){z2(sTs)� 1}(2.10)

Aplicando a transformada Z inversa, as equações são levadas ao domínio do tempo.

y[n] =1

Gy[n� 2]{4 + !0

Q2Ts + !2

0T2s }(�1)

� y[n� 1]{2!20T

2s � 8}+ x[n� 2]{2Ts}� x[n]

G = [4� 2

!0

QTs + !2

0T2s ]

(2.11)

As variáveis que ajustam o filtro às necessidades estão indicadas em 2.12 em que Q é o factor de qualidadeque vai definir o quão abrupto vai ser o descaimento da função de transferência nas frequências de cortee fs representa a frequência de amostragem.

!d=2⇡fd

fd=69kHz

Q=10

Ts=1

fs(2.12)

26

O diagrama de amplitude em 2.3a mostra um pico à volta da frequência desejada e uma atenuação comdescaimento rápido para ambos os lados. A largura de banda a 3dB é de aproximadamente 490Hz,superiores aos 200Hz previamente identificados no modelo do sinal. O diagrama de fase 2.3b apresentacaracterísticas não lineares como é de esperar numa implementação IIR.

(a) Resposta do filtro em amplitude com largura de

banda a 3dB marcada a vermelho (b) Diagrama de fase do filtro

Figura 2.3: Diagramas de Bode do filtro IIR

A implementação do filtro IIR pode ser visualizada no diagrama de blocos 2.4. O bloco Z�1 significaum atraso por uma amostra. Os vários factores a e b devem ser escolhidos de forma a que se ajustem àequação 2.11. Com esta arquitectura é possível fazer um processamento amostra a amostra o que tornaa implementação simples.

Figura 2.4: Diagrama de blocos da implementação de um IIR de segunda ordem em estrutura biquad [12]

27

2.0.3 Determinação do Time Difference of Arrival (TDOA)

28

29

30

2.0.4 Solução adoptada

31

32

2.0.5 Descodificação do ID e Checksum

33

34

2.0.6 Estimação da direcção de chegada do sinal

O problema da estimação da origem de um sinal incidente num array de receptores tem sido alvo de estudodurante várias décadas e continua a ser objecto de investigação activa. Existem variadas abordagens comoDevide and Conquer ou baseadas em Series de Taylor, no entanto exigem a existência de informaçãoprévia ou são algoritmos de natureza iterativa [13]. Um calculo baseado apenas em TDOA geralmenterecorre à intersecção de um conjunto de hyperboloides, o que pode exigir um processamento pesado [14].No entanto, um trabalho muito citado, chamado Spherical Interpolation method [15] recorre apenas aTDOAs, mas convence por ter uma solução de forma fechada com baixa complexidade computacionalque o torna um método adequado para o processamento em tempo real. Para uma estrutura USBLporem, mesmo que seja teoricamente possível determinar a posição do emissor, o desempenho em relaçãoà determinação da distância costuma ser fraco devido à baixa relação sinal-ruído normalmente presentee às dimensões reduzidas da estrutura de receptores em comparação à distância [16]. Uma abordagempassiva pelo método dos mínimos quadrados tomada em [17] e em [16] resolve o problema de estimação daDirection-of-arrival (DOA), ou seja da direcção de chegada, de um sinal a partir das TDOAs referentesaos elementos de um array de recepção.

Seguindo o raciocínio em [16] e por motivos de simplicidade assume-se um meio com velocidade depropagação constante. Pode também ser aplicada a aproximação de onda plana. Considerando umalonga distancia entre o emissor e o receptor, em comparação às dimensões da estrutura USBL, é razoávelmodelar a frente esférica do sinal acústico incidente no array de receptores por uma onda plana, comodescrito na secção 1.3. Um número de receptores igual ao número de dimensões do espaço é o suficientepara definir a posição do emissor (abordagens com menos sensores existem, embora requeiram informaçãoadicional [18]). No entanto, para minimizar os efeitos dos erros de medida na presença de ruído, recorre-sea mais receptores, introduzindo assim redundância. Na figura 2.5 é delineada uma onda plana incidenteno plano XY da estrutura USBL em dois dos N receptores (i e j) nos instantes ti e tj . O vector unitáriod = [dx dy dz]T aponta na direcção do emissor e deve ser estimado.

35

Figura 2.5: Incidência da onda acústica em dois receptores projectados no plano XY [16]

Sejam as posições dos receptores i e j dadas por ri = [xi yi zi]T e por rj = [xj yj zj ]T e a velocidade depropagação do som no meio por vp. Pode-se escrever

vp(ti � tj) = �dT (ri � rj) (2.13)

ou seja, a distância percorrida pela frente de onda entre os instantes ti e tj é igual à componente dovector (ri � rj) coincidente com a direcção de propagação da onda �d. Sem recorrer à notação vectorial,a eq. 2.13 escreve-se

vp(ti � tj) = �(dx(xi � xj) + dy(yi � yj) + dz(zi � zj)) (2.14)

sendo as incógnitas dx, dy e dz, e tendo-se escolhido representar a direcção do emissor por um vectorunitário tem-se ainda

kdk = 1 ,q

d2x + d2y + d2z = 1 (2.15)

Para calcular a direcção do emissor é necessária uma terceira condição para além daquelas expressas pelaseq. 2.14 e 2.15. Do ponto de vista físico isso significa que é necessário mais um receptor que permita mediroutro tempo de chegada e obter uma condição semelhante a 2.14. Três será portanto o número mínimode receptores necessários para calcular a direcção do emissor. No entanto, para melhorar o desempenhodo sistema é aconselhável a utilização de uma estrutura USBL com mais receptores. Considere-se entãoa generalização do problema para N receptores, considere-se a eq. 2.14 com

{i = 1, 2, ...N ; j = 1, 2, ...N ; i 6= j} (2.16)

Sejam �1 = t1 � t2, �2 = t2 � t3, ... �M = tN�1 � tN todas as combinações possíveis para as diferençasdo tempo de chegada do sinal aos N receptores, sendo portanto M =

N2 C. Se definirmos o vector �

como

� = [�1 �2 ... �M ]

T , (2.17)

que pode ser gerado recorrendo a uma matriz de combinação C 2 RM⇥N e ao vector dos tempos dechegada tm, tem-se:

� = Ctm. (2.18)

36

De forma idêntica para as posições dos N receptores define-se os vectores

x = [x1 � x2 x2 � x3 . . . xN�1 � xN ]

T , (2.19)y = [y1 � y2 y2 � y3 . . . yN�1 � yN ]

T , (2.20)z = [z1 � z2 z2 � z3 . . . zN�1 � zN ]

T . (2.21)

Então a generalização do problema pode escrever-se, com base na eq. 2.14, como

vp� = �(dxx+ dyy + dzz). (2.22)

O objectivo é encontrar a direcção do emissor, d = [dx dy dz]T , que minimiza o erro quadrático totaldado por

J =

MX

k=1

(dxxk + dyyk + vp�k)2. (2.23)

Tomando S = [x y z], a eq. 2.22 pode também ser escrita como

vp� = �Sd (2.24)

A solução dos mínimos quadrados é obtida por aplicação da pseudo inversa S# de S:

�ST vp = (STS)d� (2.25)

d = �(STS)�1ST vp� (2.26)

d = �vpS#Ctm (2.27)

com S#= (STS)�1ST .

O vector d é portanto a estimativa usada para descrever a direcção de chegada do sinal aos receptores doUSBL.

37

Capítulo 3

Filtro de Kalman

O Filtro de Kalman, inventado por R.E. Kalman e publicado em 1960, serve para estimar o estado de umsistema de forma óptima quando as medições são corrompidas por ruído. Existem muitas variantes dofiltro de Kalman. As diferentes versões podem ser associadas às categorias de filtros para uso em tempocontínuo ou discreto e filtros para problemas lineares ou não-lineares. A lista das aplicações do Filtrode Kalman é extensa e engloba vários ramos industriais, sendo um dos mais importantes, a indústria danavegação.A utilidade do filtro baseia-se na estimação óptima do estado do sistema com a utilização de medidas:

- Incompletas: relacionadas com algumas mas não com todas as variáveis de interesse

- Indirectas: relacionadas apenas indirectamente com as quantidades de interesse

- Intermitentes: disponíveis apenas em instantes espaçados no tempo de forma irregular

- Inexactas: corrompidas por várias fontes de erro

A ideia central do filtro de Kalman consiste em modelar o sistema como um sistema linear dinâmicoque é corrigido por medidas de sensores. Tanto o modelo, como também os sensores são sujeitos a ruídocausando incerteza sobre o estado real do sistema. Com o conhecimento da natureza do ruído e dassuas estatísticas, é possível obter uma estimativa óptima em sentido estatístico para o estado do sistema,mesmo no caso de sensores pouco exactos. Deve ser possível modelar o ruído como sendo ruído aditivo,branco e gaussiano (Additive White Gaussian Noise (AWGN)) e as medições são assumidas mutualmenteindependentes e também independentes da dinâmica do sistema.Um modelo linear é aplicável para descrever o comportamento de uma grande variedade de sistemas.Quando existem não-linearidades há a possibilidade de linearizar ou de usar modelos havendo assimdiversas variantes e modificações do filtro de Kalman (Extended Kalman Filter (EKF), Unscented KalmanFilter (UKF), etc.).

3.1 Exemplo de funcionamento do Filtro de Kalman

Para melhor entendimento do funcionamento básico do Filtro de Kalman serve o seguinte exemplo simplesque segue o exemplo em [1]. Pretende-se estimar uma posição unidimensional. Em dado instante t1 édeterminada a posição z1. No entanto, por falta de precisão do método de medição, o valor da medidatem associada alguma incerteza representada pelo desvio padrão �z1 ou pela variância �2

z1 . Assim sendo,pode ser estabelecida uma função de densidade de probabilidade que descreve a probabilidade de x(t1),da posição no instante t1, condicionada ao facto, de o valor medido ter sido z1. Esta função encontra-seilustrada na figura 3.1.

38

Figura 3.1: Função de densidade de probabilidade para uma estimativa da posição x, dada a medida z1e conhecendo o desvio padrão da medida �z1 [1]

A figura 3.1 mostra a curva fx(t1)|z(t1)(z|z1), uma função da posição x: Ela indica a probabilidade daposição real, dada a medida z1. Note-se que �z1 é uma medida directa da incerteza existente. Quantomaior for a incerteza, mais largo vai ser o pico da função de probabilidade estendendo-se por um maiornúmero de valores de x.Com base nesta função de densidade de probabilidade e sem mais informação, a melhor estimativa paraa posição verdadeira é

x̂(t1) = z1 (3.1)

e a variância do erro na estimativa é

�2x(t1) = �2

z1 (3.2)

Sem que a posição verdadeira se tenha alterada é tirada uma segunda medida z2 no instante t2 com umsensor melhor ou seja com maior precisão. Consequentemente, o desvio padrão �z2 também é menor,como ilustrado na figura 3.2. A figura mostra a curva que descreve a probabilidade da posição verdadeirabaseada apenas no valor da medida em instante t2 e da sua incerteza. Repare-se que esta tem uma formamais estreita em comparação á curva da medição anterior, dada a menor incerteza. Neste momentoexistem duas medidas para estimar a posição real, ambas com incerteza associada. Objectivo é então,combinar as duas medidas de forma a obter uma nova estimativa sobre a posição real com maior certezado que as duas medidas por si só.

39

Figura 3.2: Função de densidade de probabilidade para uma estimativa x, dada a medida z2 em compa-ração à função de densidade de probabilidade no instante t1 representada a tracejado [1]

Figura 3.3: Função de densidade de probabilidade que descreve a estimativa x com valor esperado µ edesvio padrão �, resultante da combinação das medidas z1 e z2 representadas a tracejado [1]

40

Se as funções de densidade de probabilidade forem descritas como distribuições normais, então a distri-buição de probabilidade condicionada por ambas medidas, z1 e z2, resulta também numa distribuiçãonormal mas com valor esperado µ e variância �2:

µ =

�2z2

�2z1 + �2

z2

z1 +�2z1

�2z1 + �2

z2

z2 (3.3)

1

�2=

1

�2z1

+

1

�2z2

(3.4)

Repare-se que � é menor que �z1 ou �z2 . Condicionado ao histórico das medidas z1 e z2, a melhorestimativa possível no instante t2 é

x̂(t2) = µ (3.5)

com variância associada �2.Olhando para a equação 3.3 de µ, é fácil entender que no caso de �z1 e �z2 terem o mesmo valor, ou seja, asduas medidas terem a mesma incerteza, a estimativa óptima para a posição apenas ser a média das duas.Por outro lado, se �z1 for maior que �z2 , havendo então mais confiança na medida z2, então a equação 3.3dá mais peso a z2 que a z1. Finalmente a variância da nova estimativa é menor que �z1 , mesmo que �z2fosse muito grande. Mesmo informação de pouca qualidade melhora a precisão da estimativa. A equaçãopara x̂(t2) pode ser rescrita como

x̂(t2) =�2z2

�2z1 + �2

z2

z1 +�2z1

�2z1 + �2

z2

z2

= z1 +�2z1

�2z1 + �2

z2

[z2 � z1]

(3.6)

Posto de outra forma, mais próxima à notação usada para filtros de Kalman, sendo que x̂(t1) = z1:

x̂(t2) = x̂(t1) +K(t2)[z2 � x̂(t1)] (3.7)

K(t2) =�2z1

�2z1 + �2

z2

(3.8)

Estas equações indicam que a estimativa óptima x̂(t2) no instante t2, corrige a estimativa óptima anteriorx̂(t1), adicionando um termo de correcção com um factor de peso óptimo K(t2). O termo de correcção éa diferença entre a medida z2 e a sua melhor estimativa no momento anterior x̂(t1).Usando K(t2) da equação 3.8, a variância dada por equação 3.4 pode ser rescrita como

�2x(t2) = �2

x(t1)�K(t2)�2x(t1) (3.9)

Os valores de x̂(t2) e de �2x(t2) contêm toda a informação de fx(t2)|z(t1),z(t2)(x|z1, z2). Ou seja, com

estas duas variáveis, a probabilidade condicional é completamente descrita para o instante t2, dada ainformação de z1 e de z2.Para aumentar a utilidade deste exemplo, é acrescentado movimento. Supõe-se então, que a posiçãoactual se altera antes de ser feita uma nova medida. Alem disso, imagine-se que a melhor forma dedescrever esse movimento é através da forma simples

dx

dt= u+ w (3.10)

onde u é a velocidade nominal e w é um termo que descreve ruído que representa a incerteza sobreo conhecimento da velocidade real. Desprezando efeitos não lineares de primeira ordem, o ruído w é

41

modelado como ruído branco gaussiano com média nula e variância �2w. A figura 3.4 mostra o que

acontece à função de densidade de probabilidade da posição no instante t2, dadas as medidas z1 e z2. Aopassar o tempo, o pico da função de densidade desloca-se de acordo com o modelo simples dinâmico aolongo do eixo x. Ao mesmo tempo a incerteza aumenta, dando uma forma mais larga ao longo do tempopois existe menor certeza sobre a posição exacta devido à adição de incerteza através do ruído w.

Figura 3.4: Devido ao modelo dinámico do sistema, o pico da função de densidade de probabilidadedesloca-se em direcção do eixo x, mas com aumento da incerteza o desvio padrão da distribuição tambémaumenta [1]

No instante t�3 , mesmo antes da medida no instante t3, a função de densidade fx(t3)|z(t1),z(t2)(x|z1, z2),representada na figura 3.4, pode ser expressa matematicamente como uma distribuição gaussiana commédia e variância dadas por

x̂(t�3 ) = x̂(t2) + u[t3 � t2] (3.11)

�2x(t

�3 ) = �2

x(t2) + �2w[t3 � t2] (3.12)

Por isso, x(t�3 ) é a predição óptima do valor de x no instante t�3 , antes da medida no instante t3, e �2(t3)

é a variância esperada nessa predição. De seguida é tirada uma medida com valor z3, e a sua variânciaé assumida como �2

z3 . Como anteriormente, existem agora duas funções de densidade de probabilidadegaussianas com informação sobre a posição, Uma que resulta de todas as informações de medidas eestimativas anteriores, e a outra contendo informação da última medida. Através do mesmo processocomo anteriormente descrito, a função de densidade de probabilidade com media x̂(t�3 ) e variância �x2(t

�3 )

é combinada com a função de media z3 e variância �2z3 resultando numa nova função de densidade de

probabilidade gaussiana com media

x̂(t3) = x̂(t�3 ) +K(t3)[z3 � x̂(t�3 )] (3.13)

�2x(t3) = �2

x(t�3 )�K(t3)�

2x(t

�3 ) (3.14)

K(t3) =�2x(t

�3 )

�2x(t

�3 ) + �2

z3

. (3.15)

A melhor previsão do valor desta estimativa, antes de ser tomada a medida z3, é corrigida por um valorde peso óptimo, multiplicado com a diferença entre z3 e a previsão do seu valor. As equações da variânciae do factor de ganho K são semelhantes às equações 3.8 e 3.9.Analisando a equação para K(t3), nota-se que o seu valor diminui quando a variância �z3 aumenta. Istosignifica que é dada pouca confiança a uma medida com muito ruído e consequentemente o seu peso ébaixo.

42

No limite quando �2 ! 1,K(t3) tende para zero, e x̂(t3) é igual a x̂(t�3 ). Uma medida z3 infinitamenteruidosa seria completamente ignorada.Se a variância do ruído do sistema dinâmico �w2 tiver um valor elevado, então �x2(t

�3 ) é elevado (de

acordo com a equação XY) e K(t3) também é elevado. Nesse caso não há muita confiança no modelodinâmico do sistema e a medida terá maior peso. No limite em que �2

w ! 1,�x2(t�3 ) ! 1 e K(t3) ! 1.

A equação (1-13) resulta então em

x̂(t3) = x̂(t�3 ) + 1[z3 � x̂(t�3 )] = z3 (3.16)

No limite em que não ha nenhuma confiança no modelo do sistema, a política óptima é de ignorar a saídado modelo e usar apenas a nova medida como estimativa óptima. Caso �2

x(t�3 ) seja zero, então K(t3)

também tem valor nulo. Isto faz sentido, pois, se �2x(t

�3 ) = 0, existe certeza sobre a estimativa anterior à

medida z3 e esta pode ser desprezada.Mesmo que as equações de Kalman não foram derivadas matematicamente na totalidade e para o casogeral, o exemplo anterior serve apenas para demonstrar o princípio de funcionamento do Filtro de Kalmane a sua estrutura.

3.2 O Filtro de Kalman

O comportamento do sistema pode ser previsto. Para isso o sistema é modelado com um modelo lineardinâmico no espaço de estado. Em caso de tempo contínuo as equações correspondentes são:

˙x̄ = Fx̄+Gw̄ Modelo de estado do processo (3.17)

z̄ = Hx̄+ v̄ Modelo de observação (3.18)

No caso de o sistema alterar os valores das variáveis de estado em instantes discretos no tempo, asequações que descrevem o modelo linear dinâmico do sistema podem ser escritas da forma

x̄k+1 = �kx̄k + �kw̄k Modelo de estado do processo (3.19)

z̄k = Hkx̄k + v̄k Modelo de observação (3.20)

O Filtro de Kalman propaga tanto o estado do sistema como também a covariância do estado no tempo.A estimativa do estado anterior à incorporação de qualquer nova medição será designada por x̂�, em queo circunflexo designa uma estimativa e o sinal de menos elevado designa o momento anterior (ou a priori)à incorporação da medição durante uma iteração das equações do filtro.As equações que descrevem o Filtro de Kalman são as seguintes:

Kk = P�k HT

k [HkP�k HT

k +Rk]�1 Ganho de Kalman (3.21)

x̂k = x̂�k +Kk[zk �Hkx̂

�k ] Actualização da estimativa de estado (3.22)

Pk = [I �KkHk]P�k Actualização da covariância do estado (3.23)

x̂�k+1 = �kx̂k Previsão do estado (3.24)

P�k+1 = �kPk�

Tk + �kQk�

Tk (3.25)

43

As equações não são avaliadas todas ao mesmo tempo. As últimas duas equações são executadas a altafrequência de forma repetida, enquanto que as primeiras três, apenas são usadas para momentos em quehaja uma observação do estado através de medições.O processo é iniciado com a introdução da estimativa a priori x̂�

k e da sua covariância P�k . Para cada

iteração do modelo do sistema deve-se conhecer a matriz de transição de estado �k e a covariância deruído do sistema Qk. Para cada iteração das equações do Filtro de Kalman, a matriz do modelo deobservação Hk e a covariância das medidas Rk devem ser conhecidas a priori ou computadas a partir dasmedidas efectuadas.As equações apresentadas podem ser divididas em dois grupos, as equações de previsão e as equações decorrecção. Conceptualmente o modelo do estado do processo em 3.20 serve como previsão do estado dosistema baseando-se no modelo do sistema avançando passo a passo no tempo. É possível, quando nãohouver medidas, manter uma estimativa do estado utilisando apenas este modelo. No entanto, o errodo modelo vai-se acumulando, aumentando a incerteza sobre a estimativa a cada passo. Em termos decontrolo a execução destas equações corresponde a uma operação em malha aberta.Por outro lado, as medidas directas ou indirectas, parciais ou completas do estado, descritas pelo modelode observação em 3.20, são consideradas pelas equações de actualização de estado quando essas medidasestiverem disponíveis. Desta forma a estimativa é corrigida pelo Filtro de Kalman, combinando a estima-tiva a priori com novas medidas ruidosas para obter uma estimativa melhorada a posteriori. O Ganho deKalman pondera a incerteza da estimativa a priori e a incerteza existente nas medidas com o objectivode diminuir a incerteza da estimativa do estado do sistema. Este passo representa em termos de controlouma retroacção em malha fechada. O processo está esquematizado na figura 3.5.O primeiro passo durante a actualização do estado após uma medida consiste em calcular o ganho deKalman Kk através da equação 3.22. O próximo passo consiste em obter a medida zk e calcular umaestimativa a posteriori do estado incorporando a medida como na equação 3.23. O último passo é obteruma estimativa da covariância a posteriori do erro através da equação 3.24.Depois de cada medição e actualização o processo de previsão é repetido usando a estimativa a poste-riori para prever uma nova estimativa a priori. Desta forma a estimativa do estado está condicionadarecursivamente por todas as medidas obtidas anteriormente.

Figura 3.5: Esquema do funcionamento do Filtro de Kalman e da aplicação das suas equações [19]

44

Se o modelo do sistema fosse perfeito não haveria necessidade de realimentação e o modelo do estadopoderia correr em malha aberta. Os sensores são necessários pois não existem modelos perfeitos. Destaforma as equações de actualização fornecem a realimentação necessária na prática. No entanto, existecompleta independência entre o modelo e a realimentação através das medidas.É importante distinguir as diferentes matrizes de covariância nas equações. A matriz Qk representa aincerteza que corrompe o modelo do sistema. A matriz Rk modela a incerteza associada às medidas. Osseus elementos devem ser expressos nas unidades e coordenadas usadas pelos sensores, não nas coordena-das dos estados que medem. Finalmente a matriz Pk é gerida pelo próprio filtro e representa a incertezaglobal existente na estimativa ao longo do tempo.

3.3 Aplicação do Filtro de Kalman

A aplicação do Filtro de Kalman para o seguimento do movimento de um animal detectado implicaprimeiramente o desenvolvimento de um modelo linear de movimento.O peixe marcado pode-se movimentar livremente podendo alterar a direcção ou velocidade em qualquermomento. O movimento do animal marinho está no entanto sujeito às leis físicas do movimento. Istosignificando que só é possível uma alteração suave e contínua do movimento. Por exemplo, um peixe nãose pode deslocar numa direcção e abruptamente, de um momento para o outro, alterar essa direcção paraa direcção oposta, pois o animal está sujeito à inércia do seu movimento. Não sendo possível prever asalterações que o peixe possa efectuar durante o seu percurso, a postura mais justa será de assumir que, deforma geral, o animal manterá a tendência do seu movimento. Um modelo adequado será então o modelode movimento com aceleração constante. A incerteza sobre a mudança de percurso ou de aceleração éconsiderada com a introdução de ruído no modelo linear do sistema.A equação 3.26 descreve de forma discreta a posição p no instante k + 1 com base na posição anteriorpk, na velocidade anterior vk e na aceleração anterior ak. De forma semelhante as equações 3.27 e 3.28permitem obter o valor da velocidade vk+1 e da aceleração ak+1 para o instante k+1, sendo �t o tempoque decorre entre duas iterações.

pk+1 = pk + vk�t+1

2

ak�t2 (3.26)

vk+1 = vk + ak�t (3.27)

ak+1 = ak (3.28)

O estado pode ser identificado como sendo o conjunto das variáveis que compõem as equações nas trêsdimensões do espaço. O vector que define o estado x é dado através de

x̄ =

⇥px py pz vx vy vz ax ay az

⇤T (3.29)

Então a equação linear que descreve o movimento discreto de forma matricial é dada por

x̄(k+1) =

2

6666666666664

1 0 0 �t 0 0

12�t2 0 0

0 1 0 0 �t 0 0

12�t2 0

0 0 1 0 0 �t 0 0

12�t2

0 0 0 1 0 0 �t 0 0

0 0 0 0 1 0 0 �t 0

0 0 0 0 0 1 0 0 �t0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 0 1

3

7777777777775

x̄k + w̄ (3.30)

Para introduzir a incerteza relativa à mudança de percurso do animal soma-se o ruído w ao modelo dosistema. Este ruído é ruído branco gaussiano que incide apenas sobre as componentes de aceleração.

45

As variáveis que descrevem as componentes de aceleração têm impacto não só nas acelerações futurasmas também têm efeito sobre as velocidades e finalmente no cálculo das posições futuras devido à formacomo as equações de sistema estão interligadas. Assim sendo, o ruído usado neste trabalho apenas temcomponentes não nulas nas componentes de aceleração e a matriz de covariância do ruído do modelo desistema pode então ser descrita por:

Qk =

2

6666666666664

14�t4 0 0

12�t3 0 0

12�t2 0 0

0

14�t4 0 0

12�t3 0 0

12�t2 0

0 0

14�t4 0 0

12�t3 0 0

12�t2

12�t30 0 0 �t2 0 0 �t 0

0

12�t3 0 0 �t2 0 0 �t 0

0 0

12�t3 0 0 �t2 0 0 �t

12�t2 0 0 �t 0 0 1 0 0

0

12�t2 0 0 �t 0 0 1 0

0 0

12�t2 0 0 �t 0 0 1

3

7777777777775

q2mag (3.31)

O parâmetro qmag descreve a magnitude do ruído e numa aplicação real terá que ser estimado. Isto poderáser complicado, pois envolve conhecimento do movimento do animal. Um valor que funcione bem paraum animal, poderá não ser adequado para outro animal. Da mesma forma se deverá adaptar a matriz decovariância dos sensores e os valores de rij aos sensores usados na práctica. Geralmente pode-se partirdo principio de que não existe correlação entre os ruídos incidentes sobre os diferentes eixos. A matriz decovariância das medidas é dada por:

Rk =

2

4rxx 0 0

0 ryy 0

0 0 rzz

3

5 (3.32)

A matriz da covariância das medidas Rk apenas tem as dimensões de 3 por 3, porque o método usadoapenas chega a uma medida da posição e não de todas as variáveis de estado que compõem o vector x.Este método encontra-se descrito com mais pormenores no capitulo 4.1.2.A matriz H que descreve o modelo de observação tem por isso a forma

H =

2

41 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0

3

5 (3.33)

Assim ficam definidas as principais matrizes necessárias para o funcionamento do Filtro de Kalman.

46

Capítulo 4

Simulação do Portable Tool em

funcionamento com aplicação do Filtro

de Kalman

4.1 Cenários de simulação

No projecto MAST/AM propôs-se desenvolver um sistema que pudesse servir para a detecção e o segui-mento de animais marinhos previamente marcados [20]. Para isso serão desenvolvidos várias ferramentaslistadas de seguida.

Surface Robotic Tool - Equipamento fixado numa bóia para flutuar à superfície a largo prazo comarray de sensores em estrutura USBL. Detecção e identificação de uma ou várias tags acústicasobtendo estimativa na direcção de chegada do sinal.

Portable Tool - Equipamento portátil para uso subaquático por um mergulhador com array de sensoresem estrutura USBL. Detecção e identificação de uma ou várias tags acústicas obtendo estimativana direcção de chegada do sinal.

Operation Console - Consola de operação a ser instalada num barco de apoio ou num laboratóriocientífico situado na costa. Junção e processamento offline dos dados recolhidos.

Estes equipamentos poderão ser usados nos seguintes cenários:

- Só se usa o Surface Robotic Tool ou o Portable Tool, ou seja apenas se usa um deles. Este cenáriopermite a detecção e identificação de uma ou de várias tags acústicas do mesmo tipo com ou seminformação de sensores de pressão. Se estiver disponível a informação de sensores de pressão, estavai permitir uma estimativa de localização do animal marcado ou então em caso que esta informaçãonão esteja disponível, obtem-se apenas a direcção de chegada do sinal.

- Se se usar o Surface Robotic Tool em conjunto com o Portable Tool é possível obter uma estimativade localização do animal marcado mesmo sem informação adicional de sensores de pressão. Numprocessamento posterior é possível obter uma estimativa sobre a trajectória percorrida pelo animalmarcado.

Com este equipamento robótico pretende-se obter períodos de operação de (semanas para o sistema desuperfície, 3 a 4 horas para o equipamento subaquático, detecção a distâncias de 600 a 1000 metros. Umaprecisão de 1 a 2 metros na estimativa da posição do animal marinho e estimativas sobre a trajectóriadistanciadas de 1 a 5 segundos (com ambos equipamentos) [20].O Portable Tool tem a possibilidade de detectar um animal marinho previamente marcado e consegueestimar a direcção de chegada do sinal. Existem vários cenários para o funcionamento do Portable Tool

47

4.1.1 Cenário apenas com o Portable Tool

No cenário em que se usa apenas o Portable Tool é possível indicar a posição do animal marcado, sehouver informação disponível sobre a pressão no local do animal. Esta informação é disponiblizada porum sensor de pressão dentro da tag. A pressão está relacionada com a profundidade, no entanto, comojá se viu na secção 1.3, o perfil de pressão pode ter grandes variações que comprometam a precisão e afiabilidade da estimativa da posição.Se se ignorar as irregularidades no perfil de pressão usando o modelo de uma coluna aquática estacionária,a pressão p ficará relacionada com a profundidade h da seguinte forma, em que ⇢ é a densidade da águae g a aceleração da terra:

p = h⇢g (4.1)

O cálculo da posição limita-se então à intersecção de uma recta, que é definida pelo próprio USBL epela sua estimativa da direcção de chegada de sinal, com um plano paralelo à superfície da água que seencontra situado a uma profundidade h.

4.1.2 Cenário com Portable Tool e Surface Robotic Tool

No cenário de utilização de dois USBLs como demonstrado na figura 4.1, é possível localizar o animalcom maior precisão do que na variante com apenas um USBL e uma tag com medidor de pressão. Nestecaso basta ter uma tag simples que apenas identifique o animal. Existe uma estação base localizada numnavio em que um dos USBLs é colocado numa bóia (SUrface Robotic Tool) na proximidade do naviocom o qual comunica via canal rádio. O outro USBL encontra-se no Portable Tool nas mãos de ummergulhador como no cenário anterior. Deve-se distinguir entre o funcionamento online e offline

Funcionamento online é o funcionamento visível para o mergulhador quando este está a seguir oanimal. O funcionamento é semelhante ao cenário anterior, sem a estimativa da profundidade.

Funcionamento offline é uma fase posterior ao mergulho em que se trabalha sobre dados gravados.Para isso é necessário sincronizar os dados de ambos os USBLs. Depois de estar estabelecida a sin-cronização, poder-se-á calcular uma estimativa da posição através da intersecção das rectas que sãoformadas pelas posições conhecidas dos USBLs e das estimativas das direcções de chegada em cadaum. Pode-se usar um Filtro de Kalman para obter uma estimativa relativamente à rota tomadapelo animal entre pontos adjacentes dado que os pontos poderão estar afastados uns dos outrosporque as tags comerciais apenas emitem aproximadamente sinal de minuto a minuto. A junçãodos dados só pode acontecer offline porque se assim não fosse, em online os USBLs teriam quecomunicar entre si em tempo real. Neste caso a capacidade de processamento seria desnecessaria-mente sobrecarregada porque a informação adquirida a mais seria de pouco ou nenhum interessepara o mergulhador que apenas quer seguir o animal.

48

Figura 4.1: Cenário típico de missão para o projecto MAST/AM com Surface Robotic Tool e PortableTool [20]

O cálculo da posição do animal marcado num cenário offline exige o cálculo da intersecção entre duasrectas. Neste caso há que considerar as seguintes possibilidades:

- Duas rectas com um ponto de intersecção

- Duas rectas paralelas sem ponto de intersecção

- Duas rectas reversas, não paralelas sem intersecção mas com um ponto de maior aproximação

No primeiro caso das rectas com um ponto de intersecção, a determinação da localização do animal ésimples. No caso de duas rectas paralelas não é possível determinar um ponto único. Este caso nãose verifica frequentemente e deve ser ignorado, não dando origem a nenhuma estimativa de posição. Ocaso mais frequente deverá ser o de duas linhas reversas, não paralelas, sem ponto de intersecção. Emseguida apresenta-se uma forma de cálculo do ponto de maior aproximação como melhor estimativa paraa posição do animal.

Ponto de maior aproximação entre duas rectas reversas

Para o cálculo da posição do animal no caso de existirem duas rectas reversas sem ponto de intersecção,a melhor estimativa para a posição do animal deve ser o ponto centrado entre os pontos de maioraproximação em cada uma das rectas.Considere-se portanto duas rectas L1 e L2:

L1 : P (s) = P0 + s(P1 � P0) = P0 + su (4.2)

L2 : Q(t) = Q0 + t(Q1 �Q0) = Q0 + tv (4.3)

Considera-se w(s, t) = Ps � Qt como o vector entre os dois pontos situados nas duas rectas. Então oobjectivo é encontrar o ponto equidistante e mais próximo possível das duas rectas.Seguindo o raciocínio em [21], é preciso, num primeiro passo, encontrar w(s, t) com tamanho mínimopara qualquer que seja s e t. Em qualquer espaço n-dimensional as duas rectas L1 e L2 minimizam a

49

sua distância nos dois pontos Pc = P(sc) e Qc = Q(tc) para os quais w(sc, tc) é o mínimo para w(s, t).O problema encontra-se esquematizado na figura 4.2. No caso de L1 e L2 não serem paralelos e não seintersectarem, então o segmento PCQC que une esses dois pontos é simultaneamente perpendicular àsduas rectas e não existe mais nenhum segmento com as mesmas propriedades. O vector wc = w(sc, tc) éperpendicular aos dois vectores de direcção das rectas u e v, o que pode ser expresso nas duas equações:u ·wc = 0 e v ·wc = 0.

Figura 4.2: Visualização dos pontos de maior aproximação para duas rectas sem intersecção [21]

É possível resolver essas duas equações substituindo wc = P (sc)�Q(tc) = w0+scu�tcv, com w0 = P0�Q0,em ambas, obtendo assim as seguintes duas equações lineares:

(u · u)sc � (u · v)tc = �u · w0

(v · u)sc � (v · v)tc = �v · w0

Admitindo que a = u · u, b = u · v, c = v · v, d = u ·w0, e e = v ·w0 e resolvendo por sc e tc, a solução é:

sc =be� cd

ac� b2e tc =

ae� bd

ac� b2

com o denominador ac�b2 diferente de zero. Note-se que ac�b2 = |u|2|v|2�(|u||v|cos✓)2 = (|u||v|sin✓)2 �0. Quando ac � b2 = 0, as duas equações são dependentes, as duas linhas são paralelas, e a distânciaentre as duas linhas é constante, um caso sem interesse para o propósito deste trabalho.Tendo resolvido por sc e tc, é possível determinar os dois pontos Pc e Qc nas duas rectas L1 e L2 quemais se aproximam um do outro. Resta determinar o ponto no centro do segmento que une os pontos Pc

e Qc. Para esse efeito basta o cálculo P = (Qc � Pc)/2.A distância entre o ponto Pc e Qc é de certa forma uma medida sobre a certeza ou a fiabilidade da medidaobtida por este método. Assim sendo, a norma do vector

���!PcQc pode ser usada para estimar a incerteza

existente sobre a medida, que é representada pela matriz de covariância das medidas Rk apresentada nasecção 3.3.

Transformação entre sistemas de coordenadas

O objectivo é a localização de um animal marcado. Para esse efeito é necessário definir um referencialpara poder descrever essa localização, a qual deverá ser definida de forma inequívoca. Para o conseguir,existem referenciais como o referencial geodésico, que descreve a posição em relação à terra usando valoresde latitude �, longitude ' e altitude h, ou o referencial ECEF (com origem Oe no centro da terra, eixo Ze

a apontar para o pólo norte, eixo Ye a intersectar a latitude de 0 � no plano equatorial e eixo Ze definidosegundo a regra da mão direita).No caso em que um sistema tenha uma variação de localização relativamente pequena em relação a estesreferenciais ligados à terra, é conveniente definir um referencial local e descrever as diferentes localizações

50

no referencial local. Quando se desejar, as coordenadas locais poderão ser convertidas do referencial localpara o referencial geodésico ou para outro referencial que seja conveniente. Para uma descrição maisdetalhada dos referenciais de coordenadas e da conversão entre eles veja-se o literatura indicada por [22].Cada um dos USBLs, tanto no Surface Robotic Tool como no Portable Tool, tem um referencial decoordenadas local associado a cada um deles e todas as medições de direcção de chegada de sinal sãofeitas nesse referencial local.O referencial local está definido como referencial de objecto e segue o movimento do USBL de forma aque as coordenadas estejam sempre definidas da mesma forma em relação ao objecto. A sua origem e oseixos que formam o referencial são definidos da seguinte forma:

1. A origem Ob está localizada no centro de gravidade do objecto associado.

2. O eixo Xb aponta para diante na direcção principal de deslocação do objecto.

3. O eixo Yb aponta para estibordo (lado direito do veiculo olhando em direcção do eixo Xb).

4. O eixo Zb aponta para baixo para seguir a regra da mão direita.

Enquanto que no Surface Robotic Tool a definição da frente poderá ser algo ambígua, pois não existeuma direcção principal de movimento (movimento da bóia apenas influenciada pelas correntes e pelaondulação), no caso do Portable Tool a direcção de movimento principal é determinada pelo movimentodo mergulhador com a estrutura de sensores a apontar para a frente na mesma direcção que o eixo Xb.Devido à semelhança geométrica do Surface Robotic Tool com o Portable Tool, o seu referencial serádefinido correspondentemente.Para juntar a informação sobre a direcção de chegada de sinal facultada por um USBL com a informaçãoobtida no outro USBL a fim de calcular o ponto de maior aproximação como já foi descrito anteriormente,é primeiramente necessário definir um referencial principal para o qual se possam converter as coordenadasde cada um dos USBLs.Por motivos de conveniência define-se este referencial principal como o referencial NED centrado noSurface Robotic Tool, até porque faz sentido definir a origem do referencial à superfície da água para seobter de forma imediata o valor da profundidade na leitura das coordenadas. Na realidade, a ondulaçãoe as correntes causam algum movimento no Surface Robotic Tool e com ele, no referencial nele centrado.Para situações de simulação pode-se modelar este referencial como sendo estático servindo como referenciaprincipal local.O referencial NED fixado ao objecto está ligado ao Surface Robotic Tool e define-se da seguinte forma(compare-se com a figura 4.3):

1. A origem On é fixada ao centro de gravidade do objecto associado.

2. O eixo Xn aponta para o norte geodésico (pólo norte).

3. O eixo Yn aponta para o leste geodésico.

4. O eixo Zn aponta para baixo na direcção da normal da elipsóide que modela a forma da terra.

51

Figura 4.3: Visualização do referencial geodésico, o referencial ECEF e do referencial NED local [22]

A problemática relativa à localização de cada um dos USBLs e da estimativa das suas atitudes não éobjecto deste trabalho dado que este tema já foi abordado em outros trabalhos (ver [23]).Uma vez que se conheça a posição e a atitude de cada um dos USBLs, a conversão de coordenadas doreferencial local para o referencial principal é a representação das coordenadas locais em coordenadas doreferencial principal (referencial NED centrado no Surface Robotic Tool).Consideremos o referencial principal aqui denominado {A}, neste caso, referencial cartesiano com basecanónica, então esta base A pode ser descrita de forma matricial (com os vectores unitários que a definemnas colunas):

A =

2

41 0 0

0 1 0

0 0 1

3

5 (4.4)

Em coordenadas homogéneas o referencial {A} será descrito por uma matriz 4 ⇥ 4. O quarto vector decoluna é dado pelas coordenadas da origem do referencial {A} descrito em relação à base do referencial{A} (neste caso [0001]

T ):

{A} =

2

664

1 0 0 0

0 1 0 0

0 0 1 0

0 0 0 1

3

775 (4.5)

Defina-se um referencial de objecto {B}, também cartesiano mas diferente de {A}, que por sua vez estádefinido pelos três vectores unitários ~Bx, ~By, ~Bz, que constituem a sua base, e a origem Bw. Descrevendoo referencial {B} também de forma matricial, mas com as suas componentes descritas em relação àscoordenadas do referencial {A} obtém-se assim a matriz A

BM de conversão de coordenadas do referencial{B} para o referencial {A}.

ABM = (

~Bx, ~By, ~Bz, Bw) (4.6)

52

Figura 4.4: O ponto p em relação aos referenciais A e B

Qualquer ponto p é representável tanto no referencial {A} como no referencial {B}. A matriz BAM

representa a transformação necessária para se poder descrever um ponto p em relação ao referencial {B}(denominado Bp), dada uma descrição em relação ao referencial {A} (denominado Ap). A figura 4.4mostra o ponto p em relação aos referenciais {A} e {B}.Para se poder especificar um ponto Bp em relação ao referencial {A} utiliza-se a matriz inversa a B

AM :

BAM ·A p =

B p (4.7)

BAM

�1 ·B p =

AB M ·B p =

A p (4.8)

A matriz BAM representa as rotações e a translação necessárias para representar o referencial {B} em

cima do referencial {A}, translação e rotações essas que podem ser subdivididas em passos sequenciais aformarem referenciais intermédios.Em primeiro lugar, o referencial principal {A} deve ser deslocado para coincidir com a sua origem naorigem do referencial do objecto {B}.Em coordenadas homogéneas uma translação T para o ponto definido por [tx, ty, tz]T pode ser apresentadapela matriz de translação

B0

A T =

2

664

1 0 0 tx0 1 0 ty0 0 1 tz0 0 0 1

3

775 . (4.9)

É possível efectuar uma operação de translação e de rotação tanto sobre os eixos do referencial (AliasTransformation), como também sobre os pontos descritos pelo referencial (Alibi Transformation). Umaoperação de translação ou de rotação efectuada sobre os eixos, corresponde à operação inversa sobre ospontos descritos pelo mesmo referencial.Assim, no caso de se querer representar o referencial {A} em cima do referencial {B}, essa translaçãocorresponde para os pontos a uma operação T (�Bw) com Bw expresso em coordenadas em relação aoreferencial {A}.Após a translação o referencial assim criado {A0} é rodado para que a orientação dos eixos coincida coma orientação relativa dos eixos do referencial {B}. Geralmente a atitude de um veículo em relação aum referencial NED nele centrado é descrito pelos ângulos denominados ângulos de Euler (Yaw, Pitch eRoll), que descrevem a atitude como a sequência de rotações intrínsecas em torno dos eixos na ordem Z-Y-X, em que a orientação dos eixos é alterada após cada rotação (compare-se a figura 4.5). Desta formaa orientação do referencial {A0} em relação ao referencial {B} pode ser decomposto em três rotaçõeselementares [22]:

Yaw é o angulo entre a projecção do eixo Xb do referencial do objecto no plano Xn - Yn do referencialNED centrado no objecto com o eixo Xn. A rotação no sentido dos ponteiros do relógio em tornodo eixo Zn dá origem a um referencial intermédio {A1}.

53

Pitch é o angulo entre o eixo XA1 do referencial intermédio {A1} para com o eixo Xb do referencialdo objecto. A rotação no sentido dos ponteiros do relógio em torno do eixo YA1 do referencialintermédio {A1} dá origem a um segundo referencial intermédio {A2} cujo eixo XA2 coincide como eixo Xb do objecto.

Roll é o angulo entre o eixo YA2 (ou ZA2) do referencial intermédio {A2} para com o eixo correspondentedo referencial do objecto. A rotação no sentido dos ponteiros do relógio em torno do eixo XA2 dosegundo referencial intermédio {A2} dá origem ao referencial do objecto {B}.

Figura 4.5: Ângulos de Euler e sequência de rotação yaw-pitch-roll [23]

Assim, a transformação (Alibi) do referencial A0 para B, dado os ângulos de Euler, pode ser expressacomo o produto de três matrizes de rotação, em que c ⌘ cos e s ⌘ sin e �, ✓ e representam os ângulosYaw, Pitch e Roll respectivamente:

BA0R =

BA2

R ·A2A1

R ·A1A0 R =

2

664

1 0 0 0

0 c� s� 0

0 �s� c� 0

0 0 0 1

3

775

2

664

c✓ 0 �s✓ 0

0 1 0 0

s✓ 0 c✓ 0

0 0 0 1

3

775

2

664

c s 0 0

�s c 0 0

0 0 1 0

0 0 0 1

3

775 (4.10)

BA0R =

2

664

c✓c c✓s �s✓ 0

s�s✓c � c�s s�s✓s + c�c s�c✓ 0

c�s✓c + s�s c�s✓s � s�c c�c✓ 0

0 0 0 1

3

775 (4.11)

Para representar um ponto, descrito em coordenadas do referencial {A}, em relação ao referencial {B},a matriz de conversão é dada pelo conjunto das matrizes de translação e de rotação:

BAM =

BA0 R ·A

0

A T (4.12)

A matriz BAM serve portanto para representar um ponto no referencial local de um USBL, dada uma

representação em coordenadas do referencial principal. Por outro lado, a matriz inversa BAM

�1=

AB

M serve para representar um ponto em coordenadas do referencial principal, dada uma descrição emcoordenadas do referencial local de um dos USBLs.

4.2 Estrutura da simulação

A simulação mostra o principio de funcionamento das estratégias de detecção e de seguimento previamentedescritas. De acordo com os dois cenários identificados no capitulo 4.1, a simulação funciona de dois modosdistintos, com um ou dois USBLs para permitir a identificação e localização de um animal marcado. Omovimento deste animal marinho pode ser descrito com uma trajectória pré-definida por uma série de

54

pontos ou pode ser simulada com um modelo de movimento do animal (ver capítulo 3.3). O modelo demovimento é dependente do estado inicial do animal e do nível de ruído a que está sujeito. Ao final dasimulação é gerado um ficheiro “output.txt” que resume os valores das estimativas de direcção, ânguloscorrespondentes em coordenadas esféricas e os pontos medidos.A simulação está implementada seguindo a metodologia de programação orientada a objectos. Destaforma, o código fonte é organizado em blocos lógicos facilitando assim a leitura e a manutenção domesmo. Existem vários tipos de objectos definidos e usados que estão brevemente descritos como segue:

Peixe contém um ID que o identifica e um sinal que é emitido, através do qual o peixe pode ser detectadoe identificado.

Canal adiciona dois tipos de ruído branco gaussiano ao sinal emitido pelos peixes. O primeiro tiposimula o ruído originado pelo canal aquático através do qual o sinal se propagou. Além disso,simula-se um ruído intrínseco aos hidrofones que formam os sensores do sistema.

Evento simboliza o momento de emissão de sinal por parte de um peixe e a sua possível detecção pelosUSBLs. Um objecto deste tipo contém informação sobre o Peixe que o originou: ID, sinal, posiçãoe instante de tempo

Lista de Eventos trata da gestão dos Eventos gerados e organiza-os segundo uma ordem cronológica.Quando dois Eventos coincidirem em momentos muito próximos a uma distância de menos que 4segundos, a Lista de Eventos junta-os para formarem apenas um. Esta operação corresponde auma colisão de dois conjuntos de sinal. Isto terá como consequência uma não detecção dos eventosoriginais por parte dos USBLs. Os USBLs recebem os Eventos da Lista de Eventos.

USBL detecta o sinal corrompido por ruído e descodifica-o. O USBL devolve o ID do peixe identificadobem como a direcção de origem determinada pelo método descrito na secção 2.0.6.

Sistema representa o núcleo da simulação. É neste objecto que são adicionados à simulação dois objectosdo tipo USBL e os peixes que devem ser simulados. Este objecto trata da gestão e da criação deeventos bem como da detecção dos mesmos, obtendo a localização do peixe através do métododescrito na secção 4.1.2. O sistema organiza a interacção de todos os objectos necessários paraa simulação sendo também responsável por acções como a apresentação gráfica dos resultados desimulação e o registo de Peixes para validação do Checksum.

Detecção combina o evento com a respectiva detecção, mantendo também métricas para a qualidade dadetecção. A precisão mede a distância entre as duas rectas reversas criadas pelo prolongamento dasdirecções de origem indicadas pelos USBLs. Ela dá uma ideia sobre a confiança a dar ao método deintersecção para cada medição, sem dar, no entanto, informação sobre a qualidade real da mediçãoem si. Quanto maior for o valor, menor será a confiança a dar à medição. O erro mede a distânciaentre a localização real do Peixe durante um Evento e a posição medida pelo Sistema indicandoassim a qualidade real da medição.

Modelo de Peixe é um modelo interno que o sistema mantém por forma a estimar o estado de umpeixe. Depois de detectado um peixe, o sistema cria ou actualiza um modelo já existente do peixedetectado. Este modelo faz uso do Filtro de Kalman havendo a possibilidade de efectuar um passode previsão ou um passo de actualização como descrito na secção 3.3.

A figura 4.6 pretende visualizar as tarefas de cada um dos objectos e da sua interligação durante asimulação.

55

Figura 4.6: Funcionamento esquematizado do programa de simulação

4.3 Resultados de Simulação

A simulação mostra o princípio de funcionamento das estratégias de detecção e de seguimento previamentedescritas. Apesar da simulação não conseguir simular todos os problemas relacionados com a propagaçãoacústica no meio subaquático, pode-se estimar o alcance de uma tag para modelos de propagação esféricaou cilíndrica. É também possível avaliar o desempenho do algoritmo usado na estimativa da direcção dechegada do sinal ou a problemática de colisão de sinal. A partir da simulação é também possível tirarconclusões sobre a precisão das medidas e sobre o desempenho do sistema para diferentes configuraçõesdo Filtro de Kalman.

Detecção de sinal

56

57

58

Seguimento do animal marcado

O algoritmo, descrito no capitulo 2.0.6, apresenta um bom desempenho na estimação da direcção de che-gada do sinal com uma precisão na ordem de 0,1 a 0,2 graus, quando descrito em ângulos de representaçãoem coordenadas esféricas.Para um cenário de apenas um USBL, a tag transmite medidas da pressão local. As variações no perfilde pressão com a profundidade, abordadas no capitulo 1.3, introduzem alguma incerteza, na relaçãoentre a pressão e a profundidade, factor principal que influencia o desempenho do método de mediçãode posição, descrito na secção 4.1.1. Dependendo das tags usadas, a precisão do método de localizaçãoestará principalmente determinada pela precisão da medida de profundidade obtida pela tag, que paraas tags comerciais varia entre ± 1,7m e ± 34 m [24]. Entende-se facilmente que o método terá maudesempenho ou que não dará mesmo resultado nenhum para situações em que o USBL se encontrar nomesmo plano ou muito próximo do plano definido pela profundidade da tag.Considerando o caso de haver dois USBLs, aplica-se o método de intersecção de duas rectas, descrito nasecção 4.1.2. A partir da simulação é possível mostrar que este método consegue atingir uma precisãopor vezes inferior a 5 metros mas quase sempre inferior a 10 metros em distâncias de 600 a 1000 metros,como visível na figura 4.7 e na figura 4.8. Para distâncias inferiores é possível até uma precisão de 1 a2 metros. Este método de localização terá pela sua natureza mau desempenho para situações em que osdois USBLs formem uma recta com o animal marcado.

59

(a) Visualização gráfica das previsões e actualizações das estimativas de posição em comparação com a trajectória

real do peixe.

(b) A precisão das medições e da estimativa conseguem estar menor que 10 metros.

Figura 4.7: Detecção de um peixe marcado. (Posições:Surface Robotic Tool - [0;0;0], Portable Tool -[800;100;150], Peixe inicialmente - [-350;100;150])

60

(a) Visualização gráfica das previsões e actualizações das estimativas de posição em comparação com a trajectória

real do peixe.

(b) A precisão das medições e da estimativa conseguem estar menor que 5 metros.

Figura 4.8: Detecção de um peixe marcado. (Posições:Surface Robotic Tool - [0;0;0], Portable Tool -[800;100;150], Peixe inicialmente - [400;400;400])

A introdução de um Filtro de Kalman não aumenta significativamente a precisão, por vezes até prejudicaa desempenho do sistema na estimativa da posição do animal. Isto deve-se também à difícil configuraçãodos valores adequados para as matrizes de covariância que fazem parte do modelo de movimento. Autilidade do Filtro de Kalman consiste antes em permitir uma estimativa a cada a um ritmo superior queas medidas e também em ter uma estimativa da posição quando não há medida, caso que se pode darnuma situação de colisão de sinal. A figura 4.9 mostra a ocorrência de colisões durante o deslocamento

61

do peixe.

Figura 4.9: Os pontos da trajectória do peixe que deram origem a uma collisão são marcados com umcirculo vermelho

62

Capítulo 5

Conclusões e Trabalho Futuro

Em qualquer estudo de marcação e de seguimento de animais marinhos deve-se em primeiro lugar colocara questão sobre a sua utilidade e o seu modo de realização. Mesmo que os produtores de equipamento detagging visem minimizar o impacto sobre a vida diária do animal, qualquer estudo que envolva uma inter-ferência na rotina dos animais, terá sempre como consequência algum transtorno ou stress no indivíduo.Um possível estudo deverá por isso ter sempre uma motivação forte. Um planeamento detalhado comuma preparação adequada é essencial para ajudar a diminuir o impacto nos animais e evitar problemasdurante a execução do estudo. A selecção de equipamento depende das condições impostas pelo ambi-ente e do animal. Desta forma uma tag que emita á frequência usada pelas tags comerciais estará maisapropriada para espécies como o salmão atlântico ou o american shad, não sensíveis a esta frequência, doque por exemplo para golfinhos (comparar audiogramas na secção 1.1).Durante o planeamento dever-se-á ter em consideração também problemas relacionados com a acústicasubaquática. Dependendo do ambiente e das condições externas do local onde se pretende efectuar umestudo, pode haver problemas na propagação de som que dificultem ou possam até inviabilizar a detecçãode animais marcados. Zonas de sombra podem surgir tanto no oceano como em zonas costeiras. Factorescomo a turbulência na água devido a correntes, à salinidade da água ou à temperatura e pressão afectama propagação do som na água. Mesmo que o sinal emitido por um animal possa ser detectado, a premissada propagação rectilínea de frente de onda plana poderá estar comprometida. A ondulação, bolhas de arà superfície, o tipo do solo e objectos na água, como os cais em zonas portuárias, podem dar origem àcriação de ecos que podem eventualmente causar sérios problemas para a detecção de um sinal acústicoou para a estimativa de direcção de chegada, uma vez que estes se baseiam na detecção de picos e nasTDOAs.Existe equipamento variado para diferentes tipos de utilização. Assim, também as tags comerciais apre-sentam apenas um compromisso entre os diversos factores. Os sensores devem ter um tamanho mínimo,mas uma longevidade máxima. Uma detectabilidade elevada é desejável para melhor poder seguir o per-curso tomado pelo animal. Por outro lado, uma emissão muito frequente de sinal tem impacto na bateriae pode dar origem a colisões de sinal que inviabilizam a detecção quando muitos animais marcados seencontram juntos uns dos outros como poderia acontecer num cardume.O espaçamento entre duas transmissões consecutivas por uma tag comercial (com cerca de um a poucosminutos) faz com que uma aplicação do Filtro de Kalman não faça muito sentido para peixes que façammuitas mudanças de direcção nesse intervalo de tempo ou que se mantenham num local muito restritocomo pode ser o caso de um peixe à procura de alimento num coral. A aplicação do Filtro de Kalman estáantes destinada para cenários de migração. A simulação mostra também a importância da parametrizaçãocorrecta do Filtro de Kalman, bem como a importância de evitar o posicionamento dos USBLs em linhacom o animal marcado ou na mesma profundidade, caso seja usado apenas um USBL. Em caso de havercolisões de sinal, a consequente falta de medida pode ser compensada pelo Filtro de Kalman.A simulação serve para mostrar o funcionamento básico do sistema de localização e de seguimento deanimais marcados, composto pelo Surface Robotic Tool e pelo Portable Tool em dois cenários de trackingactivo. Embora o sistema não consiga cumprir com os objectivos declarados pelo projecto MAST/AM,que pretende estimativas de localização com precisão de 1 a 2 metros para distâncias de 600 a 1000metros, consegue-se estimativas com precisão na ordem dos 5 a 10 metros para essa gama de distâncias.

63

Mostrou-se também que o sistema é capaz de detectar uma tag a essas mesmas distâncias para condiçõesidealizadas na propagação típica em oceanos - uma tag conseguiu ser detectada a uma distância de 780metros (comparar secção 4.3).Assim o sistema é comparável na sua precisão e no seu alcance de detecção com o sistema de localizaçãocomercial, compare-se as secções 1.2.2 e 1.2.1. No entanto o sistema apresentado neste trabalho traz avantagem de apresentar características típicas de tracking passívo, como a identificação da tag e a gravaçãode dados, em cenários de tracking activo. Alem disso é possível obter uma estimativa da posição do animalmarcado em cada 1 a 5 segundos. Desta forma, o sistema composto pelo Portable Tool e pelo SurfaceRobotic Tool, tem as suas limitações mas distingue-se da solução existente comercial, essencialmente pelasua facilidade de manuseamento e pela sua versatilidade de uso.Em trabalhos futuros a simulação poderá ser acrescentada por mais factores de propagação acústica,principalmente pela propagação por trajecto múltiplo ou outros factores que a aproximem a cenáriosmais realistas. Mesmo assim será imprescindível submeter o Surface Robotic Tool e o Portable Toola testes práticos para se poder obter mais conclusões sobre a viabilidade e as limitações reais destasferramentas. Assim, será possível identificar o tipo de cenários óptimos para a utilização do sistema dedetecção e seguimento aqui discutido.

64

65

Apêndice A

Modelo de absorção acústica do meio

aquático de Francois e Garrison

O modelo de absorção de Francois e Garrison[7] acrescenta perdas adicionais às perdas devidas à disperçãodo sinal no espaço. Seguindo este modelo, o coeficiente de absorção é decomposto em três parcelas,correspondentes à influência do acido bórico, do sulfato de magnésio e da viscosidade da água:

↵ = A1P1f1f2

f21 + f2

+A2P2f2f2

f22 + f2

+A3P3f2 (A.1)

onde ↵ é a atenuação em dB/km. z é a profundidade em m; S é a salinidade em p.s.u. (practical salinityunits); T é a temperatura em �C; e f é a frequência em kHz.A contribuição devida ao acido bórico B(OH)3 é dada por:

8>>>>>><

>>>>>>:

A1 =

8, 86

c10

(0,78pH�S))

P1 = 1

f1 = 2, 8p(

S

35

)10

(4�1245/(T+273))

c = 1412 + 3, 21T + 1, 19S + 0, 0167z

(A.2)

A contribuição devida ao sulfato de magnésio Mg(SO)4 é dada por:

8>>>>><

>>>>>:

A2 = 21, 44S

c(1 + 0, 025T )

P2 = 1� 1, 37⇥ 10

�4z

f2 =

8, 17⇥ 10

(8�1990/(T+273))

1 + 0, 0018(S � 35)

(A.3)

A contribuição devida à viscosidade da água é dada por:The contribution of pure water viscosity reads:

8><

>:

P3 = 1� 3, 83⇥ 10

�5z + 4, 9⇥ 10

�10z2

T < 20

�C ) A3 = 4, 937⇥ 10

�4 � 2, 59⇥ 10

�5T + 9, 11⇥ 10

�7T 2 � 1, 5⇥ 10

�8T 3

T > 20

�C ) A3 = 3, 964⇥ 10

�4 � 1, 146⇥ 10

�5T + 1, 45⇥ 10

�7T 2 � 6, 5⇥ 10

�10T 3

(A.4)

66

Apêndice B

Tags comerciais

Esta secção tem a intenção de apresentar os modelos das tags disponíveis, dando uma ideia dos tamanhose dos sensores disponíveis.

B.1 Tags Tipo I

B.2 Tags Tipo II

67

Acelerómetros

Temperatura / Profundidade

68

B.3 Contínuos

B.4 Sinais gravados

69

70

71

72

Apêndice C

Código fonte da simulação em MatLab

C.1 start_simulation.m

1

c l e a r ;3 c l e a r a l l ;

c l o s e a l l ;5 c l c ;

7 v_sistem = s i s tem ( ) ;

9 % ju s t needed i f non standard value i s usedv_sistem = v_sistem . setVusblPrim ( [ 0 ; 0 ; 5 0 ] , [ 0 ; � ( p i /2) ; 0 ] ) ;

11

v_sistem = v_sistem . setVusblSec ( [ 1 0 0 ; 1 2 00 ; 1 50 ] , [ � ( p i /2) ; 0 ; 0 ] ) ;13

15 v_p = [ 3 0 ; 3 0 ; 3 0 ] ;v_p_2 = [ �700 ; 600 ; 200 ] ;

17

19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

21 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fish

23

% ava i l a b l e f i s h e s25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fish ID | Checksum | f i l e name %27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 21524 | 246 | ’21524 � 400ms�500mV. i s f ’ %29 % 33517 | 81 | ’33517 � 400ms�500mV. i s f ’ %

% 11894 | 221 | ’11894 � 400ms�500mV. i s f ’ %31 % 11905 | 169 | ’11905 � 400ms�50mV. i s f ’ %

% 33518 | 102 | ’33518 � 400ms�500mV. i s f ’ %33 % 33521 | 205 | ’33521 � 400ms�500mV. i s f ’ %

% 33522 | 209 | ’33522 � 400ms�500mV. i s f ’ %35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

37

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%39 % Moving mode = 1 � prede f i n ed s t ep s

41 v_sistem = v_sistem . addFish ( ’ 21524 � 400ms�500mV. i s f ’ , 21524 , v_p , 1) ; % moving modeprede f i n ed s t ep sv_sistem = v_sistem . r e g i s t e r F i s h (21524 , 246) ;

43

45 s t ep s = [30 , 30 .15 , 30 . 5 , 30 . 85 , 31 , 31 .15 , 31 . 5 , 31 . 85 , 32 , 32 .15 , 32 . 5 , 32 .8533 , 33 .15 , 33 . 5 , 33 . 85 , 34 , 34 . 15 , 34 . 5 , 34 . 85 , 35 , 35 , 35 , 35 , 35

, 35 , 35 , 35 , 35 , 35 , 35 ,35 ,35 ,35 ,35 ,35 ,35 ,35 ,35 ,35,35 , 3 5 ;

73

30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 3030 , 30 , 30 , 30 , 30 . 5 , 30 .85 , 31 , 31 .15 , 31 . 5 , 32 , 32 . 5 , 33 ,33 . 5 , 34 , 34 . 5 , 35 , 35 . 5 , 36 , 36 .5 ,37 , 37 . 5 , 38 ,38 ,38 ,38 ,38 ,38 ,38

,38 , 3 8 ;47 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30

30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30 , 30, 30 , 30 , 30 , 30 , 30 . 5 , 31 , 31 . 5 , 32

, 3 2 . 5 , 3 3 , 3 3 . 5 , 3 4 , 3 4 . 5 , 3 5 , 3 5 . 5 , 3 6 . 5 , 3 7 ] ;

49 s t ep s = s t ep s . ∗ 1 0 ;%s t ep s ( 3 , : ) = s t ep s ( 3 , : ) + 20 ;

51 %step s ( 2 , : ) = s t ep s ( 2 , : ) + 100 ;%s t ep s ( 1 , : ) = s t ep s ( 1 , : ) + 100 ;

53

55 v_sistem = v_sistem . i n i t F i s h (21524 , s t ep s ) ;

57

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%59 % Moving mode = 2 � by f i s h model

61 %v_sistem = v_sistem . addFish ( ’33517 � 400ms�500mV. i s f ’ , 33517 , v_p_2, 2) ; % movingmode by f i s h model%v_sistem = v_sistem . r e g i s t e r F i s h (33517 , 81) ;

63

% in i tF i s hS t a t e ( p_state , p_id )65 % p_state = [v_p_2 ; 1 ; 0 ; 0 ; 0 ; 0 ; 0 ] ;

% v_sistem = v_sistem . i n i tF i s hS t a t e ( p_state , 33517) ;67

% in i tF i s hNo i s e ( p_fish_noise_mag , p_id )69 % v_sistem = v_sistem . i n i tF i s hNo i s e (0 . 000002 , 33517) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%71

% setFishPingRepeat ( p_repeat_interval , p_repeat_std_dev , p_id )73 % v_sistem = v_sistem . setFishPingRepeat (45 , 5 , 33517)

75

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%77 %% System

79 % setSoundSpeed (p_vp)% Sound speed i s around 1484 [m/ s ] . But as we want to s imulate p r e c i s i o n

81 % of s i s tem with f s =250[kHz ] but are us ing f s =2.5 [kHz ] , soundspeed must% be s e t down because o f de lay r e s o l u t i o n .

83 v_sistem = v_sistem . setSoundSpeed ( 14 . 8 4 ) ;

85 % se tF i l t e r S p e c s ( p_usbl_nr , p_f0 , p_fs , p_Q)% v_sistem = v_sistem . s e tF i l t e r S p e c s (2 , 66500 , 2500 , 10) ;

87

% standard i s ShowSignalDetect ion Off89 % v_sistem = v_sistem . setShowSigna lDetect ionOf f ( ) ;

% v_sistem = v_sistem . setShowSignalDetect ionOn ( ) ;91

% standard i s a t t enuat ion on93 v_sistem = v_sistem . se tAttenuat ionOf f ( ) ;

% v_sistem = v_sistem . setAttenuationOn ( ) ;95

% v_sistem = v_sistem . setAttenuat ionTypeSpher ic ( ) ;97 % v_sistem = v_sistem . se tAttenuat ionTypeCyl indr i ca l (20) ; % argument i s channel he ight

99

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%101 %% USBL and Channel

103 % setUsb lSenso r s ( p_usbl_nr , p_sensors )% p_sensors = [ �0 .25 , �0 .25 ,0 . 25 ,0 . 25 ;

105 % 1 ,�1 ,0 ,0 ;% 0 ,0 ,1 ,�1 ] ;

107 % v_sistem = v_sistem . s e tUsb lSenso r s (2 , p_sensors ) ;

109 % setUsb lSenso r sNo i s e ( p_usbl_nr , p_sensors_noise )% p_sensors_noise = [ 0 . 0 1 , 0 . 01 , 0 . 01 , 0 . 0 1 ] ;

74

111 % v_sistem = v_sistem . se tUsb lSenso r sNo i s e (2 , p_sensors_noise ) ;

113 % setChannelNoise ( p_noise_level )v_sistem = v_sistem . setChannelNoise ( 0 . 0 01 ) ;

115

117

119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Kalman F i l t e r

121

% initSystemModelNoise ( p_sys_noise ) 0 . 2123 v_sistem = v_sistem . initSystemModelNoise ( 0 . 0002 ) ;

125 % initMeasurementNoise ( p_meas_noise ) 0 .01v_sistem = v_sistem . initMeasurementNoise ( 0 . 2 ) ;

127

129

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%131 %%

133 % f i r s t p o s i t i o n s o f f i s h e sv_sistem = v_sistem . i n i tEven tL i s t ( ) ;

135

f o r s tep = 1 : 30137

s tep139

[ v_sistem , v_event ] = v_sistem . de t e c t ( ) ;141

143 end ;

145

147

v_sistem . show ( ) ;149 %v_sistem . p lotUsb l ( ) ;

v_sistem . p r i n t_ f i l e ( ) ;

./matlab/start_simulation.m

C.2 sistem.m

1 c l a s s d e f s i s t emp r op e r t i e s ( GetAccess = ’ p r i va t e ’ , SetAccess = ’ p r i va t e ’ )

3 s t a tu s % s ta tu s o f the systemtime % time c l o ck

5 vusbl_prim % primary usb l ob j e c tvusbl_sec % secondary usb l ob j e c t

7 no i s e % no i s e in the channel � t y p i c a l ( 0 . 1 )f i sh_vec % vecto r o f e x i s t i n g f i s h e s

9 f i s h_ id s % vecto r with i d s o f e x i s t i n g f i s h e sf ish_checksums % r e g i s t e r e d checksums

11 fish_model_vec % vec to r o f e x i s t i n g f i s h modelsf ish_model_ids % vecto r with i d s o f e x i s t i n g f i s h models

13 de t e c t i on s % vecto r conta in ing de t e c t i on ob j e c t s f o r each detec ted eventc o l l i s i o n s % vec to r conta in ing event ob j e c t s that gave no de t e c t i on r e s u l t dueto s i g n a l c o l l i s i o n

15 no_detect ions % vec to r conta in ing event ob j e c t s that gave no de t e c t i on r e s u l t

17 even t_ l i s t % l i s t o f events to be proce s s ed

19 pos % standard : [ 0 ; 0 ; 0 ]%o r i e n t a t i o n % standard : [ 0 , 0 , 0 ]

21

channel23

75

meas_noise25 sys_noise

27 dens i ty % constant o f water dens i ty used f o r measurement when j u s t primary USBLe x i s t sg rav i ty % constant o f earth g rav i ty used f o r measurement when j u s t primary USBLe x i s t s

29 vp % constant f o r propagat ion v e l o c i t y o f sound in water

31 % fo r p r i n t to f i l e , used in f � p r i n t_ f i l e ( )f i l e_heade r

33 f i l e_outputf i l e_format_header

35

exec_no % execut ion number , count ing up37

% 0 � at tenuat ion i s not cons ide r ed39 % 1 � at tenuat ion i s cons ide r ed

attenuation_on41

% 1 � at tenuat ion i s cons ide r ed s ph e r i c a l43 % 2 � at tenuat ion i s cons ide r ed c y l i n d r i c a l

attenuat ion_type45

% used f o r c y l i n d r i c a l a t t enuat ion47 channel_height

49 % parameters f o r a t t enuat ion c a l c u l a t i o ns a l i n i t y

51 temperatureatt_frequency

53

f s % sampling f requency55

end % pr op e r t i e s57

methods59

% Constructor , Getter s and S e t t e r s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%61

f unc t i on obj = s i s tem ( ) % cons t ruc to r63 obj . pos = [ 0 ; 0 ; 0 ] ;

%obj . o r i e n t a t i o n = [ 0 ; 0 ; 1 ] ;65 obj . vusbl_prim = usb l ( [ 0 ; 0 ; 0 ] , [ 0 ; � ( p i /2) ; 0 ] ) ; % i n i t usb l ob j e c t ( standard )

%obj . vusbl_prim = usb l ( [ 0 ; 0 ; 0 ] , [ 0 ; 0 ; 0 ] ) ; % i n i t usb l ob j e c t ( standard )67

obj . ev en t_ l i s t = eventL i s t ( ) ;69

% i n i t i a l i z e channel71 obj . channel = channel ( ) ;

in i t_ar ray = ze ro s (4 , 10000) ;73 v_noisy_signal = obj . channel . addNoiseToSignal ( in i t_ar ray ) ;

75 [ obj . vusbl_prim , v_threshold ] = obj . vusbl_prim . i n i t F i l t e r ( v_noisy_signal ) ;

77 obj . f i sh_vec = [ ] ; % vec to r o f e x i s t i n g f i s h e sobj . f i s h_ id s = [ ] ; % vec to r with i d s o f e x i s t i n g f i s h e s

79 obj . f ish_checksums = [ ] ; % vec to r with checksums o f r e g i s t e r e d f i s h e sobj . fish_model_vec = [ ] ; % vecto r o f e x i s t i n g f i s h models

81 obj . f ish_model_ids = [ ] ; % vec to r with i d s o f e x i s t i n g f i s h modelsobj . d e t e c t i on s = [ ] ; % vec to r with d e t e c t i on s f o r detec ted events

83

obj . meas_noise = 0 . 0 1 ;85 obj . sys_noise = 0 . 2 ;

87 obj . dens i ty = 1030 ; % kgm^�3obj . g rav i ty = 9 . 8 1 ; % m/ s^2

89 %obj . vp = 1450 ; % m/ sobj . vp = 1484 ; % m/ s^2

91

obj . vusbl_prim = obj . vusbl_prim . setSoundSpeed ( obj . vp ) ;93 obj . f i l e_heade r = [ ] ;

76

obj . f i l e_output = ’ ’ ;95 obj . f i l e_format_header = ’ ’ ;

97 obj . exec_no = 0 ;

99 obj . attenuation_on = 1 ; % standard , c on s id e r a t t enuat ionobj . attenuation_type = 1 ; % standard , c on s id e r s p h e r i c a l a t t enuat ion

101

% fo r at t enuat ion c a l c u l a t i o n103 obj . channel_height = 30 ; % f o r c y l i n d r i c a l decay law

obj . s a l i n i t y = 35 ; % in [ p . s . u . ]105 obj . temperature = 15 ; % in degree [C]

obj . att_frequency = 69 % in [ kHz ]107

end % cons t ruc to r109

111

% ge t t e r s and s e t t e r s113

f unc t i on s t a tu s = getStatus ( obj )115 s t a tu s = obj . s t a tu s ;

end % f � getStatus117

f unc t i on s e tS ta tu s ( obj , s t a tu s )119 obj . s t a tu s = s ta tu s ;

end % f � s e tS ta tu s121

f unc t i on time = getTime ( obj )123 time = obj . time ;

end % f � getTime125

f unc t i on setTime ( obj , time )127 obj . time = time ;

end % f � setTime129

f unc t i on usb l = getVusblPrim ( obj )131 usb l = obj . vusbl_prim ;

end % f � getVusblPrim133

f unc t i on obj = setVusblPrim ( obj , p , o r i en ta t i on_ang l e s )135 obj . vusbl_prim = usb l (p , o r i en ta t i on_ang l e s ) ;

137 obj . vusbl_prim = obj . vusbl_prim . setSoundSpeed ( obj . vp ) ;

139 v_noisy_signal = obj . channel . addNoiseToSignal ( z e r o s (4 , 10000) ) ;[ obj . vusbl_prim , v_threshold ] = obj . vusbl_prim . i n i t F i l t e r ( v_noisy_signal ) ;

141

end % f � setVusblPrim143

145 f unc t i on usb l = getVusblSec ( obj )usb l = obj . vusbl_sec ;

147 end % f � getVusblSec

149 f unc t i on obj = setVusblSec ( obj , p , o r i en ta t i on_ang l e s )obj . vusbl_sec = usb l (p , o r i en ta t i on_ang l e s ) ;

151

obj . vusbl_sec = obj . vusbl_sec . setSoundSpeed ( obj . vp ) ;153

v_noisy_signal = obj . channel . addNoiseToSignal ( z e r o s (4 , 10000) ) ;155 [ obj . vusbl_sec , v_threshold ] = obj . vusbl_sec . i n i t F i l t e r ( v_noisy_signal ) ;

157 end % f � setVusblSec

159 f unc t i on f i s h = getFi sh ( obj , index )i f index <= obj . getNumFishes && index > 0

161 f i s h = obj . f i sh_vec ( index ) ;e l s e

163 e r r o r ( ’ Just <<obj . getNumFishes ( )>> Fishes e x i s t . ’ ) ;end

165 end % f � getF i sh

77

167 f unc t i on num = getNumFishes ( obj )[m, n ] = s i z e ( obj . f i sh_vec ) ;

169 num = n ;end % f � getF i sh

171

f unc t i on num = getNumFishModels ( obj )173 num = length ( obj . f ish_model_ids ) ;

end % f � getNumFisheModels175

f unc t i on dets = ge tDetec t i on s ( obj )177 dets = obj . d e t e c t i on s ;

end % f � getVusblSec179

181 f unc t i on obj = setAttenuationOn ( obj )obj . attenuation_on = 1 ;

183

obj . vusbl_prim = obj . vusbl_prim . se tAttenuat ion (1 ) ;185

% Re � i n i t i a l i z e USBL f i l t e r187 i n i t_ar ray = ze ro s (4 , 10000) ;

v_noisy_signal = obj . channel . addNoiseToSignal ( in i t_ar ray ) ;189

[ obj . vusbl_prim , v_threshold ] = obj . vusbl_prim . i n i t F i l t e r ( v_noisy_signal ) ;191

i f i s o b j e c t ( obj . vusbl_sec ) % i f secondary USBL i s loaded193 obj . vusbl_sec = obj . vusbl_sec . s e tAttenuat ion (1 ) ;

% Re � i n i t i a l i z e USBL f i l t e r195 [ obj . vusbl_sec , v_threshold ] = obj . vusbl_sec . i n i t F i l t e r ( v_noisy_signal ) ;

197 end

199 end % f � setAttenuationOn

201 f unc t i on obj = setAttenuat ionOf f ( obj )obj . attenuation_on = 0 ;

203

obj . vusbl_prim = obj . vusbl_prim . se tAttenuat ion (0 ) ;205

% Re � i n i t i a l i z e USBL f i l t e r207 i n i t_ar ray = ze ro s (4 , 10000) ;

v_noisy_signal = obj . channel . addNoiseToSignal ( in i t_ar ray ) ;209

[ obj . vusbl_prim , v_threshold ] = obj . vusbl_prim . i n i t F i l t e r ( v_noisy_signal ) ;211

213 i f i s o b j e c t ( obj . vusbl_sec ) % i f secondary USBL i s loadedobj . vusbl_sec = obj . vusbl_sec . s e tAttenuat ion (0 ) ;

215 [ obj . vusbl_sec , v_threshold ] = obj . vusbl_sec . i n i t F i l t e r ( v_noisy_signal ) ;

217 end

219 end % f � s e tAttenuat ionOf f

221

f unc t i on obj = setAttenuat ionTypeSpher ic ( obj )223

obj . channel = obj . channel . setAttenuationType (1 ) ;225

end % f � setAttenuat ionTypeSpher ic227

229

f unc t i on obj = setAttenuat ionTypeCyl indr i ca l ( obj , p_channel_height )231

obj . channel = obj . channel . setAttenuationType (2 ) ;233 obj . channel_height = p_channel_height ;

235 end % f � s e tAttenuat ionTypeCyl indr i ca l

237

78

f unc t i on obj = setShowSignalDetect ionOn ( obj )239

h1 = f i g u r e ( ) ;241 obj . vusbl_prim = obj . vusbl_prim . setShowSigna lDetect ion (1 , h1 ) ;

243 i f i s o b j e c t ( obj . vusbl_sec ) % i f secondary USBL i s loadedh2 = f i g u r e ( ) ;

245 obj . vusbl_sec = obj . vusbl_sec . setShowSigna lDetect ion (1 , h2 ) ;end

247

end % f � setShowSignalDetect ionOn249

f unc t i on obj = setShowSigna lDetect ionOf f ( obj )251 obj . vusbl_prim = obj . vusbl_prim . setShowSigna lDetect ion (0 ) ;

253 i f i s o b j e c t ( obj . vusbl_sec ) % i f secondary USBL i s loadedobj . vusbl_sec = obj . vusbl_sec . setShowSigna lDetect ion (0 ) ;

255 end

257 end % f � se tShowSigna lDetect ionOf f

259

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%261

f unc t i on obj = se tUsb lSenso r s ( obj , p_usbl_nr , p_sensors )263 [m, n ] = s i z e ( p_sensors )

i f or (m~=3,n~=4)265 e r r o r ( ’ s enso r matrix must be o f dimensions 3x4 . ’ ) ;

e l s e267 i f p_usbl_nr == 1

obj . vusbl_prim = obj . vusbl_prim . s e tSenso r1 ( p_sensors ( : , 1 ) ) ;269 obj . vusbl_prim = obj . vusbl_prim . s e tSenso r2 ( p_sensors ( : , 2 ) ) ;

obj . vusbl_prim = obj . vusbl_prim . s e tSenso r3 ( p_sensors ( : , 3 ) ) ;271 obj . vusbl_prim = obj . vusbl_prim . s e tSenso r4 ( p_sensors ( : , 4 ) ) ;

e l s e273 obj . vusbl_sec = obj . vusbl_sec . s e tSenso r1 ( p_sensors ( : , 1 ) ) ;

obj . vusbl_sec = obj . vusbl_sec . s e tSenso r2 ( p_sensors ( : , 2 ) ) ;275 obj . vusbl_sec = obj . vusbl_sec . s e tSenso r3 ( p_sensors ( : , 3 ) ) ;

obj . vusbl_sec = obj . vusbl_sec . s e tSenso r4 ( p_sensors ( : , 4 ) ) ;277 end

end279 end % f � s e tUsb lSenso r s

281

f unc t i on obj = se tUsb lSenso r sNo i s e ( obj , p_usbl_nr , p_sensors_noise )283 [m, n ] = s i z e ( p_sensors_noise )

i f or (m~=1,n~=4)285 e r r o r ( ’ s enso r matrix must be o f dimensions 1x4 . ’ ) ;

e l s e287 i f p_usbl_nr == 1

obj . vusbl_prim = obj . vusbl_prim . se tNo i s eSensor1 ( p_sensors_noise (1 , 1 ) ) ;289 obj . vusbl_prim = obj . vusbl_prim . se tNo i s eSensor2 ( p_sensors_noise (1 , 2 ) ) ;

obj . vusbl_prim = obj . vusbl_prim . se tNo i s eSensor3 ( p_sensors_noise (1 , 3 ) ) ;291 obj . vusbl_prim = obj . vusbl_prim . se tNo i s eSensor4 ( p_sensors_noise (1 , 4 ) ) ;

e l s e293 obj . vusbl_sec = obj . vusbl_sec . s e tNo i s eSensor1 ( p_sensors_noise (1 , 1 ) ) ;

obj . vusbl_sec = obj . vusbl_sec . s e tNo i s eSensor2 ( p_sensors_noise (1 , 2 ) ) ;295 obj . vusbl_sec = obj . vusbl_sec . s e tNo i s eSensor3 ( p_sensors_noise (1 , 3 ) ) ;

obj . vusbl_sec = obj . vusbl_sec . s e tNo i s eSensor4 ( p_sensors_noise (1 , 4 ) ) ;297 end

end299 end % f � s e tUsb lSenso r sNo i s e

301 f unc t i on obj = setChannelNoise ( obj , p_noise_level )

303 obj . channel = obj . channel . s e tNo i s eLeve l ( p_noise_level ) ;

305 end % f � setChannelNoise

307

f unc t i on obj = i n i tF i s hS t a t e ( obj , p_state , p_id )309

79

v_idx = obj . f i ndF i sh ( p_id ) ;311 v_fish = obj . f i sh_vec ( v_idx ) ;

v_fish = v_fish . i n i t S t a t e ( p_state ) ;313 obj . f i sh_vec (1 , v_idx ) = v_fish ;

315 end % f � i n i tF i s hS t a t e

317 f unc t i on obj = in i tF i s hNo i s e ( obj , p_fish_noise_mag , p_id )v_idx = obj . f i ndF i sh ( p_id ) ;

319 v_fish = obj . f i sh_vec (1 , v_idx ) ;v_fish = v_fish . i n i tF i s hNo i s e ( p_fish_noise_mag ) ;

321 obj . f i sh_vec (1 , v_idx ) = v_fish ;end % f � i n i tF i s hNo i s e

323

f unc t i on obj = setFishPingRepeat ( obj , p_repeat_interval , p_repeat_std_dev , p_id )325 v_idx = obj . f i ndF i sh ( p_id ) ;

v_fish = obj . f i sh_vec ( v_idx ) ;327 v_fish = v_fish . s e tRepea t In t e rva l ( p_repeat_interval ) ;

v_fish = v_fish . setRepeatStdDev ( p_repeat_std_dev ) ;329 obj . f i sh_vec ( v_idx ) = v_fish ;

end % f � setFishPingRepeat331

333 f unc t i on obj = setSoundSpeed ( obj , p_vp)obj . vusbl_prim = obj . vusbl_prim . setSoundSpeed (p_vp) ;

335 i f i s o b j e c t ( obj . vusbl_sec )obj . vusbl_sec = obj . vusbl_sec . setSoundSpeed (p_vp) ;

337 endend % f � setSoundSpeed

339

341 f unc t i on obj = s e tF i l t e r S p e c s ( obj , p_usbl_nr , p_f0 , p_fs , p_Q)i f p_usbl_nr == 1

343 obj . vusbl_prim = obj . vusbl_prim . s e tF i l t e r S p e c s ( p_f0 , p_fs , p_Q) ;e l s e

345 obj . vusbl_sec = obj . vusbl_sec . s e tF i l t e r S p e c s ( p_f0 , p_fs , p_Q) ;end

347 end % f � s e tF i l t e r S p e c s

349 f unc t i on obj = initSystemModelNoise ( obj , p_sys_noise )obj . sys_noise = p_sys_noise ;

351 end % f � in itSystemModelNoise

353 f unc t i on obj = initMeasurementNoise ( obj , p_meas_noise )obj . meas_noise = p_meas_noise ;

355 end % f � in itMeasurementNoise

357

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%359

% read a f i l e , ex : 21524 � ’21524 � 400ms�500mV. i s f ’361 f unc t i on obj = in i tEven tL i s t ( obj )

363 [m, n ] = s i z e ( obj . f i sh_vec ) ;f o r i = 1 : n

365

% f i r s t i n i t i a l i z a t i o n , time i s 0367 [ v_fish , v_event ] = obj . f i sh_vec ( i ) . getNextEvent (0 ) ;

obj . f i sh_vec ( i ) = v_fish ;369

obj . ev en t_ l i s t = obj . ev en t_ l i s t . addEvent ( v_event ) ;371

end373

end % f � i n i tEven tL i s t375

377 % f i nd s index in f i s h vec to r f o r f i s h i d e n t i f i e d by idfunc t i on idx = f indF i sh ( obj , id )

379 [m, n ] = s i z e ( obj . f i sh_vec ) ;f o r i = 1 : n ;

381 v_fish = obj . f i sh_vec ( i ) ;

80

i f v_fish . ge t Id ( ) == id383 idx = i ;

r e turn385 end

end387

% case o f f i s h id not found in f i s h vec to r389 e r r o r ( ’ Fish ID does not e x i s t in system . ’ ) ;

%idx = �1;391 end % f � f i ndF i sh

393

% i n i t s f i s h f o r movement type 1 with prede f i n ed p o s i t i o n s395 f unc t i on obj = i n i t F i s h ( obj , id , s t ep s )

397 idx = obj . f i ndF i sh ( id ) ;

399 obj . f i sh_vec ( idx ) = obj . f i sh_vec ( idx ) . i n i t S t e p s ( s t ep s ) ;

401 end % f � i n i t F i s h

403

% read a f i l e , ex : 21524 � ’21524 � 400ms�500mV. i s f ’405 f unc t i on s i g n a l = openFi l e ( obj , f i le_name )

s i g = dlmread ( file_name , ’ , ’ ) ;407 [m, n ] = s i z e ( s i g ) ;

s i g n a l = s i g ( 1 :m, 2 ) ;409 s i g n a l = transpose ( s i g n a l ) ;

411 end % f � openFi l e

413 % Add a f i s h to the system . A f i s h i s i d e n t i f i e d by a code read from a% f i l e , ex : 21524 � ’21524 � 400ms�500mV. i s f ’

415 f unc t i on obj = addFish ( obj , fi le_name , id , pos , v_moving_mode)s i g n a l = obj . openFi l e ( f i le_name ) ;

417 % moving mode , l a s t argument o f f i s h con s t ruc to r% 1 � prede f i n ed p o s i t i o n s

419 % 2 � movement modelobj . f i sh_vec = cat (2 , obj . f i sh_vec , f i s h ( id , s i gna l , pos , 45 , 5 , v_moving_mode ,

obj . dens i ty , obj . g r av i ty ) ) ;421 obj . f i s h_ id s = [ obj . f i sh_ids , id ] ;

423 end % f � getUsbl

425

f unc t i on obj = r e g i s t e r F i s h ( obj , p_id , p_checksum)427

v_col = [ p_id ; p_checksum ] ;429 obj . f ish_checksums = [ obj . f ish_checksums , v_col ] ;

431 end % f � getUsbl

433 f unc t i on r e s = checkChecksum ( obj , p_id , p_checksum)r e s = 0 ;

435 [m, n ] = s i z e ( obj . f ish_checksums ) ;f o r i = 1 : n

437 i f obj . f ish_checksums (1 , i ) == p_idi f obj . f ish_checksums (2 , i ) == p_checksum

439 r e s = 1 ;end

441

r e turn ;443 end

end445

end % f � getUsbl447

449 % other f un c t i on s

451 % crea t e a s i g n a l o f input s i gna l_ in on a l l f our channe l s with de lay ( de lay arraydepending on po s i t i o n o f f i s h )

81

f unc t i on s ignal_out = crea t eP ings ( obj , s igna l_in , de lay )453 s ignal_out = ze ro s ( l ength ( de lay ) , l ength ( s i gna l_in ) ) ;

455 end % f � getUSBL

457

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%459 % ca l c u l a t e po int out o f two usb l d i r e c t i o n measurements

% in : usbl_posit ion_1 , direct ion_1 , usbl_posit ion_2 , d i rec t ion_2461 % Points must be in same coord inate�system

% out : Point in space463 % Two c l o s e s t po in t s on the l i n e s are c a l c u l a t ed . Output i s the po int

% inbetween them .465 f unc t i on [ point , p r e c i s i o n ] = measureTwoUsbl ( obj , p_1 , v_1 , p_2 , v_2)

467 w = (p_1 � p_2) ;a = dot (v_1 , v_1) ;

469 b = dot (v_1 , v_2) ;c = dot (v_2 , v_2) ;

471 d = dot (v_1 , w) ;e = dot (v_2 , w) ;

473

i f ( ( a∗c )�(b∗b) ) > 0475 s = ( ( b∗e )�(c∗d) ) / ( ( a∗c )�(b∗b) ) ;

t = ( ( a∗e )�(b∗d) ) / ( ( a∗c )�(b∗b) ) ;477

%s = (�d + t ∗b) /a479

point_1 = p_1 + s ∗ v_1 ;481 point_2 = p_2 + t ∗ v_2 ;

483 e l s epoint_1 = p_1 ;

485 point_2 = p_2 ;end

487

point = point_2 + ( ( point_1 � point_2 ) / 2) ;489

p r e c i s i o n = norm ( ( point_1 � point_2 ) ) ; % The lower the number , the h igher thep r e c i s i o n

491

end % f � measureTwoUsbl493

495

f unc t i on [ pos , p r e s s i s i o n ] = measureOneUsbl ( obj , p_pos , p_ud , p_pressure )497

% ca l c u l a t e a l l in g l oba l c oo rd ina t e s499 v_p0z = p_pressure / ( obj . dens i ty ∗ obj . g rav i ty ) ;

501 v_p0x = p_pos (1 )+(v_p0z�p_pos (3 ) ) /p_ud(3) ∗p_ud(1) ;v_p0y = p_pos (2 )+(v_p0z�p_pos (3 ) ) /p_ud(3) ∗p_ud(2) ;

503

pos = [ v_p0x ; v_p0y ; v_p0z ] ;505

p r e s s i s i o n = 0 ;507 end % f � measureOneUsbl

509

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%511 % main rou t ine

func t i on [ obj , v_event ] = det ec t ( obj )513

% i n i t l o c a l v a r i a b l e s515

v_id = NaN;517 v_checksum = NaN;

v_detected = 0 ;519

521 obj . exec_no = obj . exec_no + 1 ;

82

523 v_attenuation_1 = 0 ;v_attenuation_2 = 0 ;

525 %%%%%%%%%%%%%%%%

527 [ obj . event_l i s t , v_event ] = obj . ev en t_ l i s t . getNextEvent ( ) ;

529 % current eventsv_times = v_event . getTimes ( ) ;

531 v_f i shes = v_event . g e tF i she s ( ) ;

533

%%%%%%%%%%%%%%%535 % USBL_1

537 v_noise_s1 = obj . vusbl_prim . getNoi seSensor1 ( ) ;obj . channel = obj . channel . s e tNo i s eSensor1 ( v_noise_s1 ) ;

539

v_noise_s2 = obj . vusbl_prim . getNoi seSensor2 ( ) ;541 obj . channel = obj . channel . s e tNo i s eSensor2 ( v_noise_s2 ) ;

543 v_noise_s3 = obj . vusbl_prim . getNoi seSensor3 ( ) ;obj . channel = obj . channel . s e tNo i s eSensor3 ( v_noise_s3 ) ;

545

v_noise_s4 = obj . vusbl_prim . getNoi seSensor4 ( ) ;547 obj . channel = obj . channel . s e tNo i s eSensor4 ( v_noise_s4 ) ;

549

i f obj . attenuation_on551 v_usbl_pos = obj . vusbl_prim . getPos ( ) ;

v_depth = ( v_f i shes (1 ) . getDepth ( ) + v_usbl_pos (3 ) ) /2 ; %depth i s mean betweenf i s h e s depth and usb l s depth

553 [ obj . channel , v_attenuation_1 ] = obj . channel . s e tAttenuat ion ( v_depth , obj .att_frequency , obj . s a l i n i t y , obj . temperature , v_f i shes (1 ) . getPos ( ) , obj . vusbl_prim .getPos ( ) , obj . channel_height ) ;

555 end

557 v_noisy_signal = obj . c r e a t eNo i syS i gna l ( obj . vusbl_prim , v_event ) ;[ v_dets , v_angles1 , v_d1 ] = obj . vusbl_prim . de t e c t ( v_noisy_signal ) ;

559 [ v_ids_1 , v_checksums_1 ] = obj . vusbl_prim . ca lcFi shID ( v_dets ) ;

561

v_pos_meas = [NaN, NaN, NaN ] ; % to prevent e r r o r when no de t e c t i on i s a v a i l a b l e563

i f i s o b j e c t ( obj . vusbl_sec ) % i f secondary USBL i s loaded565

567 %%%%%%%%%%%%%%%% USBL_2

569

v_noise_s1 = obj . vusbl_sec . getNo i seSensor1 ( ) ;571 obj . channel = obj . channel . s e tNo i s eSensor1 ( v_noise_s1 ) ;

573 v_noise_s2 = obj . vusbl_sec . getNo i seSensor2 ( ) ;obj . channel = obj . channel . s e tNo i s eSensor2 ( v_noise_s2 ) ;

575

v_noise_s3 = obj . vusbl_sec . getNo i seSensor3 ( ) ;577 obj . channel = obj . channel . s e tNo i s eSensor3 ( v_noise_s3 ) ;

579 v_noise_s4 = obj . vusbl_sec . getNo i seSensor4 ( ) ;obj . channel = obj . channel . s e tNo i s eSensor4 ( v_noise_s4 ) ;

581

583 i f obj . attenuation_onv_usbl_pos = obj . vusbl_sec . getPos ( ) ;

585 v_depth = ( v_f i shes (1 ) . getDepth ( ) + v_usbl_pos (3 ) ) /2 ; %depth i s meanbetween f i s h e s depth and usb l s depth

[ obj . channel , v_attenuation_2 ] = obj . channel . s e tAttenuat ion ( v_depth , 69 ,35 , 15 , v_f i shes (1 ) . getPos ( ) , obj . vusbl_sec . getPos ( ) , obj . channel_height ) ;

587 end

589

83

591 v_noisy_signal = obj . c r e a t eNo i syS i gna l ( obj . vusbl_sec , v_event ) ;[ v_dets , v_angles2 , v_d2 ] = obj . vusbl_sec . de t e c t ( v_noisy_signal ) ;

593 [ v_ids_2 , v_checksums_2 ] = obj . vusbl_sec . ca lcFi shID ( v_dets ) ;

595 % In case o f d e t e c t i on � no c o l i s i o ni f ~or ( isempty ( v_ids_1 ) , isempty ( v_ids_2 ) )

597

id_tmp = v_ids_1 ( 1 , : ) == v_ids_2 ( 1 , : ) ;599

% point o f most aproximation to both usb l d i r e c t i o n l i n e s601 [ v_pos_meas , v_prec i s ion ] = obj . measureTwoUsbl ( obj . vusbl_prim . getPos ( ) ,

v_d1 , obj . vusbl_sec . getPos ( ) , v_d2) ;

603

605 f o r j = 1 : min ( l ength ( v_ids_1 ) , l ength ( v_ids_2 ) )i f and ( id_tmp( j ) == 1 , ( v_ids_1 (1 , j ) ~= 0) )

607

v_id = v_ids_1 (1 , j ) ;609 v_checksum = v_checksums_1 (1 , j ) ;

611 v_detected = 1 ; % a f i s h was detec ted

613 endend

615 end

617

619 e l s e % i f j u s t primary USBL i s loaded

621 % In case o f d e t e c t i on � no c o l i s i o ni f ~isempty ( v_ids_1 )

623

v_pressure = v_event . ge tPre s su re ( ) ;625

[ v_pos_meas , v_prec i s ion ] = obj . measureOneUsbl ( obj . vusbl_prim . getPos ( ) ,v_d1 , v_pressure ) ;

627

f o r j = 1 : l ength ( v_ids_1 )629 i f ( v_ids_1 (1 , j ) ~= 0)

631 v_id = v_ids_1 (1 , j ) ;v_checksum = v_checksums_1 (1 , j ) ;

633

v_detected = 1 ; % a f i s h was detec ted635

end637 end

end639

end641

% i f a f i s h was detec ted643 i f a l l ( [ v_detected == 1 ,~ i snan (v_pos_meas (1 ) ) , ~ i s i n f (v_pos_meas (1 ) ) , ~ i snan

(v_pos_meas (2 ) ) , ~ i s i n f (v_pos_meas (2 ) ) , ~ i snan (v_pos_meas (3 ) ) , ~ i s i n f (v_pos_meas (3 )) ] )

% in case o f c o l i s i o n , time o f c o l i s i o n event event i s time645 % of f i r s t event

647 [ real_ang_1 , rea l_loca l_d i rec t i on_1 ] = obj . vusbl_prim .getDirectionFromPosGeo ( v_event . g e tPo s i t i o n s ( ) ) ;

real_glob_dir_1 = obj . vusbl_prim . l o ca lToGloba lD i r e c t i on (rea l_loca l_d i rec t i on_1 ) ;

649

i f i s o b j e c t ( obj . vusbl_sec )651 [ real_ang_2 , rea l_loca l_d i rec t i on_2 ] = obj . vusbl_sec .

getDirectionFromPosGeo ( v_event . g e tPo s i t i o n s ( ) ) ;real_glob_dir_2 = obj . vusbl_sec . l o ca lToGloba lD i r e c t i on (

rea l_loca l_d i rec t i on_2 ) ;653 end

84

655 % compare detec ted checksum with r e g i s t e r e di f obj . checkChecksum (v_id , v_checksum) == 1

657

[ obj , v_pos ] = obj . updateFishModel ( v_id , v_pos_meas , v_times (1 ) ,v_prec i s ion ) ;

659

% save de t e c t i on ob j e c t661 obj . d e t e c t i on s = [ obj . de t e c t i on s , d e t e c t i on ( v_event , v_pos_meas ,

v_prec is ion , v_pos , v_id ) ] ;

663 v_real_pos = v_f i shes (1 ) . getPos ( ) ;

665 % crea t e new l i n e f o r output to f i l ei f i s o b j e c t ( obj . vusbl_sec )

667 obj . f i l e_output = [ obj . f i l e_output , num2str ( obj . exec_no ) , ’ |’ , num2str ( v_times (1 ) ) , ’ | ’ , num2str ( v_id ) , ’ | ’ , num2str (v_pos_meas (1 ) ) , ’ , ’ ,num2str (v_pos_meas (2 ) ) , ’ , ’ , num2str (v_pos_meas (3 ) ) , ’ | ’ , num2str ( v_real_pos (1 ) ) , ’ , ’, num2str ( v_real_pos (2 ) ) , ’ , ’ , num2str ( v_real_pos (3 ) ) , ’ | ’ , num2str (v_d1(1) ) , ’ , ’ ,num2str (v_d1(2) ) , ’ , ’ , num2str (v_d1(3) ) , ’ | ’ , num2str ( real_glob_dir_1 (1) ) , ’ , ’ ,num2str ( real_glob_dir_1 (2) ) , ’ , ’ , num2str ( real_glob_dir_1 (3 ) ) , ’ | ’ , num2str (v_angles1 (1 ) ) , ’ , ’ , num2str ( v_angles1 (2 ) ) , ’ | ’ , num2str ( real_ang_1 (1) ) , ’ , ’ , num2str (real_ang_1 (2) ) , ’ | ’ , num2str ( v_attenuation_1 ) , ’ | ’ , num2str (v_d2(1 ) ) , ’ , ’ , num2str (v_d2(2) ) , ’ , ’ , num2str (v_d2(3) ) , ’ | ’ , num2str ( real_glob_dir_2 (1 ) ) , ’ , ’ , num2str (real_glob_dir_2 (2) ) , ’ , ’ , num2str ( real_glob_dir_2 (3 ) ) , ’ | ’ , num2str ( v_angles2 (1 ) ) , ’ ,

’ , num2str ( v_angles2 (2 ) ) , ’ | ’ , num2str ( real_ang_2 (1) ) , ’ , ’ , num2str ( real_ang_2 (2) ) , ’| ’ , num2str ( v_attenuation_2 ) , ’ | ’ , ’ \n ’ ] ;

e l s e669 obj . f i l e_output = [ obj . f i l e_output , num2str ( obj . exec_no ) , ’ |

’ , num2str ( v_times (1 ) ) , ’ | ’ , num2str ( v_id ) , ’ | ’ , num2str (v_pos_meas (1 ) ) , ’ , ’ ,num2str (v_pos_meas (2 ) ) , ’ , ’ , num2str (v_pos_meas (3 ) ) , ’ | ’ , num2str ( v_real_pos (1 ) ) , ’ , ’, num2str ( v_real_pos (2 ) ) , ’ , ’ , num2str ( v_real_pos (3 ) ) , ’ | ’ , num2str (v_d1(1) ) , ’ , ’ ,num2str (v_d1(2) ) , ’ , ’ , num2str (v_d1(3) ) , ’ | ’ , num2str ( real_glob_dir_1 (1) ) , ’ , ’ ,num2str ( real_glob_dir_1 (2) ) , ’ , ’ , num2str ( real_glob_dir_1 (3 ) ) , ’ | ’ , num2str (v_angles1 (1 ) ) , ’ , ’ , num2str ( v_angles1 (2 ) ) , ’ | ’ , num2str ( real_ang_1 (1) ) , ’ , ’ , num2str (real_ang_1 (2) ) , ’ | ’ , num2str ( v_attenuation_1 ) , ’ | ’ , ’ \n ’ ] ;

671 end

673 end

675 e l s e

677 i f v_event . getType ( ) == 1 % c o l l i s i o nobj . c o l l i s i o n s = [ obj . c o l l i s i o n s , v_event . g e tPo s i t i o n s ( ) ] ;

679 obj . f i l e_output = [ obj . f i l e_output , ’SIGNAL COLLISION \n ’ ] ;e l s e

681 obj . no_detect ions = [ obj . no_detect ions , v_event . g e tPo s i t i o n s ( ) ] ;obj . f i l e_output = [ obj . f i l e_output , ’NO DETECTION \n ’ ] ;

683 end

685 end

687 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

689

691 f o r i = 1 : l ength ( v_f i shes )

693 % next events[ v_f i shes ( i ) , v_event ] = v_f i shes ( i ) . getNextEvent ( v_times ( i ) ) ;

695 obj . ev en t_ l i s t = obj . ev en t_ l i s t . addEvent ( v_event ) ;

697 % obj . ev en t_ l i s t . p r i n t ( ) ;

699

obj = obj . updateFish ( v_f i shes ( i ) ) ;701

end703

end % f � det e c t705

85

707 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

709 f unc t i on [ obj , v_pos ] = updateFishModel ( obj , p_id , p_pos_meas , p_t , p_prec i s ion )

711 v_pos = p_pos_meas ;

713

[ obj , v_fish_model , v_idx ] = obj . getFishModel ( p_id , p_pos_meas ) ;715

v_t = v_fish_model . getTime ( ) ;717 v_dt = p_t � v_t ;

719

f o r t_step = 1 : 1 0 : ( f l o o r (v_dt ) ) % in s t ead o f us ing p r ed i c t type events , j u s tp r ed i c t i n g a c e r t a i n number o f t imes

721 [ v_fish_model , v_pos ] = v_fish_model . p r ed i c tNextPos i t i on (10) ; % TODO: Achtung, 10 a l s argument ko r r e c t ?

end723

v_fish_model = v_fish_model . s e tUsb lNo i se ( p_prec i s ion ) ;725

[ v_fish_model , v_pos ] = v_fish_model . updateNextPosit ion (p_pos_meas , p_t) ;727

obj . fish_model_vec ( v_idx ) = v_fish_model ;729

end % f � updateFishModel731

733

735 % get f i s h model from l i s t by i t s id or add a new f i s h�model i f not% found

737 f unc t i on [ obj , v_fish_model , v_idx ] = getFishModel ( obj , p_id , p_pos )

739 f o r i = 1 : l ength ( obj . f ish_model_ids )

741 i f obj . f ish_model_ids ( i ) == p_idv_fish_model = obj . fish_model_vec ( i ) ;

743 v_idx = i ;r e turn ;

745 end ;

747 end ;

749 v_fish_model = fish_model ( p_id , p_pos , obj . sys_noise , obj . meas_noise ) ;

751 % f i sh�model was not found , add new f i s h�modelobj . fish_model_vec = [ obj . fish_model_vec , v_fish_model ] ;

753 obj . f ish_model_ids = [ obj . f ish_model_ids , p_id ] ;

755 v_idx = length ( obj . f ish_model_ids ) ;

757 end % f � getFishModel

759

761 f unc t i on obj = updateFish ( obj , v_fish )

763 v_id = v_fish . ge t Id ( ) ;

765 f o r i = 1 : l ength ( obj . f i s h_ id s )i f obj . f i s h_ id s ( i ) == v_id ; % i f f i s h a l r eady e x i s t s in l i s t o f f i s h e s

767

% th i s i s nece s sa ry as each f i s h keeps t rack o f i t s p o s i t i o n s .769 % the ob j e c t in the f i s h vec to r needs to be updated with

% the f i s h from the l a s t event .771 obj . f i sh_vec ( i ) = v_fish ;

end773 end

775 end % f � updateFish

86

777 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

779 f unc t i on v_noisy_signal = c r ea t eNo i syS i gna l ( obj , v_usbl , v_event )

781 t imes = v_event . getTimes ( ) ;p o s i t i o n s = v_event . g e tPo s i t i o n s ( ) ;

783 s i g n a l s = v_event . g e t S i gna l s ( ) ;

785 i f v_event . getType ( ) == 0v_delay = v_usbl . getDelayFromPos ( p o s i t i o n s ) ;

787 v_noisy_signal = obj . channel . no i s yS i gna l ( s i gna l s , v_delay ) ;

789 e l s e % in case o f c o l i s i o n[m, n ] = s i z e ( p o s i t i o n s ) ;

791 f o r i = 1 : n

793 v_delay = v_usbl . getDelayFromPos ( p o s i t i o n s ( : , i ) ) ;v_delay = v_delay + times ( i ) � t imes (1 ) ;

795

i f i ==1797 v_noisy_signal = obj . channel . no i s yS i gna l ( s i g n a l s ( i ) , v_delay ) ;

e l s e799 v_noisy_signal = v_noisy_signal + obj . channel . no i s yS i gna l ( s i g n a l s ( i ) ,

v_delay ) ;end

801

end803

end805

end % f � c r ea t eNo i syS i gna l807

809 f unc t i on obj = show( obj )

811 % plo t p o s i t i o n h i s t o r i e s o f f i s h e s and f i shmode l s

813 co l ou r = l i n e s ( l ength ( obj . f i sh_vec ) ) ;

815 % example :% p lo t ( a , ’ marker ’ , ’ x ’ , ’ co lo r ’ , ’ r ’ , ’ l i n e s t y l e ’ , ’ : ’ )

817

% show a l l f i s h e s toge the r with c o l l i s i o n s819 h1 = f i g u r e ( ) ;

hold on ;821

x l ab e l ( ’X ’ ) ;823 y l ab e l ( ’Y ’ ) ;

z l a b e l ( ’Z ’ ) ;825

g r id on ;827

t i t l e ( ’ A l l f i s h e s and c o l l i s i o n s ’ ) ;829

f i sh_legend = {} ;831

handles = [ ] ;833 l egends = {} ;

835 f o r i = 1 : l ength ( obj . f i sh_vec )

837 v_id = obj . f i s h_ id s ( i ) ;

839 v_pos_h = obj . f i sh_vec ( i ) . getPosH ( ) ;

841 s e t (0 , ’ CurrentFigure ’ , h1 )hold on ;

843 f i sh_legend = [ f i sh_legend ; cat (2 , ’ Fish ID : ’ , i n t 2 s t r ( v_id ) ) ] ;

845 % re a l f i s h po s i t i o n

87

p_fish = plo t3 (v_pos_h ( 1 , : ) ,v_pos_h ( 2 , : ) ,v_pos_h ( 3 , : ) , ’ marker ’ , ’ x ’ , ’ c o l o r ’, c o l ou r ( i , : ) ) ;

847

849 h2 = f i g u r e ( ) ;hold on ;

851

x l ab e l ( ’X ’ ) ;853 y l ab e l ( ’Y ’ ) ;

z l a b e l ( ’Z ’ ) ;855

g r id on ;857

% re a l f i s h po s i t i o n859 p_fish = plo t3 (v_pos_h ( 1 , : ) ,v_pos_h ( 2 , : ) ,v_pos_h ( 3 , : ) , ’ marker ’ , ’ x ’ , ’ c o l o r ’

, ’ b ’ ) ;l = ’ r e a l Fish pos ’ ;

861 l egends = [ l egends ; l ] ;

863 h3 = f i g u r e ( v_id ) ;h = [ h2 ; h3 ] ;

865

handles = [ handles , h ] ;867

869 end

871

s e t (0 , ’ CurrentFigure ’ , h1 )873 hold on ;

f i sh_legend = c e l l s t r ( f i sh_legend ) ;875 l egend ( f i sh_legend ) ;

877 [m, n ] = s i z e ( obj . c o l l i s i o n s ) ;f o r i = 1 : n

879 % re a l f i s h po s i t i o np_fish = plo t3 ( obj . c o l l i s i o n s ( 1 , : ) , obj . c o l l i s i o n s ( 2 , : ) , obj . c o l l i s i o n s ( 3 , : ) , ’

marker ’ , ’O ’ , ’ c o l o r ’ , ’ r ’ , ’ l i n e s t y l e ’ , ’ none ’ ) ;881

end883

885 f o r i = 1 : l ength ( obj . fish_model_vec )

887 v_id = obj . f ish_model_ids ( i ) ;

889 idx = obj . f i ndF i sh ( v_id ) ;h = handles (1 , idx ) ;

891

l = l egends ( idx ) ;893

s e t (0 , ’ CurrentFigure ’ ,h )895

897 x l ab e l ( ’X ’ ) ;y l ab e l ( ’Y ’ ) ;

899 z l a b e l ( ’Z ’ ) ;

901 g r id on ;

903 t i t l e ( cat (2 , ’ Fish ID : ’ , i n t 2 s t r ( v_id ) ) ) ;

905 hold on ;

907 v_pos_m_h = obj . fish_model_vec ( i ) . getPosH ( ) ;

909 v_pos_meas_h = obj . fish_model_vec ( i ) . getMeasPosH ( ) ;

911 v_pos_pred_h = obj . fish_model_vec ( i ) . getPredPosH ( ) ;

913 v_pred_P_h = obj . fish_model_vec ( i ) . getPredPH ( ) ;

88

915 v_upd_P_h = obj . fish_model_vec ( i ) . getUpdPH( ) ;

917 i f l ength (v_pos_m_h)>0% est imated po s i t i o n by update s tep

919 p_model = p lo t3 (v_pos_m_h( 1 , : ) ,v_pos_m_h( 2 , : ) ,v_pos_m_h( 3 , : ) , ’ marker ’ , ’. ’ , ’ c o l o r ’ , ’ g ’ , ’ l i n e s t y l e ’ , ’ : ’ ) ;

a = ’ est imated pos ’ ;921 l = [ l ; a ] ;

end923

i f l ength (v_pos_meas_h)>0925 % measured po s i t i o n by USBLs

p_meas = plo t3 (v_pos_meas_h ( 1 , : ) ,v_pos_meas_h ( 2 , : ) ,v_pos_meas_h ( 3 , : ) , ’marker ’ , ’ . ’ , ’ c o l o r ’ , ’ r ’ , ’ l i n e s t y l e ’ , ’ : ’ ) ;

927 a = ’ measured pos ’l = [ l ; a ] ;

929 end

931 i f l ength (v_pos_pred_h)>0% pred i c t ed p o s i t i o n s by p r ed i c t i on step

933 p_model_pred = plo t3 (v_pos_pred_h ( 1 , : ) , v_pos_pred_h ( 2 , : ) , v_pos_pred_h( 3 , : ) , ’ marker ’ , ’+ ’ , ’ c o l o r ’ , ’ k ’ , ’ l i n e s t y l e ’ , ’ : ’ ) ;

a = ’ p r ed i c t ed pos ’935 l = [ l ; a ] ;

end937

% show unce r ta in ty by drawing e l l i p s o i d s939 %fo r j = 1 : l ength (v_pos_pred_h)

% obj . showP(h , v_pos_pred_h ( : , j ) , v_pred_P_h ( : , j ) ) ;941 %end

943 l = c e l l s t r ( l ) ;l egend ( l ) ;

945

947 x l ab e l ( ’X ’ ) ;y l ab e l ( ’Y ’ ) ;

949 z l a b e l ( ’Z ’ ) ;

951 g r id on ;

953 obj . vusbl_prim . show (h) ;v_p_usbl = obj . vusbl_prim . getPos ( ) ;% + [ 1 0 ; 1 0 ; 0 ] ;

955 t ext ( v_p_usbl (1 ) , v_p_usbl (2 ) , v_p_usbl (3 ) , ’ Sur face Robotic Tool � 1 ’ ) ;

957 i f i s o b j e c t ( obj . vusbl_sec )obj . vusbl_sec . show (h) ;

959 v_p_usbl = obj . vusbl_sec . getPos ( ) ;% + [ 1 0 ; 1 0 ; 0 ] ;t ex t ( v_p_usbl (1 ) , v_p_usbl (2 ) , v_p_usbl (3 ) , ’ Portab le Tool � 2 ’ ) ;

961 endend

963

965 obj . showStats ( handles ) ;

967 end % f � p lo t

969

f unc t i on obj = showP( obj , h , c , r )971 % xc , yc , zc are c en te r and xr , yr , z r r ep r e s en t unce r ta in ty in each

% d i r e c t i o n973 s e t (0 , ’ CurrentFigure ’ ,h ) ;

975 norm( r )

977 i f norm( r ) < 200 % in case o f too high unce r ta in ty drawing w i l l be a messhold on ;

979 [ x , y , z ] = e l l i p s o i d ( c (1 ) , c (2 ) , c (3 ) , r (1 ) , r (2 ) , r (3 ) ) ;mesh (x , y , z ) ;

981 endend

983

89

985

f unc t i on obj = showStats ( obj , handles )987

989 f o r i = 1 : l ength ( obj . f ish_model_ids )

991 v_times = [ ] ;v_prec i s i ons = [ ] ;

993 v_meas_errors = [ ] ;v_ids = [ ] ;

995 v_upd_errors = [ ] ;

997 f o r j = 1 : l ength ( obj . d e t e c t i on s )

999

v_id = obj . d e t e c t i on s ( j ) . getDetID ( ) ;1001

i f v_id == obj . f ish_model_ids ( i )1003

v_event = obj . d e t e c t i on s ( j ) . getEvent ( ) ;1005

v_times = [ v_times , v_event . getTimes ( ) ] ;1007 v_prec i s i ons = [ v_prec i s ions , obj . d e t e c t i on s ( j ) . g e tP r e c i s i o n ( ) ] ;

v_meas_errors = [ v_meas_errors , obj . d e t e c t i on s ( j ) . getMeasError ( ) ] ;1009 v_upd_errors = [ v_upd_errors , obj . d e t e c t i on s ( j ) . getUpdError ( ) ] ;

v_ids = [ v_ids , obj . d e t e c t i o n s ( j ) . getDetID ( ) ] ;1011

end1013

1015 end

1017 idx = obj . f i ndF i sh ( obj . f ish_model_ids ( i ) ) ;h = handles (2 , idx ) ;

1019 s e t (0 , ’ CurrentFigure ’ ,h )

1021 hold on ;t i t l e ( cat (2 , ’ S t t i s t i c s f o r f i s h ID : ’ , i n t 2 s t r ( obj . f ish_model_ids ( i ) ) ) ) ;

1023

1025 s i z e ( v_prec i s i ons )s i z e ( v_meas_errors )

1027 s i z e ( v_upd_errors )

1029 v_prec i s i ons = v_prec i s ions ’ ;v_meas_errors = v_meas_errors ’ ;

1031 v_upd_errors = v_upd_errors ’ ;

1033 v_stats = cat (2 , v_prec i s ions , v_meas_errors , v_upd_errors ) ;

1035 x l ab e l ( ’X ’ ) ;y l ab e l ( ’Y ’ ) ;

1037 z l a b e l ( ’Z ’ ) ;

1039 bar ( v_stats )

1041 g r id on ;l egend ( ’ P r e c i s i on ’ , ’ Measurement e r r o r ’ , ’ Est imation e r r o r ’ )

1043

%subplot ( 3 , 1 , 1 ) ;1045 %i f l ength ( v_prec i s i ons )>0

% %stem ( v_prec i s i ons )1047 % bar ( v_prec i s i ons ) ;

% x l ab e l ( ’m’ ) ;1049 % y lab e l ( ’ i t e r a t i o n step ’ ) ;

%end1051 %t i t l e ’ Measurement Prec i s i on ’ ;

1053 %subplot ( 3 , 1 , 2 ) ;%i f l ength ( v_meas_errors )>0

1055 % %stem ( v_meas_errors ) ;

90

% bar ( v_meas_errors ) ;1057 % x lab e l ( ’m’ ) ;

% y l ab e l ( ’ i t e r a t i o n step ’ ) ;1059 %end

%t i t l e ’ Measurement Error ’ ;1061

%subplot ( 3 , 1 , 3 ) ;1063 %i f l ength ( v_upd_errors )>0

% %stem ( v_upd_errors ) ;1065 % bar ( v_upd_errors ) ;

% x l ab e l ( ’m’ ) ;1067 % y lab e l ( ’ i t e r a t i o n step ’ ) ;

%end1069 %t i t l e ’ Est imation Error ’ ;

1071

end1073

end1075

1077

1079 f unc t i on p r i n t_ f i l e ( obj )

1081 % crea t e a matrix y , with two rowsx = 0 : 0 . 1 : 1 ;

1083 y = [ x ; exp (x ) ] ;

1085 % open a f i l e f o r wr i t i ngf i d = fopen ( ’ output . txt ’ , ’w ’ ) ;

1087

% pr in t a t i t l e , f o l l owed by a blank l i n e1089 f p r i n t f ( f i d , ’ Detec t ions : \ n\n ’ ) ;

1091

obj . f i l e_format_header = cat (2 , obj . f i le_format_header , ’%s %s\n ’ ) ;1093 obj . f i l e_heade r = [ obj . f i l e_heade r ; ’ exec no | ’ , ’ time [ s ] | ’ , ’ Fish ID | ’ , ’ measured

Pos [m] | ’ , ’ r e a l Pos [m] | ’ , . . .’ d i r e c t i o n 1 [m] | ’ , ’ r e a l d i r 1 [m] | ’ , ’ ang l e s 1 [ ] | ’ , ’ r e a l ang l e s 1 [ ] | ’ , ’

a t t enuat ion 1 [dB ] | ’ , . . .1095 ’ d i r e c t i o n 2 [m] | ’ , ’ r e a l d i r 2 [m] | ’ , ’ ang l e s 2 [ ] | ’ , ’ r e a l ang l e s 2 [ ] | ’ , ’

a t t enuat ion 2 [dB ] | ’ ] ;

1097 %obj . f i le_format_output = cat (2 , obj . f i le_format_output , ’%u %f %f %f \n ’ ) ;

1099 % pr in t va lue s in column order% x va lues appear on each row o f the f i l e

1101

1103 f p r i n t f ( f i d , obj . f i le_format_header , obj . f i l e_heade r ) ;%f p r i n t f ( f i d , obj . f i le_format_output , obj . f i l e_output ) ;

1105

f p r i n t f ( f i d , ’ \n ’ ) ;1107 f p r i n t f ( f i d , obj . f i l e_output ) ;

1109

f c l o s e ( f i d ) ;1111

1113 %dlmwrite ( ’ output . txt ’ , obj . f i l e_output , ’�append ’ , ’ d e l im i t e r ’ , ’\ t ’ ) ;% ’ p r e c i s i on ’ , 6)

1115

1117

1119 end

1121 end % methodsend

./matlab/sistem.m

91

C.3 usbl.m

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

C.4 channel.m

1 c l a s s d e f channel%UNTITLED Summary o f t h i s c l a s s goes here

3 % Adding no i s e to o r i g i n a l s i g n a l . This happens on one dimension [1∗ x ]% matrix and not on [4∗ x ] dimensions as four channe l s are to c l o s e

5 % fo r no i s e not to be c o r r e l a t e d . Small no i s e i s added l o c a l l y on% each channel at USBL.

7

p r op e r t i e s ( GetAccess = ’ pub l i c ’ , SetAccess = ’ p r i va t e ’ )9 no i s e_ l ev e l % l e v e l o f no i s e to be added to s i g n a l

noise_sensor_111 noise_sensor_2

noise_sensor_313 noise_sensor_4

array_length % length o f array o f no i sy s i gna l , standard i s 1000015 s i g n a l % s e r i e s o f p ings bu i l d i ng a s i g n a l that codes the i d e n t i t y

% o f the f i s h17 no i sy_s igna l % o r i g i n a l s i g n a l with WAGNoise and a l e a t o ry po s i t i o n

% in a f i x array o f x samples � TODO: lenght o f array19

at tenuat ion % attenuat ion c a l c u l a t ed by Francois�Garr ison Model21 % 1 � at tenuat ion i s cons ide r ed s ph e r i c a l

% 2 � at tenuat ion i s cons ide r ed c y l i n d r i c a l23 attenuat ion_type

opening_angle % angle f o r c y l i n d r i c a l decay law25

end27

methods29

% Creator , Getter s and S e t t e r s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%31

f unc t i on obj = channel ( )33

obj . no i s e_ l ev e l = 0 . 0 0 2 ;35 obj . noise_sensor_1 = 0 . 0 0 1 ;

obj . noise_sensor_2 = 0 . 0 0 1 ;37 obj . noise_sensor_3 = 0 . 0 0 1 ;

obj . noise_sensor_4 = 0 . 0 0 1 ;39 obj . array_length = 10000 ;

41 obj . a t t enuat ion = 1 ;obj . attenuation_type = 1 ; % sph e r i c a l i s standard

43 obj . opening_angle = 45 ; % angle f o r c y l i n d r i c a l decay lawend % f � Constructor

45

f unc t i on obj = in i tChanne l ( no i s e_leve l , array_length )47 obj . no i s e_ l ev e l = no i s e_ l ev e l ;

obj . array_length = array_length ;49 end % f � in i tChanne l

51 f unc t i on s i g n a l = ge tS i gna l ( obj )

108

s i g n a l = obj . s i g n a l ;53 end % f � ge tS i gna l

55 f unc t i on obj = s e t S i g n a l ( obj , s i g n a l )obj . s i g n a l = s i g n a l ;

57 end % f � s e t S i g n a l

59 % A ge t t e r that i s not suposed to be used , as no i sy s i g n a l i s% created by f � no i s yS i gna l

61 f unc t i on no i sy_s igna l = getNo i syS igna l ( obj )no i sy_s igna l = obj . no i sy_s igna l ;

63 end % f � ge tNo i syS igna l

65 % A s e t t e r that i s not supposed to be used , as no i sy s i g n a l i s% created from o r i g i n a l s i g n a l in f � no i s yS i gna l

67 f unc t i on obj = se tNo i syS i gna l ( obj , no i sy_s igna l )obj . s i g n a l = no i sy_s igna l ;

69 end % f � s e tNo i syS i gna l

71 f unc t i on no i s e_ l eve l = getNo i s eLeve l ( obj )no i s e_ l ev e l = obj . no i s e_ l eve l ;

73 end % f � getNo i s eLeve l

75 f unc t i on obj = se tNo i s eLeve l ( obj , no i s e_ l ev e l )obj . no i s e_ l ev e l = no i s e_ l ev e l ;

77 end % f � s e tNo i s eLeve l

79

f unc t i on noise_sensor_1 = getNoi seSensor1 ( obj )81 noise_sensor_1 = obj . noise_sensor_1 ;

end % f � getNoi seSensor183

f unc t i on obj = se tNo i s eSenso r1 ( obj , noise_sensor_1 )85 obj . noise_sensor_1 = noise_sensor_1 ;

end % f � s e tNo i s eSensor187

f unc t i on noise_sensor_2 = getNoi seSensor2 ( obj )89 noise_sensor_2 = obj . noise_sensor_2 ;

end % f � getNoi seSensor291

f unc t i on obj = se tNo i s eSenso r2 ( obj , noise_sensor_2 )93 obj . noise_sensor_2 = noise_sensor_2 ;

end % f � s e tNo i s eSensor295

f unc t i on noise_sensor_3 = getNoi seSensor3 ( obj )97 noise_sensor_3 = obj . noise_sensor_3 ;

end % f � getNoi seSensor399

f unc t i on obj = se tNo i s eSenso r3 ( obj , noise_sensor_3 )101 obj . noise_sensor_3 = noise_sensor_3 ;

end % f � s e tNo i s eSensor3103

f unc t i on noise_sensor_4 = getNoi seSensor4 ( obj )105 noise_sensor_4 = obj . noise_sensor_4 ;

end % f � getNoi seSensor4107

f unc t i on obj = se tNo i s eSenso r4 ( obj , noise_sensor_4 )109 obj . noise_sensor_4 = noise_sensor_4 ;

end % f � s e tNo i s eSensor4111

113 f unc t i on array_length = getArrayLength ( obj )array_length = obj . array_length ;

115 end % f � getArrayLength

117 f unc t i on obj = setArrayLength ( obj , array_length )obj . array_length = array_length ;

119 end % f � setArrayLength

121

f unc t i on obj = setAttenuationType ( obj , p_type )123 obj . attenuation_type = p_type ;

109

end % f � setArrayLength125

127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Function to add no i s e and <<hide>> o r i g i n a l s i g n a l in l onge r

129 % s i g n a l . Sub s t i t u t e s g e t t e r and s e t t e r A g e t t e r that i s not suposed to be used ,as no i sy s i g n a l i s

% created by f � no i s yS i gna l131 f unc t i on no i sy_s igna l = no i s yS i gna l ( obj , s i gna l , v_delay )

133 s i g n a l = s i g n a l .∗ obj . a t t enuat ion ;

135

no i sy_s igna l = cat (2 , z e r o s (1 , round ( r e a l ( v_delay (1 ) ) ) ) , s i g n a l ) ;137 no i sy_s igna l = cat (2 , no i sy_s igna l , z e r o s (1 , obj . array_length � round ( r e a l (

v_delay (1 ) ) ) ) ) ;

139 % noi sy_s igna l = in t e rp ( no i sy_s igna l , 1 00 ) ; % i n t e r p o l a t i o n

141

f o r i = 2 : l ength ( v_delay )143 noisy_signal_tmp = cat (2 , z e r o s (1 , round ( r e a l ( v_delay ( i ) ) ) ) , s i g n a l ) ;

noisy_signal_tmp = cat (2 , noisy_signal_tmp , z e r o s ( 1 , ( obj . array_length �round ( r e a l ( v_delay ( i ) ) ) ) ) ) ;

145

% noisy_signal_tmp = in t e rp ( noisy_signal_tmp ,100 ) ; % i n t e r p o l a t i o n147

no i sy_s igna l=cat (1 , no i sy_s igna l , noisy_signal_tmp ) ;149

end151

[m, n ] = s i z e ( no i sy_s igna l ) ;153 no i sy_s igna l = no i sy_s igna l + obj . no i s e_ l ev e l .∗ randn (m, n) ;

no i sy_s igna l ( 1 , : ) = no i sy_s igna l ( 1 , : ) + obj . noise_sensor_1 .∗ randn (1 , n) ;155 no i sy_s igna l ( 2 , : ) = no i sy_s igna l ( 2 , : ) + obj . noise_sensor_2 .∗ randn (1 , n) ;

no i sy_s igna l ( 3 , : ) = no i sy_s igna l ( 3 , : ) + obj . noise_sensor_3 .∗ randn (1 , n) ;157 no i sy_s igna l ( 4 , : ) = no i sy_s igna l ( 4 , : ) + obj . noise_sensor_4 .∗ randn (1 , n) ;

end % f � no i s yS i gna l159

161 % TODO: complete . . . l ength c o r r e c t ? i n t e rp or not ?% Function to add no i s e and <<hide>> o r i g i n a l s i g n a l in l onge r

163 % s i g n a l . Sub s t i t u t e s g e t t e r and s e t t e r A g e t t e r that i s not suposed to be used ,as no i sy s i g n a l i s

% created by f � addNoiseToSignal165 f unc t i on no i sy_s igna l = addNoiseToSignal ( obj , s i g n a l )

167 [m, n ] = s i z e ( s i g n a l ) ;

169 no i sy_s igna l = s i g n a l + obj . no i s e_ l ev e l .∗ randn (m, n) ;

171 no i sy_s igna l ( 1 , : ) = no i sy_s igna l ( 1 , : ) + obj . noise_sensor_1 .∗ randn (1 , n) ;no i sy_s igna l ( 2 , : ) = no i sy_s igna l ( 2 , : ) + obj . noise_sensor_2 .∗ randn (1 , n) ;

173 no i sy_s igna l ( 3 , : ) = no i sy_s igna l ( 3 , : ) + obj . noise_sensor_3 .∗ randn (1 , n) ;no i sy_s igna l ( 4 , : ) = no i sy_s igna l ( 4 , : ) + obj . noise_sensor_4 .∗ randn (1 , n) ;

175

end % f � addNoiseToSignal177

179 % Calcu la te a t t enuat ion at t in [ dB/km] by s ph e r i c a l% propagat ion model and Model o f Francois�Garr ison

181 % from pos1 to pos2% type dec ide s f o r sphe r i c or c y l i n d r i c a l decay

183 % type = 1 � sphe r i c a l , 2 � c y l i n d r i c a l% z = depth in [m]

185 % S = s a l i n i t y in [ p . s . u . ]% T = Temperature in Ce l s i u s [C]

187 % f = frequency o f the s i g n a l in [ kHz ]func t i on [ obj , a t t_tota l ] = setAttenuat ion ( obj , z , f , S , T, pos1 , pos2 ,

channel_height )189

% assume neut ra l pH value in water

110

191 pH = 7 ;

193 % cont r i bu t i on o f bo r i c ac id B(OH)3 :c = 1412 + 3.21∗T + 1.19∗S + 0.0167∗ z ;

195 A1 = 8.86/ c ∗ 10^(0.78 ∗ pH � 5) ;P1 = 1 ;

197 f 1 = 2 .8∗ s q r t (S/35) ∗10^(4�1245/(T + 273) ) ;

199

% cont r i bu t i on o f magnesium s u l f a t e Mg(SO) 4 :201 A2 = 21.44 ∗ S/c ∗(1+0.025∗T) ;

P2 = 1�1.37∗10^(�4)∗z + 6 .2 ∗ 10^(�9)∗z ^2;203 f 2 = 8 .17 ∗ 10^(8�1990/(T+273) ) / (1 + 0.0018 ∗ (S�35) ) ;

205 % cont r i bu t i on due to pure water v i s c o s i t yP3 = 1 � 3 .83 ∗ 10^(�5)∗z + 4.9∗10^(�10) ∗ z ^2;

207 i f T < 20A3 = 4.937 ∗ 10^(�4)�2.59 ∗ 10^(�5)∗T + 9.11 ∗ 10^(�7) ∗ T^2 � 1 .5 ∗

10^(�8) ∗ T^3;209 e l s e

A3 = 3.964 ∗ 10^(�4)�1.146 ∗ 10^(�5)∗T + 1.45 ∗ 10^(�7) ∗ T^2 � 6 .5 ∗10^(�10) ∗ T^3;

211 end

213 %%%%%%%% Attenuation due to absorbt ion in [ dB/km]

215 att_absorbt ion = A1 ∗ P1 ∗ f 1 ∗ f ^2 / ( f1 ^2 + f ^2) + A2 ∗ P2 ∗ f 2 ∗ f ^2 / ( f2 ^2 +f ^2) + A3 ∗ P3 ∗ f ^2;

217 % from [dB/km] to [ dB/m]att_absorbt ion = att_absorbt ion / 1000 ;

219 %%%%%%%

221 i f obj . attenuat ion_type == 1 % sph e r i c a l% sph e r i c a l decay

223 a t t_sphe r i c a l = 20 ∗ l og10 (norm( pos2�pos1 ) ) ;

225 at t_tota l = att_absorbt ion + at t_sphe r i c a l

227 e l s e% c y l i n d r i c a l decay

229 % assuming tag on mid he ight o f channel% assuming graz ing ang le at 45 degree s

231

r_0 = channel_height /(2∗ tan ( obj . opening_angle ) ) ;233

i f norm( pos2�pos1 ) <= r_0235 a t t_cy l i n d r i c a l = 20 ∗ l og10 (norm( pos2�pos1 ) ) ;

e l s e237 a t t_cy l i n d r i c a l = 10∗ l og10 ( ( norm( pos2 ( 1 : 2 )�pos1 ( 1 : 2 ) ) ) ∗ r_0) ;

end239

at t_tota l = att_absorbt ion + a t t_cy l i n d r i c a l ;241

end243

obj . a t t enuat ion = (1/ 10^(( a t t_tota l ) /20) ) ;245

end % f � se tAttenuat ion247

% retu rn s ac tua l a t t enuat ion in [dB ]249 f unc t i on v_attenuation = getAttenuat ion ( obj )

251 v_attenuation = �20∗ l og10 ( obj . a t t enuat ion ) ;

253 end % f � getAttenuat ion

255

f unc t i on obj = setAttenuat ionOf f ( obj )257 obj . a t t enuat ion = 1 ;

259 end

111

261 end

263 end

./matlab/channel.m

C.5 eventList.m

c l a s s d e f even tL i s t2 %Li s t o f Events

4 p r op e r t i e s ( GetAccess = ’ p r i va t e ’ , SetAccess = ’ p r i va t e ’ )events_to_process % L i s t o f events that are to be proce s sed

6 event_times % Times o f events f o r s o r t i n g o f queue events_to_processprocessed_events % L i s t o f events that are to be proce s s ed

8

end10

methods12

f unc t i on obj = eventL i s t ( )14 obj . events_to_process = [ ] ;

obj . event_times = [ ] ;16 obj . processed_events = [ ] ;

end18

% Adds a new event to events_to_process mainta in ing order in time .20 f unc t i on obj = addEvent ( obj , v_event )

22 t imes = v_event . getTimes ( ) ;t = times (1 ) ;

24

i f isempty ( obj . event_times )26 obj . events_to_process = [ v_event ] ;

temp_t = v_event . getTimes ( ) ;28 obj . event_times = temp_t (1) ;

30 r e turn ; % event l i s t should be a so r t ed l i s t . I f event was so r t ed in ,r ou t in e can be ended

32 end

34

i f obj . event_times (1 ) > t36

i f obj . event_times (1 ) > t + 4 % no c o l i s i o n38 obj . event_times = [ t , obj . event_times ] ;

obj . events_to_process = [ v_event , obj . events_to_process ] ;40 e l s e % case o f c o l i s i o n , s e v e r a l events form one , jo inEvents

obj . event_times (1 ) = t ;42 obj . events_to_process (1 ) = v_event . jo inEvents ( obj . events_to_process

(1 ) ) ;end

44

r e turn ; % event l i s t should be a so r t ed l i s t . I f event was so r t ed in ,r ou t in e can be ended

46 end

48 f o r i = 1 : l ength ( obj . event_times )i f obj . event_times ( i ) > t

50 i f obj . event_times ( i ) > t + 4 % no c o l i s i o n

52 temp_t = obj . event_times ( i : l ength ( obj . event_times ) ) ;temp_e = obj . events_to_process ( i : l ength ( obj . event_times ) ) ;

54

obj . event_times = obj . event_times ( 1 : ( i ) ) ;56 obj . events_to_process = obj . events_to_process ( 1 : ( i ) ) ;

112

58 obj . event_times ( i ) = t ;obj . events_to_process ( i ) = v_event ;

60

obj . event_times = [ obj . event_times , temp_t ] ;62 obj . events_to_process = [ obj . events_to_process , temp_e ] ;

64 e l s e

66

68

i70

72

74 obj . event_times ( i ) = t ;

76 obj . events_to_process ( i ) = v_event . jo inEvents ( obj .events_to_process ( i ) ) ;

end78

r e turn ; % event l i s t should be a so r t ed l i s t . I f event was so r t ed in ,r ou t in e can be ended

80

end82

end84

i f t > obj . event_times ( l ength ( obj . event_times ) ) + 486 obj . event_times = [ obj . event_times , t ] ;

obj . events_to_process = [ obj . events_to_process , v_event ] ;88 e l s e

obj . events_to_process ( l ength ( obj . events_to_process ) ) = v_event . jo inEvents( obj . events_to_process ( l ength ( obj . events_to_process ) ) ) ;

90

end92

94 end % f � addEvent

96

% Returns the event at index index from events_to_process .98 f unc t i on event = getEvent ( obj , index )

event = obj . events_to_process ( index ) ;100 end % f � getEvent

102

% Returns the event at index index from events_to_process .104 f unc t i on n = numberEventsToProcess ( obj )

n = length ( obj . events_to_process ) ;106 end % f � numberEventsToProcess

108

% Removes the event at index index from events_to_process .110 f unc t i on obj = removeEvent ( obj , index )

obj . events_to_process ( index ) = [ ] ;112 obj . event_times ( index ) = [ ] ;

end % f � removeEvent114

% Removes the f i r s t event from events_to_process and adds i t to processed_events .116 f unc t i on [ obj , v_event ] = getNextEvent ( obj )

118 v_event = obj . events_to_process (1 ) ;obj . processed_events = [ obj . processed_events , obj . events_to_process (1 ) ] ;

120

i f obj . numberEventsToProcess ( ) > 1122 obj . events_to_process = obj . events_to_process ( 2 : l ength ( obj .

events_to_process ) ) ;obj . event_times = obj . event_times ( 2 : l ength ( obj . event_times ) ) ;

124

e l s e

113

126 obj . events_to_process = [ ] ;obj . event_times = [ ] ;

128 end

130 end % f � getNextEvent

132

134 f unc t i on obj = pr in t ( obj )’ Processed Events : ’

136

f o r i = 1 : l ength ( obj . processed_events )138 obj . processed_events ( i ) . p r i n t ( ) ;

end140

’ Events to Process : ’142

f o r i = 1 : l ength ( obj . events_to_process )144 obj . events_to_process ( i ) . p r i n t ( ) ;

end146

148 end % f � pr in t

150

end152

154 end

./matlab/eventList.m

C.6 event

c l a s s d e f event2 % A s i n g l e event ( a f i s h sends a s i g n a l )

% In case o f c o l i s i o n an event may have two f i s h e s , s i g n a l s .4

6 p r op e r t i e s ( GetAccess = ’ p r i va t e ’ , SetAccess = ’ p r i va t e ’ )t imes % Time the event happens . In case o f c o l i s i o n , the re w i l l be s t i l l ,

j u s t one i n i t i a l time8 f i s h e s % Fish that emits the s i g n a l

p o s i t i o n s % Pos i t i on at which the f i s h i s when emit t ing the s i g n a l10 s i g n a l s % emitted s i g n a l by f i s h , in case o f c o l i s i o n , combination o f both

s i g n a l stype % Type d e f i n e s wether i t i s a p r ed i c t i on or update event

12 % 0 � s i n g l e ; 1 � c o l i s i o n

14 standard_duration % standard durat ion o f a s i n g l e s i g n a l

16 pr e s su r e % Pressure measure � i s only used i f measurement method with only oneUSBL i s used

18 end

20 methods

22 f unc t i on obj = event (p_time , p_fish , p_posit ion , p_signal , p_pressure )i f ( narg in > 0)

24 obj . t imes = p_time ;obj . f i s h e s = p_fish ;

26 obj . p o s i t i o n s = p_posit ion ;obj . s i g n a l s = p_signal ;

28 obj . type = 0 ; % s i n g l e event

30 obj . p r e s su r e = p_pressure ;

32 obj . standard_duration = 4 ; % standard durat ion o f a s i n g l e s i g n a l i s 4 s

114

end34 end % f � Constructor

36

38 f unc t i on t imes = getTimes ( obj )t imes = obj . t imes ;

40 end % f � getTime

42 f unc t i on f i s h e s = ge tF i she s ( obj )f i s h e s = obj . f i s h e s ;

44 end % f � getF i sh

46 f unc t i on p o s i t i o n s = ge tPo s i t i o n s ( obj )p o s i t i o n s = obj . p o s i t i o n s ;

48 end % f � ge tPo s i t i on

50 f unc t i on s i g n a l s = ge tS i gna l s ( obj )s i g n a l s = obj . s i g n a l s ;

52 end % f � ge tS i gna l

54 f unc t i on type = getType ( obj )type = obj . type ;

56 end % f � getType

58 f unc t i on setType ( obj , p_type )obj . type = p_type ;

60 end % f � setType

62 f unc t i on depth = getPre s su re ( obj )depth = obj . p r e s su r e ;

64 end % f � ge tPre s su re

66 % return durat ion o f event in [ s ] . A s i n g l e events s i g n a l has a% standard durat ion o f 4 s .

68 f unc t i on durat ion = getDurat ion ( obj )

70 durat ion = obj . t imes ( l ength ( obj . t imes ) )�obj . t imes (1 ) + obj . standard_duration ;

72 obj . type = v_type ;end % f � getDurat ion

74

76 f unc t i on obj = addEvent ( obj , p_event )

78 % in case event to add i s a c o l i s i o n , j u s e jo inEvent ( ) i n s t eadi f p_event . getType ( ) == 0 % i f added event i s a s i n g l e event and not a

c o l i s i o n80

t = p_event . getTimes ( ) ;82 p = p_event . g e tPo s i t i o n s ( ) ;

s = p_event . g e t S i gna l s ( ) ;84 f = p_event . g e tF i she s ( ) ;

86

i f obj . t imes (1 ) > t88

obj . t imes = [ t , obj . t imes ] ;90 obj . p o s i t i o n s = [ p , obj . p o s i t i o n s ] ;

obj . s i g n a l s = [ s , obj . s i g n a l s ] ;92 obj . f i s h e s = [ f , obj . f i s h e s ] ;

94 r e turn ; % event l i s t should be a so r t ed l i s t . I f event was ordered in, r ou t in e can be ended

end96

f o r i = 1 : l ength ( obj . t imes )98 i f obj . t imes ( i ) > t

100

temp_t = obj . t imes ( i : l ength ( obj . t imes ) ) ;102 temp_p = obj . p o s i t i o n s ( : , i : n ) ;

115

temp_s = obj . s i g n a l s ( i : l ength ( obj . s i g n a l s ) ) ;104 temp_f = obj . f i s h e s ( i : l ength ( obj . f i s h e s ) ) ;

106

obj . t imes = obj . t imes ( 1 : ( i ) ) ;108 obj . p o s i t i o n s = obj . p o s i t i o n s ( : , 1 : i ) ;

obj . s i g n a l s = obj . s i g n a l s ( 1 : ( i ) ) ;110 obj . f i s h e s = obj . f i s h e s ( 1 : ( i ) ) ;

112

obj . t imes ( i ) = t ;114 obj . p o s i t i o n s ( : , i ) = p ;

obj . s i g n a l s ( i ) = s ;116 obj . f i s h e s ( i ) = f ;

118 obj . t imes = [ obj . times , temp_t ] ;obj . p o s i t i o n s = [ obj . po s i t i on s , temp_p ] ;

120 obj . s i g n a l s = [ obj . s i gna l s , temp_s ] ;obj . f i s h e s = [ obj . f i s h e s , temp_f ] ;

122

124 r e turn ; % event l i s t should be a so r t ed l i s t . I f event wasordered in , r ou t in e can be ended

end126

end128

obj . t imes = [ obj . times , t ] ;130 obj . p o s i t i o n s = [ obj . po s i t i on s , p ] ;

obj . s i g n a l s = [ obj . s i gna l s , s ] ;132 obj . f i s h e s = [ obj . f i s h e s , f ] ;

134 obj . type = 1 ; % c o l i s i o n type

136 e l s ee r r o r ( ’Not a l lowed to add an event o f type c o l i s i o n . Use jo inElements ( )

i n s t ead ’ ) ;138 end

140 end % f � addEvent

142

f unc t i on obj = jo inEvents ( obj , p_event )144

i f p_event . getType ( ) == 0 % type s i n g l e event146

obj = obj . addEvent ( p_event ) ; % add s i n g l e event148

e l s e % type c o l i s i o n150

t = p_event . getTimes ( ) ;152 p = p_event . g e tPo s i t i o n s ( ) ;

s = p_event . g e t S i gna l s ( ) ;154 f = p_event . g e tF i she s ( ) ;

156 f o r i = 1 : l ength ( v_event . getTimes )obj = obj . addEvent ( event ( t ( i ) , f ( i ) , p ( : , i ) , s ( i ) ) ) ; % add s i n g l e

events158 end

160 end

162 obj . type = 1 ; % c o l i s i o n type

164 end % f � setType

166

168 f unc t i on p r i n t ( obj )’ Event : ’

170 i f obj . type == 0’ _Single Event ’

116

172 e l s e’ _Col i s ion ’

174 end’_Times : ’

176 obj . t imes’ _Pos it ions : ’

178 obj . p o s i t i o n s

180 end % f � pr in t

182

end % � methods184 end % � c l a s s d e f

./matlab/event.m

C.7 fish

1 c l a s s d e f f i s h%UNTITLED Summary o f t h i s c l a s s goes here

3 % Deta i l ed exp lanat ion goes here

5 p r op e r t i e s ( GetAccess = ’ p r i va t e ’ , SetAccess = ’ p r i va t e ’ )id % I d e n t i f i c a t i o n number o f f i s h ex . 33517

7 s i g n a l % s e r i e s o f p ings bu i l d i ng a s i g n a l that codes the i d e n t i t y o f the f i s hpos % cur rent p o s i t i o n o f the f i s h

9 pos_h % h i s t o r y o f the past f i s h p o s i t i o n s

11 r epea t_ in t e rva l % mean time between two pings ( events )repeat_std_dev % standard dev i a t i on from mean time between two pings ( events )

13

% There are two models o f movement :15 % 1 � By prede f i n ed po in t s

% 2 � By dynamical model17

% moving_mode d e f i n e s the model o f movement19 moving_mode

21 % 1 � prede f i n ed s t ep ss tep % actua l s tep

23 s t ep s % a l l s t ep s o f movement 3 x n matrix

25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2 � Kalman

27 A % Matrix f o r dynamic model o f f i s hB % Matrix f o r dynamic model o f f i s h

29 C % Matrix f o r dynamic model o f f i s hu % meassure f o r a c c e l e r a t i o n

31 fish_noise_magf i sh_no i s e

33 Q_loc % s t a t e va r i a b l e with po s i t i on , v e l o c i t y and a c e l e r a t i o n

35 dt % time i n t e r v a l between s t ep s

37 pres sure_no i se % no i s e that adds to p r e s su r e measurements f � ge tPre s su re ( )depth_noise % no i s e that adds to depth measurements f � getDepth ( )

39 dens i ty % dens i ty o f water , used f o r p r e s su r e measurements f �ge tPre s su re ( )

g rav i ty % earth grav i ty , used f o r p r e s su r e measurements f � ge tPre s su re ( )41 end

43 methods

45 % Constructor , Getter s and S e t t e r s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

47 f unc t i on obj = f i s h (p_id , p_signal , p_posit ion , p_repeat_interval ,p_repeat_std_dev , p_moving_mode , p_density , p_gravity )

49 obj . id = p_id ; % i d e n t i f i c a t i o n o f f i s h ex . 33517

117

51 % r e a j u s t s i g n a l s t r ength to 1 mVPP (VPP � Volt Peak to Peak )obj . s i g n a l = ( p_signal . /max( p_signal ) ) . / 2 ;

53 obj . pos = p_posit ion ;obj . pos_h = [ ] ;

55 obj . dt = p_repeat_interval ; % time i n t e r v a l between s t ep s

57 obj . r epea t_ in t e rva l = p_repeat_interval ;obj . repeat_std_dev = p_repeat_std_dev ;

59

obj . moving_mode = p_moving_mode ;61

i f p_moving_mode == 1 % prede f i n ed p o s i t i o n s63 % For movement model 1 i n i t S t e p s ( s t ep s ) needs to be c a l l e d

obj . s tep = 1 ;65 e l s e

% by i n e r t i a67 % c a l l i n g in i tS i s temMode l i s op t i ona l to d e f i n e i n i t i a l

% ve l o c i t y , a c e l e r a t i o n and no i s e l e v e l69

obj = obj . in i tMode l ( ) ;71 obj . fish_noise_mag = 0 .000002 ; %proce s s no i s e : the v a r i a b i l i t y in how

f a s t the Fish i s speed ing up ( stdv o f a c c e l e r a t i o n : meters / sec ^2)obj . Q_loc = [ obj . pos ; 0 . 1 ; 0 ; 0 ; 0 ; 0 ; 0 ] ; % the f i s h r e a l s t a t e ( d e f au l t

i n i t i a l va lue )73

end75

obj . p re s sure_no i se = 0 . 0 0 1 ; % w i l l only be used when j u s t one USBL i sa v a i l a b l e

77 obj . depth_noise = 0 . 0 0 1 ;obj . dens i ty = p_density ;

79 obj . g rav i ty = p_gravity ;

81 end % f � Constructor

83 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% With movement model 1 , t h i s f unc t i on should be c a l l e d from out s id e

85 % d i r e c t l y a f t e r f i s h con s t ruc to r% in : 3xn matrix with a l l p o s i t i o n s

87 f unc t i on obj = i n i t S t e p s ( obj , s t ep s )obj . s t ep s = s t ep s ;

89 end % f � i n i t S t e p s

91 f unc t i on obj = i n i t S t a t e ( obj , s t a t e )obj . Q_loc = s t a t e ;

93 end % f � i n i t S t a t e

95

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%97 % Receive as argument i n i t i a l v e l o c i t y vector , a c e l e r a t i o n vector , no i s e

amplitude% For movement model 0 , by i n e r t i a

99 f unc t i on obj = initSystemModel ( obj , p_vel , p_acel , p_fish_noise_mag )

101 obj . fish_noise_mag = p_fish_noise_mag ; %proce s s no i s e : the v a r i a b i l i t y in howf a s t the Fish i s speed ing up ( stdv o f a c c e l e r a t i o n : meters / sec ^2)

103 obj . Q_loc = [ obj . pos ; p_vel ; p_acel ] ; % the f i s h r e a l s t a t e

105 end % f � in i tS i s temMode l

107

f unc t i on obj = in i tF i s hNo i s e ( obj , p_fish_noise_mag )109 obj . fish_noise_mag = p_fish_noise_mag ;

end % f � i n i t S t e p s111

f unc t i on v_id = get Id ( obj )113 v_id = obj . id ;

end % f � get Id115

f unc t i on obj = s e t Id ( obj , p_id )

118

117 obj . id = p_id ;end % f � s e t I d

119

f unc t i on s i g n a l = ge tS i gna l ( obj )121 s i g n a l = obj . s i g n a l ;

end % f � ge tS i gna l123

f unc t i on obj = s e t S i g n a l ( obj , p_signal )125 obj . s i g n a l = p_signal ;

end % f � s e t S i g n a l127

f unc t i on pos = getPos ( obj )129 pos = obj . pos ;

end % f � getPos131

f unc t i on obj = setPos ( obj , p o s i t i o n )133 obj . pos = po s i t i o n ;

end % f � setPos135

f unc t i on obj = c l ea rPosH i s to ry ( obj )137 obj . pos_h = [ ] ;

end % f � c l ea rPosH i s to ry139

f unc t i on pos_h = getPosH ( obj )141 pos_h = obj . pos_h ;

end % f � c l ea rPosH i s to ry143

f unc t i on z = getDepth ( obj )145

z = obj . pos (3 ) + obj . depth_noise ∗ randn ( ) ;147 end % f � getDepth

149 f unc t i on i n t e r v a l = getRepeat Inte rva l ( obj )i n t e r v a l = obj . r epea t_ in t e rva l ;

151 end % f � ge tRepeat In te rva l

153 f unc t i on obj = se tRepea t In t e rva l ( obj , p_repeat_interval )obj . r epea t_ in t e rva l = p_repeat_interval ;

155 end % f � s e tRepea t In t e rva l

157

f unc t i on std_dev = getRepeatStdDev ( obj )159 std_dev = obj . repeat_std_dev ;

end % f � getRepeatStdDev161

f unc t i on obj = setRepeatStdDev ( obj , p_repeat_std_dev )163 obj . repeat_std_dev = p_repeat_std_dev ;

end % f � setRepeatStdDev165

f unc t i on pr e s su r e = getPre s su re ( obj )167 pre s su r e = obj . pos (3 ) ∗ obj . dens i ty ∗ obj . g rav i ty + obj . pre s sure_no i se ∗

randn ( ) ;end % f � ge tPre s su re

169

171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Calcu la te p o s i t i o n ( , v e l o c i t y and a c e l e r a t i o n ) f o r next s tep o f

173 % f i s h% return po s i t i o n

175 f unc t i on [ obj , pos ] = nextPos i t i on ( obj )

177 % There are two d i f e r e n t moving modelsi f obj . moving_mode == 1 % by prede f i n ed p o s i t i o n s

179

i f obj . s t ep < length ( obj . s t ep s )181 obj . s tep = obj . s tep + 1 ;

end183 pos = obj . s t ep s ( 1 : 3 , obj . s t ep ) ;

185 e l s e % by i n e r t i a

187 obj . f i s h_no i s e = obj . fish_noise_mag .∗ . . .

119

[ 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; randn ; randn ; randn ; ] ;189

obj . Q_loc = obj .A ∗ obj . Q_loc + obj . f i s h_no i s e ;191 pos = obj . Q_loc ( 1 : 3 ) ;

193 end

195 obj . pos = pos ;

197 % r e g i s t e r p o s i t i o n in h i s t o r yobj . pos_h = [ obj . pos_h , pos ] ;

199

end % f � nextPos i t i on201

203 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Actua l i z e matrix A o f s i s t em model depending on lapsed time

205 % return Fishfunc t i on obj = in i tMode l ( obj )

207

%% Def ine update equat ions ( Coe f f i c e n t matr i ce s ) : A phys i c s based209 %% model f o r where we expect the Fish to be

obj .A = [1 0 0 obj . dt 0 0 (1/2∗ obj . dt ^2) 0 0;

211 0 1 0 0 obj . dt 0 0 (1/2∗ obj . dt ^2) 0;

0 0 1 0 0 obj . dt 0 0 (1/2∗ obj . dt^2) ;

213 0 0 0 1 0 0 obj . dt 0 0;

0 0 0 0 1 0 0 obj . dt 0;

215 0 0 0 0 0 1 0 0 obj . dt;

0 0 0 0 0 0 1 0 0;

217 0 0 0 0 0 0 0 1 0;

0 0 0 0 0 0 0 0 1] ;

219

end % f � in i tMode l221

223 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Returns the next event ( a f i s h sends a s i g n a l ) by adding

225 %repea t_ in t e r va l l p lus a random time with in repeat_std_dev to%the provided time and s e t t i n g a new po s i t i o n f o r the f i s h

227 f unc t i on [ obj , v_event ] = getNextEvent ( obj , p_current_time )

229 obj . dt = obj . r epea t_ in t e rva l + obj . repeat_std_dev ∗ randn ( ) ; % time s i n c el a s t s tep

231 obj = obj . in i tMode l ( ) ;

233 v_t = p_current_time + obj . dt ; % time the next event w i l l occur

235

[ obj , v_p ] = obj . nextPos i t i on ( ) ; % po s i t i o n o f the f i s h in the next event atfu tu r e time t

237 obj . pos = v_p ;

239 v_event = event (v_t , obj , v_p , obj . s i gna l , obj . g e tPre s su re ( ) ) ; % s i n g l e eventtype

end % f � ca lcu lateNextEvent241

end243

end

./matlab/fish.m

120

C.8 fish_model

1 c l a s s d e f f ish_model%UNTITLED Summary o f t h i s c l a s s goes here

3 % Deta i l ed exp lanat ion goes here

5 p r op e r t i e s ( GetAccess = ’ p r i va t e ’ , SetAccess = ’ p r i va t e ’ )pos % cur rent po s i t i o n o f the f i s h

7 pos_h % h i s t o r y o f est imated f i s h p o s i t i o n spos_meas_h % h i s t o r y o f measured p o s i t i o n s f o r f i s h

9 id % Detcted and a s s o c i a t ed ID

11 pos_pred_h % pred i c t ed p o s i t i o n s ( vs updated p o s i t i o n s )

13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Kalman

15 A % Matrix f o r dynamic model o f f i s hB % Matrix f o r dynamic model o f f i s h

17 C % Matrix f o r dynamic model o f f i s hu % meassure f o r a c c e l e r a t i o n

19 Q_estimateP

21 Ezfish_noise_mag

23 usbl_noise_magtkn_x %usbl_noise_mag_x

25 tkn_y %usbl_noise_mag_ytkn_z %usbl_noise_mag_z

27

pred_P_h % h i s t o r y o f a s s o c i a t ed unce r ta in ty matrix f o r p r ed i c t i o n s t ep s29 upd_P_h % h i s t o r y o f a s s o c i a t ed unce r ta in ty matrix f o r update s t ep s

31

dt % time i n t e r v a l between s t ep s33 last_update_pos % po s i t i o n at l a s t update s tep

35 time

37 end

39 methods

41 % Constructor , Getter s and S e t t e r s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

43

f unc t i on obj = fish_model ( p_id , p_posit ion , p_sys_noise , p_meas_noise )45

% r e a j u s t s i g n a l s t r ength to 1 mVPP (VPP � Volt Peak to Peak )47 obj . id = p_id ;

obj . pos = p_posit ion ;49 obj . pos_h = [ ] ;

51 obj . last_update_pos = p_posit ion ;

53 obj . pos_pred_h = [ ] ;

55 obj . Q_estimate = [ p_posit ion ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ] ; % the f i s h r e a l s t a t e . No movementcan i n i t i a l l y be asumed

57 %s t a r t with same no i s e in measurements and f i s h ( model )obj . fish_noise_mag = p_sys_noise ; % proce s s no i s e ( o f model ) : the

v a r i a b i l i t y in how f a s t the Fish i s speed ing up ( stdv o f a c c e l e r a t i o n : meters / sec ^2)59 obj . usbl_noise_mag = p_meas_noise ; % measurement_noise o f usb l

61 obj . tkn_x = 1 ; % s t a r t with same no i s e in measurements and f i s h ( model )obj . tkn_y = 1 ; % s t a r t with same no i s e in measurements and f i s h ( model )

63 obj . tkn_z = 1 ; % s t a r t with same no i s e in measurements and f i s h ( model )

65 obj . dt = 1 ; % Todo : � a l t e r to argument with r e a l time i n t e r v a l between step

67

121

%% Def ine update equat ions ( Coe f f i c e n t matr i ce s ) : A phys i c s based69 %% model f o r where we expect the Fish to be

obj .A = [1 0 0 obj . dt 0 0 (1/2∗ obj . dt ^2) 0 0;

71 0 1 0 0 obj . dt 0 0 (1/2∗ obj . dt ^2) 0;

0 0 1 0 0 obj . dt 0 0 (1/2∗ obj . dt^2) ;

73 0 0 0 1 0 0 obj . dt 0 0;

0 0 0 0 1 0 0 obj . dt 0;

75 0 0 0 0 0 1 0 0 obj . dt;

0 0 0 0 0 0 1 0 0;

77 0 0 0 0 0 0 0 1 0;

0 0 0 0 0 0 0 0 1] ;

79

obj .C = [1 0 0 0 0 0 0 0 0 ;81 0 1 0 0 0 0 0 0 0 ;

0 0 1 0 0 0 0 0 0 ] ;83

85 Ex = [0 . 0 100 0 0 0 .0020 0 0 0 .00020 0

0 0 .0100 0 0 0 .0020 0 00 .0002 0

87 0 0 0 .0100 0 0 0 .0020 00 0 .0002

0 .0020 0 0 0 .0004 0 0 0 .00010 0

89 0 0 .0020 0 0 0 .0004 0 00 .0001 0

0 0 0 .0020 0 0 0 .0004 00 0 .0001

91 0 .0002 0 0 0 .0001 0 0 0 .00010 0

0 0 .0002 0 0 0 .0001 0 00 .0001 0

93 0 0 0 .0002 0 0 0 .0001 00 0 . 0 0 0 1 ] ;

95

obj .P = Ex ; % est imate o f i n i t i a l Fish po s i t i o n var i ance ( covar iance matrix )97

99 obj . Ez = [ obj . dt^4/4 0 0 ; . . .0 obj . dt^4/4 0 ; . . .

101 0 0 obj . dt^4/4 ] . ∗ obj . usbl_noise_mag ^2;

103 end % f � Constructor

105 f unc t i on s i g n a l = get Id ( obj )s i g n a l = obj . id ;

107 end % f � get Id

109 f unc t i on obj = s e t Id ( obj , id )obj . id = id ;

111 end % f � s e t I d

113 f unc t i on time = getTime ( obj )time = obj . time ;

115

end % f � getTime117

f unc t i on s i g n a l = getPos ( obj )119 s i g n a l = obj . pos ;

end % f � getPos121

122

f unc t i on obj = setPos ( obj , p o s i t i o n )123 obj . pos = po s i t i o n ;

obj . Q_loc_estimate ( 1 : 3 ) = po s i t i o n ;125 end % f � setPos

127 f unc t i on obj = c l ea rPosH i s to ry ( obj )obj . pos_h = [ ] ;

129 obj . pos_meas_h = [ ] ;end % f � c l ea rPosH i s to ry

131

f unc t i on z = getDepth ( obj )133 z = obj . pos (3 ) ;

end % f � c l ea rPosH i s to ry135

f unc t i on pos_h = getPosH ( obj )137 pos_h = obj . pos_h ;

end % f � c l ea rPosH i s to ry139

f unc t i on pos_h = getMeasPosH ( obj )141 pos_h = obj . pos_meas_h ;

end % f � c l ea rPosH i s to ry143

f unc t i on pos_h = getPredPosH ( obj )145 pos_h = obj . pos_pred_h ;

end % f � getPredPosH147

f unc t i on pos_h = getPredPH ( obj )149 pos_h = obj . pred_P_h ;

end % f � getPredPosH151

f unc t i on pos_h = getUpdPH( obj )153 pos_h = obj .upd_P_h;

end % f � getPredPosH155

f unc t i on obj = se tF i shNo i s e ( obj , fish_noise_mag )157 obj . fish_noise_mag = fish_noise_mag ;

end159

f unc t i on obj = setUsb lNo i se ( obj , usbl_noise_mag )161 obj . usbl_noise_mag = usbl_noise_mag ;

end163

165 f unc t i on obj = in i tS i s t emMatr ix ( obj )% A phys i c s based

167 % model f o r how we expect the Fish to behave% needs to beupdated because o f p o s s i b l y va r i a b l e dt

169

obj .A = [1 0 0 obj . dt 0 0 (1/2∗ obj . dt ^2) 0 0;

171 0 1 0 0 obj . dt 0 0 (1/2∗ obj . dt ^2) 0;

0 0 1 0 0 obj . dt 0 0 (1/2∗ obj . dt^2) ;

173 0 0 0 1 0 0 obj . dt 0 0;

0 0 0 0 1 0 0 obj . dt 0;

175 0 0 0 0 0 1 0 0 obj . dt;

0 0 0 0 0 0 1 0 0;

177 0 0 0 0 0 0 0 1 0;

0 0 0 0 0 0 0 0 1] ;

179

end % f � i n i tS i s t emMatr ix181

183

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

123

185

f unc t i on [ obj , pos ] = pred i c tNextPos i t i on ( obj , v_dt )187

obj . dt = v_dt ;189 obj = obj . in i tS i s t emMatr ix ( ) ;

191 % Pred i c t next s t a t e o f the f i s h with the l a s t s t a t e and pred i c t ed motion .obj . Q_estimate = obj .A ∗ obj . Q_estimate ;

193 %pred i c t next covar iance

195

% Ex convert the proce s s no i s e ( stdv ) in to covar iance matrix197 Ex = [ obj . dt^4/4 0 0 obj . dt^3/2 0 0 obj . dt^2/2 0 0

; . . .0 obj . dt^4/4 0 0 obj . dt^3/2 0 0 obj . dt^2/2 0

; . . .199 0 0 obj . dt^4/4 0 0 obj . dt^3/2 0 0 obj . dt

^2/2; . . .obj . dt^3/2 0 0 obj . dt^2 0 0 obj . dt 0 0

; . . .201 0 obj . dt^3/2 0 0 obj . dt^2 0 0 obj . dt 0

; . . .0 0 obj . dt^3/2 0 0 obj . dt^2 0 0 obj . dt

; . . .203 obj . dt^2/2 0 0 obj . dt 0 0 1 0 0 ;

. . .0 obj . dt^2/2 0 0 obj . dt 0 0 1 0 ;

. . .205 0 0 obj . dt^2/2 0 0 obj . dt 0 0 1

] . ∗ obj . fish_noise_mag ^2;

207

obj .P = obj .A ∗ obj .P ∗ obj .A’ + Ex ;209

pos = obj . Q_estimate ( 1 : 3 ) ;211

% r e g i s t e r est imated po s i t i o n213 %obj . pos_h = [ obj . pos_h , pos ] ;

215 obj . pos = pos ;

217 obj . pos_pred_h = [ obj . pos_pred_h , pos ] ;

219 pred_P = [ sq r t ( obj .P(1 , 1 ) ) ; s q r t ( obj .P(2 , 2 ) ) ; s q r t ( obj .P(3 , 3 ) ) ] ;

221 obj . pred_P_h = [ obj . pred_P_h , pred_P ] ;

223

end % f � e s t imateNextPos i t i on225

227 f unc t i on [ obj , pos ] = updateNextPosit ion ( obj , pos_meas , v_time )

229

obj . time = v_time ;231

% r e g i s t e r measured p o s i t i o n s233 obj . pos_meas_h = [ obj . pos_meas_h , pos_meas ] ;

235 Q_meas = pos_meas ;

237 % pred i c t ed USBL measurement covar iance% Ez convert the measurement no i s e ( stdv ) in to covar iance

239 % matrixobj . Ez = [1 0 0 ; . . .

241 0 1 0 ; . . .0 0 1 ] . ∗ obj . usbl_noise_mag ^2;

243

% Kalman Gain245 K = obj .P∗ obj .C’∗ inv ( obj .C∗ obj .P∗ obj .C’+ obj . Ez ) ;

247 % Update the s t a t e e s t imate .

124

obj . Q_estimate = obj . Q_estimate + K ∗ (Q_meas � obj .C ∗ obj . Q_estimate ) ;249

% update covar iance e s t imat i on .251 obj .P = ( eye (9 )�K∗ obj .C) ∗ obj .P ;

253 pos = obj . Q_estimate ( 1 : 3 ) ;

255 % r e g i s t e r est imated po s i t i o nobj . pos_h = [ obj . pos_h , pos ] ;

257

obj . pos = pos ;259

obj . last_update_pos = pos ; % to keep track o f v e l o c i t y at beg inning o f t h i sprocedure

261

263 upd_P = [ sq r t ( obj .P(1 , 1 ) ) ; s q r t ( obj .P(2 , 2 ) ) ; s q r t ( obj .P(3 , 3 ) ) ] ;obj .upd_P_h = [ obj .upd_P_h, upd_P ] ;

265

end % f � updateNextPos it ion267

269 end

271 end

./matlab/fish_model.m

C.9 detection

c l a s s d e f d e t e c t i on2

p r op e r t i e s ( GetAccess = ’ p r i va t e ’ , SetAccess = ’ p r i va t e ’ )4 event % event that got detec ted

measurement % measured po s i t i o n6 p r e c i s i o n % p r e c i s i o n o f measurement . The lower the be t t e r

upd_pos % po s i t i o n es t imate a f t e r p r ed i c t i on and update s tep8 det_id % detec ted id

10 end

12 methods

14 f unc t i on obj = de t e c t i on ( p_event , p_meas , p_prec is ion , p_upd_pos , p_id )

16 obj . event = p_event ;obj . measurement = p_meas ;

18 obj . p r e c i s i o n = p_prec i s ion ;obj . upd_pos = p_upd_pos ;

20 obj . det_id = p_id ;

22 end % f � Constructor

24

f unc t i on e r r = getMeasError ( obj )26

i f obj . event . getType ( ) == 0 % i f event o f type s i n g l e event28

p = obj . event . g e tPo s i t i o n s ( ) ;30 e r r = norm ( ( p ( : , 1 ) ) � obj . measurement ) ; % d i s t ance between measured

po s i t i o n and r e a l f i s h e s p o s i t i o n during event

32 e l s ee r r = �1; % in case o f c o l i s i o n and no detec t i on , measurement i s not

usab le34

end36

end % f � getMeasError

125

38

40 f unc t i on e r r = getUpdError ( obj )

42 i f obj . event . getType ( ) == 0 % i f event o f type s i n g l e event

44 p = obj . event . g e tPo s i t i o n s ( ) ;e r r = norm ( ( p ( : , 1 ) ) � obj . upd_pos ) ; % d i s t ance between measured po s i t i o n

and r e a l f i s h e s p o s i t i o n during event46

e l s e48 % TODO: try d i s t anc e va lue c a l c u l a t i o n in case o f a c o l i s i o n

e r r = �1; % in case o f c o l i s i o n and no detec t i on , measurement i s notusab le

50

end52

end % f � getUpdError54

56 f unc t i on r = getEvent ( obj )r = obj . event ;

58 end % f � getEvent

60 f unc t i on r = getMeasurement ( obj )r = obj . measurement ;

62 end % f � getMeasurement

64 f unc t i on r = ge tP r e c i s i o n ( obj )r = obj . p r e c i s i o n ;

66 end % f � g e tP r e c i s i o n

68 f unc t i on r = getDetID ( obj )r = obj . det_id ;

70 end % f � getDetID

72 f unc t i on r = getUpdPos ( obj )r = obj . upd_pos ;

74 end % f � getUpdPos

76 end % � methodsend % � c l a s s d e f

./matlab/detection.m

126

Bibliografia

[1] Maybeck, Peter S. Stochastic Models, Estimation, and Control, volume Vol. 1. Outubro 2014:http://www.cs.unc.edu/~welch/kalman/media/pdf/maybeck_ch1.pdf.

[2] Kenneth C. Balcomb III and Diane E. Claridge. A Mass Stranding of Cetaceans caused by NavalSonar in the Bahamas. Outubro 2014: http://www.bahamaswhales.org/resources/Stranding_Article.pdf.

[3] Committee on Potential Impacts of Ambient Noise in the Ocean on Marine Mammals, NationalResearch Council. Ocean Noise and Marine Mammals. The National Academies Press, 2003. Outubro2014: http://users.ece.utexas.edu/~ling/2A_US1.pdf.

[4] The Society for Marine Mammalogy. Guidlines for the Treatment of Marine Mammals in Field Re-search. Outubro 2014: http://www.marinemammalscience.org/images/stories/file/Ethics/Ethics%20Guidelines.pdf.

[5] Hwaii Institute of Marine Biology. Shark & Reef Fish Research. Outubro 2014: http://www.hawaii.edu/himb/ReefPredator/Tools.htm.

[6]

[7] Pacific Salmon Foundation. Outubro 2014: http://www.psf.ca/.

[8]

[9] Pacific Ocean Shelf Tracking. Outubro 2014: http://kintama.com/wp-content/uploads/2010/09/BPA-2005-Report.pdf.

[10] Ocean Tracking Network. Outubro 2014: http://oceantrackingnetwork.org/.

[11] Australian Acoustic Tracking And Monitoring System. Outubro 2014: http://imos.org.au/aatams.html.

[12] Fish Aggregating Devices as Instrumented Observatories of pelagic ecosystems. Outubro 2014:http://www.fadio.ird.fr/.

[13]

[14] Zaihan Jiang. Underwater Acoustic Networks – Issues and Solutions. International Journal ofIntelligent Ccontrol and Systems, Vol. 13(3), Setembro 2008. Outubro 2014: http://www.ezconf.net/IJICS/files/134/01-Zaihan-Jiang-IJCS-10pages.pdf.

[15]

[16]

[17]

[18] Lurton,X. An Introduction to Underwater Acoustics, Principles and Applications. Springer, 2010.Outubro 2014: http://www.beck-shop.de/fachbuch/leseprobe/9783540784807_Excerpt_001.pdf.

128

[19] Fundamentals of underwater sound. The International Association of Oil and Gas Producers, OGPPublications, Maio 2008. Report No: 406. Outubro 2014: http://www.ogp.org.uk/pubs/406.pdf.

[20] Yao Wang. Sampling and Interpolation. EE3414, Multimedia Communication Systems - 1, Fevereiro2006. Outubro 2014: http://eeweb.poly.edu/~yao/EE3414/sampling.pdf.

[21] Saeed V. Vaseghi. Advanced Digital Signal Processing and Noise Reduction, Chapter 2. WILEY,3rd edition edition, 2006. Outubro 2014: http://dea.brunel.ac.uk/cmsp/Home_Saeed_Vaseghi/Chapter03-Fourier.pdf.

[22] Charan Langton. Linear Time Invariant (LTI) Systems and Matched Filter. ComplexToReal.com.Outubro 2014: http://complextoreal.com/wp-content/uploads/2013/01/mft.pdf.

[23] Matthias Viertler Michael Hraschan. FIR vs. IIR Filter.

[24] Digital biquad filter. Artigo Wikipedia.com. Outubro 2014: http://en.wikipedia.org/wiki/Digital_biquad_filter.

[25] Gaetano Scarano Giovanni Jacovitti. Discrete Time Techniques for Time Delay Estima-tion. 1993. Outubro 2014: ftp://p79.apl.washington.edu/incoming/williams/InSAS/Time%20Delay%20Estimation/giovani_discrete-time-tde.pdf.

[26] João Nuno Silva Picão Oliveira. Desenvolvimento de um Sistema de Posicionamento Acústico USBLe Validação com Testes de Mar. Master’s thesis, Instituto Superior Técnico, Universidade Técnicade Lisboa, Outubro 2009. Outubro 2014: https://dspace.ist.utl.pt/bitstream/2295/571912/1/tese_final.pdf.

[27] Thomas C. Austin. The Application of Spread Spectrum Signaling Techniques to Underwater Acous-tic Navigation. 1994. Outubro 2014: http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=00518658.

[28] Virtium. Data Coding And Error Checking Techniques. Whitepaper. Outubro2014: http://www.virtium.com/Whitepaper/Virtium%20-%20Data%20Coding%20and%20Error%20Checking%20Techniques%20-%20VT-WP-001.pdf.

[29] Y. T. Chan, K. C. Ho. A Simple and Efficient Estimator for Hyperbolic Location. IEEE Transactionson Signal Processing, Vol. 42(8), Agosto 1994. Outubro 2014: http://www.vtvt.ece.vt.edu/research/references/uwb/ranging_mobile_location/estimator_hyperbolic_location.pdf.

[30] Julius O. Smith, Jonathan S. Abel. Closed-Form Least-Squares Source Location Estimation fromRange-Difference Measurements. IEEE Transactions on Acoustics, Vol. ASSP-35(12), Dezembro1987.

[31] Julius O. Smith, Jonathan S. Abel. The Spherical Interpolation Method of Source Localization. IEEEJournal Of Oceanic Engineering, Vol. OE-12(1), Janeiro 1987. Outubro 2014: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1145217&tag=1.

[32] Carlos Jorge Silvestre Marco Martins Morgado, Paulo Jorge Oliveira. Sistema De Navegação Inercialcom Ajuda USBL, Relatório Técnico, Março 2006. Outubro 2014: http://users.isr.ist.utl.pt/~pjcro/RUMOS/RT1.pdf.

[33] Jonathan S. Abel, Khosrow Lashkari. Track Parameter Estimation from Multipath Delay Informa-tion. IEEE Journal of Oceanic Engineering, Vol. OE-12(1), Janeiro 1987.

[34] Alonzo Kelly. A 3D State Space Formulation of a Navigation Kalman Filter for Autonomous Vehicles.Fevereiro 2006. Outubro 2014: http://www.ri.cmu.edu/pub_files/pub1/kelly_alonzo_1994_6/kelly_alonzo_1994_6.pdf.

[35] Morgado M. Erzini K. Bentes L. Afonso A. Hazin F. Block B. Tuna Research Oliveira P., Silvestre C.and Conservation Center. Preliminary Results from Project MAST/AM- Advanced Tracking andTelemetry Methodologies to Study Marine Animals.

[36] GeomAlgorithms.com. Outubro 2014: http://geomalgorithms.com/a07-_distance.html.

129

[37] Cai G., Chen B.M., Lee T.H. Unmanned Rotorcraft Systems. Springer, 2011. Chapter 2. CoordinateSystems and Transformations.

[38] Bruno Miguel Simões Carvalho Cardeira. Arquitecturas para Navegação Inercial/GPS com Aplicaçãoa Veículos Autónomos. Fevereiro 2009.

[39]

130