( PDF ) Estimativa do expoente de hurst, por meio da transformada wavelet, de séries temporais de precipitação de chuvas das regiões climáticas do Estado de São Paulo no período de 1978 a 1997

Download PDF

ads:

UNIVERSIDADE ESTADUAL PAULISTA

Instituto de Geociências e Ciências Exatas

Departamento de Física

Campus de Rio Claro

Estimativa do Expoente de Hurst, por meio da Transformada Wavelet, de

Séries Temporais de Precipitação de Chuvas das Regiões Climáticas do

Estado de São Paulo no período de 1978 a 1997.

Luiz Roberto Salomão

Orientador: Prof. Dr. José Roberto Campanha

Dissertação de Mestrado elaborada

junto ao Programa de Pós-Graduação

em Física - Área de concentração em

Física Aplicada, para obtenção do

título de Mestre em Física.

Rio Claro (SP)

2006

ads:

Livros Grátis

http://www.livrosgratis.com.br

Milhares de livros grátis para download.

517.355 Salomão, Luiz Roberto

S173e Estimativa do expoente de Hurst, por meio da transforma-

da wavelet, de séries temporais de precipitação de chuvas das

regiões climáticas do estado de São Paulo no período de l978

a 1997 / Luiz Roberto Salomão. – Rio Claro : [s.n.], 2006

110 p. : il., figs., gráfs., tabs., mapas

Dissertação (mestrado) – Universidade Estadual Paulista,

Instituto de Geociências e Ciências Exatas

Orientador: José Roberto Campanha

1. Fourier, Transformações de. 2. Sistema complexo. 3. Di-

mensão fractal. 4. Persistência e antipersistência. 5. Movimen-

to browniano fracionário. 6. Espectro de potência wavelet.

7. Fourier, Transformada de. I. Título.

Ficha Catalográfica elaborada pela STATI – Biblioteca da UNESP

Campus de Rio Claro/SP

ads:

UNIVERSIDADE ESTADUAL PAULISTA

Instituto de Geociências e Ciências Exatas

Departamento de Física

Campus de Rio Claro

Estimativa do Expoente de Hurst, por meio da Transformada Wavelet, de

Séries Temporais de Precipitação de Chuvas das Regiões Climáticas do

Estado de São Paulo no período de 1978 a 1997.

Aluno: Luiz Roberto Salomão

Orientador: Prof. Dr. José Roberto Campanha

Comissão Examinadora

Prof. Dr. José Roberto Campanha - Orientador

IGCE/UNESP/Rio Claro

Prof. Dr. Hari Mohan Gupta

IGCE/UNESP/Rio Claro

Prof. Dr. Roberto Nicolau Onody

IFSC/USP/São Carlos

Rio Claro, 26 de outubro de 2006.

Resultado: Aprovado.

iii

A minha esposa Sueli e a meus filhos Rafael e Maria Carolina.

Agradecimentos

Ao Prof. Dr. José Roberto Campanha pela atenção, dedicação e paciência que

sempre marcou a sua orientação e pela amizade que permitiu que este trabalho fosse realizado

de forma agradável.

Aos Prof. Dr. Hari Mohan Gupta e Prof. Dr. Roberto Nicolau Onody pelas

sugestões que permitiram o aprimoramento deste trabalho.

Aos professores e funcionários do Departamento de Física do Instituto de

Geociências pela amizade e carinho com que me receberam.

Aos amigos da Academia da Força Aérea pelo incentivo.

Aos companheiros de estudo Rosana, Natale, Luciano e Shinaider pela colaboração

e apoio.

A minha esposa por suportar com amor e carinho todas as minhas manias.

A todos muito obrigado.

SUMÁRIO

Índice --------------------------------------------------------------------------------------------- vi

Índice de figuras--------------------------------------------------------------------------------viii

Índice de tabelas-------------------------------------------------------------------------------- xii

Resumo-------------------------------------------------------------------------------------------xiii

Abstract ------------------------------------------------------------------------------------------xiv

ÍNDICE

CAPÍTULO 1 -------------------------------------------------------------------------------------- 01

Objetivo--------------------------------------------------------------------------------------- 01

1.1 - Introdução ------------------------------------------------------------------------------ 02

1.2 - Características climáticas do estado de São Paulo -------------------------------- 04

CAPÍTULO 2 -------------------------------------------------------------------------------------- 10

2.1 - Série de Fourier------------------------------------------------------------------------ 10

2.2 - Integral de Fourier--------------------------------------------------------------------- 11

2.3 - Transformada de Fourier ------------------------------------------------------------- 12

2.4 - Transformada wavelet ---------------------------------------------------------------- 17

2.5 - Escala e Deslocamento --------------------------------------------------------------- 18

2.6 - Tipos de wavelet (mais comuns) ---------------------------------------------------- 20

2.7 - Aplicação da transformada wavelet------------------------------------------------- 23

CAPÍTULO 3 -------------------------------------------------------------------------------------- 28

3.1 - Fractais---------------------------------------------------------------------------------- 28

3.2 - Dimensão ------------------------------------------------------------------------------- 30

3.3 - Poeira de Cantor ----------------------------------------------------------------------- 31

3.4 - Triângulo de Sierpinsky -------------------------------------------------------------- 33

3.5 - Contagem de caixas ------------------------------------------------------------------- 33

3.6 - Expoente de Hurst--------------------------------------------------------------------- 34

3.7 - Auto-Afinidade ------------------------------------------------------------------------ 40

3.8 - Estimativa de H para série temporal auto-afim------------------------------------ 42

vii

CAPÍTULO 4 -------------------------------------------------------------------------------------- 44

4.1 - Validade do expoente de Hurst estimado por wavelet---------------------------- 45

4.2 - Expoente de Hurst e o tamanho da série temporal -------------------------------- 52

4.3 - Ausência de medições em séries temporais de precipitação de chuvas -------- 54

4.4 - Escolha da wavelet mais adequada para estimar o expoente de Hurst---------- 65

CAPÍTULO 5 -------------------------------------------------------------------------------------- 67

5.1 - Obtenção e tratamento dos dados pluviométricos--------------------------------- 67

5.2 - Estimativa do expoente de Hurst para cada série de precipitação--------------- 68

5.3 - Séries com os períodos seco e úmido bem definidos ----------------------------- 79

5.4 - Espectro de potência wavelet para as séries temporais de chuva---------------- 81

CAPÍTULO 6 -------------------------------------------------------------------------------------- 89

Referências Bibliográficas----------------------------------------------------------------------- 93

Apêndice I ------------------------------------------------------------------------------------------ 100

Apêndice II ----------------------------------------------------------------------------------------- 101

Apêndice III ---------------------------------------------------------------------------------------- 102

Apêndice IV ---------------------------------------------------------------------------------------- 103

Apêndice V------------------------------------------------------------------------------------------ 105

Apêndice VI ---------------------------------------------------------------------------------------- 106

Apêndice VII --------------------------------------------------------------------------------------- 108

viii

ÍNDICE DE FIGURAS

Figura 1.1 - Mapa do estado de São Paulo com as regiões climáticas

(MONTEIRO, 1973).

Figura 2.1 - Gráfico da função f

(t) para valores de t próximos de 5. 13

Figura 2.2 - Gráfico da função f

(t) para valores de t próximos de 5. 14

Figura 2.3 - a) Gráfico da transformada de Fourier da função f

b) Gráfico da transformada de Fourier da função f

Figura 2.4 - Esquema simplificado de janelamento. 16

Figura 2.5 - Deslocamento da wavelet "chapéu mexicano". 19

Figura 2.6 - Mudanças de escala para a wavelet "chapéu mexicano". 19

Figura 2.7 - Escalogramas das transformadas wavelet das funções f

e f

função da escala.

Figura 2.8 - Escalogramas das transformadas wavelet das funções f

e f

função da freqüência.

Figura 2.9 - a) Espectro de potência em função do período para a função f

b) Espectro de potência em função do período para a função f

Figura 3.1 - Possível linha costeira de um país em três escalas diferentes. 29

Figura 3.2 - Representação da escala para um segmento de reta (a), para uma

figura plana (b) e para um sólido (c).

Figura 3.3 - Representação do conjunto poeira de Cantor 32

Figura 3.4 - Cálculo da dimensão da figura poeira de Cantor. 32

Figura 3.5 - Representação do triângulo de Sierpinsky. 33

Figura 3.6 - Quadriculado com caixas de lado δ cobrindo a figura que se quer

calcular a dimensão fractal

Figura 3.7 - Esquema de um reservatório de água (FEDER, 1989). 36

Figura 3.8 - Exemplos de movimento browniano fracionário gerado pelo

algoritmo de Joachim Johansson (Apêndice III).

Figura 3.9 - Variação da correlação de um movimento browniano fracionário

em função do expoente de Hurst.

Figura 3.10 - Quadriculado usado para cálculo da dimensão fractal de processo

auto-afim.

Figura 4.1 - Gráfico com os primeiros passos do algoritimo mpdwalk. 46

Figura 4.2 - Gráficos de chuva cumulativa dos desvios da média e série fbm

gerada com o mesmo H que a série de chuva.

Figura 4.3 - a) Desvios da média e histograma para uma série de chuva com

média 0, e variância 0,010.

b) Desvios da média e histograma para uma série de fbm com

média 0 e variância 0,058.

Figura 4.4 - Gráficos mostrando o valor médio estimado do expoente de Hurst

para séries de movimento browniano obtidas por meio de quatro

tipos de geradores de movimento browniano fracionário.

Figura 4.5 - a) Gráfico do valor médio de H para cem séries com número de

elementos entre 10 e 500.

b) Gráfico do valor médio de H para cem séries com número de

elementos entre 500 e 8192.

Figura 4.6 - a) Gráfico de precipitação de chuvas no município de Pilar do

Sul, posto Usina Batista no período de 01/01/1978 a

31/12/1997 (E4028).

b) Gráfico da soma cumulativa dos dados de precipitação de

chuvas da série E4028 (corrigidos pela média e pela reta que

passa no primeiro e último ponto).

Figura 4.7 - a) Medidas substituídas em posições aleatórias.

b) Medidas substituídas seguidamente a partir da posição 2796

no sentido do final da série.

c) Medidas substituídas seguidamente a partir da posição 5973

no sentido do início da série.

Figura 4.8 - a) Gráfico de precipitação de chuvas do município de Cananéia,

posto Morro Redondo (G5002).

b) Gráfico de precipitação de chuvas do município de Eldorado,

posto Ouro Leve (F5020).

c) Gráfico de H com o número de elementos usados para o

cálculo da média (G5002).

d) Gráfico de H com o número de elementos usados para o

cálculo da média (F5020).

Figura 4.9 - a) Medidas retiradas em posições aleatórias

b) Medidas retiradas seguidamente a partir da posição 5973 em

direção ao início da série

Figura 4.10 - Valor médio do expoente de Hurst em função da wavelet usada

para estimá-lo.

Figura 5.1 - a) Gráfico de uma série temporal de precipitação de chuvas.

b) Gráfico de uma série temporal de precipitação de chuvas

corrigida pela média. Para esta série a média das precipitações

considerando os 20 anos é de 3,6260 mm e temos 5491 dias

com precipitação zero.

Figura 5.2 - a) Gráfico da soma cumulativa de uma série temporal de

precipitação de chuvas.

b) Gráfico da soma cumulativa de uma série temporal de

precipitação de chuvas corrigida pela reta que passa pelo

primeiro e último ponto.

Figura 5.3 - Valor médio de H para cada região e sub-região climática do

estado de São Paulo.

Figura 5.4 - a) gráfico de precipitação de chuvas da cidade de Assaré.

b) gráfico de precipitação de chuvas da cidade de Cratos.

Figura 5.5 - Espectro de potência de séries de precipitação de chuvas da

região climática II (a, b e c).

Figura 5.6 - Espectro de potência de série de precipitação de chuvas da região

climática III.

Figura 5.7 - Espectro de potência de séries de precipitação de chuvas da

região climática V.

Figura 5.8 - Espectro de potência de série de precipitação de chuvas da região

climática VII.

Figura 5.9 - Espectro de potência de séries de precipitação de chuvas da

região climática VIII.

Figura 5.10 - Espectro de potência de séries de precipitação de chuvas da

região climática IX.

Figura VII-1 Iterações da função de dilatação. 109

Figura VII-2 Wavelet Daubechies de ordem 8 construída após oito iterações da

função de dilatação.

110

xii

ÍNDICE DE TABELAS

abela 1.1 - As regiões climáticas segundo as unidades geomorfológicas

(MONTEIRO, 1973)

abela 1.2 - Resumo do clima das regiões do estado de São Paulo e as

cidades onde foram coletados os dados de chuva.

abela 4.1 - Média dos valores do expoente de Hurst estimados por wavelet

das séries de movimento browniano fracionário. Com cada

algoritmo utilizado para simular fbm foram geradas 100 séries

para cada valor de H. Epm é o erro padrão da média (VIEIRA,

1999).

abela 4.2 - Valores de H usando a média das medidas anteriores para as

ausências de medições.

abela 4.3 - Valores de H usando a união dos trechos sem ausências de

medições.

abela 5.1 - Valores de H para as séries da região climática I 71

abela 5.2- Valores de H para as séries da região climática II 72

abela 5.3 - Valores de H para as séries da região climática III 73

abela 5.4- Valores de H para as séries da região climática IV 73

abela 5.5 - Valores de H para as séries da região climática V 74

abela 5.6 - Valores de H para as séries da região climática VI 75

abela 5.7 - Valores de H para as séries da região climática VII 75

abela 5.8 - Valores de H para as séries da região climática VIII 76

abela 5.9 - Valores de H para as séries da região climática IX 77

abela 5.10 - Valores de H para as séries da região nordeste. 80

xiii

RESUMO

Este trabalho pode ser separado em duas partes. Na primeira apresentamos,

de forma resumida, a natureza fractal das séries de precipitação de chuvas, a estimativa do

expoente de Hurst (H), a transformada wavelet, o seu uso na estimativa de H e a

comprovação, por meio de testes com séries de movimento browniano fracionário geradas

com H estabelecido a priori, que temos um bom método para estimar H com a transformada

wavelet. Na segunda parte é feita uma análise das séries temporais de precipitação de chuvas

(séries de chuva). As séries de chuvas foram obtidas junto ao Departamento de Águas e

Energia Elétrica (DAEE) do estado de São Paulo e são constituídas de medições diárias, em

milímetros, da quantidade de chuva em postos meteorológicos distribuídos em todo o estado.

Algumas séries apresentam ausências de medições que acabam transformando uma série

longa em duas ou mais séries menores. Em função disso foi estabelecido um procedimento

que permite aproveitar essas séries com ausências (dias sem medições) para estimar o

expoente de Hurst. A partir de considerações climáticas o estado foi dividido em nove

regiões. Foram selecionadas séries com 20 anos de registros de chuva (de janeiro de 1978 a

dezembro de 1997). Dos valores de H encontrados foi feita uma média para cada região

climática do estado de São Paulo.

Palavras chaves: sistema complexo, série temporal, chuva, movimento browniano fracionário,

wavelet, expoente de Hurst.

xiv

ABSTRACT

This work can be separated into two parts. In the first part is briefly

presented the fractal nature of the rain precipitation series, the estimate of the Hurst exponent

(H), the wavelet transform, its use in the estimate of H and the evidence, by tests with

generated series of fractional brownian motion with H established a priori, that we have a

good method to estimate H with wavelet transform. In the second part an analysis of rain

precipitation series is made. The rain series had been obtained from the Departmento de

Águas e Energia Elétrica (DAEE) of the state of São Paulo and contain daily measurements,

in milimeters, of the amount of rain in distributed meteorological stations all over the state.

Some series present absences of measurements, transforming a long series into two or more

lesser series. A procedure was established, allowing to use these series with absences (days

with no measurements) to estimate the Hurst exponent. From climatic considerations the state

was divided into nine regions. The regions with rain data over a period of 20 years have been

selected as rain series (from january 1978 to december 1997). From the several values of H,

we computed an averaged value of H for each climatic regions of the state of São Paulo.

Keyswords: complex system, temporal series, rain, fractional Brownian motion, wavelet,

Hurst exponent.

CAPÍTULO 1

Objetivo

A análise de séries temporais de precipitação de chuvas tem importância

fundamental em termos de previsão meteorológica e também em termos de estabelecer as

características do clima de uma determinada região.

Neste trabalho, iremos estimar o expoente de Hurst das séries temporais de

precipitação de chuvas das regiões climáticas do estado de São Paulo no período de janeiro de

1978 a dezembro de 1997.

Para a estimativa do expoente de Hurst usaremos o método da transformada

wavelet aplicado às séries de precipitação de chuvas. Além disso, discutiremos a precisão do

método ao estimar o expoente de Hurst (H), a adequação deste método para a análise de séries

de chuvas e o aproveitamento das séries de chuvas com interrupções de suas medidas para a

estimativa do expoente de Hurst.

1.1 - Introdução

Na leitura da previsão do tempo em um jornal normalmente estamos

interessados na possibilidade de chuva durante o dia. Costumamos ignorar os outros dados

presentes, tais como, pressão, temperatura, umidade, ventos, deslocamento de frentes polares

e ou continentais entre outras. Na previsão é levada em conta uma série de variáveis

atmosféricas, posição geográfica e relevo.

A precipitação de chuvas faz parte de um sistema chamado clima que pela

quantidade de variáveis envolvidas e pelo relacionamento não linear entre elas

(GLEICK,1990), podemos dizer que se trata de um sistema complexo. De um lado temos uma

parte cíclica que apresenta as estações climáticas com o clima se repetindo, com mais ou

menos intensidade, todo ano e por outro lado a incerteza do momento em que acontecerão os

fenômenos climáticos como, por exemplo, a chuva.

Apesar da aparente aleatoriedade da chuva estudos recentes têm encontrado

certa 'ordem' sob a forma de uma invariância de escala, certos padrões de precipitação se

repetem, em escala maior ou menor, ao longo de uma série temporal de chuvas. Esta

invariância de escala também é encontrada na Geometria Fractal como um dos seus

fundamentos (LOVEJOY e MANDELBROT, 1985).

As séries temporais de precipitação de chuvas foram estudadas por

PETERS, HERTLEIN e CHRISTENSEN, 2002, OLSSON, NIEMCZYNOWICZ e

BERNDTSSON, 1993, RODRIGUEZ, ITURBE e DE POWER, 1989, SVENSSON,

OLSSON, BERNDTSSON, 1996, que as caracterizaram como um sistema complexo por

possuírem propriedades fractais e ou multifractais.

Na natureza encontramos várias estruturas fractais e podemos citar como

exemplos as árvores, a folha da samambaia, a linha da costa erodida pela água e pelos ventos

e o contorno das nuvens. Os processos temporais também apresentam essa similaridade

estatística em várias escalas e temos como exemplo a ocorrência de terremotos, a variação de

caudais de rios, o aparecimento de manchas solares e as séries temporais de precipitação de

chuvas.

Harold E. Hurst, engenheiro inglês, trabalhou no projeto de construção de

uma represa no rio Nilo, a partir de 1906 e um dos seus trabalhos foi determinar como seria a

liberação de água do reservatório para não secar e também não transbordar. Na análise de

dados das cheias do rio Nilo, desenvolveu um método estatístico chamado Análise de R/S

(HURST, BLACK, SIMAIKA, 1965) que deu origem ao expoente de Hurst.

A estimativa do expoente de Hurst (H) nos fornece informações sobre

correlação e persistência (PETERS, 1996), que é um excelente indicador para o caso de séries,

como a de precipitação de chuvas, que se situam entre uma série periódica e uma aleatória.

O expoente de Hurst de séries temporais de chuvas descreve a irregularidade

ou a aleatoriedade presente na série e também se há ou não efeitos de longo alcance presente,

contribuindo, assim, para conhecimento de como serão as séries no futuro e de como o clima

está ou não se modificando (BRESLIN e BELWARD, 1999).

Existem vários estimadores do expoente de Hurst, a análise de R/S, método

da transformada de Fourier, método variacional, método do periodograma, método da

transformada wavelet, entre outros, e cada tem os seus pontos fortes e fracos. Alguns

subestimam e outros superestimam o valor de H e existem, também, aqueles que apresentam

um erro muito grande. Neste trabalho optamos pelo método da transformada wavelet por

apresentar, nos testes feitos (capítulo quatro) com séries de movimento browniano fracionário

com o valor de H estabelecido a priori, um bom resultado para todos os valores de H.

As séries temporais analisadas neste trabalho são séries de precipitação de

chuvas, resultado de medições em postos meteorológicos distribuídos por todo o estado de

São Paulo, fornecidas pelo Departamento de Águas e Energia Elétrica do estado de São Paulo.

As séries apresentam medições em períodos diferentes e algumas apresentam ausências de

medições por períodos de alguns dias a anos. No capítulo quatro apresentamos um teste para

verificar o possível aproveitamento destas séries, com ausências de medições, para a

estimativa do expoente de Hurst.

Na escolha das séries de chuva estabelecemos que todas as séries deveriam

ser analisadas para um mesmo período de tempo (capítulo cinco), assim todas as séries

estariam sob a influência dos mesmos fatores determinantes do clima do estado de São Paulo.

Para determinar o tamanho das séries de chuvas a serem analisadas levamos em conta a

abundância de séries longas em cada região climática do estado de São Paulo (se o período de

medição fosse muito grande o número de séries diminuía rapidamente) e a variação de H em

função do número de elementos da série (capítulo quatro).

Ainda no capítulo cinco foi feita uma comparação das séries de chuvas do

estado de São Paulo com séries de chuvas do estado do Ceará que apresentam períodos seco e

úmido bem definidos e uma análise do espectro de potência wavelet que mostra além do

período anual um segundo e terceiro períodos em algumas séries.

1.2 - Características climáticas do estado de São Paulo

O estado de São Paulo localiza-se, em termos de clima, entre a zona

intertropical, caracterizada pela predominância de massas tropicais e equatoriais, e a zona

subtropical controlada pelas massas tropicais e polares. As variações das grandes áreas de

circulação de ar, criadas pela rotação da Terra e pela transferência de calor do Equador para

os pólos, afetam diretamente o clima no estado.

Esta localização permite segundo MONTEIRO, 1973, separar três climas

regionais determinados, principalmente, pela quantidade anual de chuvas: o sudeste e noroeste

permanentemente úmidos e o norte com período seco bem definido. Ainda em função da

distribuição diária das chuvas as regiões foram subdivididas classificando em nove regiões o

clima do estado.

Na figura 1.1 podemos ver a separação em três climas regionais, A

, A

e B,

separados pela linha vermelha e as nove regiões, I, II, III, IV, V, VI, VII, VIII e IX separadas

pela linha preta. As regiões aparecem separadas por uma faixa branca, uma vez que, é

impossível estabelecer uma linha para a separação do clima. Na tabela 1.1 temos a

caracterização das regiões climáticas segundo as unidades geomorfológicas.

Devido ao alto índice pluviométrico existente no território paulista como um

todo, as diferenças entre as regiões se reduzem praticamente às variações que ocorrem no

período outono-inverno, quando as variações de abastecimento e atividades das massas

polares são maiores. Nos anos em que há grande atividade polar o estado de São Paulo recebe

grande quantidade de chuva e quando há pouca atividade polar os sistemas intertropicais se

sobressaem e os períodos secos são maiores.

A Frente Polar Atlântica é a principal causadora das chuvas no território

paulista e a sua oscilação do rio do Prata até próximo ao Equador é que determina o tamanho

dos períodos de seca e chuva (PEREIRA JUNIOR e CHRISTOFOLETTI, 2003).

Na tabela 1.2 apresentamos um resumo do clima de cada região e as cidades

onde foram coletados os dados de chuva usados neste trabalho.

Tabela 1.1 As regiões climáticas segundo as unidades geomorfológicas (MONTEIRO,

1973)

Figura 1.1 Mapa do estado de São Paulo com as regiões climáticas (MONTEIRO, 1973).

TABELA 1.2

Região

Clima Precipitação

anual média

(1978 a 1997)

LITORAL E PLANALTO ATLÂNTICO NORTE

Clima úmido da costa, controlado por massas equatoriais e

tropicais, expostas às massas Tropicais Atlânticas.

Sub-região I:Cunha (Ferraz), São Luiz do Paraitinga

(Catucaba), Natividade da Serra (Laranjal), Natividade da

Serra, Paraibuna (Alegre), Natividade da Serra (B. Alto), São

Luiz do Paraitinga (Briet) e Cunha (Vargem do Tanque)

Sub-região Ia:Ubatuba (Picinguaba), Ilha Bela, São Sebastião

(São Francisco), Caraguatatuba (Porto Novo), Ubatuba

(Maranduba) e Cunha (Rio Manso).

1620 mm

LITORAL E PLANALTO ATLÂNTICO SUL

Clima úmido da face oriental e subtropical do continente

dominado por massa Tropical e controlado por massas polares.

Sub-região IIa: Cananéia (Morro Redondo), Iguape

(Momuna), Itariri (B. Igrejinha), Miracatu (Pedra do Largo),

Juquiá (Escalvado), Registro (Jurumirim), Registro, Eldorado,

Eldorado (Ouro Leve) e Cajati (Barra do Azeite).

Sub-região IIb: Monguaguá, São Vicente, Guarujá (Pereque) e

Santos (Caeté).

Sub-região IIc: Pilar do Sul (U. Batista), Coronel Macedo,

Itapeva, Itararé (B. Boa Vista), Barra do Chapéu, Ribeirão

Branco (Pinara), Capão Bonito (F. Santa Inês), São Miguel

Arcanjo (Taquaral), Piedade (B. Ribeirão Bonito), Guapiara (B.

Pinheiro), Itapeva (Guarizinho), Itapetininga (Eng. Hermilio) e

Itapetininga (Gramadinho).

Sub-região IId: Mogi das Cruzes (taiacupeba), Mauá, Santo

André (Paranapiacaba), São Bernardo do Campo (Rio Acima) e

São Paulo (Santo Amaro).

1764 mm

III

O VALE DO PARAÍBA

Clima úmido da costa, controlado por massas equatoriais e

tropicais, expostas às massas Tropicais Atlânticas.

Guararema (Paratei), Jacarei, São José dos Campos

(Pararangaba), Caçapava (Quirino), Taubaté (Macuco),

Guaratinguetá (Brumado), Silveiras e Lorena

1442 mm

TABELA 1.2 (continuação)

Região

Clima Precipitação

anual média

(1978 a 1997)

A MANTIQUEIRA

Clima úmido da costa, controlado por massas equatoriais e

tropicais, expostas às massas Tropicais Atlânticas.

Sub-região IVa: Piracaia (Crioulos), São José dos Campos

(Guirra), São José dos Campos (S. F. Xavier), Monteiro Lobato

(S. Benedito) e Campos do Jordão (Capivari).

Sub-região IVb: Atibaia, Pinhalzinho, Socorro, Lindóia,

Aguas de Lindóia, Espirito Santo do Pinhal (U. Pinhal), São

João da Boa Vista (F. Paraíso), Aguas da Prata (Fartura), São

José do Rio Pardo e Caconde.

1707 mm

O CENTRO-NORTE

Clima tropical alternadamente seco e úmido, controlado por

massas equatoriais e tropicais.

Sub-região Va: Mococa (Sitio Esplanada), Santa Rosa do

Viterbo (Bom Sucesso), Tambau, Santa Rita do Passa Quatro

(Usina), Porto Ferreira, Casa Branca (Venda Branca),

Pirassununga (Fernando Costa), Santa Cruz da Conceição

(S.S.Geraldo), Leme (Cresciumal), Mogi Guaçu (Capetinga),

Rio Claro (F. São José), Santa Gertrudes e Araras (Santa Cruz).

Sub-região Vb: Matão, Tabatinga (Araruba), Santa Lucia,

Araraquara (Chibarro), Bocaina, Ribeirão Bonito, Brotas

(Campo Redondo), Itirapina e Torrinha (U. Três Saltos).

Sub-região Vc: Icem, Colombia (F. Campo Grande), Guaira

(F. Antas), Ituverava, Pedregulho (Alto do Lajeado), Itirapuá,

Sales de Oliveira (F. Conquista), Morro Agudo, Colina,

Olímpia (F. Cruz Alta), Novais, Monte Alto, Sertãozinho (U.

Sta Elisa) e Altinópolis (Aguas Virtuosas).

1525 mm

A PERCÈE

(1)

DO TIETE

Clima tropical alternadamente seco e úmido, controlado por

massas equatoriais e tropicais.

Anhembi, Itu (Pirapitingui), Mairinque (D. Catarina), Boituva,

Laranjal Paulista, Capivari, Rio das Pedras e Campinas.

1361 mm

(1) Caminho natural (FREIRE, 1900, CORRÊA & STEINBERG, 1980).

TABELA 1.2 (continuação)

Região

Clima Precipitação

anual média

(1978 a 1997)

VII

A SERRA DE BOTUCATU

Clima tropical alternadamente seco e úmido, controlado por

massas equatoriais e tropicais.

Aguas de Santa Bárbara, Itatinga (Lobo), Bauru, Iacanga (F.

Barreirinho), Borborema (F. Ano Bom do Tiete), Borborema

(F. Laranjal), Santa Cruz do Rio Pardo (Sodrelia), Avaré (F. Sta

Gabriela), Paulistânia, Pardinho, Lençóis Paulista (Gleba Rio

Claro), Novo Horizonte (F. Uirapuru), Avaré, Agudos e

Reginópolis.

1422 mm

Região

Clima Precipitação

anual média

(1978 a 1997)

VIII

O OESTE

Clima tropical alternadamente seco e úmido, controlado por

massas equatoriais e tropicais.

Andradina, Alto Alegre, Andradina (Três Irmãos), Auriflama,

Brauna (F. N.S. Aparecida), Floreal, Flórida Paulista (F. Sto

André), Gastão Vidigal, Itapura, Jales, José Bonifácio, Marília

(Dirceu), Marília (Amadeu Amaral), Mirandópolis (Yuba),

Muritinga do Sul (F. Boa Vista), Osvaldo Cruz, Panorama,

Pereira Barreto (Ideal), Planalto (São Jerônimo), Pontalinda,

Rubiacéia (Caramuru), Santa Fé do Sul, Sud Mennucci, Tupã

(Varpa), Tupã, Turiuba, Ubarana (F. cataco) e Valparaíso

(Perpetuo Socorro).

1317 mm

O SUDOESTE

Clima úmido da face oriental e subtropical do continente

dominado por massa Tropical e controlado por massas polares.

Presidente Epitácio (Sul Mineira), Maracai (B. Agua do Mato),

Mirante do Paranapanema, Marabá Paulista, Santo Anastácio

(F. Santa Isabel), Narandiba, Presidente Venceslau (F.

Clotilde), Iepe (Jaguarete), Rosana (Itaporã), Teodoro Sampaio

(Santaida), Teodoro Sampaio, Rosana (Nova do Pontal),

Presidente Epitácio (Sucurita), Florinea eCandido Mota (F. Sta

Rosa).

1332 mm

Tabela 1.2 - Resumo do clima das regiões do estado de São Paulo e as cidades onde foram

coletados os dados de chuva..

CAPÍTULO 2

ANÁLISE DE SÉRIES TEMPORAIS

Neste capítulo pretende-se fazer uma breve introdução à análise de séries

temporais baseando-se em transformada de Fourier e transformada wavelet.

2.1 - Série de Fourier

Joseph Fourier, em 1807, estabeleceu que se uma função f(x) for periódica

de período 2L, definida entre -L e +L, com x ∈ R, ela pode ser expressa por meio de uma

soma de funções trigonométricas senos e cossenos. Define-se série de Fourier de f(x),

segundo TOLSTOV (1962) e MORETTIN (1999), como:

)]

(senb)

cos(a[

)x(f

∑

∞

(2.1)

onde

∫

L-

dx )

(sen f(x)

dx )

( cos f(x)

Os valores de a

e b

determinam a contribuição de cada seno e cosseno da

série na representação de f(x).

2.2 - Integral de Fourier

Na representação de f(t) por meio da série de Fourier (equação 2.1), se o

período da função tende a infinito, a série de Fourier se torna a integral de Fourier

(SPIEGEL, 1976).

Considerando que f(t) e f'(t) são seccionalmente contínuas e f(t) é

absolutamente integrável em (-∞,∞), o teorema da integral de Fourier afirma que uma função

f(x) pode se expandida em integrais de senos e cossenos:

∫

∞

ααα+αα=

d ] ) x sen( ) B( ) x cos( ) A( [ ) x ( f

onde

x d ) x cos( ) x f(

) ( A

∫

∞

=α

x d ) x sen( ) x f(

) ( B

∫

∞

=α

Aqui A e B são funções que medem a contribuição dos cossenos e senos

para a representação de f(x).

A integral de Fourier pode ser colocada em outras formas

αα

∫∫

∞

=α

∞

−∞=

ddu ]) u- x(cos[ )u f(

)x(f

ou usando a fórmula de Euler )x(seni)xcos(e

∫∫

∞

)u - x(i

d du e )u f(

) x (f

∫∫

∞

∞−

αα

ui-xi

due )u f( de

) x (f (2.2)

A segunda integral da equação 2.2 é uma função de α e será integrada

juntamente com a f(x) de menos a mais infinito. Os valores de α é que definem o argumento

das funções seno e cosseno.

2.3 - Transformada de Fourier

A transformada de Fourier é definida (TOLSTOV, 1999) por:

∫

∞

=α

u i-

du e )u f( ) (F (2.3)

e a transformada inversa por:

αα

∞

−

∫

d e ) F(

) x (f

x i

(2.4)

A transformada de Fourier separa o sinal em suas componentes (senos e

cossenos) de diferentes freqüências. Podemos também entender a transformada como uma

ferramenta matemática para transportar o sinal observado no tempo ou espaço para o domínio

da freqüência.

Em muitas situações não conhecemos a função matemática que representa o

resultado de uma série de medidas. O que temos, na realidade, são registros discretos de um

determinado fenômeno f(t) do tipo f

, f

, .....f

N-1

, com N observações amostradas em

intervalos de tempo ∆t, durante um período T. Define-se a transformada de Fourier discreta

como

∑

−

ω−

(2.5)

e a transformada inversa por

∑

−

eFf (2.6)

com

∆

=ω .

O número de operações computacionais para o cálculo da transformada de

Fourier discreta ou sua inversa, a partir das equações 2.5 e 2.6, é da ordem de N

. Isto se torna

um problema quando o N é muito grande. COOLEY e TUKEY, 1965, desenvolveram um

algoritmo chamado de FFT (Fast Fourier Transform) que reduz o número de operações para

Nlog

Como um exemplo, consideremos duas funções f

(t) e f

(t) dadas pelas

expressões abaixo;

= sen(2π10t) para t < 5 e f

= sen(2π30t) para t ≥ 5

Figura 2.1 - Gráfico da função f

(t) para valores de t próximos de 5.

4.6 4.7 4.8 4.9 5 5.1 5.2 5.3 5.45.1

-1

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

tempo

f1 (t)

= sen(2π10t) + sen(2π30t) para qualquer t

Figura 2.2 - Gráfico da função f

(t) para valores de t próximos de 5.

Aplicando a definição da transformada rápida de Fourier (equação 2.5) para

as funções f

e f

obtemos os gráficos

a) b)

Figura 2.3 - a) Gráfico da transformada de Fourier da função f

b) Gráfico da transformada de Fourier da função f

4.6 4.8 5 5.2 5.4

-2

-1.5

-1

-0.5

0.5

1.5

tempo

f2 (t)

0 10 20 30 40 50 60 70

Freqüência

Coeficiente

Transformada de Fourier de f1

0 10 20 30 40 50 60 70

100

150

200

250

Freqüência

Coeficiente

Transformada de Fourier de f2

Praticamente não há diferença entre os gráficos da transformada de Fourier

das funções f

e f

Nos dois gráficos os picos estão em 10 e 30, exatamente os valores usados

para a construção das funções f

e f

. Para algumas situações, por exemplo, em sinais

estacionários, a transformada de Fourier representa uma grande informação que é a freqüência

da onda. As funções f

e f

são compostas de combinações diferentes de senos com

freqüências numéricamente iguais a 10 e 30 e por meio dos gráficos presentes na

figura 2.3,

não há como identificar onde estão ocorrendo as variações da freqüência ao longo do tempo.

Aplicando-se a transformada de Fourier a um sinal encontramos as freqüências presentes, mas

sem relação com o tempo ou espaço. Este é o maior problema da transformada de Fourier.

Na transposição do sinal do domínio do tempo (espaço) para o domínio da

freqüência perde-se a localização temporal (espacial), ou seja, olhando os coeficientes da

transformada de Fourier de um sinal é impossível saber onde determinado evento ocorre.

Para minimizar esse efeito Dennis Gabor (1946) adaptou a Transformada de

Fourier para analisar uma pequena porção do sinal. Essa adaptação é chamada de

Transformada de Fourier 'janelada' (do inglês short time Fourier transform) - STFT e é

definida por KUMAR et All, 1994, como:

∫

∞

∞−

= due )t -u g( )u f( ) t w,f( G

u wi-

∫

∞

∞−

= due )u (g )u f( ) t w,f( G

u wi-

(2.7)

A transformada inversa é dada por

∫∫

∞

u wi

dw due )t -u g( )t w,f(G

) t (f (2.8)

Na equação da transformada de Fourier 'janelada' - STFT - a função g é

chamada de função janela e é normalmente uma função real, com g(0) ≠ 0, que decai

rápidamente a zero quando nos afastamos do zero. A função f multiplicada pela função g dá

como resultado uma espécie de modulação da função f pela função g e como resultado só

temos a onda original nas posições em que g é diferente de zero, como podemos ver na

figura

2.4

Após a formação da janela, convolução de f(x) com g(x), aplica-se a

transformada de Fourier e o resultado, as freqüências encontradas, se referem ao intervalo

onde a função g(x) é diferente de zero.

Figura 2.4 - Esquema simplificado de janelamento.

Continuando o processo a função g(x) é deslocada e os passos representados

figura 2.4 são repetidos até que a função janela percorra toda a extensão da função f(t).

Neste processo a mesma janela g(x) percorre todo o sinal independente das

variações que ocorrem na função f(t), ou seja, a janela é a mesma para todas as freqüências

presentes em f(t). Dependendo do tamanho da janela pode-se analisar um espectro de

freqüência maior ou menor e conseqüentemente localizar com menor ou maior precisão uma

determinada freqüência.

Quanto maior o tamanho da janela maior a informação sobre a freqüência e

menor a informação sobre o tempo e vice-versa.

O sucesso deste método depende da escolha da janela g(x). Em alguns casos

a escolha de g(x) permite a visualização clara de dado evento e uma escolha diferente da

janela faz o evento não aparecer.

O procedimento correto seria fazer várias análises com janelas diferentes o

que resultaria em uma quantidade muito grande de gráficos para analisar, além de gastar um

tempo enorme em cálculos computacionais.

O próximo passo é trabalhar com janelas de tamanho variável e ao mesmo

tempo fazer o deslocamento, o que nos leva à transformada WAVELET.

2.4 - Transformada wavelet

Segundo CHIERECE, 2003, os trabalhos de Alfred Haar, Paul Levey,

Alberto Calderon, J. Morlet, Yves Meyer, Stephane Mallat e Ingrid Daubechies formalizaram

a transformada wavelet e a sua aplicabilidade, nos mais diversos ramos da ciência, não pára

de crescer.

Na análise de Fourier é feita uma decomposição do sinal em funções senos e

cossenos com freqüências diversas e na análise usando a transformada wavelet a

decomposição do sinal é feita por meio da função wavelet que será modificada por meio de

dilatações e compressões (uso da escala) e será transladada ao longo de todo sinal. Com isto

teremos indicação das freqüências presentes no sinal, em que tempo elas ocorrem (informação

não fornecida com a transformada de Fourier) e, também, evita a possibilidade de escolha

errada da janela na transformada de Gabor.

Segundo DAUBECHIES, 1992, LIMA, 2002, MORETTIN, 1979, podemos

definir a transformada wavelet como:

dt)t()t(f)b;a(C

b,a

∫

∞

∞−

ψ= (2.9)

)

(

)t(

b,a

−

ψ=ψ (2.10)

onde ψ(t) é a função (wavelet) escolhida para se fazer a análise do sinal f(t). O valor

a é o

parâmetro de escala e

b o parâmetro de localização com a e b ∈ R e a≠0. Mudando o valor de

a tem-se o efeito de dilatação (a>1) ou de contração (a<1), enquanto que mudanças no

parâmetro

b tem o efeito de analisar a função f(t) em torno deste ponto.

A função ψ (t) deve possuir as seguintes características (DAUBECHIES,

1992):

a) satisfazer as condições de normalização

1dt)t(

=ψ

∫

b) decair suficientemente rápido para se obter uma boa localização

∞<ψ

∫

∞

−

dt)t(

c) satisfazer a condição de admissibilidade

∫

∞

−

∞<= dwC

)w(

onde Ψ(w) é a transformada de Fourier de ψ (t). Isso é o que garante a transformada wavelet

inversa

dadb

)

(

)b;a(C

)t(f

−

Ψ=

∫∫

∞

∞−

∞

∞−

d) a média ser igual a zero, ψ (t) comporta-se como uma onda (daí o nome wavelet).

∫

∞

−

=ψ 0dt)t(

2.5 - Escala e Deslocamento

O deslocamento de uma wavelet tem o mesmo sentido que o janelamento na

transformada de Gabor, o deslocamento da função ao longo do sinal.

Matemáticamente ψ(t) e ψ(t-b) estão separados por uma distância b.

Podemos ver isso na

figura 2.5 onde foi usado b=1.5 e a wavelet "Chapéu Mexicano".

Figura 2.5 - Deslocamento da wavelet "chapéu mexicano".

A escala faz uma dilatação ou uma compressão na wavelet e, como o fator

de escala aparece na equação dividindo a variável, fatores de escala menores do que 1 fazem a

compressão da wavelet e fatores maiores que 1 fazem a dilatação (

figura 2.6).

Figura 2.6 - Mudanças de escala para a wavelet "chapéu mexicano".

-5 -4 -3 -2 -1 0 1 2 3 4 5

-0.4

-0.2

0.2

0.4

0.6

0.8

Deslocamento da wavelet para b=1,5

Wavelet "chapéu mexicano"

-5 0 5

-0.5

0.5

a=1

-5 0 5

-0.5

0.5

a=1/2

-5 0 5

-0.5

0.5

a=1/4

-5 0 5

-0.5

0.5

a=1

-5 0 5

-0.5

0.5

a=2

-5 0 5

-0.5

0.5

a=4

Com escalas pequenas, wavelets mais comprimidas, estamos comparando

com detalhes do sinal que variam rapidamente, ou seja, estamos procurando alta freqüência no

sinal. Por outro lado, quando trabalhamos com escalas grandes, wavelets mais dilatadas,

estamos comparando com detalhes do sinal que variam lentamente, ou seja, estamos

procurando baixa freqüência.

Os deslocamentos e as mudanças de escala geram novas funções chamadas

de wavelet filhas. Na verdade estamos trabalhando com várias funções ortogonais

(determinadas pela wavelet "mãe", pela escala e pelo deslocamento) com cada uma cobrindo

todo espaço de definição do sinal.

2.6 - Tipos de wavelet (mais comuns)

Wavelet Haar

g(t) = 1 para 0≤ t ≤ 0,5

g(t) = -1 para 0,5 ≤ t ≤ 1

g(t) = 0 para t > 1 e t < 0

0 0.5 1 1.5

-1.5

-1

-0.5

0.5

1.5

g(t)

wavelet de Haar

Wavelet Chapéu Mexicano

2/t2

e)t1()t(g

−

−=

Wavelet Morlet

2/t

tiw

ee)t(g

−

-5 0 5

-0.4

-0.2

0.2

0.4

0.6

0.8

g(t)

Wavelet chapéu mexicano

-5 0 5

-3

-2

-1

g(t)

Wavelet Morlet

Wavelet Daubechies

A wavelet daubechies é obtida por meio de iterações de uma função de

dilatação e no apêndice VII apresentamos, resumidamente, o desenvolvimento matemático e o

resultado da aplicação de funções codificadas no software MATLAB para a construção da

wavelet Daubechie de ordem 8.

0 0.5 1 1.5 2 2.5 3

-1.5

-1

-0.5

0.5

1.5

Wavelet Daubechies de ordem 4

Psi(x)

0 1 2 3 4 5 6 7

-1

-0.5

0.5

1.5

Wavelet daubechies de ordem 8

Psi(x)

2.7 - Aplicação da transformada wavelet

As equações 2.9 e 2.10 nem sempre podem ser utilizadas diretamente. Em

alguns casos conhecemos a função f(t), mas, em outros, não conhecemos a f(t) e temos que

fazer a integração numérica, que é o caso de sinais medidos em intervalos de tempo, de

qualquer maneira um programa computacional deve ser usado para a aplicação da

transformada wavelet.

O programa MATLAB oferece duas possibilidades para a aplicação da

transformada wavelet: DWT (do inglês Discret Wavelet Transform) - transformada wavelet

discreta e CWT (do inglês Continuos Wavelet Transform) - transformada wavelet contínua. A

denominação contínua se refere ao fato da CWT operar com qualquer escala e também pelo

fato da wavelet ser deslocada vagarosamente sobre todo o domínio da função (MISIT, 1996).

O algoritmo DWT foi desenvolvido por MALLAT, 1988 e consiste no uso

de filtros complementares, um passa-alta e outro passa-baixa. Após a passagem pelos filtros o

sinal dobra de tamanho e uma parte é eliminada alternadamente. A parte do sinal resultante do

filtro passa-alta é a alta freqüência correspondendo a escalas pequenas e é chamada

usualmente de detalhe. A parte do sinal que passou pelo filtro passa-baixa é a freqüência

baixa correspondendo a escalas grandes e é chamada de aproximação. Para muitos tipos de

sinais a parte de baixa-freqüência contém a maior quantidade de informação A parte do

sinal correspondendo à aproximação é novamente submetida ao conjunto de filtros repetindo

o processo até a aproximação fique com um só elemento. O conjunto formado pela última

aproximação e todos os detalhes forma o conjunto dos coeficientes da transformada wavelet

discreta. (MISIT, 1996, MALLAT, 1988).

A transformada contínua wavelet é a soma do sinal multiplicado pela

wavelet-mãe e filhas (obtidas por meio de deslocamentos e mudança da escala) ao longo de

todo o sinal.

A receita para o cálculo da CWT é apresentada no "Wavelet ToolBox"

(Matlab).

1º Passo:

Tome a wavelet e compare com a seção inicial do sinal original.

2º Passo:

Calcule os coeficientes C, que representam o quanto a wavelet está correlacionada

com a seção do sinal. Quanto maior

C maior a similaridade entre a wavelet e a seção do sinal.

O valor de

C depende da forma da wavelet escolhida.

3º Passo:

Desloque a wavelet para a direita e repita os passos 1 e 2 até cobrir o sinal todo.

4º Passo

Mude a escala da wavelet e repita os passos 1, 2 e 3.

5º Passo

Repita os passos de 1 a 4 para todas as escalas.

Como resultado da aplicação dos cinco passos obtemos os coeficientes da

CWT para diferentes escalas para cada seção do sinal.

A matriz dos coeficientes

C tem M linhas por N colunas, onde N é o número

de pontos usados no sinal e M o número de escalas usadas. Normalmente a matriz dos

coeficientes é representada num gráfico xy e cada ponto da matriz dos coeficientes recebe

uma tonalidade de uma cor que quanto maior o coeficiente mais clara é a tonalidade. Este tipo

de gráfico é chamado de 'escalograma' (esta palavra não tem no Aurélio).

Retornando às funções f

e f

(figuras 2.1 e 2.2) e aplicando a transformada

wavelet, com função de Morlet, obtemos os escalogramas representados na

figura 2.7.

Figura 2.7 - Escalogramas das transformadas wavelet das funções f

(esquerda) e f

(direita) em função da escala.

Os escalogramas mostram que temos um máximo dos coeficientes em torno

da escala 55 e em torno da escala 165 e que no caso da função f

acontece ao longo de todo o

sinal e no caso de f

há uma indicação de onde acontece uma e outra escala.

Ainda estamos sem conhecer o valor das freqüências presentes nos sinais.

Quando se trabalha com escalas pequenas e encontramos um grande valor dos coeficientes

significa que encontramos uma seção do sinal com alta freqüência. Desta maneira podemos

Valor absoluto dos coeficientes para a = 1, 2, 3, ...

tempo (b)

escala (a)

500 1000 1500 2000

101

111

121

131

141

151

161

171

181

191

Valor absoluto dos coeficientes para a = 1, 2, 3, ...

tempo (b)

escala (a)

500 1000 1500 2000

101

111

121

131

141

151

161

171

181

191

relacionar escala com freqüência (na verdade o correto seria chamar de pseudo-freqüência

uma vez que os coeficientes da CWT dependem da escala e do tempo ou espaço ) pela relação

a seguir:

∆

(2.11)

Onde

a é a escala, ∆ é o período do sinal analisado, F

é a freqüência central

da wavelet usada e

é a pseudo-freqüência correspondente à escala a (MISITI et all, 1996) .

Utilizando a equação 2.11 podemos colocar o gráfico da

figura 2.7 com o

eixo Y em freqüências.

Figura 2.8 - Escalogramas das transformadas wavelet das funções f

(esquerda) e f

(direita) em função da freqüência.

Nos gráficos da figura 2.8 temos a indicação das freqüências presentes no

sinal (10 e 30) e onde elas ocorrem.

Na análise usando transformadas wavelet ainda temos uma outra opção que

é trabalhar com o espectro de potência da wavelet (do inglês "global wavelet spectra")

semelhante ao que é feito em termos do espectro de energia obtido com a transformada rápida

de Fourier.

Segundo FARGE, 1992, ZANANDREA et All, 2004, BOLZAN, 2004,

STOEV et All, 2004, pode-se definir uma densidade de potência espectral, como

Valor absoluto dos coeficientes

tempo ( b)

freqüência (a)

500 1000 1500 2000

1664

151.3

79.2

53.7

40.6

32.6

27.3

23.4

20.5

18.3

16.5

13.8

12.7

11.8

10.3

9.7

9.2

8.7

Valor absoluto dos coeficientes

tempo (b)

freqüência (a)

500 1000 1500 2000

1664

151.3

79.2

53.7

40.6

32.6

27.3

23.4

20.5

18.3

16.5

13.8

12.7

11.8

10.3

9.7

9.2

8.7

b)C(a,

)b,a(D

= (2.12)

e o espectro de potência wavelet por

D(a,b) =

)b,a(D

∑

(2.12)

A análise de potência espectral será usada para identificar os períodos

presentes nas séries e a sua respectiva energia.

Voltando às funções f

e f

, que usamos para exemplificar as transformadas

de Fourier e wavelet, aplicando a equação 2.12 para os coeficientes da transformada wavelet e

representando em um gráfico monolog o espectro de potência wavelet (D) pelo período (no

eixo logarítimo), encontramos novamente os períodos usados na construção das funções. Nos

gráficos da

figura 2.9 encontramos os períodos presentes em f

e f

e não temos mais a

indicação de onde ocorrem os períodos.

Figura 2.9 Espectro de potência em função do período para a função f

(esquerda) e

espectro de potência em função do período para a função f

(direita).

-4

-3

-2

-1

0.05

0.1

0.15

0.2

0.25

Período

Espectro de potência wavelet

Espectro de potência wavelet para a função f1

-4

-3

-2

-1

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Período

Espectro de potencia wavelet

Espectro de potência wavelet para a função f2

CAPÍTULO 3

Fractais e expoente de Hurst de séries auto-afins

Com o intuito de analisar as séries temporais pretende-se neste capítulo

introduzir o conceito de dimensão fractal, caracterizar o expoente de Hurst, conceituar auto-

afinidade e estimar o expoente de Hurst usando a transformada wavelet.

3.1 - Fractais

Os mapas de alguns países apresentam uma linha de contorno com uma

série de reentrâncias e saliências que poderão ser vistas ou não dependendo da escala que o

mapa foi construído. À medida que aumentamos a escala do mapa verificamos que a

rugosidade da linha costeira se torna cada vez mais visível. Na

figura 3.1 temos uma linha

sinuosa, em escalas diferentes, que pode representar uma parte de da linha costeira de um

país. Se utilizarmos um padrão de medida grande (por exemplo 1km) para a medir a costa o

valor encontrado se aproxima do comprimento da linha tracejada e se utilizarmos um padrão

de medida pequeno (por exemplo 50m) encontraremos uma valor maior.

Podemos dizer que o comprimento da costa desenhada na

figura 3.1

depende do padrão de medida usado.

Figura 3.1 - Possível linha costeira de um país em três escalas diferentes.

Nos três desenhos da figura 3.1 encontramos certos desenhos da costa que

se repetem em tamanho maior ou menor, são padrões de linha costeira similares. Esta

similaridade acontece com muita freqüência na natureza. Além do contorno da costa,

encontramos também nas plantas, como na folha da samambaia, na couve-flor, a erosão feita

pelos rios, as fraturas no solo após um terremoto, no caminho percorrido pelas cargas elétricas

durante os relâmpagos, nas ondas do mar e, assim, poderíamos relacionar uma série muito

grande de exemplos.

O conceito de similaridade vem do estudo da Geometria. Dois objetos são

similares quando tem a mesma forma embora de tamanhos diferentes. Os ângulos devem ser

iguais e os segmentos de reta devem ter um único fator de escala independente da direção.

Um triângulo equilátero é similar a qualquer outro triângulo equilátero.

Podemos diminuir os lados indefinidamente e ainda assim teremos um triângulo equilátero.

Quando se trabalha com objetos físicos, reais, aumentando-se a escala para ver pedaços cada

vez menores do todo, a partir de certo ponto chegamos à dimensão de moléculas e átomos e

não teremos mais a similaridade.

Peitgen at all, 1992, utiliza, nos casos em que há detalhes reconhecíveis em

todas as escalas macroscópicas enquanto se tratar de objetos reais, o termo auto-similaridade e

quando há escalas diferentes para direções diferentes o termo auto-afim.

As figuras em que há similaridade entre o todo e a parte foram chamadas de

fractais por Benoit Mandelbrot em 1967 e definiu um fractal com a frase

"Um fractal é uma forma feita de partes similares de algum modo ao todo."

3.2 - Dimensão

Uma característica fundamental dos fractais, segundo MANDELBROT,

1983, é ter uma dimensão fracionária. A dimensão fractal de um objeto foge um pouco do

conhecimento comum em termos de dimensão. A geometria euclidiana nos dá uma dimensão

zero para o ponto, um para a curva, dois para a superfície e três para um volume. Uma outra

maneira de definir a dimensão de um objeto é por meio da medida da quantidade de objetos,

em escalas menores, necessários para formar o original.

Como ponto de partida tomemos um segmento de reta de comprimento L.

Podemos dividir o segmento em 3 segmentos de comprimento L/3 e cada um deles, se

multiplicado 3

, se assemelha ao original (figura 3.2 a).

Tomando-se um quadrado de lado L e dividindo-o em quadrados menores

de lado L/3, cada um deles, quando multiplicados por 3

, se assemelham ao original (figura

3.2 b).

Do mesmo modo um cubo de lado L pode ser dividido em cubos menores

com lado L/3 e cada um deles, quando multiplicado por 3

, se assemelham ao original (figura

3.2 c).

Em cada caso o número N de figuras menores que constituem a figura maior

é igual ao fator de escala S elevado a potencia D que nos indica a dimensão da figura

geométrica.

S log

N log

Dou S N

== (3.1)

Figura 3.2 - Representação da escala para um segmento de reta (a), para uma figura plana

(b) e para um sólido (c).

Esse processo de determinação da dimensão de uma figura pode ser

estendido para os fractais e podemos entender melhor o significado de dimensão fractal por

meio de exemplos.

3.3 - Poeira de Cantor

O conjunto denominado poeira de Cantor foi apresentada inicialmente por

Georg Cantor em 1883. Inicialmente consta de um segmento de reta, de comprimento L, que é

dividido em três partes iguais e a do meio é retirada. Como primeiro resultado ficamos com

dois segmentos de reta iguais e com comprimento igual a um terço do original. Cada um dos

segmentos é novamente dividido em três partes e a parte do meio é retirada, resultando, em

quatro segmentos iguais de comprimento um sexto de L. Este procedimento é repetido

resultando (

figura 3.3) em uma estrutura semelhante a poeira, daí o nome poeira de Cantor

(GLEICK, 1990).

Figura 3.3 - Representação do conjunto poeira de Cantor

Do mesmo modo como foi calculada a dimensão para a linha, plano e

sólido, a cada interação do conjunto de Cantor, a partir do segmento de reta inicial, temos:

1ª interação S=3 N=2 nível=1

2ª interação S=9 N=4 nível=2

3ª interação S=27 N=8 nível=3

4ª interação S=81 N=16 nível=4

Usando um gráfico log-log encontramos o valor de D a partir do coeficiente

angular da reta D log S = log N.

0.6309

log3

log2

D ==

Figura 3.4 - Cálculo da dimensão da figura poeira de Cantor.

log ( S )

log ( N ) = D log ( S )

log ( N )

3.4 - Triângulo de Sierpinsky

Para a construção do triângulo de Sierpinsky começamos com um triângulo

equilátero. Ligam-se os pontos médios dos lados do triângulo e retira-se o triângulo assim

formado (

figura 3.5). Para cada triângulo que restou aplica-se novamente a regra acima e

assim continuamente. No limite temos o triângulo de Sierpinsky (FEDER, 1989).

Figura 3.5 - Representação do triângulo de Sierpinsky.

Quando retiramos de um triângulo a figura formada pelos segmentos de reta

que ligam os pontos médios dos lados, restam três triângulos internos cujos lados precisam ser

multiplicados por dois para igualar ao tamanho do lado do triângulo inicial. Comparando com

o procedimento para determinar a dimensão de figuras geométricas simples,

equação 3.1,

temos:

1.585

2 log

3 log

S log

N log

D 2 S e 3 N =====

3.5 - Contagem de caixas

A figura da qual se quer calcular a dimensão fractal, nem sempre tem um

procedimento conhecido de geração e nestes casos deve-se usar o método de contagem de

caixas (box-couting). Consiste em cobrir a figura com um quadriculado, com caixas de lado δ,

e contar a quantidade N de caixas que são necessárias para cobrir totalmente a figura. Na

segunda etapa o tamanho das caixas é diminuído e a contagem se repete. Fazendo o tamanho

das caixas tender a zero a dimensão será dada por

→δ

log

N log

lim D

(3.2)

O valor do limite é encontrado por meio do coeficiente angular da reta log

N(δ) =D log 1/ δ para valores de δ cada vez menores.

O processo de contagem de caixas não pode ser aplicado diretamente em

casos de processos em que as escalas são diferentes em direções diferentes (processos auto-

afins). Para os processos auto-afins a utilização do método de contagem de caixas só pode ser

aplicado se as caixas sofrerem a mesma mudança de escala que as variáveis são submetidas

(CAMPANHA, 2004).

Figura 3.6 - Quadriculado com caixas de lado δ cobrindo a figura que se quer calcular a

dimensão fractal

3.6 - Expoente de Hurst

Harold Edwin Hurst, engenheiro inglês trabalhando no Egito, dedicou uma

grande parte de sua vida ao estudo do fluxo das águas do rio Nilo e seus afluentes.

Inicialmente contratado pelo governo inglês, em 1906, para a construção de uma barragem no

Egito, deparou com o problema da determinação da quantidade ideal de água no reservatório,

de modo que não faltasse nem houvesse um excesso de água o que colocaria a barragem em

perigo. Considerando as chuvas que aconteciam nas nascentes dos afluentes e do próprio Nilo,

qual deverá ser a vazão da represa para ficar em um nível aceitável? Para responder a esta

pergunta, Hurst estudou os mais de 800 anos de registros de observações das cheias anuais do

rio Nilo e a partir de seus estudos propôs um modelo para o controle da vazão do reservatório

da represa. Os seus estudos e modelos estão relatados no livro Long-term storage: an

experimental study (Hurst e all, 1965).

Com o intuito de estabelecer uma regra para a descarga de água e manter o

nível do reservatório dentro do esperado, Hurst inventou um método estatístico "análise de

reescalonamento" (Análise de R/S) para analisar os dados de entrada e saída de água do

reservatório em função do tempo. A água que entra no reservatório depende da vazão dos rios

e seus afluentes que, por sua vez, está fortemente relacionada com a quantidade de chuvas na

região e nas nascentes dos rios. A dependência com a quantidade de chuva indica uma

componente com comportamento aleatório que será levada em conta na política de descargas

do reservatório para manter seu nível dentro do esperado.

Considerando um determinado ano t, o reservatório recebe uma quantidade

de água

(t) e deve liberar uma quantidade para manter o nível. A entrada de água média

durante n anos é dada por:

∑

ξ=ξ

)t(

(3.3)

Essa deve ser a quantidade média liberada por ano para manter o nível do

reservatório. A água acumulada no reservatório em função do valor médio liberado e da

quantidade recebida é

∑

ξ−ξ=

])m([)n,t(X (3.4)

A diferença entre o valor máximo acumulado e o valor mínimo acumulado é

chamado de amplitude ("range") e identificado pela letra R (

equação 3.5), que representa a

capacidade necessária para manter uma saída média

ξ em função de uma entrada anual

ξ (t). Retornando ao início, para um reservatório que nunca transborda e nem fica vazio, R

representa a diferença entre a máxima e a mínima quantidade de água presente no reservatório

(

figura 3.7).

(

)

(

)

n,tXminn,tXmax)n(R

nt1

≤≤

−

= (3.5)

Figura 3.7- Esquema de um reservatório de água (FEDER, 1989).

Hurst usou a relação de R dividido pelo desvio padrão (equação 3.6) para

poder trabalhar com uma quantidade adimensional.

∑

ξ−ξ=

])t([

S (3.6)

Hurst aplicou a sua estatística para vazão de rios, quantidade de sedimentos

após cada cheia dos rios, no tamanho dos anéis dos troncos de árvores e encontrou

empiricamente que R/S era proporcional ao tempo de observação elevado a um expoente que

posteriormente foi chamado por Benoit Mandelbrot de expoente de Hurst e representado pela

letra H.

tetanconsCnC

== (3.7)

Na análise dos dados de cheias do rio Nilo, Hurst notou que haviam

períodos em que uma cheia era seguida de outras cheias e anos de seca eram seguidos de

outros anos de seca. Esta ocorrência mostrava haver uma componente não aleatória no

histórico de cheias e secas. Havia um efeito memória em que os anos de seca e cheias eram

dependentes do que havia acontecido nos anos anteriores. MANDELBROT e WALLIS, 1969,

usaram o termo "Efeito José" se referindo à passagem bíblica em que os sonhos do Faraó

foram interpretados por José e alertaram sobre sete anos de fartura seguidos de sete anos de

dificuldades. As cheias do Nilo acabam depositando uma grande quantidade de sedimentos

nas áreas de plantio localizadas às margens do rio, renovando a fertilidade do solo a cada ano.

Uma grande cheia acarretava uma grande fertilização do solo e conseqüentemente uma grande

colheita. Anos com pequenas cheias (anos de seca) eram indicativos de pequenas colheitas.

O "Efeito José" se refere a essa memória encontrada por Hurst, anos de

cheias eram seguidos de anos de cheias e anos de seca eram seguidos de anos de seca. A

explicação encontrada por Hurst foi que no início da estação de chuvas, antes de começar

elevar o nível das águas dos rios, a chuva tem que umedecer toda a região e só depois as

águas começam a correr em direção aos leitos dos rios. Existe um tempo entre o começo das

chuvas e as inundações que vai depender do quanto a região está seca. Se a região está muito

seca, grande parte das águas da chuva será usada para umedecer a terra, sobrando pouca água

para ocasionar as inundações e assim fertilizar o solo.

Hurst chegou à conclusão que havia um efeito memória. A cheia do ano x

n+1

dependia da cheia do ano anterior x

. Aplicando o método R/S para as cheias do rio Nilo,

Hurst encontrou um valor de H aproximadamente igual a 0.7 (HURST, 1951).

FELLER, 1951, demonstrou que para uma série composta de variáveis

aleatórias que sejam identicamente distribuídas e com variância finita o valor de R/S é igual a

2/1

(S/R

= (3.8)

Para este caso, H=1/2, não há correlação entre os elementos, a série tem as

características dos deslocamentos encontrados no movimento descrito por Robert Brown em

1928, que posteriormente foi chamado de movimento browniano.

Na análise dos dados de nível do rio Nilo, Hurst encontrou um valor de

H=0,7 mostrando que há uma correlação positiva entre os dados. A cheia de um determinado

ano influenciava a cheia do ano seguinte.

Os processos em que H≠1/2 foram chamados por Mandelbrot e Van Ness de

movimento browniano fracionário (fbm). Na

figura 3.8 temos a representação de um

movimento browniano fracionário para diferentes valores do expoente de Hurst.

Figura 3.8 - Exemplos de movimento browniano fracionário gerado pelo algoritmo de

Joachim Johansson (Apêndice III).

0 1000 2000 3000 4000 5000

-5

Y(n)

H = 0.2

0 1000 2000 3000 4000 5000

-60

-50

-40

-30

-20

-10

Y(n)

H = 0.5

0 1000 2000 3000 4000 5000

-50

100

150

200

250

300

350

Y(n)

H = 0.8

Segundo PEART, 2005, PEITGEN, 1992 e CAMPANHA, 2004, a

correlação entre os incrementos do movimento browniano fracionário, B

(t), é dada por

1H2

1 2

])t(B[E

)]t(B)t(B[E

)t(

−

−−

=ρ

−

(3.9)

onde E[ ] é o valor médio e V

é uma constante.

equação 3.9 mostra que a correlação é determinada em função do

parâmetro H e não depende do tempo. A

figura 3.9 mostra a variação da correlação com a

variação do expoente de Hurst.

Figura 3.9 - Variação da correlação de um movimento browniano fracionário em função

do expoente de Hurst.

Para H=1/2 a correlação é zero, os valores da série temporal são

independentes uns dos outros e a série é chamada de movimento browniano.

A série em que H>1/2,

)t(

>0, é chamada de persistente, ou seja, os

incrementos futuros são similares aos incrementos do passado. Situação encontrada por Hurst

no caso das cheias do rio Nilo.

0 0.2 0.4 0.6 0.8 1

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Expoente de Hurst

Correlação

Variação da correlação de um fbm com o expoente de Hurst

Quando H<1/2,

)t(

<0, a série é chamada de anti-persistente e os

incrementos positivos no passado indicam incrementos negativos no futuro, há uma inversão

nos incrementos com maior freqüência.

3.7 - Auto-Afinidade

Na figura 3.9, os movimentos brownianos fracionários Y(n) são funções de

n que na maioria das vezes está diretamente relacionado com o tempo e a diferença entre os

comportamentos dos fbm está caracterizado pelo expoente de Hurst.

Segundo Peitgen, 1992, os incrementos do movimento browniano

fracionário são proporcionais aos incrementos na variável n elevado ao parâmetro de Hurst.

)n(Y ∆∝∆

Esta relação mostra que se a variável n sofrer um escalonamento, passando

para an, o valor da variável Y(n) será escalonada para Y(n) a

, o que caracteriza um processo

auto-afim, escalas diferentes para diferentes variáveis.

)n(Ya)n(Yann

→→ (3.10)

Para este caso a dimensão fractal pode ser determinada pelo processo de

contagem de caixas desde que as dimensões das caixas sejam escalonadas diferentemente.

Como exemplo tomemos a representação do fbm Y(n) com N pontos,

conforme mostra a

figura 3.10.

O quadriculado foi construído de forma que cada caixa na direção n tem

lado 1/N e na direção Y(n) tem lado (1/N)

. Desta forma estamos garantindo a diferença de

escalonamento da

equação 3.10.

Cada caixa tem área 1/N

H+1

que corresponde N

1-H

vezes a área da caixa

utilizada no cálculo da dimensão fractal de um processo auto-similar representado na

figura

3.6.

Figura 3.10 - Quadriculado usado para cálculo da dimensão fractal de processo auto-afim.

Usando a equação 3.2 o valor da dimensão fractal, para o processo auto-

afim, será

Nlog

N log

lim

Nlog

N.N log

lim D

H-2

H-1

N ∞→∞→

H2D

−

(3.11)

No caso da

equação 3.11 a dimensão fractal pode assumir valores não

inteiros. Quanto menor for H, mais a dimensão fractal se aproxima de 2, e o perfil se torna tão

rugoso que é como se preenchesse o espaço como um plano. Por outro lado, quando H cresce,

aproximando de 1, o perfil fica suave e a dimensão fractal se aproxima de 1 (

figura 3.9).

Quando temos o movimento browniano H=0,5 a dimensão fractal será 1,5 (PETERS, 1996).

3.8 - Estimativa de H para série temporal auto-afim

Entre os diversos estimadores para o expoente de Hurst, foi escolhido o

método que usa a transformada de wavelet, que, segundo SIMONSEN et All, 1998, apresenta

resultados mais consistentes em todas as faixas de valores de H.

Seja f(t) uma série temporal auto-afim e escalonando a variável t por λ a

função f(t) passará a λ

f(t), conforme a equação 3.10.

t)f( )t(f

λλ=

−

(3.12)

onde o sinal de igual é entendido como igualdade estatística.

Aplicando a transformada wavelet nos dois lados da

equação 3.12, temos

b)t)](a,f( C[ )b,a)](t(f[C

-H

λλ=

ou usando a definição da transformada wavelet da

equação 2.9

dt )

b-t

( t)f(

)b,a)](t(f[C

H-

ψλλ=

∫

∞

Chamando λt de y teremos

∫

∞

)

b-y/

( f(y)

)b,a)](t(f[C

-H

∫

∞

)

b-y

( f(y)

)b,a)](t(f[C

-H

dy )

b-y

( f(y)

)b,a)](t(f[C

1/2--H

∫

∞

b)a,C[f(t)]( )b,a)](t(f[C

1/2--H

λλλ= ou

)b,a)](t(f[C b)a,C[f(t)](

1/2H+

λ=λλ (3.13)

Ao escalonarmos a variável t por λ, a função f(t) é escalonada para λ

-H

f(λt) e

os coeficientes da transformada wavelet são escalonados por λ

H+1/2

SIMONSEN ET All, 1988, no desenvolvimento do método AWC -

"AVERAGE WAVELET COEFFICIENT", propõem que para evitar a dependência do fator

de deslocamento na

equação 3.13, seja calculada a média do módulo dos coeficientes da

transformada wavelet para cada escala, uma vez que quando se trata de processos auto-afim

estamos muito mais interessados no fator de escala do que no de deslocamento.

|)b,a)](t(f[C| C[f](a)

〉

〈

(3.14)

O operador

〉〈 determina o valor da média aritmética com respeito à

variável b. Com a inserção do valor absoluto dos coeficientes temos a representação da

energia da wavelet para cada escala.

Usando a notação da

equação 3.14, a equação 3.13 passa a

)a](f[C a)C[f](

1/2H+

λ=λ

Fazendo-se o gráfico log C[f](a) por log a, se encontrarmos uma linha reta,

significa que a função f é auto-afim e que o coeficiente angular da reta é igual a H+1/2.

Dessa forma podemos sintetizar os passos do método AWC :

1- Calcular os coeficientes da transformada wavelet para a série temporal;

2- Calcular a média do valor absoluto dos coeficientes da transformada wavelet;

3- Fazer o gráfico do log C[f](a) por log a, determinar a reta pelo método dos mínimos

quadrados e encontrar o coeficiente angular da reta, e,

4- Encontrar o valor do expoente de Hurst por coeficiente angular = H+1/2.

CAPÍTULO 4

Wavelet, Expoente de Hurst e Série Temporal de Chuva

Neste capítulo será feita uma discussão de tópicos que permitam a

estimativa do expoente de Hurst das séries temporais de precipitação de chuvas.

Quando se usa um estimador para encontrar o expoente de Hurst, sempre

fica uma dúvida: qual a validade do estimador que estamos usando? A resposta, pretendemos

dar estimando o valor de H, por meio da transformada wavelet, de séries que foram geradas

com o expoente de Hurst conhecido e assim, podemos comparar os resultados encontrados.

O segundo tópico a ser analisado se refere à dependência do expoente de

Hurst, estimado por meio da transformada wavelet, com o tamanho da série. Esta análise

serve de base para definir um tamanho da série temporal de chuva para estimar o expoente de

Hurst para as regiões climáticas do estado de São Paulo.

Ao aplicar a transformada wavelet temos que escolher a função wavelet e

nem sempre temos uma base segura para isto. Novamente usando séries geradas com o valor

de H conhecido estimaremos para uma mesma série o valor de H com várias wavelet e

teremos assim um meio de comparação.

Uma outra dificuldade que encontramos ao trabalhar com séries temporais

de chuvas é o fato das séries conterem interrupções de medições por um ou mais dias. Como o

método de estimar o expoente de Hurst por meio de transformadas wavelet considera a série

toda, a falta de medições, impede a estimativa de H. Neste capítulo proporemos alternativas

para solucionar esse impasse.

4.1 - Validade do expoente de Hurst estimado por wavelet

Existem vários métodos para a estimativa do expoente de Hurst e entre eles

podemos destacar o método R/S, o periodograma, o variacional, o da transformada de Fourier

e o da transformada wavelet.

CLEGG, 2005, SIMONSEN e HANSEN, 1998, BRESLIN E BELWARD,

1999 E MIRANDA, 1997, apresentaram em seus trabalhos comparações entre os métodos

estimadores fazendo a estimativa de H em séries de fbn e fbm em que se conhece a priori o

expoente de Hurst.

O método R/S superestima H para valores do expoente de Hurst próximos

de zero e subestima para valores próximos de um. O método variacional subestima o valor de

H de um modo geral. O método da transformada de Fourier apresenta um erro, na

determinação da inclinação da reta (H é determinado pelo coeficiente angular da reta), maior

que no caso do método da transformada wavelet. De um modo geral todos os estimadores

apresentam pontos altos e baixos o que dificulta a comparação entre eles. Em função disto

optamos por testar o método da transformada wavelet por meio de séries de fbm geradas para

um determinado H.

O movimento browniano (BROWN, 1827, EINSTEIN, 1905, WIENER,

1923) e movimento browniano fracionário (MANDELBROT, 1977) já foram amplamente

estudados e podem ser tomados como referência para a estimativa do expoente de Hurst.

Existem vários algoritmos que geram um fbm com um valor de H

determinado e, entre eles, o que produz um movimento browniano de forma rápida e direta é

o chamado de deslocamento do ponto médio, descrito por VOSS, PEITGEN et all, 1988, que

é descrito resumidamente a seguir.

Para um processo X(t) que será gerado para o tempo t entre 0 e 1,

começamos atribuindo para X(0)=0 e para X(1) uma amostra de um número randômico

gaussiano.

O próximo ponto será a média entre X(0) e X(1) acrescido de uma

compensação D que é um número randômico gaussiano multiplicado por um fator de escala

igual a 1/2.

Na próxima etapa serão dois pontos médios entre X(0) e X(1/2) e X(1/2) e

X(1) acrescidos de quantidades D

e D

(números randômicos gaussianos multiplicados pela

escala anterior reduzida por 2

0,5

). A figura 4.1 mostra as primeiras etapas para o processo

X(t).

Figura 4.1- Gráfico com os primeiros passos do algoritimo mpdwalk.

Os próximos serão os pontos médios entre X(0) e X(1/4), X(1/4) e X(1/2),

X(1/2) e X(3/4), X(3/4) e X(1), acrescidos da compensação D

multiplicada pela escala, da

etapa anterior, reduzida pelo fator 2

0,5

O processo se repete e novos pontos vão sendo colocados nos pontos

médios entre os pontos já existentes e acrescidos do número randômico gaussiano corrigido

pela escala diminuída pelo fator 2

0,5

Quando estamos interessados em gerar movimento browniano fracionário

(H≠1/2) a compensação D será multiplicada, na primeira vez, pelo número de elementos que

terá a série menos um, elevado a H. As escalas das etapas seguintes serão sempre a escala da

etapa anterior reduzida por um fator 2

(MIDDLETON,1995).

Pretende-se usar estas séries, assim geradas, para estimar o expoente de

Hurst por meio da transformada wavelet e comparar os valores de H estimados com os valores

usados na geração das séries.

Foram escolhidos quatro algoritmos geradores de fbm:

fbmfeder: utiliza a equação 9.25 do livro FRACTALS de Jens Feder, 1988 e

foi codificado para o

software MATLAB por Juha Merikoski (Apêndice I);

mpdwalk: utiliza o algoritmo do deslocamento do ponto médio descrito por

VOSS, PEITGEN et all, 1988 e o programa fornecido por MIDDLETON, 1995(Apêndice II);

fbmsim: baseia-se no artigo "Simulation of stacionary gaussian processes

[0,1]

" de WOOD e CHAN, 1994 e a codificação para o software MATLAB foi apresentada

por JOHANSSON, 2000 (Apêndice III), e,

fbmlevinson: utiliza a fatoração de Cholesky/Levinson e o programa escrito

por GONÇALVES, 1997 faz parte do FRACTAL ANALYSIS SOFTWARE (Apêndice IV).

O fato de termos escolhido o movimento browniano fracionário para testes

do método de estimativa de H é por que podemos gerar fbm com valor de H determinado e

também pelo fato das séries temporais de precipitação de chuvas (soma cumulativa dos

desvios da média) serem similares (BRESLIN & BELWARD, 1999) com o movimento

browniano fracionário (fbm) como podemos ver na

figura 4.2.

Com as séries da

figura 4.2, foi calculada a correlação e encontramos um

coeficiente de correlação de 0,19 com uma probabilidade de não correlação igual a zero.

A medida da correlação foi feita pelo método de Spearman (VIEIRA, 2003)

na função

corr, do MATLAB, que determina o coeficiente de correlação r (-1 ≤ r ≤ 1) e nos

dá, ainda, a probabilidade do teste da hipótese da não correlação.

Figura 4.2- Gráficos de chuva cumulativa dos desvios da média e série fbm gerada com o

mesmo H que a série de chuva.

Figura 4.3 - a) Desvios da média e histograma para uma série de chuva com média 0, e

variância 0,010.

b) Desvios da média e histograma para uma série de fbm com média 0 e

variância 0,058.

0 1000 2000 3000 4000 5000 6000 7000 8000

-1000

-500

500

1000

chuva

0 1000 2000 3000 4000 5000 6000 7000 8000

100

150

fbm

0 5000 10000

-0.5

0.5

Chuva

-0.5 0 0.5 1

2000

4000

6000

Chuva

0 5000 10000

-1

-0.5

0.5

fbm

-1 0 1

200

400

600

fbm

A partir das séries da

figura 4.2, foram construídas as séries dos respectivos

deslocamentos, foi feito uma normalização deixando o maior deslocamento igual a um,

construído o histograma e calculada a média e a variância. O resultado pode ser visto na

figura 4.3.

Pela semelhança entre as séries da

figura 4.1 podemos dizer que as séries de

movimento browniano fracionário são adequadas para testar o método de estimativa do

expoente de Hurst, mas, tendo em vista a

figura 4.2, há necessidade de um estudo maior para

usar as séries de fbm como modelo para a precipitação de chuvas.

Com estes programas geradores de fbm foram elaboradas cem séries para

cada valor de H (0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9), estimados os valores de H por meio da

transformada wavelet para cada série, feita uma média dos valores de H e comparados com o

expoente de Hurst segundo a qual a série foi gerada. Cada série foi gerada com 7.305

elementos que corresponde, em tamanho, a uma série de precipitação diária de chuva de 20

anos que iremos analisar no capítulo 5.

Para cada gerador foi feito um gráfico com o valor de H utilizado na geração

do fbm no eixo das abscissas e o valor médio do H estimado por wavelet no eixo das

ordenadas, que pode ser visto na

figura 4.4.

figura 4.4 a)

0 0.2 0.4 0.6 0.8 1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Valor de H usado na geração das séries

Valor médio dos H estimados nas 100 séries

Séries geradas usando o algorítimo fbmfeder

figura 4.4 b)

figura 4.4 c)

0 0.2 0.4 0.6 0.8 1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Valor de H usado na geração das séries

Valor médio dos H estimados nas 100 séries

Séries geradas usando o algorítimo mpdwalk

0 0.2 0.4 0.6 0.8 1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Valor de H usado na geração das séries

Valor médio dos H estimados nas 100 séries

Séries geradas usando o algorítimo fbmsim

figura 4.4 d)

Figura 4.4 - Gráficos mostrando o valor médio do expoente de Hurst para séries de fbm

obtidas por meio de quatro tipos de geradores de movimento browniano

fracionário. As barras de erro nos gráficos se referem ao desvio padrão.

As análises dos gráficos da

figura 4.4 e da tabela 4.1 mostram que os

valores estimados do expoente de Hurst, por meio da transformada wavelet, de séries de fbm

estão bem próximos dos valores esperados e que apresenta uma regularidade na estimativa de

H para todos os valores do expoente de Hurst. O erro padrão da média que aparece na

tabela

4.1

indica que os valores estimados de H nas cem séries estão bem próximos uns dos outros e

que existe uma coerência nos valores.

0 0.2 0.4 0.6 0.8 1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Valor de H usado na geração das séries

Valor médio dos H estimados nas 100 séries

Séries geradas usando o algorítimo fbmlevinson

Tabela 4.1

Simuladores de movimento browniano fracionário

H mpdwalk fbmsim fbmfeder fbmlevinson

médio

Epm H

médio

Epm H

médio

Epm H

médio

Epm

0,1 0,1345 0,0211 0,0320 0,0179 0,1345 0,0211 0,0396 0,0181

0,2 0,2146 0,0194 0,1652 0,0165 0,2146 0,0194 0,1645 0,0154

0,3 0,3072 0,0165 0,2773 0,0157 0,3072 0,0165 0,2821 0,0154

0,4 0,3954 0,0164 0,3847 0,0149 0,3954 0,0164 0,3874 0,0159

0,5 0,4872 0,0170 0,4857 0,0165 0,4872 0,0170 0,4843 0,0151

0,6 0,5784 0,0163 0,5806 0,0174 0,5784 0,0163 0,5877 0,0163

0,7 0,6695 0,0178 0,6846 0,0178 0,6695 0,0178 0,6863 0,0176

0,8 0,7578 0,0189 0,7765 0,0180 0,7578 0,0189 0,7768 0,0178

0,9 0,8420 0,0212 0,8685 0,0221 0,8420 0,0212 0,8672 0,0219

Tabela 4.1 - Média dos valores do expoente de Hurst estimado por wavelet das séries de

movimento browniano fracionário. Com cada algoritmo utilizado para simular

fbm foram geradas 100 séries para cada valor de H. Epm é o erro padrão da

média (VIEIRA, 1999).

4.2 - Expoente de Hurst e o tamanho da série temporal

Para estabelecer a dependência do expoente de Hurst com o tamanho da

série, vamos recorrer novamente, às séries de movimento browniano fracionário. A partir de

uma série temporal de fbm, gerada pelo algoritmo

mpdwalk (VOSS, PEITGEN et all, 1988),

com tamanho aproximado das séries de chuvas diárias que vamos estudar, estimamos o valor

de H e depois diminuímos o tamanho da série e estimamos novamente o valor de H.

Foram geradas cem séries com H=0,5 (movimento browniano) e com 8192

) elementos e para cada série iniciamos diminuindo o tamanho de 500 em 500 e conforme

foi diminuindo a série tomamos um passo menor.

figura 4.5 temos dois gráficos que mostram a variação de H com o

tamanho da série.

Figura 4.5 a)gráfico do valor médio de H para cem séries com tamanho entre 10 e 500.

b)gráfico do valor médio de H para cem séries com tamanho entre 500 e 8192.

Observação: A barra de erro representada se refere à média dos erros cometidos na

estimativa de H por meio da transformada wavelet.

0 100 200 300 400 500

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Número de elementos da série

Valor de H (média de 100 séries)

Valor de H em função do tamanho da série

0 1000 2000 3000 4000 5000 6000 7000 8000 9000

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Número de elementos da série

Valor de H (média de 100 séries)

Valor de H em função do tamanho da série

Verifica-se nos gráficos que o valor de H apresenta uma tendência de

convergir, para o valor com que a série foi gerada, quando aumenta o número de elementos e

também que, para o número de elementos correspondente a vinte anos, o valor estimado de H

é bom.

Outro dado importante é que o erro ao estimar H vai diminuindo, na média,

conforme aumenta o número de elementos.

4.3 - Ausência de medições em séries temporais de precipitação de chuvas

Quando se trabalha com séries numéricas longas, construídas a partir de

dados de precipitação diária de chuvas, é comum encontrar intervalos de um ou mais dias sem

medições. As razões que motivaram a falta de medições são as mais variadas,

impossibilidades de acesso ao posto meteorológico, dificuldades com o pessoal que faz a

coleta diária e até a desativação por falta de equipamento.

A falta de medições numa série (mais ou menos 20 anos) transforma a série

longa em um conjunto de séries menores e ao estimarmos o expoente de Hurst (H) para séries

pequenas o valor de H varia muito em função do tamanho da série, como foi mostrado no

item 4.2, fazendo que o H das séries pequenas seja diferente da série como um todo (longa).

Como o expoente de Hurst está relacionado diretamente com a memória da

série, os valores futuros de precipitação de chuva dependem dos registros do presente,

propomos, inicialmente, que a medida em falta seja substituída pela média de medidas

anteriores e como o expoente de Hurst, também é uma característica da série inteira, uma

segunda proposta, é unir os trechos da série com medições e assim estimar o valor do

expoente de Hurst.

Para as duas propostas fica a pergunta: até quantas medições ausentes

podem ser substituídas pela média ou retiradas sem que o valor do expoente de Hurst se altere

demasiadamente.

Com uma série de precipitação de chuvas sem interrupções de medidas foi

estimado o valor do expoente de Hurst, em seguida retiradas algumas medidas, aplicadas as

propostas e novamente estimado o expoente de Hurst para cada caso. Comparando os

resultados podemos verificar até que ponto a falta de uma ou mais medidas afetam a

estimativa do expoente de Hurst.

Os dados das séries temporais de chuva foram obtidos no Departamento de

Águas e Energia Elétrica do estado de São Paulo (DAEE) e são registros diários de

quantidade de chuva, medidas em milímetros, durante períodos variados dependendo de cada

posto de coleta de dados. Analisando os registros diários de cada posto, no período de

01/01/1978 a 31/12/1997, encontramos algumas séries sem ausência de medidas, algumas

com ausências de poucos dias e outras com falta de registros por meses seguidos.

A série escolhida como padrão foi a do posto Usina Batista, do município de

Pilar do Sul, com 7305 registros diários no período de 01/01/1978 a 31/12/1997 (série E4020)

e para a estimativa do expoente de Hurst foi feito inicialmente uma correção pelo valor médio

da série. Esta correção consiste em retirar de cada valor da série o valor médio das

precipitações durante todo o período, fazendo que a série passe a ter um valor médio igual a

zero.

0 1000 2000 3000 4000 5000 6000 7000 8000

100

120

140

Número de dias a partir de 01/01/1978

Precipitação em mm

E4028

Figura 4.6 a) Gráfico de precipitação de chuvas no município de Pilar do Sul, posto

Usina Batista no período de 01/01/1978 a 31/12/1997 (E4028).

b) Gráfico da soma cumulativa dos dados da série E4028 (corrigidos pela

média e pela reta que passa no primeiro e último ponto).

Com essa nova série é feita a soma cumulativa dos dados e nova correção é

feita com base na reta que passa pelo primeiro e último ponto da série (

figura 4.6 b). Com esta

última correção pretende-se retirar qualquer tendência originada pelo fato dos valores de

precipitação serem todos maior ou igual a zero. Esta última correção é feita retirando de cada

ponto da soma cumulativa dos desvios da média, o valor correspondente da reta que passa

pelo primeiro e último ponto da série.

A substituição de medições ausentes pela média de medições anteriores foi

feita de quatro maneiras diferentes:

a) a medida que será substituída é sorteada dentro da série, é feita a substituição e o expoente

de Hurst é estimado. Uma nova medida é substituída e o processo se repete;

b) a substituição ocorre em posições consecutivas, a partir de uma posição sorteada, em

direção ao final da série;

c) semelhante ao anterior só que em direção ao início da série, e

d) a substituição é feita baseada nas ausências presentes em séries de precipitação de postos

de coletas próximos ao padrão escolhido.

0 1000 2000 3000 4000 5000 6000 7000 8000

-1400

-1200

-1000

-800

-600

-400

-200

200

400

E4028

Número de dias a partir de 01/01/1978

Precipitação acumulada

Para a estimativa do expoente de Hurst foi utilizada a função cwt do

software Matlab versão 7 (R14), a wavelet de Daubechies do tipo 8 e a escala baseada em

potência de dois com o expoente variando de 0 a n, onde n é o maior inteiro menor que o

logaritmo, da base dois, do número de elementos da série (CHIERICE, 2003, DAUBECHIES,

1992, MISITI ET All, 1996).

Uma vez calculado os coeficientes da transformada wavelet foi feito o

gráfico do logaritmo, na base dois, da média do valor absoluto dos coeficientes pelo

logaritmo, na base dois, da escala. O valor estimado do expoente de Hurst é dado pelo

coeficiente angular da reta resultante menos 0,5 (CAMPANHA, 2004 e SIMONSEN et All,

1998).

Segundo TAYLOR, 1982, a incerteza ao se calcular o coeficiente angular da

reta é dada por:

onde

logaritmo, na base 2, da escala

logaritmo, na base 2, da média do valor absoluto dos coeficientes

N número de pontos usados para determinar a reta

valores calculados sobre a reta para cada x

Nos procedimentos a, b e c, acima, a substituição ocorre até que a diferença

entre o valor estimado da série original e o da série substituída seja igual ao erro presente na

estimativa do expoente de Hurst. Como a substituição é feita uma por vez foi escolhida, de

forma arbitrária, fazer a média das dez medições anteriores à medida que será substituída.

Na utilização do padrão de ausências de medições (item d) foi feito um teste

para verificar qual o número de medições anteriores que deve ser usado no cálculo da média.

O teste consiste em variar o número de medições anteriores no cálculo da média e ver qual

permite estimar o expoente de Hurst mais próximo do valor da série padrão.

Com a série padrão (E4028) foi estimado o expoente de Hurst e

encontramos o valor de 0,6339 ±0,0337.

figura 4.7 a)

figura 4.7 b)

0 500 1000 1500

0.63

0.635

0.64

0.645

0.65

0.655

0.66

0.665

Número de medidas trocadas

Valor de H

Medidas trocadas aleatóriamente

0 200 400 600 800 1000 1200

0.62

0.64

0.66

0.68

0.7

0.72

0.74

Número de medidas trocadas

Valor de H

Medidas trocadas seguidamente de 2796 para o fim da série

figura 4.7 c)

Figura 4.7 a) Medidas substituídas em posições aleatórias.

b) Medidas substituídas a partir da posição 2796 no sentido do final da série.

c) Medidas substituídas a partir da posição 5973 no sentido do início da série.

Com o primeiro tipo de substituição a diferença entre o valor de H da série

padrão e o estimado após a substituição ainda é menor que o erro (±0,0338) depois de 1424

substituições como podemos ver no gráfico

a da figura 4.7.

No segundo tipo de substituição o início das trocas foi feito a partir da

posição 2796 (sorteada) e continuou-se, seguidamente, no sentido do final da série, atingindo

o erro após 174 substituições. Da mesma forma a terceira maneira de se fazer as substituições

foi iniciada a partir da posição 5973 e seguidamente em sentido ao início da série. Os

resultados podem ser vistos nos gráficos

b e c da figura 4.7.

Para o último tipo de substituição foram escolhidas nove séries da mesma

região que a série E4028 que apresentam diversos padrões de ausências de medições. Na

figura 4.8 temos os gráficos das séries G5002 e F5020 em que as posições em que não há

registros de precipitação de chuvas foram marcadas com o número -20 e os gráficos com a

variação do expoente de Hurst em função do número de elementos usados para o cálculo da

média. Na

tabela 4.2 há a relação das nove séries, com identificação, o número de ausências

0 500 1000 1500 2000

0.63

0.64

0.65

0.66

0.67

0.68

0.69

Número de medidas trocadas

Valor de H

Medidas trocadas seguidamente de 5973 para o início da série

de medidas, número de medições usadas no cálculo da média, o valor estimado do expoente

de Hurst mais próximo do padrão e a diferença encontrada.

figura 4.8 a)

figura 4.8 b)

0 1000 2000 3000 4000 5000 6000 7000 8000

-50

100

150

200

250

300

350

400

450

Número de dias a partir de 01/01/1978

Precipitação em mm

G5002

0 200 400 600 800 1000 1200 1400

0.65

0.66

0.67

0.68

0.69

0.7

0.71

0.72

0.73

0.74

0.75

g5002

Nº de elementos usados na média

figura 4.8 c)

figura 4.8 d)

Figura 4.8 a) Precipitação de chuvas de Cananéia, posto Morro Redondo (G5002).

b) Variação de H com o número de elementos usados para o cálculo da média

(G5002).

c) Precipitação de chuvas de Eldorado, posto Ouro Leve (F5020).

d) Variação de H com o número de elementos usados para o cálculo da média

(F5020).

0 1000 2000 3000 4000 5000 6000 7000 8000

-20

100

120

140

160

180

200

Número de dias a partir de 01/01/1978

Precipitação em mm

F5020

0 100 200 300 400 500

0.637

0.6375

0.638

0.6385

0.639

0.6395

0.64

0.6405

0.641

0.6415

F5020

Nº de elementos usados na média

TABELA 4.2

Série de referência E4028 - H=0,633947±0,0337344 0,6002126<H<0,6676814

Posto Cidade

Número

ausências

Número de

elementos

da média

Valor de

Diferença

G5002

Cananéia (Morro

Redondo)

1037 852 0,6549 0,0215

E3043 Guarujá (Perequê) 365 54 0,6469 0,0135

E5018 Itapetininga (E. Hermilio) 266 26 0,6384 0,0050

F5020 Eldorado (Ouro Leve) 139 41 0,6373 0,0039

F4040 Iguape (Momuna) 52 59 0,6334 0,0

F4024 Registro (Jurumirim) 36 82 0,6350 0,0016

E5045 Itapeva 8 3 0,6335 0,0001

Tabela 4.2 - Valores de H usando a média das medidas anteriores para as ausências de

medições.

Nos quatro casos de substituição das ausências pela média das medidas

anteriores, encontramos que um grande número de medidas pode ser substituído pela média

de medidas anteriores sem alterar significativamente o expoente de Hurst. No primeiro caso

substituímos 1456 medidas (4 anos) e ainda não atingimos o limite do erro (

figura 4.7 a).

No caso de substituição em seqüência (de 2976 no sentido do fim da série),

o erro é atingido após 174 substituições, bem diferente dos demais casos. Provavelmente isto

se dá pelo fato de que após 10 substituições começamos a usar na média valores já

substituídos anteriormente (

figura 4.7 b).

No terceiro caso chegamos ao ponto de 1572 substituições o que representa

mais de quatro anos de medidas (

figura 4.7c) lembrando que para estes três casos não foi feito

um teste para ver qual o número ideal de elementos que devem ser usados na média.

O principal resultado fica por conta da substituição utilizando um padrão de

ausência de medidas de um posto meteorológico da mesma região climática. Verificamos que

o número ideal de medições a ser usado para calcular a média é diferente para cada série

usada (

figura 4.8) e os valores estimados para o expoente de Hurst são praticamente iguais ao

estimado na série padrão (

Tabela 4.2).

A segunda proposta nos leva a juntar as partes sem interrupções de

medições e tratar como se fosse uma série completa.

Para esta proposta usamos a mesma série padrão, E4028, e fizemos a

retirada de medidas de três maneiras: aleatória, seguidamente a partir de uma posição e

usando padrões de ausências de séries da mesma região climática que a padrão (as mesmas

usadas na substituição pela média).

Nos dois primeiros casos a cada retirada de uma medida, os dois trechos são

unidos, o novo expoente de Hurst é estimado e comparado com o valor da série padrão.

No terceiro caso as posições das ausências das séries escolhidas são

transferidas para a série padrão, os trechos com medidas são unidos e o expoente de Hurst é

estimado.

figura 4.9 a mostra o gráfico da variação do expoente de Hurst com a

retirada aleatória de medições e vemos que a retirada de mais de 1300 medidas o expoente de

Hurst não atingiu o limite do erro encontrado na estimativa de H.

figura 4.9 a)

0 200 400 600 800 1000 1200 1400

0.605

0.61

0.615

0.62

0.625

0.63

0.635

0.64

Número de medidas retiradas aleatóriamente

Valor de H

Série E4028 - Retirada aleatória de medidas

X: 1318

Y: 0.6062

figura 4.9 b)

Figura 4.9

a) Medidas retiradas em posições aleatórias.

b) Medidas retiradas seguidamente a partir da posição 5973 em direção ao início

da série.

No caso da retirada seguida de medições a

figura 4.9 b mostra que mesmo

após a retirada de 2000 medidas o expoente de Hurst não atingiu o limite do erro.

tabela 4.3 podemos ver que a diminuição da série pela retirada das

medidas não altera significativamente o valor de H

TABELA 4.3

Série de referência E4028 - H=0,633947±0,0337344 0,6002126<H<0,6676814

Posto Cidade

Número

ausências

Número de

elementos

da série

Valor de

Diferença

G5002

Cananéia (Morro

Redondo)

1037 6268 0,6422 0,0083

E3043 Guarujá (Perequê) 365 6940 0,6347 0,0008

E5018 Itapetininga (E. Hermilio) 266 7039 0,6417 0,0078

F5020 Eldorado (Ouro Leve) 139 7166 0,6367 0,0028

F4040 Iguape (Momuna) 52 7253 0,6329 -0,0011

F4024 Registro (Jurumirim) 36 7269 0,6354 0,0015

E5045 Itapeva 8 7297 0,6322 -0,0018

Tabela 4.3 - Valores de H usando a união dos trechos sem ausências de medições.

0 500 1000 1500 2000

0.615

0.62

0.625

0.63

0.635

0.64

0.645

0.65

0.655

0.66

Número de medidas retiradas da posição 5973 para o início da série

Valor de H

Série E4028 - Medidas retiradas seguidamente

X: 1947

Y: 0.6196

Nesta segunda proposta, união dos trechos com medições, apresenta

resultados melhores que os resultados da substituição pela média das medidas anteriores e

será utilizada no nosso trabalho.

Os gráficos da

figura 4.9 e a tabela 4.3 nos mostram que em nenhum dos

casos analisados os valores estimados do expoente de Hurst ficou fora dos limites impostos

pelo erro cometido ao estimar o valor de H.

4.4 - Escolha da wavelet mais adequada para estimar o expoente de Hurst

Não existe um método científico para a escolha da wavelet a ser usada na

análise de uma série usando a transformada wavelet e em função disto, mais uma vez

recorreremos às séries de movimento browniano fracionário com o valor de H conhecido a

priori.

Figura 4.10 Valor médio do expoente de Hurst em função da wavelet usada para estimá-lo.

Db2 Db4 Db6 Db8 Db10 Db14 Morl Mexh

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tipo de wavelet

Valor de H (média de 100 séries)

Séries geradas com o algorítimo fbmfeder

Utilizando o algoritmo do deslocamento do ponto médio descrito por VOSS,

PEITGEN et all, 1988, foram geradas 100 séries, com H=0.7 (valor de H próximo aos valores

encontrados para as série de precipitação de chuvas) e 8192 elementos. Para cada série o valor

de H foi estimado usando oito wavelets diferentes e por fim um valor médio de H foi

calculado para cada wavelet.

As wavelets usadas foram as da família Daubechies, (ordem 2, 4, 6, 8,

10,14), Morlet e Chapéu mexicano.

O resultado do teste pode ser visto na

figura 4.10. Praticamente todas as

waveletes estimaram, na média, o mesmo valor do expoente de Hurst e em função disto no

nosso trabalho usamos a wavelet Daubechies de ordem 8.

O fato dos valores médios de H aparecerem abaixo da linha H=0,7 é uma

característica do gerador de fbm utilizado (mpdwalk), uma vez que, na

figura 4.4 b também

aparece o valor médio abaixo de 0,7.

CAPÍTULO 5

Análise de séries temporais de chuvas

Nos trabalhos de CHIERICE, 2003 e FAVARETTO, 2004, foram

analisados dados de chuva do estado de São Paulo em séries de chuvas de tamanhos, épocas e

regiões climáticas diferentes. Neste trabalho, inicialmente, pretendia-se que todas as séries

analisadas fossem do mesmo período e tivessem o mesmo tamanho. Posteriormente

resolvemos agrupar as séries por região climática usando a divisão apresentada por

MONTEIRO, 1973, em seu atlas A DINÂMICA CLIMÁTICA E AS CHUVAS NO

ESTADO DE SÃO PAULO.

A análise das séries foi feita com o intuito de caracterizar cada região

climática do estado de São Paulo com o expoente de Hurst.

5.1 - Obtenção e tratamento dos dados pluviométricos

Em 1998 o Departamento de Águas e Energia do estado de São Paulo

disponibilizou o seu banco de dados pluviométricos, por meio de um CD('compact disc'), com

registros de todos os postos até o ano de 1997. Os dados estão classificados por cidade, postos

meteorológicos, por data e se apresentam com valores de precipitações de chuva, medidos em

milímetros, para cada dia. Apresenta ainda o valor total mensal, a chuva máxima e a

quantidade de dias com chuvas para cada mês.

Uma primeira análise foi feita no sentido de separar os postos

meteorológicos que continham medições por vinte anos ou mais e principalmente que

continham medições no período de 01/01/1978 a 31/12/1997. O período de vinte anos foi

escolhido em função dos dados disponíveis e por que séries com este tamanho apresentam um

expoente de Hurst com poucas variações em função do tamanho. Se aumentássemos o

tamanho das séries o número de séries em cada região diminuiría drasticamente.

Como já mencionado no capítulo 4 algumas séries apresentam faltas de

medições e na transformação dos dados para um formato compatível com o 'software'

MATLAB (arquivos com extensão .mat), as ausências de medições foram transformadas em

uma precipitação de -20 ficando assim mais fácil a sua identificação.

Após a transformação dos dados em uma série de precipitações em

seqüência, dia após dia, foi aplicada, para as séries com faltas de medições, a segunda

proposta discutida no item 4.3, a união dos trechos da série sem faltas de medições. A

utilização da segunda proposta se deve ao fato de que na simulação feita com o padrão de

faltas numa série completa, foi a que estimou o expoente de Hurst mais próximo do original.

O próximo passo foi separar para cada região e sub-região climática as

séries que seriam analisadas. Em cada região o número de séries usadas foi diferente devido à

disponibilidade de séries que atendiam aos nossos requisitos.

5.2 - Estimativa do expoente de Hurst para cada série de precipitação

O gráfico a) da figura 5.1 mostra uma série de precipitação de chuvas e a

primeira coisa a fazer é uma correção da série pela média. Esta correção consiste em calcular

a média das precipitações e retirar de cada valor da série. Esta correção faz que a série fique

com média zero, ou seja, cada pico representa o deslocamento da média que podemos ver no

gráfico

b) da figura 5.1.

Figura 5.1 - a) Gráfico de uma série temporal de precipitação de chuvas.

b) Gráfico de uma série temporal de precipitação de chuvas corrigida pela

média. Para esta série a média das precipitações considerando os 20 anos é

de 3,6260 mm e temos 5491 dias com precipitação zero.

0 1000 2000 3000 4000 5000 6000 7000 8000

100

120

Tempo em dias a partir de 01/01/1978

Precipitação (mm)

Série E6012 - Coronel Macedo - Região II

0 1000 2000 3000 4000 5000 6000 7000 8000

-20

100

120

Tempo em dias a partir de 01/01/1978

Precipitação (mm)

Série E6012 - corrigida pela média

O próximo passo é fazer a soma cumulativa dos elementos da série (

figura

5.2-a) e aplicar uma última correção para retirar qualquer tendência de deslocamento em uma

direção preferencial.

Figura 5.2 - a) Gráfico da soma cumulativa de uma série de precipitação de chuvas.

b) Gráfico anterior corrigido pela reta que passa pelo primeiro e último ponto.

0 1000 2000 3000 4000 5000 6000 7000 8000

-800

-600

-400

-200

200

400

600

800

Tempo em dias a partir de 01/01/1978

Precipitação (mm)

E6012 - soma cumulativa

0 1000 2000 3000 4000 5000 6000 7000 8000

-800

-600

-400

-200

200

400

600

800

Série E6012 - corrigida pela média e pela reta 1° e último ponto

Tempo em dias a partir de 01/01/1978

Precipitação (mm)

Esta correção consiste em traçar uma reta unindo o primeiro e o último

ponto, e retirar de cada elemento da série, na posição t, o valor calculado sobre a reta na

posição t (

figura 5.2-b).

Neste ponto temos uma série da soma cumulativa dos desvios da média de

precipitação de chuvas, sem interrupções, corrigida pela reta que passa pelo primeiro e último

ponto pronta para ter o seu expoente de Hurst estimado por meio da transformada wavelet.

Usando o método ACW descrito no item 3.8, com a wavelet Daubechies de

ordem 8 e a escala em potência de dois, foram estimados o expoente de Hurst para todas as

séries selecionadas das regiões climáticas.

O resultado para as nove regiões climáticas do estado de São Paulo é

apresentado nas

tabelas 5.1 a 5.9

REGIÃO I

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Sub-região I

Cunha (Ferraz) e2134 1 209 0.7243 0.0290

São Luiz do Paraitinga (Catucaba) e2055 2 4 0.6050 0.0476

Natividade da Serra (Laranjal) e2116 3 8 0.6294 0.0489

Natividade da Serra e2008 4 25 0.6173 0.0444

Paraibuna (Alegre) e2130 5 12 0.6204 0.0480

Natividade da Serra (B. Alto) e2024 6 14 0.5800 0.0519

São Luiz do Paraitinga (Briet) e2135 7 5 0.5876 0.0509

Cunha (Vargem do Tanque) d1025 8 689 0.6063 0.0543

Sub-região Ia

Ubatuba (Picinguaba) e1004 9 122 0.6760 0.0216

Ilha Bela e2012 10 532 0.6548 0.0306

São Sebastião (São Francisco) e2045 11 51 0.6536 0.0286

Caraguatatuba (Porto Novo) e2128 12 844 0.6548 0.0360

Ubatuba (Maranduba) e2122 13 145 0.6361 0.0312

Cunha (Rio Manso) e1006 14 11 0.5821 0.0572

Média dos postos 1 a 8 H=0,6213 erro médio 0,0469

Média dos postos 9 a 14 H=0,6429 erro médio 0,0342

Média geral H=0,6306 erro médio 0,0414

Tabela 5.1 - Valores de H para as séries da região climática I

REGIÃO II

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Sub-região IIa

Cananéia (Morro Redondo) g5002 1 1037 0.6598 0.0291

Iguape (Momuna) f4040 2 52 0.6281 0.0411

Itariri (B. Igrejinha) f4011 3 2 0.6090 0.0428

Miracatu (Pedra do Largo) f4012 4 236 0.6442 0.0290

Juquiá (Escalvado) f4037 5 177 0.6437 0.0297

Registro (Jurumirim) f4024 6 36 0.6497 0.0304

Registro f4005 7 1 0.6217 0.0409

Eldorado f5007 8 5 0.6382 0.0322

Eldorado (Ouro Leve) f5020 9 139 0.6512 0.0312

Cajati (Barra do Azeite) f5014 10 193 0.6598 0.0291

Sub-região IIb

Monguaguá f3002 11 65 0.6373 0.0233

São Vicente e3056 12 352 0.6357 0.0308

Guarujá (Pereque) e3043 13 365 0.5994 0.0357

Santos (Caeté) e3041 14 27 0.6079 0.0304

Sub-região IIc

Pilar do Sul (U. Batista) e4028 15 0 0.6339 0.0337

Coronel Macedo e6012 16 0 0.6386 0.0184

Itapeva e5045 17 8 0.6115 0.0286

Itararé (B. Boa Vista) e6032 18 0 0.6327 0.0219

Barra do Chapéu f6006 19 39 0.6179 0.0259

Ribeirão Branco (Pinara) f5039 20 10 0.6435 0.0227

Capão Bonito (F. Santa Inês) e5071 21 7 0.621 0.0317

São Miguel Arcanjo (Taquaral) f4031 22 1 0.6268 0.0362

Piedade (B. Ribeirão Bonito) e4132 23 1 0.6222 0.0345

Guapiara (B. Pinheiro) f5010 24 5 0.6672 0.0200

Itapeva (Guarizinho) e5032 25 2 0.6404 0.0286

Itapetininga (Eng. Hermilio) e5018 26 266 0.6163 0.0338

Itapetininga (Gramadinho) e5034 27 1 0.6177 0.0393

Sub-região IId

Mogi das Cruzes (taiacupeba) e3223 28 233 0.5964 0.0512

Mauá e3148 29 278 0.5966 0.0498

Santo André (Paranapiacaba) e3037 30 143 0.6864 0.0226

São Bernardo do Campo (Rio Acima) e3244 31 52 0.6279 0.0395

São Paulo (Santo Amaro) e30006 32 203 0.6315 0.0413

Média dos postos 1 a 10 H= 0,6385 erro médio 0,0347

Média dos postos 11 a 14 H=0,6201 erro médio 0,0300

Média dos postos 15 a 27 H=0,6300 erro médio 0,0289

Média dos postos 28 a 32 H=0,6278 erro médio 0,0409

Média geral H=0,6311 erro médio 0,0325

Tabela 5.2- Valores de H para as séries da região climática II

REGIÃO III

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Guararema (Paratei) e3054 1 26 0.6216 0.0440

Jacarei e2031 2 202 0.6092 0.0479

São José dos Campos (Pararangaba) e2099 3 32 0.6307 0.0429

Caçapava (Quirino) e2048 4 6 0.6180 0.0470

Taubaté (Macuco) e2127 5 440 0.5993 0.0535

Guaratinguetá (Brumado) d2009 6 118 0.5937 0.0561

Silveiras d1020 7 0 0.6071 0.0529

Lorena d2035 8 2 0.6337 0.0438

Média geral H= 0,6142 erro médio 0,0485

Tabela 5.3 - Valores de H para as séries da região climática III

REGIÃO IV

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Sub-região IVa

Piracaia (Crioulos) e3229 1 181 0.6681 0.0345

São José dos Campos (Guirra) d3070 2 11 0.6189 0.0484

São José dos Campos (S. F. Xavier) d2021 3 1 0.6081 0.0549

Monteiro Lobato (S. Benedito) d2026 4 323 0.6647 0.0427

Campos do Jordão (Capivari) d2001 5 5 0.6024 0.0553

Sub-região IVb

Atibaia e3074 6 139 0.6425 0.0422

Pinhalzinho d3036 7 83 0.6966 0.0378

Socorro d3030 8 16 0.671730 0.0370

Lindóia d3014 9 20 0.6423 0.0458

Aguas de Lindóia d3024 10 9 0.6270 0.0490

Espirito Santo do Pinhal (U. Pinhal) d3003 11 24 0.6241 0.0474

São João da Boa Vista (F. Paraíso) d3066 12 24 0.6447 0.0463

Aguas da Prata (Fartura) c3043 13 5 0.5682 0.0678

São José do Rio Pardo c3035 14 5 0.6140 0.0536

Caconde c3014 15 23 0.6668 0.0482

Média dos postos 1 a 5 H=0,6324 erro médio 0,0472

Média dos postos 6 a 15 H=0,6398 erro médio 0,0475

Média geral H=0,6373 erro médio 0,0474

Tabela 5.4- Valores de H para as séries da região climática IV

REGIÃO V

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Sub-região Va

Mococa (Sitio Esplanada) c4069 1 1 0.6292 0.0518

Santa Rosa do Viterbo (Bom Sucesso) c4103 2 104 0.6394 0.0483

Tambau c4093 3 10 0.6404 0.0500

Santa Rita do Passa Quatro (Usina) c4095 4 1 0.6414 0.0481

Porto Ferreira c4094 5 31 0.6113 0.0541

Casa Branca (Venda Branca) c4097 6 22 0.6191 0.0504

Pirassununga (Fernando Costa) c4085 7 37 0.5983 0.0565

Santa Cruz da Conceição (S.S.Geraldo) d4032 8 34 0.6283 0.0458

Leme (Cresciumal) d4030 9 0 0.6324 0.0465

Mogi Guaçu (Capetinga) d4105 10 0 0.6120 0.0500

Rio Claro (F. São José) d4016 11 0 0.6313 0.0446

Santa Gertrudes d4059 12 31 0.6381 0.0417

Araras (Santa Cruz) d4034 13 7 0.5962 0.0519

Sub-região Vb

Matão c5074 14 99 0.6046 0.0561

Tabatinga (Araruba) c5107 15 69 0.6351 0.0437

Santa Lucia c5011 16 1 0.6204 0.0531

Araraquara (Chibarro) c5017 17 14 0.5346 0.0531

Bocaina d5077 18 0 0.6246 0.0472

Ribeirão Bonito d5003 19 103 0.6264 0.0457

Brotas (Campo Redondo) d4098 20 5 0.6659 0.0350

Itirapina d4014 21 270 0.6551 0.0425

Torrinha (U. Três Saltos) d5006 22 3 0.6669 0.0392

Sub-região Vc

Icem b6001 23 35 0.6493 0.0513

Colombia (F. Campo Grande) b5063 24 16 0.5845 0.0640

Guaira (F. Antas) b5040 25 1 0.6013 0.0599

Ituverava b4034 26 1 0.6462 0.0522

Pedregulho (Alto do Lajeado) b4062 27 56 0.6206 0.0638

Itirapuá b4053 28 3 0.6610 0.0491

Sales de Oliveira (F. Conquista) b4012 29 12 0.6216 0.0518

Morro Agudo b5004 30 6 0.5983 0.0586

Colina b5028 31 33 0.6150 0.0564

Olímpia (F. Cruz Alta) b6053 32 30 0.5977 0.0585

Novais b5003 33 15 0.6256 0.0553

Monte Alto c5070 34 5 0.6534 0.0477

Sertãozinho (U. Sta Elisa) c5114 35 77 0.6228 0.0550

Altinópolis (Aguas Virtuosas) c4106 36 29 0.6166 0.0537

Média dos postos 1 a 13 H= 0,6244 erro medio 0,0492

Média dos postos 14 a 22 H= 0, 6259 erro medio 0,0462

Média dos postos 23 a 36 H= 0, 6224 erro medio 0,0555

Média geral H= 0,6240 erro medio 0,0509

Tabela 5.5 - Valores de H para as séries da região climática V

REGIÃO VI

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Anhembi d5037 1 41 0.6323 0.0371

Itu (Pirapitingui) e4023 2 77 0.6320 0.0351

Mairinque (D. Catarina) e4043 3 10 0.6874 0.0333

Boituva e4046 4 1 0.6273 0.0374

Laranjal Paulista e4050 5 74 0.6212 0.0398

Capivari d4069 6 70 0.6385 0.0377

Rio das Pedras d4068 7 5 0.6599 0.0349

Campinas d4044 8 0 0.6376 0.0390

Média geral H= 0.6420 erro medio 0.0368

Tabela 5.6 - Valores de H para as séries da região climática VI

REGIÃO VII

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Aguas de Santa Bárbara d6006 1 115 0.6095 0.0433

Itatinga (Lobo) d5040 2 6 0.6440 0.0345

Bauru d6036 3 43 0.6284 0.0393

Iacanga (F. Barreirinho) c5029 4 2 0.6308 0.0466

Borborema (F. Ano Bom do Tiete) c6056 5 3 0.6755 0.0359

Borborema (F. Laranjal) c5082 6 102 0.6007 0.0529

Santa Cruz do Rio Pardo (Sodrelia) d6028 7 92 0.6230 0.0411

Avaré (F. Sta Gabriela) d5080 8 56 0.6537 0.0360

Paulistânia d6091 9 1 0.6415 0.0322

Pardinho e5060 10 13 0.6441 0.0362

Lençóis Paulista (Gleba Rio Claro) d5081 11 2 0.6173 0.0385

Novo Horizonte (F. Uirapuru) c6099 12 38 0.6512 0.0394

Avaré e5014 13 8 0.6210 0.0367

Agudos d5041 14 196 0.6360 0.0394

Reginópolis c6050 15 5 0.6368 0.0424

Média geral H= 0.6342 erro medio 0.0396

Tabela 5.7 - Valores de H para as séries da região climática VII

REGIÃO VIII

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Andradina b8004 1 161 0.6023 0.0457

Alto Alegre c7011 2 3 0.5759 0.0523

Andradina (Três Irmãos) b8029 3 12 0.5956 0.0500

Auriflama b7045 4 38 0.5954 0.0544

Brauna (F. N.S. Aparecida) c7016 5 2 0.6143 0.0415

Floreal b7037 6 19 0.6090 0.0558

Flórida Paulista (F. Sto André) c8030 7 0 0.6149 0.0452

Gastão Vidigal b7043 8 39 0.6179 0.0478

Itapura b8023 9 7 0.5971 0.0504

Jales b7008 10 31 0.5570 0.0655

José Bonifácio c6088 11 90 0.6252 0.0456

Marília (Dirceu) d6098 12 0 0.6218 0.0428

Marília (Amadeu Amaral) d7067 13 1 0.6365 0.0343

Mirandópolis (Yuba) c8056 15 3 0.6278 0.0439

Muritinga do Sul (F. Boa Vista) c8051 16 17 0.6254 0.0468

Osvaldo Cruz c7067 17 21 0.6404 0.0382

Panorama c8019 18 10 0.6136 0.0420

Pereira Barreto (Ideal) b8030 19 5 0.5841 0.0508

Planalto (São Jerônimo) c6101 20 13 0.5772 0.0568

Pontalinda b7038 21 19 0.5609 0.0627

Rubiacéia (Caramuru) c7081 22 12 0.5654 0.0552

Santa Fé do Sul b7024 23 0 0.5973 0.0563

Sud Mennucci b7046 24 32 0.5538 0.0636

Tupã (Varpa) d7061 25 11 0.5810 0.0479

Tupã c7043 26 2 0.5715 0.0523

Turiuba b7050 27 81 0.6029 0.0531

Ubarana (F. cataco) c6089 28 7 0.6484 0.0407

Valparaíso (Perpetuo Socorro) c7049 29 4 0.6153 0.0444

Média geral H= 0,6010 erro medio 0,0495

Tabela 5.8 - Valores de H para as séries da região climática VIII

REGIÃO IX

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

Presidente Epitácio (Sul Mineira) c9007 1 36 0.6087 0.0339

Maracai (B. Agua do Mato) d7041 2 2 0.5935 0.0398

Mirante do Paranapanema d8008 3 2 0.6149 0.0306

Marabá Paulista d8035 4 8 0.6360 0.0265

Santo Anastácio (F. Santa Isabel) d8038 5 0 0.6287 0.0308

Narandiba d8040 6 28 0.6253 0.0302

Presidente Venceslau (F. Clotilde) d8048 7 309 0.6173 0.0265

Iepe (Jaguarete) d8053 8 10 0.5835 0.0431

Rosana (Itaporã) d9001 9 47 0.6284 0.0222

Teodoro Sampaio (Santaida) d9002 10 21 0.5776 0.0355

Teodoro Sampaio d9003 11 10 0.5872 0.0363

Rosana (Nova do Pontal) d9016 12 31 0.5765 0.0366

Presidente Epitácio (Sucurita) d9019 13 17 0.6283 0.0269

Florinea d7054 14 103 0.5969 0.0341

Candido Mota (F. Sta Rosa) d7056 15 33 0.5882 0.0385

Média geral H= 0,6060 erro medio 0,0328

Tabela 5.9 - Valores de H para as séries da região climática IX

Figura 5.3 - Valor médio de H para cada região e sub-região climática do estado de São

Paulo.

5.3 - Séries com os períodos seco e úmido bem definidos

Com o intuito de comparar o valor estimado do expoente de Hurst,

selecionamos séries de precipitação de chuvas em que o período seco e úmido são bem

definidos. Tomamos o cuidado de escolher séries do mesmo tamanho e período que as séries

analisadas no item anterior.

As séries temporais utilizadas são provenientes de cidades do Ceará e todos

os procedimentos utilizados na análise das séries do estado de São Paulo foram usados. Na

figura 5.4 apresentamos dois gráficos de precipitação, das cidades de Assaré e Cratos, onde é

marcante a presença do período seco.

tabela 5.10 temos os valores de H estimados, o erro cometido na

estimativa, o valor médio e o erro médio das séries.

0 1095 2190 3285 4380 5475 6570 7665

100

Tempo em dias a partir de 01/01/1978

Assaré

Precipitação (mm)

Figura 5.4 - a) gráfico de precipitação de chuvas da cidade de Assaré.

b) gráfico de precipitação de chuvas da cidade de Cratos.

POSTOS ANALISADOS DA REGIÃO NE

MUNICÍPIO POSTO Nº FALTAS H ERRO(H)

ASSARÉ 1 0 0.7369 0.0412

BOA VIAGEM 2 0 0.7414 0.0429

ORÓS 3 31 0,7591 0,0417

QUIXERAMOBIM 4 0 0.7689 0.0382

CRATO 5 0 0.7148 0.0476

QUIXADÁ 6 0 0.7716 0.0456

SOBRAL 7 0 0.7943 0.0482

TIANGUÁ 8 0 0.7884 0.0503

Média geral H=0,7494 erro medio 0,0445

Tabela 5.10 - Valor de H para séries da região nordeste.

0 1095 2190 3285 4380 5475 6570 7665

100

120

140

Tempo em dias a partir de 01/01/1978

Crato

Precipitação (mm)

5.4 - Espectro de potência wavelet para as séries temporais de chuva

No capítulo 2 foi mostrado a utilização da transformada wavelet para a

identificação dos períodos presentes em uma série por meio de escalogramas e pelo espectro

de potência wavelet.

TRIGO et all, 1996, estudando a variação dos caudais diários de rios

brasileiros e portugueses, BOLZAN, 2004, analisando a variação do número de manchas

solares e ZANANDREA et all, 2005, trabalhando com a detecção de sinais de radares,

utilizaram o espectro de potência wavelet para identificar períodos presentes nas séries

temporais.

Nas séries de precipitação de chuvas completas (7305 elementos) foi

calculado o espectro de potência da wavelet, com o intuito de encontrar os períodos presentes

em cada série.

O procedimento para a determinação do espectro de potência da wavelet

pode ser descrito em cinco etapas:

1- a wavelet usada foi a Daubechies de ordem 8 e a escala em potência de dois;

2- com as séries de precipitação de chuvas (soma cumulativa dos desvios da média

devidamente já corrigidas), foi calculada a transformada wavelet, usando a função cwt

("continuou wavelets transforms") do programa MATLAB, obtendo os coeficientes da

transformada wavelet em função da escala;

3- para cada escala foi feita a soma do quadrado do valor absoluto dos coeficientes e depois

dividido pelo valor da escala determinando, assim, o espectro de potência wavelet (

D);

4- usando a função scal2freq do programa MATLAB, encontramos as freqüências (pseudo-

freqüências) correspondentes a cada escala e por meio da relação período = freqüência

-1

obtemos os períodos (T), e,

5- confeccionando o gráfico monolog D x T (T no eixo logarítimo) podemos encontrar os

períodos presentes nas séries temporais de chuva.

Para cada série esse procedimento foi aplicado e os gráficos presentes nas

figuras 5.5 a 5.13.

x 10

Região II - E6012

Período em dias

Espectro de potência wavelet

0.5

1.5

2.5

3.5

x 10

Região II - E4028

Período em dias

Espectro de potência wavelet

Figura 5.5 - Espectro de potência de séries de precipitação de chuvas da região climática

II (a, b e c).

Figura 5.6 - Espectro de potência de série de precipitação de chuvas da região climática

III.

0.5

1.5

2.5

x 10

Região II - E6032

Período em dias

Espectro de potência wavelet

x 10

Região III - D1020

Período em dias

Espectro de potência wavelet

x 10

Região V - D4030

Período em dias

Espectro de potência wavelet

x 10

Região V - D4105

Período em dias

Espectro de potência wavelet

x 10

Região V - D5077

Período em dias

Espectro de potência wavelet

Figura 5.7 - Espectro de potência de séries de precipitação de chuvas da região climática V

(a, b, c, d).

x 10

Região V - D4016

Período em dias

Espectro de potência wavelet

Figura 5.8 - Espectro de potência de série de precipitação de chuvas da região climática

VII.

x 10

Região VI - D4044

Período em dias

Espectro de potência wavelet

0.5

1.5

2.5

3.5

4.5

x 10

Região VIII - C8030

Período em dias

Espectro de potência wavelet

Figura 5.9 - Espectro de potência de séries de precipitação de chuvas da região climática

VIII (a, b, c).

x 10

Região VIII - D6098

Período em dias

Espectro de potência wavelet

x 10

Região VIII - B7024

Período em dias

Espectro de potência wavelet

0.5

1.5

2.5

x 10

Região IX - D8038

Período em dias

Espectro de potência wavelet

Figura 5.10 - Espectro de potência de série de precipitação de chuvas da região climática

IX.

CAPÍTULO 6

Conclusões

1- No capítulo quatro foi discutido sobre a estimativa do expoente de Hurst por

meio da transformada wavelet e o processo foi testado utilizando-se séries de fbm geradas

com um valor de H conhecido a priori.

O resultado mostrou que o método de estimativa do expoente de Hurst por

meio da transformada wavelet é consistente ao apresentar, para os valores de H testados,

valores muito próximos dos valores de H estabelecidos a priori.

Um outro ponto a ser ressaltado é o fato do desvio padrão (calculado quando

foi feita a média do valor de H para cada conjunto de cem séries) ser pequeno demonstrando

coerência nos resultados.

Com estes resultados podemos dizer que o método da transformada wavelet

para estimar o expoente de Hurst para um movimento browniano fracionário é confiável e

abrangente.

2- Na escolha do tamanho das séries temporais de precipitação de chuvas que

foram analisadas, foi levado em conta, principalmente, a disponibilidade dessas séries para

cada região climática do estado de São Paulo. O valor escolhido foi de 7305 elementos que

corresponde a uma série diária de precipitação de chuvas no período de primeiro de janeiro de

1978 a trinta e um de dezembro de 1997.

Para testar a dependência do tamanho da série na estimativa do expoente de

Hurst recorremos novamente às séries de movimento browniano fracionário. A partir de séries

geradas com H=0,5 e com tamanho 2

, o tamanho foi sendo reduzido e o valor de H

recalculado.

Os gráficos da figura 4.2 nos mostram que para séries, com tamanhos

correspondentes a vinte anos, o valor de H está em uma região de estabilidade (sem grandes

variações) e que o valor do erro cometido ao estimar o valor de H é pequeno. Seguindo a

tendência dos gráficos se o número de elementos aumentasse o valor de H seria melhor

definido, porém, a dificuldade em encontrar séries de precipitação de chuvas seria muito

grande.

3- É comum encontrar em séries de precipitação de chuvas ausências de

medições que interrompem a série transformando em duas ou mais séries menores. Duas

sugestões para aproveitar o maior número de elementos da série foram analisadas. A primeira

trata de substituir as ausências pela média de precipitações anteriores, justificada pelo fato que

o expoente de Hurst tem características de memória e a segunda simplesmente unir as partes

da série sem interrupções baseado no fato que o valor de H é característica da série como um

todo.

As duas sugestões testadas apresentaram resultados muito próximos e com

um erro muito pequeno e no nosso trabalho optamos pela segunda sugestão por ser mais

simples, apresentar um erro menor e principalmente por ser a série de chuvas auto-afim cada

um dos trechos unidos guarda as mesmas características estatísticas da série como um todo.

4- Com os métodos testados e aprovados o próximo passo foi estimar o

expoente de Hurst para as regiões climáticas do estado de São Paulo.

Nas tabelas 5.1 a 5.9 verificamos que os valores de H dentro de cada região

e também entre as regiões não tem uma variação muito grande ficando em torno de 0,6. Na

figura 5.3 temos o valor médio de H para cada região e sub-região.

No estado de São Paulo as chuvas, em todas as regiões climáticas, são

reguladas pela atuação das correntes da circulação atmosférica da vertente atlântica da

América do Sul sendo a mais influente a Frente Polar Atlântica que atua desde o rio Prata até

as proximidades da linha do Equador. Dependendo da atividade das massas polares (tamanho

e freqüência com que são geradas) é que aparecem as diferenças nas quantidades anuais de

chuvas. Um ano de elevada atividade polar implica em elevada pluviosidade para todo o

território paulista.

Os índices pluviométricos anuais do estado de São Paulo oscilam entre 1100

a 2000 mm, existindo pequenas áreas serranas do litoral que chegam em torno de 4500 mm

(uma das maiores do país) e mesmo as regiões mais secas (pequenas áreas) os índices não são

inferiores a 950 mm.

Em função das regiões climáticas receberem praticamente a mesma

influência das frentes polares e pelo alto índice pluviométrico era de se esperar que o valor

estimado do expoente de Hurst para as regiões estivessem muito próximos uns dos outros,

uma vez que, as séries de precipitação de chuvas se referem ao mesmo período de observação.

As diferenças de valores de H, encontrados entre as regiões e sub-regiões, se

devem ao fator relevo e ao fato de receber as frentes com força total ou já aliviadas de sua

quantidade de água.

Quando foi feito, para efeito de comparação, a estimativa do expoente de

Hurst para séries temporais de chuvas de algumas cidades do estado do Ceará, encontramos

um valor médio 0,75 para cidades que tem um valor médio anual de precipitação em torno de

857 mm e os períodos seco e chuvoso bem definidos. Isto vem mostrar que os valores de H

encontrados no estado de São Paulo foram bem estimados, uma vez que, aqui os períodos

seco e chuvoso não são bem definidos.

5- A análise dos períodos presentes nas séries de precipitação de chuvas, por

meio do espectro de potência wavelet, só foi feita para as séries que apresentavam todas as

medições do período (7305 dias). O resultado do espectro de potência wavelet é um gráfico

com o tempo no eixo das abscissas e quando se usa uma série que teve os trechos com

medições unidos é equivalente dizer que trechos do eixo do tempo foram eliminados e a

análise espectral perde o sentido. Para usar séries que tiveram trechos unidos ou substituídos

pela média das medições anteriores há necessidade de se fazer um estudo mais aprofundado

para ver a influência da união dos trechos ou a substituição pela média no deslocamento dos

períodos presentes.

Dos períodos encontrados nas séries de precipitação de chuvas, o que mais

se destaca é o período em torno de 365 dias que está presente em todas as séries.

Devido à pequena quantidade de séries completas dentro de cada região a

análise fica prejudicada ficando as diferenças encontradas caracterizadas pelas especificidades

do relevo de cada região, uma vez que todas as regiões do estado de São Paulo estão sujeitas

às mesmas influências das frentes polares.

Com exceção da região IX todas as outras apresentam um segundo período

variando de 1200 a 4000 dias. A região II apresenta um segundo período com valor médio de

1400 dias. O segundo pico da única série disponível da terceira região aparece em torno de

4000 dias, o maior período encontrado.

A região VI apresenta uma única série com período de 1200 dias e a região

V aparece com um valor médio de 2200 dias.

As regiões VIII e IX recebem influência direta da Frente Polar do Atlântico

sendo que esta passa primeiro na região IX, descarregando a maior parte das chuvas, e depois

vai para a região VIII. Provavelmente por causa disso a região VIII não apresenta o segundo

período em torno de 1000 dias, como na região IX, coincidindo o terceiro período da IX com

o segundo da oitava região. Na região oitava somente as frentes polares de alta intensidade

que tem influência na precipitação de chuvas.

A análise feita da periodicidade usando o espectro de potência wavelet foi

superficial por termos usado séries com tamanho limitado a 20 anos. Sugerimos como

trabalho futuro um estudo mais aprofundado do espectro de potência wavelet com relação à

limitação do tamanho da série na identificação da periodicidade e a identificação dos períodos

em séries de precipitação de chuvas das várias regiões climáticas do Brasil fazendo uma

comparação entre séries do litoral e de regiões centrais.

Referências Bibliográficas

ABRY, P.; VEITCH, D. Wavelet Analysis of Long Range Dependent Traffic,

IEEE Trans.

on Info. Theory, v. 44, n. 1, p. 2-15, 1998.

ABRY, P.;VEITCH, D.; FLANDRIN, P.; TAQQU, M. S. Wavelets for the analysis,

estimation and synthesis of scaling data.

Self-Similar Network Traffic and

Performance Evaluation. New York: Kihong Park and Walter Willinger. Wiley, p. 39-

88, 2000.

AYOADE, J. O.

Introdução à climatologia para os trópicos. 8-ed. - Rio de Janeiro:

Bertrand Brasil, 2002. 332p.

BARTON, C. C.; LA PONTE, P. R.

Fractals in the earth sciences. New York: Plenum

Press, 1995.

BOLZAN, M. J. A. Análise da transformada em ondaletas aplicada em sinal geofísico.

Revista brasileira de ensino de física, v. 26, n. 1, p. 37-41, 2004

BRESLIN, M. C.; BELWARD, J. A. Fractal dimensions for rainfall time series.

Mathematics

and Computers in Simulation, Amsterdam, v.48, p. 437 - 446, 1999.

BURGOS, T. R.; VALDÉS, E. P.

Estimation of the fractal dimension of a rainfall times

series over a zone relevant to the agriculture in Havana. Departamento

Agrometeorologia INSMET, v. 5, n. 1, january 1999.

BURROUGHS, P. A. Fractal dimensions of landscapes and other environmental data.

Nature, London, v. 294, n.5837, p. 240 - 242, Nov. 1981.

BURRUS, C. S., GOPINATH, R. A. and GUO, H.

Introduction to wavelets and wavelet

transforms, A Primer.

N. J., Prentice Hall, 1997.

CALDERON, A. Intermediate spaces and interpolation, the complex method. Studia-math,

24, p. 113-190, 1964.

CAMPANHA, J. R.

Sistemas complexos e aplicações. 2004. - Tese de Livre Docência -

Instituto de Geociências e Ciências Exatas, Universidade Estadual Paulista, Rio Claro,

2004.

CARR, J. R. Statistical self-affinity, fractal dimension and geologic interpretation.

Engineering Geology, Amsterdam, v. 48, p. 269 - 282, 1997.

______; BENZER, W. B. On the practice of estimating fractal dimension.

Mathematical

Geology, New York, v. 23, n. 7, p. 945 - 959, 1991.

CHIERICE, R. A. F.

O uso de wavelets na determinação do expoente de Hurst de uma

série temporal diária de chuvas do município de Araras-SP de 1955-2000. 2003. –

f.71. Dissertação (Mestrado) – Instituto de Geociências e Ciências Exatas, Universidade

Estadual Paulista, Rio Claro, 2003.

CHRISTOFOLETTI, A. L. H.

Análise fractal e multifractal da estrutura de estações

chuvosas em localidades do Estado de São Paulo. 1997. -. Tese (Doutorado em

Geociências e Meio Ambiente) - Instituto de Geociências e Ciências Exatas,

Universidade Estadual Paulista, Rio Claro, 1997.

CLEGG, R., G. A Practical Guide to Measuring the Hurst Parameter,

School of Computing

Science Technical Report Series

CS-TR-916, University of Newcastle, 2005

COOLEY, J. W.; TUKEY, O. W. An algorihm for the machine calculation of complex fourier

series.

Mathematics and Computers, v 19, p. 297 - 301,1965.

DAUBECHIES, I.

Ten lectures on wavelets. Philadelphia: Rutgers University and AT&T

Bell Laboratories, 1992.

DUARTE, P. F. Q.

Análise das periodicidades das séries temporais de nível do mar e

pressão atmosférica de Cananéia, utilizando a transformada de ondaletas (wavelets).

1998. - f. 66. Dissertação (Mestrado em Ciências) – Instituto Oceanográfico,

Universidade de São Paulo, São Paulo, 1998.

EINSTEIN, A. On the movement of small particles suspended in a stationary liquid demanded

by molecular-kinetic theory of heat.

Annalen der Physik, v. 17, p. 549-560, 1905.

FARGE, M. Wavelets transforms and their applications to turbulence.

Annual Review of

Fluid Mechanics, Palo Alto, V. 24, p. 395-457, 1992.

FAVARETTO, A. B.

Estimativa do expoente de Hurst de séries temporais de chuva do

estado de São Paulo usando as transformadas de Fourier, Wavelets e Analise de

R/S. 2004. - f. 89. Dissertação de (Mestrado) - Instituto de Geociências e Ciências

Exatas, Universidade Estadual Paulista, Rio Claro, 2004.

FEDER, J.

Fractals. New York: Plenum Press, 1989.

FELLER, W. The asymptotic distribution of the range of sums of independent random

variables.

Annals of Mathematical Statistics, 22,427, 1951.

GABOR, D. Theory of communication.

Journal of Institute of Electrical Engineers,

v. 93, p. 429 - 457, 1946.

GALLANT, J. C.; MOORE, I. D.; HUTCHINSON, M. F.; GESSLER, P. Estimating fractal

dimension of profiles: a comparison of methods.

Mathematical Geology,

New York, v. 26, n. 4, p. 455 - 481, 1994.

GAN, K. C.; McMAHON, T. A.; FINLAYSON, B. L. Analysis of periodicity in streamflow

and rainfall data by colwell's indices.

Journal of Hydrology, Amsterdam, v. 123, n. 1/2,

p. 105 - 118, 1991.

GLEICK, J.

Caos: a criação de uma nova ciência. Rio de Janeiro: Campus, 1990. 310 p.

GOLDBERGER, A. L. and WEST, B. J. Applications of nonlinear dynamics to clinical

cardiology. Annals of the New York Academy of Sciences, vol 504, Issue 1, 195-213,

1987.

GOMES, S. M.; DOMINGUES, M. O.; KAIBARA, M. K. Navegando de Fourier a Wavelets:

Aplicações em Dinâmica e Controle/J. M. Balthazar...[et all]-São Carlos,

SP: ABCM: SBMAC, 2001, 501p.

GRABBE, J. O.

Chaos & Fractals in Financial Markets. Disponível em: <http://

orlingrabbe.org/chãos index.htm>

GROSSMANN, A.; MORLET, J. Decomposition of hardy functions into square-integrable

wavelets of constant shape.

SIAM, Journal. Mathematics Analysis, v. 15, p. 723 - 736,

1984.

HAAR, A. Zur theorie der ortthogonaler fumktionensysteme.

Mathematiche Annalen,

Berlin, v. 69, p. 331 - 371, 1910.

HARPEN, M. D. An introduction to wavelet theory and application for the radiological

physicist.

Med. Phys. vol 25, p 1985-1993, 1998.

HUDGINS, L.; FRIEHE, C. A.; MAYER, M.E. Wavelet transforms and atmospheric

turbulence.

Physical Review Letters, New York, v. 71, n. 20, p. 3279 - 3282, 1993.

HURST, H. E. Long-term storage of reservoirs.

Transaction of the american society of civil

engineers,

1951.

HURST, H. E.; BLACK, R. P.; SIMAIKA, Y. M.

Long-term storage: an experimental

study.

London: Constable London, 1965. 145p.

IFEACHOR, E. C.; JERVIS, B. W.

Digital signal processing - a practical approach. Addison

- Wesley, 1994.

JONES, C. L.; LONERGAN, G. T. ; MAINWARING., D. E. Wavelet packet computation of

the Hurst exponent.

Journal of Physics A: Mathematical and General, Bristol, v. 29, p.

2509 - 2527, 1996.

KORVIN, G.

Fractal models in the earth sciences. Amsterdam: Elsevier, 1992. 396 p.

KUMAR, P. F.; GEORGIOU, E.

Wavelet analysis in geophysics: an introduction. Efi

Fourfoula - Georgiou and Praveen Kumar (eds). 1994.

LAM, N. S. N.; DE COLA, L.

Fractals in Geography. New Jersey: PTR Prentice Hall, 1993.

LIMA, P. C. Wavelet: uma introdução.

Matemática Universitária, n. 33, p. 13-34,

dezembro. 2002.

LOVEJOY, B. S.; MANDELBROT, B. B. Fractal properties of rain, and a fractal model.

Tellus. Series A. Dynamic Meteorology and Oceanography, Copenhagen, v. 37,

p. 209 - 232, 1985.

MACAU, E. E. N., Introdução à análise wavelet-fractal de séries temporais.

4º Congresso

temático de dinâmica, controle e aplicações - DINCON 2005, Bauru, 2005.

MAEDER, R. E. Fractional Brownian Motion,

The Mathematica Journal, 6(1), 1995.

MALAMUD, B. D.; TURCOTTE, D. Self-affine time series: measures of weak and strong

persistence.

Journal of Statistical Planning and Inference, Amsterdam, v. 80,

p. 173 - 196, 1999.

MALLAT, S. A theory for multiresolution signal decomposition : the wavelet representation.

IEEE Transaction on Pattern Analysis and Machine Intelligence, v. 11, p. 674-693,

1989.

MALONE, D.

Solutions to Dilation Equations. 2000, Ph. D. Thesis. University of Dublin,

2000.

MANDELBROT, B. B.

The fractal geometry of nature. New York: W. H. Freeman, 1983.

MANDELBROT, B. B.; WALLIS, J. R. Water Resources Research, 5, 321-340, 1969.

MANDELBROT, B. B.; VAN NESS, J. W. Fractional Brownian motion, fractional noise and

applications.

SIAM Reviw, 10, 422-427, 1968.

MATSUSHITA, M.; OUCHI, S. On the self-affinity of various curves.

Physica D,

Non Linear Phenomena, Amsterdam, v. 38, p. 246 - 251, 1989.

MAZZARELLA, A. Multifractal dynamic rainfall processes in Italy.

Theoretical and

Applied Climatology

, Wien, v. 63, p. 73 - 78, 1999.

MIRANDA, J. G. V.

Análise fractal de reescalonamento temporal para chuvas.

1997. - f. 75. Dissertação (Mestrado em Física) – Instituto de Física, Universidade

Federal da Bahia, Salvador, 1997.

MISITI, M.; MISITI, Y.;OPPENHEIM, G.; POGGI, JM.

Wavelet Toolbox: For use with

Matlab. The Mathworks. 1996.

MONTEIRO, C. A. F.

A dinâmica climática e as chuvas no estado de São Paulo: estudo

geográfico sob a forma de atlas. São Paulo: IG/USP, 1973.

MORETTIN, P. A.

Análise harmônica de processos estocásticos. Rio de Janeiro: IMPA,

1979. 176 p.

______.

Ondas e Ondaletas: Da análise de Fourier à análise de ondaletas. São Paulo: Udusp,

1999. 272 p.

OLSSON, J.; NIEMCZYNOWICZ, J.; BERNDTSSON, R. Fractal Analysis of High-

Resolution Rainfall Time Series.

Journal off Geophysical Research, v. 98, n. D12,

p. 23, 265-23, 274. December 20, 1993.

OLSSON, J.; BERNDTSSON, R.; SVENSSON, C. Multifractal properties of daily rainfall in

two different climates.

Water Resources Research, v. 32 n. 8, p. 2463-2472,

August 1996.

PEART, C.

Modelling network traffic via a multifractal wavelet model. 2005. -

Dissertation submitted for 1 MSc in networks, non-linear dynamics and data analysis. -

Department of mathematics. University of York, 2005.

PEREIRA, D. C.

Análise dinâmica por wavelets em um sistema com fricção seca e

amortecimento. 2002. - f. 72. Dissertação (Mestrado em Física) - Instituto de

Geociências e Ciências Exatas, Universidade Estadual Paulista, Rio Claro, 2002.

PEREIRA JUNIOR, A.; CHRISTOFOLETTI, A. L. H. Análise fractal da distribuição

espacial das chuvas no Estado de São Paulo.

Geografia, Rio Claro, v. 28, n. 1, p. 97 -

133, jan./abr. 2003.

PETERS, E. E.

Chaos and order in the capital markets: a new view of cycles, prices and

market volatility. Canada: A Wiley Finance Edition, 1996.

PETERS, O.; HERTLEIN, C.; CHRISTENSEN, K. A complexity view on rain.

Physical

Review Letters, v. 88, p. 1-4, 2002.

PEITGEN, H. O.; JURGENS, H.; SAUPE, D.

Chaos and fractals: new frontiers of science.

New York: Springer Verlag, 1992.

PLOTNICK, R. E.; PRESTEGAARD, K.Fractal Analysis of Geological Time Series.

Fractals in Geography

. Nina Siu-Ngan Lam, Lee de Cola. Prentice Hall, 1993.

POLIKAR, R.

The story of wavelets. In: INTERNACIONAL JOINT CONFERENCE ON

CIRCUITS, SYSTEMS, COMMUNICATIONS AND COMPUTERS, 3., 1999,

Athens:

IEEE / IMACS, 1999.

QUINN, D.

Survey of Wavelet Analysis. 1997.- f 94. M. Sc. thesis - Department of Applie

Mathematics, National University of Irelan.

RODRIGUEZ, I. I.; DE POWER, B. F. ...{et al}. Chaos in Rainfall,

Water Resources

Research, v. 25, n.7, p. 1667-1675, 1989.

ROWE, A. C.; ABBOTT, P. C. Daubechies wavelets and mathematica.

Computers in

physics

, Australia, v. 9, n. 6, p. 635 - 648, 1995.

SAUPE, D. Algorithms for randon fractals. pages: 71-113.

The Science of Fractal Images.

Heinz-Otto Peitgen, Dietmar Saupe (eds.). Springer-Verlag, 1988.

SEVCIK, C.

A Procedure to Estimate de Fractal Dimension of Waveforms.Complexity

International, 5, 1998. Disponível em: <http://life.esu.edu.au/complex/ci/vol5/sevcik

SIMONSEN, I.; HANSEN, A.; NES, O. M. Determination of the Hurst exponent by use of

wavelet transforms.

Physical Review E, New York, v. 58, n. 3, p. 2779 - 2787, 1998.

SPIEGEL, M. R.

Análise de Fourier. Mcgraw-Hill do Brasil, 1976.

STRANG, G. Wavelets and Dilatation Equations: A Brief Introduction.

Siam Review, vol 31,

p 613-627, 1989.

STRANG, G. Wavelets.

American Scientist, vol 82, p 250-255, 1994.

STOEV, S; TAQQU, M. S.; PARK, C.; MARRON, J. S. Strengths and limitations of the

wavelet spectrum method in the analysis of internet traffic.

SAMSI technical report #

2004-8, Boston, 2004.

TAYLOR, J. R. An introduction to error analysis. California university science books,

1982.

TOLSTOV, G. P.

Fourier Series. New York: Dover Publication, INC, 1962. 336 p.

TORRENCE, C.

Wavetest Example Matlab Script For Wavelet, using NINO3 SST

DATAS. Disponível em: <http://www.paos.colorado.edu/research/wavelet/>. Acesso em:

19 mar. 2002.

TORRENCE, C.; COMPO, G. P. A pratical guide to wavelet analysis.

Bulletin of the

American Meteorological Society, Colorado, V. 79, n. 1, p. 61-78, janeiro 1998.

TURCOTTE, D. L. Fractals and Chaos in Geology and Geophysics. Cambridge University

Press, (2 edition), 1997.

TRIGO, R. M.; GALVÃO, C. O.; TRIGO, I. F. Aplicação de Wavelets ao estudo da

variabilidade de caudais diários: uma comparação entre algumas bacias portuguesas e

brasileiras.

IV SILUSBA, Coimbra, Portugal, 1999.

VIEIRA, S.

Estatística experimental. São Paulo: Editora Atlas S. A., 1999.

VOSS, R. F. Fractals in nature: from characterization to simulation. pages: 21-70.

The

Science of Fractal Images. Heinz-Otto Peitgen, Dietmar Saupe (eds.). Springer-

Verlag,1988.

WERON, R.; PRZYBYLOWICZ, B. Hurst analysis of eletricity price dynamics.

Physica A, New York, v. 283, p. 462 - 468, 2000.

YEVJEVICH, V. Stochastic Processes in Hidrology.

Water Resources Publications, 276

pp., Mcgraw-Hill do Brasil, 1976.

ZANANDREA, A.; SAOTOME, O.; FREITAS, S. V. Deteção de sinais de pulsos de radar

através de transformada contínua de wavelets.

XXVIII Congresso Nacional de

Matemática Aplicada e Computacional - CNMAC, Porto Alegre, p. 478, 2004.

100

APÊNDICE I

Programa fbmfeder: utiliza a equação 9.25 do livro FRACTALS de Jens Feder, 1988 e foi

codificado para o

software MATLAB por Juha Merikoski.

function Ball=fbmfeder(H,TN,Nruns)

% fbm.m - Fractional Brownian Motion (Juha Merikoski)

% Eq.(9.25) of 'Fractals' by Feder (with error corrected)

% Before run set e.g. as follows

%H=0.7; TN=1000; M=700; n=8; Nruns=5;

%TN numero de pontos gerados

%Nruns numero de series geradas

%dB ruido

%Ball fbm

M=700; n=8;

%disp(['H=' num2str(H) '; TN=' int2str(TN) '; M=' int2str(M) ...

% '; n=' int2str(n) '; Nruns=' int2str(Nruns)]);

Hmh=H-1/2;

normf=1/((n^H)*gamma(H+1/2));

powmax=max([n*TN+n*M]);

ipow=1:1:powmax;

ipow=(ipow/n).^Hmh;

Ball=zeros(Nruns,TN);

dBall=zeros(Nruns,TN);

for irun=1:Nruns,

nxi=n*TN+n*M;

xi=randn(1,nxi);

B=zeros(1,TN);

dB=zeros(1,TN);

for t=1:1:TN,

dB(t)=0;

for i=1:1:n

dB(t)=dB(t)+ipow(i)*xi(1+n*(M+t)-i);

end;

for i=1:n*(M-1),

dB(t)=dB(t)+(ipow(n+i)-ipow(i))*xi(1+n*(M+t-1)-i);

end;

Ball(irun,:)=cumsum(dB)*normf;

dBall(irun,:)=dB*normf;

if mod(irun,Nruns/5)==0,

cloc=clock;% disp([' : ' num2str([irun/Nruns cloc(4:6)])]);

end;

101

APÊNDICE II

Programa mpdwalk: utiliza o algoritmo do deslocamento do ponto médio descrito por VOSS,

PEITGEN et all, 1988 e foi codificado para o

software MATLAB por MIDDLETON, 1995.

function y=mpdwalk(yb,yf,n,h)

% y=mpdwalk(yb,yf,n,H)

% returns a vector y consisting of steps in a fractional random

% walk from y = yb to y = yf. The number of steps is 2^n, i.e.,

% n is the number of times the line is divided into midpoints.

% H is the Hurst coefficient: H = 0.5 for Brownian walks

% 0.5 < H < 1 implies long-range correlation (persistence)

% 0 < H 0.5 implies antipersistence. Press space-bar to start.

% Written by Gerry Middleton, December 1995, using the mid-point

% displacement algorithm.

y=zeros(2^n+1,1);

for i=n:-1:1

xl=1;

dx=2^(i-1);

jlim=2^(n-i);

for j=1:jlim

xr=xl+2*dx;

xm=xl+dx;

y(xm)=0.5*y(xl)+0.5*y(xr)+randn*dx^h;

xl=xr;

end

%plot(y);pause;close

102

APÊNDICE III

Programa fbmsim: baseia-se no artigo "Simulation of stacionary gaussian processes [0,1]

" de

WOOD e CHAN, 1994 e a codificação para o

software MATLAB foi apresentada por

JOHANSSON, 2000.

function BH = fbmsim(gridpts,H,s2)

% BH = fbmsim(gridpts,H,s2)

% fbmsim generates a realisatiofl of fractional Brownian motion, BH, with

% Hurst parameter H in (0,1) and Var(BH(1)»~s2 at the times specified iii

% the vector gridpts which is of the form 0:n where n is an integer > O.

% The simulated process can later easily be rescaled to another grid

% resolution.

% The metod of simulation is based on the article

% Joachim .Johansson

% Dept. Mathematical Statistics

% Chalmers University of Technology

% Gàteborg, Sweden

% 2000-01-05

maxg = 23; % Maximum value of g bef ore we make an

% approximation of the circulant matrix

% The names of the variables coinside vith the notation in the article

% whenever possible.

n = length(gridpts); % Nuniber of timepoints

% Find our value of g

ok = 0;

g = ceil(log(2*(n-1))/log(2));

while (ok == 0) & (g <= maxg)

m=2^g;

t = [(0:(m/2)) (m-((m/2+1):(m-1)))];

c = 0.5*(abs(t-1).^(2*H) + abs(t+1).^(2*H) - 2*(abs(t).^(2*H)));

lambda = fft(c);

if (max(abs(imag(lambda))) < 0.0001) & (min(real(lambda)) > -0.001) ok = 1;

lambda = real(lambda); lambda = max(0,lambda); else

g = g+1;

end;

% We might have to make an approximation of the circulant matrix if

% the value of g is too large.

if (g > maxg)

['No suitable g found, approximating C']

g =ceil(log(2*(n-1))/log(2));

m=2^g;

t = [(0:(m/2)) (m-((m/2+1):(m-1)))];

c = 0.5*(abs(t-1).^(2*H) + abs(t+1).^(2*H) - 2*abs(t).^(2*H));

lambda1 = fft(c);

lambda = max(0,lambda);

rho = sum(lambda)/sum(lambda1);

lambda = rho*lambda;

end;

% Geiierate our realisation of the increment prOCeBS, Y ...

randn('state',sum(100*clock));

U = randn(1,m/2);

V = randn(1,m/2);

tmpl = (1/sqrt(2*m))*sqrt(lambda(2:(m/2))).*(U(2:(m/2)) + i*V(2:(m/2)));

tmp2 = (1/sqrt(2*m))*sqrt(lambda(2:(m/2))).*(U(2:(m/2)) - i*V(2:(m/2)));

tmp2 = fliplr(tmp2);

a = [sqrt(lambda(1)/m)*U(1) tmpl(:)' sqrt(lambda(m/2+1)/m)*V(1) tmp2(:)'];

Y = real(fft(a));

Y =Y(1:n);

BH = sqrt(s2)*(cumsum(Y)-Y(1));

103

APÊNDICE IV

Programa fbmlevinson: utiliza a fatoração de Cholesky/Levinson, foi confeccionado por

GONÇALVES, 1997 e faz parte do FRACTAL ANALYSIS SOFTWARE.

function [x,y,r] = fbmlevinson2(N,H,tmax,sigma,seed) ;

% fbmlevinson

% Levinson synthesis of a fractional Brownian motion

% Paulo Goncalves

% June 6th 1997

% Generates a Fractional Brownian Motion (fBm) using Cholesky/Levinson

% factorization

% 1. Usage

% ______________________________________________________________________

% [x,y,r] = fbmlevinson(N,H,[seed])

% ______________________________________________________________________

% 1.1. Input parameters

% o N : Positive integer

% Sample size of the fBm

% o H : Real in [0,1]

% Holder exponent

% o seed : real scalar

% Random seed generator

% 1.2. Output parameters

% o x : real vector [1,N]

% Time samples of the fBm

% o y : real vector [1,N]

% Vector of N i.i.d. white standard Gaussian r.v.'s (input process of

% the generator)

% o r : real vector [1,N]

% First row of the var/cov Toeplitz matrix R of the increment

% process w[k] = x[k+1] - x[k].

% 2. See also:

% mbmlevinson

% 3. Examples

% % Synthesis of a fractional Brownian Motion with Holder exponent H

% N = 1024 ; H = 0.8 ;

% t = linspace(0,1,N) ;

% clf ;

% [x] = fbmlevinson(N,H) ;

% plot(t,x) ;

% title ('Fractional Brownian Motion - H = 0.8') ;

% xlabel ('time') ;

% INRIA holds all the ownership rights on the Software.

% The scientific community is asked to use the SOFTWARE

% in order to test and evaluate it.

% INRIA freely grants the right to use modify the Software,

% integrate it in another Software.

% Any use or reproduction of this Software to obtain profit or

% for commercial ends being subject to obtaining the prior express

% authorization of INRIA.

% INRIA authorizes any reproduction of this Software.

% - in limits defined in clauses 9 and 10 of the Berne

% agreement for the protection of literary and artistic works

% respectively specify in their paragraphs 2 and 3 authorizing

% only the reproduction and quoting of works on the condition

% that :

% - "this reproduction does not adversely affect the normal

104

% exploitation of the work or cause any unjustified prejudice

% to the legitimate interests of the author".

% - that the quotations given by way of illustration and/or

% tuition conform to the proper uses and that it mentions

% the source and name of the author if this name features

% in the source",

% - under the condition that this file is included with

% any reproduction.

% Any commercial use made without obtaining the prior express

% agreement of INRIA would therefore constitute a fraudulent

% imitation.

% The Software beeing currently developed, INRIA is assuming no

% liability, and should not be responsible, in any manner or any

% case, for any direct or indirect dammages sustained by the user.

% Any user of the software shall notify at INRIA any comments

% concerning the use of the Sofware (e-mail : FracLab@inria.fr)

% This file is part of FracLab, a Fractal Analysis Software

if all(H-H(1))

error('Time-dependant regularity not allowed yet.')

end

switch nargin

case 2

tmax = N-1 ;

sigma = 1 ;

seed = rand(1) * 1e6 ;

shift = 1 ;

case 3

sigma = 1 ;

seed = rand(1) * 1e6 ;

shift = tmax/N ;

case 4

seed = rand(1) * 1e6 ;

shift = tmax/N ;

case 5

shift = tmax/N ;

end

t = linspace(0,tmax,N) ;

s = 0 ;

alpha = 2*H ;

r = sigma^2*(abs(t+shift-s).^alpha + abs(t-s-shift).^alpha - ...

2*abs(t-s).^alpha)/2 ;

randn('seed',seed) ;

y = randn(N,1) ;

x = zeros(1,N) ;

inter1 = r ;

inter2 = [0 r(2:N) 0] ;

Y = y(1)*r ;

k = -inter2(2) ;

aa = sqrt(r(1)) ;

for j = 2:N

aa = aa*sqrt(1-k^2) ;

inter = k*inter2(j:N) + inter1(j-1:N-1) ;

inter2(j:N) = inter2(j:N) + k*inter1(j-1:N-1) ;

inter1(j:N) = inter ;

bb = y(j)/aa ; ;

x(j:N) = x(j:N) + bb*inter1(j:N) ;

k = -inter2(j+1)/(aa^2) ;

end

x = cumsum(x(:)) ;

105

APÊNDICE V

Programa para gerar séries de fbm e estimar o valor de H

% Programa para gerar séries de fbm e estimar o valor de H

% As séries são geradas com 16384 elementos e usados os 7305 centrais

% A escolha do tipo de gerador a ser usado é feita nas linhas 23 a 27

% As séries serão geradas com H= 0.3, 0.4, 0.5, 0.6, 0.7 (variável hu)

% Para cada valor de H serão geradas 100 séries e para cada uma o valor

% de H é estimado

% Programa elaborado por Luiz Roberto Salomão jun/2005

clear

close all

clc

hu=0.3;% primeiro H

nupon=16384;%tamanho da serie gerada;

nusada=7304;%tamanho da serie usada (centro da serie gerada);

nh=5; %número de Hs usados

nsg=100;%número de séries geradas para cada H;

for ih=5:nh

for i100=1:nsg

% xx=função geradora do mov brow

%xx=mpdwalk(0,1,14,hu);%********

[xx] = fbmlevinson(nupon,hu) ;%********

%xx=fbmfeder(hu,nupon,1);%*******

%xx = fbmsim(0:nupon,hu,0.5);%********

%y=xx(4540:11844)';

plot(y(1:nusada));pause;close;

% Corrige pela reta primeiro e ultimo ponto

len=length(y);

nusada=len;

a=polyfit([1,nusada],[y(1),y(nusada)],1);

for v=1:nusada

x(v)=y(v)-(a(1)*v+a(2));

end

plot(x(1:nusada));pause;close;

% Cálculo dos coeficientes da transformada wavelet

escala=floor(log2(len));

coefs=cwt(x,2.^(0:escala),'db8');

for i=1:escala

z(i)=log2(mean(abs(coefs(i,:)))); % log2 abs mean coeficient.

end

% Estimativa de H

eixox=1:escala;

coef=polyfit(eixox,z,1);

s1=num2str(coef(1)-0.5);

f1=polyval(coef,eixox);

ss = ['Estimativa de H:' s1 ];

%plot(eixox,z,'.'),xlabel('log2 scales'),ylabel('log2 mean abs coef')

%hold on

%plot(eixox, f1, 'g'), title(ss)

%hold off;pause;close;

% Cálculo do erro

reta=polyval(coef,eixox);

erroy=sqrt(sum((z-reta).^2)/(length(z)-2));

delta=length(eixox)*sum(eixox.^2)-(sum(eixox))^2;

erroh=erroy*sqrt(length(z)/(delta));

err(ih,i100)=erroh;

hmed(ih,i100)=coef(1)-0.5;

end

hu=hu+0.1;

end

%Salva os valores de H e do erro

save('hmedfbmlevinson.mat','hmed');

save('errffbmlevinson.mat','err');

for i=1:nh

hmedio(i)=mean(hmed(i,:));

erromedio(i)=mean(err(i,:));

disp(['H médio = ' num2str(hmedio(i)) ' Erro médio = ' num2str(erromedio(i))]);

end

106

APÊNDICE VI

Programa para estimar o expoente de hurst de séries temporais de chuva

% Programa para estimar o expoente de hurst de séries temporais de chuva

% utilizando transformada wavelet e fourier

% Desenvolvido por Luiz Roberto Salomão jun/2005

clear

close all

clc

fprintf('\n Estima o H por wavelets (db8) e Fourier %s\n',datestr(now))

fprintf('\n \n')

% Carregando o arquivo que contém a série

ss=input('Entre com o nome do arquivo ( nome.mat ): ','s');

disp(['Arquivo ' ss]);

load(ss);

plot(x0);pause;close;

len=length(x0);

% Correção pela média

media=mean(x0);

x1=x0-media;

plot(x1);pause;close;

% Soma cumulativa

x2=cumsum(x1);

plot(x2);pause;close;

%corrige pela reta que passa pelo primeiro e ultimo ponto

a=polyfit([1,len],[x2(1),x2(len)],1);

for v=1:len

x(v)=x2(v)-(a(1)*v+a(2));

end

plot(x);pause;close;

% Define a escala

esc=floor(log2(len));

escala=2.^(0:esc);

% Calcula os coeficientes da transformada wavelet

coefs=cwt(x,escala,'db8');

% Calcula o log da média do valor absoluto dos coeficientes

for i=1:esc

z(i)=log2(mean(abs(coefs(i,:))));

end

% Construção do gráfico log2 mean coef X log2 escala (estimativa de H)

eixox=1:esc;

coef=polyfit(eixox,z,1);

s1=num2str(coef(1)-0.5);

f1=polyval(coef,eixox);

s2 = ['Estimativa de H por wavelet: ' s1 ];

plot(eixox,z,'.'),xlabel('log2 escala'),ylabel('log2 mean abs coef')

hold on

plot(eixox, f1, 'g'), title(s2)

hold off;pause;close;

% Cálculo do erro ao estimar H

reta=polyval(coef,eixox);

erroy=sqrt(sum((z-reta).^2)/(length(z)-2));

delta=length(eixox)*sum(eixox.^2)-(sum(eixox))^2;

erroh=erroy*sqrt(length(z)/(delta));

s3=[s2 ' Erro= ' num2str(erroh)];

disp(s3)

% Estimativa de H por Fourier

% Cálculo da transformada de Fourier

y = fft(x);

% Cálculo da 'Potência' e log da primeira metade dos coeficientes

yy = y .* conj(y);

p = log(yy(fix(2:len/2)));

% Cálculo do log da freqüencia

f = log(1: len/2 - 1);

% Traçado do gráfico log Potencia X log freqüência

plot(f, p), xlabel('log frequencia'), ylabel('log potencia');

hold on;

Nmin = 16;

Nmax = floor(15*len/32);

p1 = polyfit(f, p, 1);

p2 = polyfit(f(Nmin:Nmax), p(Nmin:Nmax),1)';

s4 = num2str(-(p2(1)+1)*0.5);

s5 = ['Estimativa de H por Fourier: ' s4];

f1 = polyval(p1, f);

f2 = polyval(p2, f(Nmin:Nmax));

107

plot(f(Nmin:Nmax), f2, 'r'), title(s5);

hold off;pause;close;

% Cálculo do erro ao estimar o H por Fourier

length(p(Nmin:Nmax));

length(f2);

erroP=sqrt(sum((f2-p(Nmin:Nmax)).^2)/(length(Nmin:Nmax)-2));

eixofreq=(Nmin:Nmax);

deltafreq=length(eixofreq)*sum((f(Nmin:Nmax)).^2)-(sum(f(Nmin:Nmax)))^2;

errohf=erroP*sqrt(length(eixofreq)/deltafreq);

s6=[s5 ' Erro= ' num2str(errohf)];

disp (s6);

108

APÊNDICE VII

WAVELETS DAUBECHIES

O resumo aqui apresentado para a construção das wavelets Daubechies se

baseia nos trabalhos de DECLAN, 1997, STRANG, 1989, STRANG, 1994, HARPEN, 1998,

MALONE, 2000 e BURRUS et all, 1998.

Dada uma função de dilatação (em inglês

dilation)

+∞<<∞−Φ=Φ

∑

x- para )kx2(h )x(

1-M

(1)

com coeficientes h

constantes para k=0,1,2,3,......M-1 com M par, a função wavelet pode ser

definida por:

∑

Φ=Ψ

1-M

k-1-M

k)-(2x h (-1) )x( (2)

Os valores de h

fora do intervalo {0,1,.......M-1} são iguais a zero.

A condição de ortonormalidade estabelece que:

∫

∞

∞−

δ=−Φ−Φ dx)bx()ax(

com a e b Є Z

e teremos:

∑

−

=δ=

0q2qkk

.0,1,2,3...q 2 h h (3)

A condição que a wavelet tenha decaimento rápido, com M≥2 e M par, nos

leva a:

∑∑

−

1-M

12kk2

1 h h (4)

Existe uma grande quantidade de conjuntos de h

que satisfazem as

equações 3 e 4. Entretanto se considerarmos que a wavelet tem média igual a zero, teremos

um único conjunto de h

109

0,1.......m com 0 h k )1(

k-1-m

−===−

∑

−

(5)

Se os coeficientes satisfazem as equações

3, 4, e 5 podemos a partir da

equação 1 encontrar a equação de dilatação e calcular a equação da wavelet.

No apêndice C do livro "Introduction to wavelets and wavelet transforns"

(BURRUS et all,1998) e no site

http://www.dsp.rice.edu/software/wavebook.shtml (acesso em

18 set. 2006), encontramos a codificação para o software MATLAB, da determinação dos

coeficientes h

(função DAUB), da iteração da função de dilatação (função PSA) e da

construção da wavelet (função WAVE).

Para a wavelet de ordem 8 os coeficientes h

, determinados pela função

DAUB, são:

= 0,2303 h

= 0,7148 h

= 0,6308 h

= -0,0279

= -0,1870 h

= 0,0308 h

= 0,0328 h

= -0,0105

Partindo de uma função constante, usando a função PSA, a função de

dilatação é iterada. A

figura VII-1 mostra a função de dilatação após cada iteração e a figura

VII-2 a wavelet construída a partir da função de dilatação pela função WAVE.

Figura VII-1 - Iterações da função de dilatação

0 500 1000

-2

0 500 1000

0.1

0.2

0 500 1000

-0.5

0.5

0 500 1000

-0.5

0.5

0 500 1000

-1

0 500 1000

-2

0 500 1000

-2

0 500 1000

-2

110

Figura VII-2 - Wavelet Daubechies de ordem 8 construída após oito iterações da função de

dilatação.

0 1 2 3 4 5 6 7

-1

-0.5

0.5

1.5

Wavelet daubechies de ordem 8

Psi(x)

Livros Grátis
( http://www.livrosgratis.com.br )
 
Milhares de Livros para Download:
 
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas

Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
 
 