Localização hipocentral com método probabilı́stico
não-linear num modelo de velocidades 3–D
Ana Luı́sa Grilo Pinho
Dissertação para obtenção do Grau de Mestre em
Engenharia Fı́sica Tecnológica
Júri
Presidente: Prof. João Seixas
Orientador: Prof. João Fonseca
Vogais: Doutora Sandra Heleno
Doutora Susana Vilanova
Novembro de 2008
Resumo
Este trabalho apresenta um método de localização hipocentral para SWIberia + RA.
SWIberia compreende a zona de Portugal Continental com maior actividade sı́smica,
sendo na região compreendida entre o Banco de Gorringe e a costa oeste portuguesa
onde se localizam os epicentros de maior magnitude [27].
A modelação hipocentral foi realizada com o software NonLinLoc, que segue a inversão probabı́listica não-linear para PDF’s LS-L2 e EDT.
Utilizou-se um modelo de velocidades tridimensional de ondas P e S para SWIberia
de Grandin et.al., [1]. Este modelo é reconstruı́do através de um esquema de sólidos
geométricos para o NonLinLoc. Extendeu-se SWIberia pela RA para incluir as estações
dos Açores e Madeira que pertencem à rede IST.
Os tempos de propagação foram calculados através das diferenças finitas à equação
de Eikonal. O algoritmo é adaptado a modelos de velocidade com alto contraste, sendo
sensı́vel à variação da velocidades das ondas entre a crosta oceânica e a continental.
Para testar o método implementado no laboratório, procedeu-se a ensaios com os
sismos de 12/02/2007 e 01/07/2007.
Os resultados finais foram avaliados para três algorimos de cálculo das PDF’s: GridSearch, Metropolis-Gibbs e Oct-Tree.
Compara-se os resultados com os valores dados pelo EMSC e IM. Considerou-se
também os exemplos descritos em [6]. Os melhores resultados foram obtidos com o
método OCT-Tree e PDF EDT.
1
Abstract
The aim of this thesis is to present an hypocenter non-linear location method for SWIberia plus Atlantic area. SWIberia has a special seismic interest because it comprehends
Gorringe Bank, a geological ocean structure with a relevant seismic activity.
It was chosen a nonlinear earthquake location featured by NonLinLoc. This program
follows the probabilistic nonlinear formulation for earthquake location using or LS-L2 and
EDT PDF’s typecast.
A 3–D velocity P and S model by Grandin et.al. for SWIberia was used [1]. Taking this
model, it was extended the P and S wave ocean pattern towards Azores and Madeira’s
region. The intent of this extrapolation was to include the IST stations located at Azores
and Madeira.
The travel times were evaluated with the Eikonal finite-difference scheme developed
for high contrasted velocity models, which is the present case due to wave pattern variation between oceanic and continental crustal.
In order to test the NonLinLoc method, it was performed a seismic signal analysis of
two seismics events occurred at 12/02/2007 and 01/07/2007.
The final results were computed through three different numerical PDF evaluation
algorithms: Grid-Search, Metropolis-Gibbs and Oct-Tree.
The results are compared to EMSC and IM data. It was taken into consideration the
case studies analysed in [6] (págs. 30-46), too. Our best fit was with OCT-Tree method
applied to the EDT PDF.
2
Palavras-chave (Key words)
modelo de velocidades tridimensional; three dimensional velocity model
localização hipocentral; hypocentral locations
NonLinLoc
3
Índice
Resumo
1
Abstract (Resumo em Inglês)
2
Palavras-chave (Key words)
3
Lista de Tabelas
6
Lista de Figuras
8
Lista de Abreviaturas
9
1 Introdução
10
2 Fundamentos Teóricos
2.1 Determinação da malha dos tempos de chegada e
ângulos take-off . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Cálculo dos tempos de propagação numa malha 2–D [2]
2.1.2 Cálculo dos tempos de chegada numa malha 3–D [2] . .
2.1.3 Frentes de onda difractadas [2] . . . . . . . . . . . . . . .
2.1.4 Cálculo dos ângulos take-off [5] [7] . . . . . . . . . . . . .
2.2 Localização hipocentral . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Inversão não-linear [3] [4] [5] . . . . . . . . . . . . . . . .
2.2.2 Estimadores de Gauss [5] [7] . . . . . . . . . . . . . . . .
2.2.3 Algoritmo de cálculo da PDF: Grid-Search [3] [5] . . . . .
2.2.4 Algoritmo de cálculo da PDF: Metropolis-Gibbs [3] [5] . .
2.2.5 Algoritmo de cálculo da PDF: Oct-Tree [5] . . . . . . . . .
12
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
12
12
15
16
17
18
19
24
25
26
28
3 Metodologia do trabalho
31
3.1 Modelo de Velocidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Leitura das fases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3 Funcionamento do NonLinLoc . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4 Resultados Finais
4.1 Modelo de Velocidades e Slow Len . . . .
4.1.1 Secções XY . . . . . . . . . . . . .
4.1.2 Secções XZ . . . . . . . . . . . . .
4.2 Tempos de propagação e ângulos take-off
4.2.1 Tempos de propagação . . . . . .
4.2.2 Ângulos take-off . . . . . . . . . .
4.3 Localizações Hipocentrais . . . . . . . . .
4.3.1 Grid-Search . . . . . . . . . . . . .
4.3.2 Metropolis-Gibbs . . . . . . . . . .
4.3.3 Oct-Tree . . . . . . . . . . . . . . .
4.3.4 Análise dos Resultados . . . . . .
. . . . .
. . . . .
. . . . .
teóricos
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
39
39
39
40
43
43
44
45
47
50
54
57
5 Considerações finais
63
Referências bibliográficas
64
4
A Anexo
67
A.1 Secções XZ do modelo de velocidades . . . . . . . . . . . . . . . . . . . . . 67
A.1.1 Intervalo [32◦ ; 34◦ [ N: zona extrapolada . . . . . . . . . . . . . . . . . 67
A.1.2 Intervalo [34◦ ; 40◦ ] N . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
A.1.3 Intervalo [40◦ ; 40, 5◦ [ N: zona extrapolada . . . . . . . . . . . . . . . . 72
A.2 Tempos de propagação e ângulos take-off teóricos para as estações sismológicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
A.2.1 Secções XY e XZ dos tempos de propagação . . . . . . . . . . . . . 73
A.2.2 Secções XZ dos ângulos take-off . . . . . . . . . . . . . . . . . . . . 76
A.3 Código C do programa de visualização do modelo G RANDIN et.al. para a
plataforma gráfica Gnuplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
A.4 Sintaxe do script para o modelo de velocidades . . . . . . . . . . . . . . . . . 83
A.5 Macro para o SAC das leituras dos picos de fase (exemplo) . . . . . . . . . . 84
A.6 Informação adicional relativa à localização hipocentral dos sismos de 12/02/2007
e 01/07/2007 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
A.7 Batimetria SWIberia + região Atlântica dos Açores e da Madeira . . . . . . . 86
5
Lista de Tabelas
1
2
3
4
5
6
7
Localizações hipocentrais determinadas segundo o método Grid-Search para
o evento sı́smico de 12/02/2007. . . . . . . . . . . . . . . . . . . . . . . . . .
Localizações hipocentrais determinadas segundo o método Grid-Search para
o evento sı́smico de 01/07/2007. . . . . . . . . . . . . . . . . . . . . . . . . .
Localizações hipocentrais determinadas segundo o método Metropolis-Gibbs
com fs = 8, 0 e l=3 para o evento sı́smico de 12/02/2007. . . . . . . . . . . .
Localizações hipocentrais determinadas segundo o método Metropolis-Gibbs
com fs = 9, 0 e l=3 para o evento sı́smico de 12/02/2007. . . . . . . . . . . .
Localizações hipocentrais determinadas segundo o método Metropolis-Gibbs
com fs = 8, 0 e l=65 para o evento sı́smico de 01/07/2007. . . . . . . . . . .
Localizações hipocentrais determinadas segundo o método OCT-Tree para
o evento sı́smico de 12/02/2007. . . . . . . . . . . . . . . . . . . . . . . . . .
Localizações hipocentrais determinadas segundo o método OCT-Tree para
o evento sı́smico de 01/07/2007. . . . . . . . . . . . . . . . . . . . . . . . . .
6
47
48
50
51
52
54
56
Lista de Figuras
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
Algumas estações da rede sı́smica Navigator . . . . . . . . . . . . . . . . . .
Princı́pio de Huygens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Transmissão da frente de onda numa grelha 2–D [2] . . . . . . . . . . . . . .
Transmissão da frente de onda numa grelha 3–D. [2] . . . . . . . . . . . . . .
Zona sombra no ponto P numa malha 2–D. [2] . . . . . . . . . . . . . . . . .
Representação dos ângulos take-off numa geometria plana. . . . . . . . . .
Ilustração do algoritmo Grid-Search [5] . . . . . . . . . . . . . . . . . . . . .
Ilustração do algoritmo MetropolisGibbs [5] . . . . . . . . . . . . . . . . . . .
Ilustração do algoritmo OctTree: subdivisão octal de uma célula i nas suas
células-filhas [5] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Algoritmo OctTree: ilustração das subdivisões sistemáticas [5] . . . . . . . .
Algoritmo OctTree: representação compacta da PDF [5] . . . . . . . . . . . .
Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das
ondas P para a secção XZ, com y = 38, 77◦ . . . . . . . . . . . . . . . . . .
Modelo SWIberia de Raphaël Grandin et.al [1]: curvas de nı́vel das ondas P
para a secção XZ, com y = 38, 77◦ . . . . . . . . . . . . . . . . . . . . . . . .
Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das
ondas P para a secção XY , z = 0km . . . . . . . . . . . . . . . . . . . . . .
Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das
ondas S para a secção XY , z = 0km . . . . . . . . . . . . . . . . . . . . . .
Leitura dos tempos de chegada das ondas P e S na componente vertical (Z)
do sinal recebido pela estação de Manteigas (MTE) para o sismo de 12 de
Fevereiro de 2007 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção XY , z = −1.01km . . . . . . . . . . . . . . . .
Modelo da Slow Len para as ondas P na região SWIberia + região Atlântica
dos Açores e Madeira: secção XY , z = −1.01km . . . . . . . . . . . . . . . .
Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção (parcial) XZ com y = 38, 77◦ . . . . . . . . . .
Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção (parcial) XZ na linha da costa com y = 36, 5◦ .
Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção (parcial) XZ com y = 36, 5◦ , evidenciando o
padrão de velocidades irregular na zona do Banco de Gorringe . . . . . . . .
Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das
ondas P para a secção XZ com y = 36, 5◦ , onde se verificam as irregularidades associadas à zona do Banco de Gorringe. . . . . . . . . . . . . . . . .
Tempos de propagação calculados para a estação sı́smica de Alcochete:
secção XY , z = −1, 01km. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Tempos de propagação calculados para a estação sı́smica de Alcochete:
secção XZ, y = 38, 77◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Ângulos take-off calculados para a estação sı́smica de Alcochete: secção
XZ, y = 38, 77◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
Grid-Search e PDF do tipo EDT. . . . . . . . . . . . . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
Grid-Search e PDF do tipo LS-L2. . . . . . . . . . . . . . . . . . . . . . . . .
7
10
13
13
15
16
17
26
28
29
29
30
32
32
33
35
37
39
40
41
41
42
42
43
43
44
48
49
28
29
30
31
32
33
34
35
36
37
Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
Grid-Search e PDF do tipo EDT. . . . . . . . . . . . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
Metropolis-Gibbs, fs = 8, 0, l=3 e PDF do tipo EDT. . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
Metropolis-Gibbs, fs = 9, 0, l=3 e PDF do tipo EDT. . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
Metropolis-Gibbs, fs = 8, 0, l=65 e PDF do tipo EDT. . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
OCT-Tree e PDF do tipo LS-L2. . . . . . . . . . . . . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
OCT-Tree e PDF do tipo EDT. . . . . . . . . . . . . . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
OCT-Tree e PDF do tipo LS-L2. . . . . . . . . . . . . . . . . . . . . . . . . .
Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
OCT-Tree e PDF do tipo EDT. . . . . . . . . . . . . . . . . . . . . . . . . . .
Visualização 3–D do hipocentro do sismo de 12/02/2007 com o método OCTTree e PDF do tipo EDT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Visualização 3–D do hipocentro do sismo de 01/07/2007 com o método OCTTree e PDF do tipo EDT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
. 49
. 51
. 52
. 53
. 55
. 55
. 56
. 57
. 61
. 62
Lista de Abreviaturas
EDT
EMSC
PDF Equal-Differential-Time
European-Mediterranean Seismological Centre
LS-L2
função de norma L2 com RMS
IM
Instituto de Meteorologia, Portugal
NonLinLoc
Probabilistic, Non-Linear, Global-Search Earthquake Location in 3D Media, software desenvolvido e mantido por Anthony Lomax
PDF
Função de densidade a posteriori
RA
RMS
Região Atlântica
Root Mean Square, função com desvio médio
quadrático
Scatter Sample
conjunto de pontos associados aos valores da
PDF
Seismic Analysis Code
Sudoeste da Penı́nsula Ibérica
SAC
SWIberia
9
1
Introdução
O objectivo deste trabalho foi implementar o método de localização hipocentral mais adequado para a região SWIberia + zona atlântica dos Açores e Madeira no Laboratório de
Sismologia do IST. A dita região de estudo alberga a rede sı́smica Navigator, à qual pertencem as estações do IST1 e do IM. De momento estão em construção as estações de
ViaLonga, Porto, Chã de Alvares, Aouinet Torkoz (Marrocos), Madeira, Açores (Terceira) e
Cabo Verde.
Figura 1: Algumas estações da rede sı́smica Navigator
Na escolha do algoritmo a implementar para o desenvolvimento deste trabalho, foram
tomados em consideração três software’s distintos de localização: o HYPOINVERSE-2000
[39], o hypoDD e o NonLinLoc. Os dois primeiros programas, escritos em Fortran, adoptam
métodos de localização linear. Em particular, o hypoDD desenvolve o método dos mı́nimos
quadrados na construção da matriz de inversão das derivadas parciais do modelo linear
(ver [7], págs. 417 e 418 ou consultar site [38]). No entanto, o software “vencedor” foi o
1 As actuais estações em funcionamento são as de Alcochete (PACT: 38◦ 46.40870 N
Monsanto (PMST: 38◦ 44.08020 N 9◦ 11.09880 W).
10
8◦ 50.11790 W) e a de
NonLinLoc escrito em linguagem de programação C, desenvolvido e mantido por Anthony
Lomax (ver [5]). A versão instalada foi a v4.40.1 04Mar2008.
Contrariamente aos programas anteriores, o NonLinLoc aplica um processo de inversão não-linear que segue a formulação probabı́listica apresentada por Tarantola e Valette [4].
A escolha de um método de localização por inversão não-linear prende-se com a geometria da rede Navigator, fortemente condicionada pelas caracterı́sticas geográficas da
região. A distribuição das estações não é homogénea, concentrando-se a maioria ao longo
da linha da costa continental. Desta forma, os métodos convencionais de inversão linear
tornam-se pouco auspiciosos na obtenção de bons resultados. A existência de grandes
“gap’s” azimutais no mapa da rede torna bastante difı́cil a construção de um modelo teórico
de raytracing de equações linearmente independentes. Por outro lado, o modelo de velocidades 3–D de alto contraste para esta região reforça ainda mais este problema.
O método de inversão não-linear associado ao método de raytracing elaborado por
Podvin e Lecomte [2] vão resolver respectivamente estes dois problemas.
Na secção seguinte, descrevem-se os fundamentos teóricos que estão na base do
algoritmo implementado no NonLinLoc e evidenciam-se as caracterı́sticas mais determinantes na escolha deste método, para o nosso caso de estudo.
11
2
2.1
Fundamentos Teóricos
Determinação da malha dos tempos de chegada e
ângulos take-off
O cálculo dos tempos de chegada das ondas P e S (raytracing) é determinado com base no
algoritmo numérico de Podvin e Lecomte [2]. Para isso considera-se uma malha tridimensional cúbica, onde para cada célula (cúbica) a velocidade é considerada como constante.
Os tempos de propagação entre as estações e todos os pontos da grelha tridimensional
são calculados através do método de diferenças finitas para equação de Eikonal:
(∇t)2 = s2
(1)
onde s é a slowness do meio e t(x) representa o tempo de chegada da frente de onda no
ponto x. Cada iteração corresponde a um ponto de partida diferente na malha 3–D. Em
cada iteração e para cada ponto, determina-se o menor tempo de chegada com base na
interpolação linear dos tempos de chegada de dois pontos vizinhos e no pressuposto da
velocidade ser constante em cada célula.
Este método numérico supera um problema presente em métodos equivalentes de
resolução da (1), nomeadamente o algoritmo de Vidale (1988) [25]: o raytracing para
modelos de velocidade com zonas de alto contraste. A computação dos primeiros tempos de chegada com raytracing num meio heterogéneo é difı́cil porque qualquer ponto
num modelo recebe um conjunto grande de frentes de onda, em que só uma é “admitida”
na continuação do cálculo numérico. No entanto, muitas vezes há certas zonas da grelha slowness que nunca chegam a ser calculadas (“zonas sombra”), porque os sinais são
difractados nas interfaces de meios diferentes.
Existem alguns métodos numéricos que para uma dada fonte numa grelha, calculam
iterativamente diferentes ângulos take-off e escolhem o equivalente a um caminho de menor tempo de propagação. Estes métodos minimizam o problema acima descrito mas são
computacionalmente muito morosos; para além disso, não entram em linha de conta frentes de onda difractadas (e.g. frentes de onda da descontinuidade de Moho) que são muitas
vezes as primeiras ondas registadas ([6], pág.8).
A aplicação sistemática do princı́pio de Huygens aproximado às diferenças finitas tenta
resolver este problema. O princı́pio de Huygens reconhece que todos os pontos de uma
malha são fontes de uma nova frente de onda. Desta forma, todos os pontos vizinhos
de um dado ponto emitem um segundo impulso, isto é, actuam como segundas fontes no
momento em que são atingidos pela frente de onda. Considerando, então, que para cada
ponto existe uma frente de onda independente, o algoritmo determina as descontinuidades
do gradiente temporal de (1) a partir da intersecção local dessas frentes de onda.
As subsecções seguintes descrevem o algoritmo numérico de raytracing para uma malha 2–D, generalizando-se depois para o caso 3–D; apresentam o conceito de ângulos
take-off associado à determinação do raytracing e descrevem o seu cálculo numérico.
2.1.1
Cálculo dos tempos de propagação numa malha 2–D [2]
Considere-se a malha 2–D quadrada representada na figura 3.
A diferença temporal entre dois pontos vizinhos M e N da grelha refere-se à transmissão
da frente de onda plana através da interface MN. Se s é a slowness numa das células
adjacentes a essa interface, então:
12
Figura 2: Princı́pio de Huygens
Figura 3: Transmissão da frente de onda numa grelha 2–D [2]
13
∂2t
∂2t
+ 2
2
∂x
∂y
(∇t)2 =
∂t
tN − tM
=
∂x
h
onde h representa o tamanho da aresta da malha. Inserindo (3) em (2), tem-se que:
r
(tN − tM )2
∂t
= ± s2 −
∂y
h2
(2)
(3)
(4)
∂t
As duas soluções de ∂y
tornam ambı́gua a direcção de propagação, determinando se
a onda atinge ou não o ponto P. Duas condições são necessárias para que a frente de
onda alcance o ponto P:
i) A distância máxima que a frente de onda pode
percorrer dentro de uma célula é o
√
valor da hipotenusa do quadrado, isto é, 2h. Logo, pela equação (2), o tempo de
propagação máximo da frente de onda dentro da célula é:
2(hs)2 = (tN − tM )2 + (tP − tM )2
(5)
Como a malha é quadrada e s é constante dentro da célula, então (tN − tM )2 =
(tP − tM )2 , pelo que:
(tN − tM )2 =
(hs)2
2
⇒
hs
0 6 (tN − tM ) 6 √
2
Por outro lado e atendendo à equação (4),
s2 −
(tN − tM )2
> 0
h2
⇔
0 6 (tN − tM ) 6 hs
Intersectando as condições anteriores, tem-se a primeira condição:
⇒
hs
0 6 (tN − tM ) 6 √
2
(6)
Esta condição assegura que a onda, que passa por MN, possa alcançar o ponto P.
Mesmo que a onda passe para a célula adjacente NPQR, s0 pode ser suficientemente
baixo para que ocorra uma reflexão total, gerando uma frente de onda na interface
hs
NP (ver lei de Snell). Caso (tN − tM ) > √
, a frente de onda passa para a célula
2
acima sem passar por P (ver Figura 3).
ii) A segunda condição é:
∂t
>0
∂y
(7)
Esta condição assegura que ~s está orientado para P.
No entanto, recorrendo ao princı́pio de Fermat podemos verificar que somente a primeira condição é estritamente necessária. O princı́pio de Fermat postula que o caminho de
uma onda corresponderá sempre ao caminho com o menor tempo de propagação possı́vel
em relação a caminhos adjacentes. Logo, o tempo de chegada de dada frente de onda ao
14
ponto P será sempre maior ou igual ao primeiro tempo de chegada para este ponto. Por
∂t
outro lado, se ∂y
6 0 significaria que o tempo de chegada para o ponto P seria inferior em
relação ao do ponto N, resultando na violação da equação de Eikonal. Desta forma, só a
condição (6) é que garante a “iluminação” do ponto P.
Assim, e pela equação (4), o tempo de chegada ao ponto P é:
p
tP = tN + (hs)2 − (tN − tM )2
(8)
2.1.2
Cálculo dos tempos de chegada numa malha 3–D [2]
O cálculo dos tempos de propagação numa grelha 3–D segue os mesmo princı́pios teóricos
descritos para o caso 2–D. Consider-se a seguinte figura 4.
Figura 4: Transmissão da frente de onda numa grelha 3–D. [2]
Quatro frentes diferentes frentes de onda são consideradas na transmissão através do
plano MNPQ; cada frente de onda envolve o cálculo dos tempos de chegada em três
cantos da interface.
Desta forma, quatro topologias são consideradas em cada célula para a determinação
dos tempos de chegada ao ponto R, onde cada topologia relaciona a transmissão local
duma frente de onda através de metade da interface MNPQ. Tal como em 2–D, a condição
imposta é assegurar a ligação pela frente de onda da interface MNPQ ao ponto R.
Por analogia à condição de tN − tM para o caso 2–D expressa por (6), para assegurarmos que a frente de onda que atravessa a célula passa pelo ponto R, então, para a
primeira topologia, tem-se que:
tM ≤ tN ,
tM ≤ tP ,
p
(hs)2 − (tN − tM )2
√
(tP − tM ) ≤
2
p
2
(hs) − (tP − tM )2
√
(tN − tM ) ≤
2
indicando que o ângulo de projecção da slowness não pode exceder os 45 ◦ .
15
(9)
(10)
Pode-se, então, deduzir que para cada topologia o tempo de chegada ao ponto R é:
MNP → R
tM ≤ tN ,
tM ≤ tP ,
p
tR = tN + tP − tM + (hs)2 − (tN − tM )2 − (tP − tM )2
(11)
QN P → R
tN ≤ tQ ,
tR = tQ +
tP ≤ tQ ,
q
(hs)2 − (tQ − tN )2 − (tQ − tP )2
(12)
NMQ → R
0 ≤ (tN − tM ) ≤ (tQ − tN ),
q
tR = tQ + (hs)2 − (tQ − tN )2 − (tN − tM )2
(13)
PMQ → R
0 ≤ (tN − tM ) ≤ (tQ − tN ),
q
tR = tQ + (hs)2 − (tQ − tP )2 − (tP − tM )2
(14)
Na estrutura 3–D, são consideradas 24 interfaces para cada ponto, sendo, por isso,
calculada 96 transmissões possı́veis de cada onda plana. (As interfaces são aquelas que
pertencem às células que dado ponto pertence mas às quais esse ponto não pertence.)
2.1.3
Frentes de onda difractadas [2]
Figura 5: Zona sombra no ponto P numa malha 2–D. [2]
O ponto M actua como segunda fonte pelo princı́pio de Huygens.
Se um dado ponto P estiver localizado numa zona sombra para um determinado padrão
local de frentes de onda que passam por M, então, pelo princı́pio de Huygens esse ponto
M vai agir como uma segunda fonte, emitindo uma onda difractada para o ponto P (ver
Figura 2.1.3). O tempo de chegada será:
16
√
tP = tM + hs 2
(15)
Analogamente para 3–D, tem-se (ver Figura 4):
√
tR = tM + hs 3
(16)
Este cálculo pode ser efectuado sistematicamente pois o seu valor será sempre igual
ou superior ao primeiro tempo de chegada do ponto a determinar.
2.1.4
Cálculo dos ângulos take-off [5] [7]
Figura 6: Representação dos ângulos take-off numa geometria plana.
Os dois raios sı́smicos que atinjem a interface mostram a relação existente entre o ângulo
de incidência i e o declive função do tempo de propagação T (x). (ver [7], pág. 134)
Observe a figura 6 que representa a propagação de ondas sı́smicas numa geometria
plana. Quando a onda chega a uma interface, é possı́vel determinar o ângulo de incidência
i que o raio sı́smico constrói entre si e a normal a essa interface. Note-se que, como se
está a considerar a geometria plana, isso significa que as interfaces são vistas como planos
que separam diferentes meios por onde a onda sı́smica se propaga.
Da figura 6 extrai-se facilmente que:
sin i =
vdT
dx
(17)
logo,
dT
sin i
=
(18)
dx
v
Substituindo v pela velocidade da onda ao sair da fonte, o ângulo take-off será proporcional
ao declive da função T (x).
Para o caso de estudo presente neste trabalho, considera-se que quando o raio está
num ponto exactamente vertical ao ponto da interface seguinte, então i = 180◦ .
Assim, os ângulos take-off relativos aos primeiros tempos de chegada para cada ponto
são estimados a partir dos gradientes dos tempos de propagação. Dois gradientes são
estimados para cada direcção x, y e z: Gbaixo entre um dado ponto e o respectivo ponto
vizinho precedente ao longo do eixo e Galto entre esse ponto e o ponto seguinte ao longo
17
desse mesmo eixo. O gradiente total Geixo ao longo do eixo é a média desses dois gradientes. O gradiente total ao longo dos três eixos determina o vector gradiente da função
tempo de propagação.
É também determinado um factor de qualidade Qeixo entre 0 e 10 pela seguinte expressão:
20Gbaixo Galto
(19)
Qeixo =
G2baixo G2alto
Se Qeixo 6 0 (i.e. os gradientes têm sinais opostos), Qeixo = 0; se Qeixo = 10, então os
dois gradientes têm a mesma magnitude e o mesmo sinal. É calculado um factor de qualidade final para os ângulos take-off a partir da média ponderada dos factores de qualidade
ao longo de cada eixo, onde o peso de cada factor é dado pela média do gradiente ao
longo do eixo:
|Gx |Qx + |Gy |Qy + |Gz |Qz
Q=
(20)
|Gx | + |Gy | + |Gz |
2.2
Localização hipocentral
A localização hipocentral baseia-se essencialmente na inversão da equação que relaciona
os dados recebidos com as coordenadas do hipocentro a determinar:
(21)
d = M (p)
Os métodos de inversão linear convencionais passam por determinar iterativamente a
matriz G, correspondente às derivadas parciais em primeira ordem da expansão em série
de Taylor da equação (21) em torno de um valor de p determinado na iteração anterior.
O valor hipotético de p para a primeira iteração é calculado segundo o modelo escolhido
para a região de estudo. Desta forma, tem-se que:
(22)
∆d = G∆p
Invertendo G, pelo ao método dos mı́nimos quadrados, obtém-se a inversa G−g . (Notese que não é possı́vel inverter trivialmente G pois trata-se certamente de uma matriz n × 4
com n > 4, sendo, por isso, uma matriz não-quadrada.) A equação numérica a calcular é:
∆p = G−g ∆d
(23)
iterando até que o resı́duo ∆d seja mı́nimo.
G−g é definido por:
G−g = (GT G)−1 GT
(24)
T
Para modelos teóricos de equações linearmente não-independentes, a matriz G G será
singular pois terá um ou mais valores próprios nulos; consequentemente, o seu determinante é zero e portanto esta matriz não poderá ser invertida. A decomposição de Lanczos permite-nos determinar a inversa generalizada de matrizes singulares, no entanto o
espaço das soluções será a dos vectores próprios correspondentes a valores próprios nãonulos. Assim, o espaço das soluções da inversão linear será inferior ao espaço dos dados
observados, havendo porções do modelo que não poderão ser detectadas. Para além
disso, os resultados finais terão um erro associado à imprecisão do modelo, porque não
serão incluı́dos na inversão os vectores próprios associados aos valores próprios nulos.
([7], páginas 416-429)
18
A formulação probabilı́stica não-linear possibilita a utilização de modelos de velocidades tridimensionais heterogéneos, permitindo a determinação de resultados fidedignos
com uma boa estimativa da incerteza associada. O cálculo não-linear hipocentral envolve
a construção de uma função de densidade de probabilidade para os parâmetros a estimar. Esses parâmetros são as três coordenadas espaciais de localização do hipocentro
e o tempo de origem do sismo. A representação hiper-volúmica desta localização inclui
multiplas soluções óptimas para modelos de elevado contraste, tendo, portanto, uma forma
bastante irregular. Os métodos numéricos lineares de localização hipocentral mais comuns
produzem uma única solução com uma incerteza associada, que segue a distribuição normal. Tal solução só será credı́vel, caso o modelo tenha uma única solução óptima e cuja
incerteza seja aproximadamente elipsoidal. Desta forma, o método de localização nãolinear permite a identificação de uma função de densidade de probabilidade com uma
forma volumétrica irregular, sendo por isso mais adequado ao nosso caso de estudo.
Note-se que o problema de localização hipocentral é inerentemente um problema nãolinear. Genericamente, pode-se definir o tempo de chegada tobs como:
Z
tobs = t0 +
u(r0 )ds
(25)
r0 (s)
onde u representa a slowness num espaço não-homogéneo e r0 (s) representa um ponto
à distância s ao longo do caminho r0 entre a fonte sı́smica e as estações receptoras. A
equação (25) é não-linear, porque uma modificação na localização da fonte irá provocar
uma modificação no caminho que é integrado
Nas subsecções seguintes, descreve-se o método probabilı́stico de inversão não-linear
e os algoritmos computacionais Grid-Search, Metropolis-Gibbs e Oct-Tree utilizados para
o cálculo da função de densidade de probabilidade a posteriori (PDF) com norma LS-L2
ou norma EDT.
2.2.1
Inversão não-linear [3] [4] [5]
Função Densidade de Probabilidade
Considere-se um sistema fı́sico discreto L no sentido amplo do termo, caracterizado pelas suas variáveis, pela famı́lia de parâmetros associados às medidas instrumentais e respectivos outputs. Assume-se que L pode ser definido por uma famı́lia de parâmetros
X = {X1 , ..., Xm }. Qualquer conjunto de valores possı́veis para esses parâmetros é representado por x = {x1 , ..., xm }, caracterizando um possı́vel modelo para L no espaço
m-dimensional E m . Note-se que a parametrização de um sistema não é única, podendo
haver duas parametrizações equivalentes. Por exemplo, a velocidade v e a slowness s
são dois parâmetros equivalentes. A escolha entre dois parâmetros equivalentes pode ser
vista também como diferentes escolhas do sistema de coordenadas em E m .
O grau de conhecimento que temos acerca dos valores dos parâmetros do nosso sistema podem variar entre um estado de total ignorância e um estado de total conhecimento.
Assim, a descrição de qualquer estado de conhecimento sobre os valores de X é postulada por:
Z
m(A) =
f (x)d(x) (A ∈ E m )
(26)
A
que é absolutamente contı́nua na medida de Lesbegue sobre m(E m ). A quantidade m(A)
é denominada por “medida” de A. Se m(E m ) é finita, então f (x) pode ser normalizada, tal
19
que m(E m ) = 1; neste caso, f (x) é denominada por “função de densidade de probabilidade” para A.
Função de densidade para o estado de informação nula de um sistema
Introduz-se, agora, o conceito de informação nula de um dado sistema. Para isso, considerese o problema da localização de um dado hipocentro para um conjunto de dados recolhidos. As coordenadas do sistema são cartesianas (X, Y, Z. A pergunta é: qual será a
função de densidade µ(x, y, z) que contém a menor informação da localização do hipocentro? A resposta é: aquela que tem a mesma probabilidade dP para todas as regiões de
igual volume dV :
dP = const.dV
(27)
Em coordenadas cartesianas dV = dx dy dz, logo a equação (26) dá a solução:
µ(x, y, z) = const.
(28)
Conjunção de dois estados independentes
Seja si e sj dois estados de informação (obtidos independentemente um do outro), representados pelas funções fi , fj . Considere-se também µ expresso pela equação (28).
Por definição, a conjunção dos estados si e sj é um estado de informação representado
pela função de densidade f (x) dado por:
f (x) =
fi (x)fj (x)
µ(x)
(29)
A função de densidade f (x) não é necessariamente normalizável, mas para o presente
caso considera-se fi (x), fj (x) e f (x) normalizáveis. Desta forma, (29) é, de facto, uma
função de densidade de probabilidade.
Função de densidade a priori
Para a famı́lia de parâmetros X que descreve L , pode-se definir uma subclasse genérica
de parâmetros D = (D1 , ..., Dr ) referentes aos tempos de chegada medidos nas estações.
As coordenadas do hipocentro desconhecidas podem ser escritas como P = (P1 , ..., P4 ).
Então, define-se a função de probabilidade de densidade para determinado conjunto de
valores d e p como:
ρ(x) = ρ(d, p)
(30)
A equação (30) é a função de densidade a priori que representa os resultados medidos
para um conjunto de parâmetros p a determinar.
20
Função de densidade do modelo teórico
A equação análoga à equação (21) na formulação probabilı́stica é:
θ(x) = θ(d, p)
(31)
Tendo em conta a definição de probabilidade condicional para acontecimentos independentes, a função de densidade de d dado p é:
θ(d|p) =
θ(d, p)
θp (p)
(32)
onde θp (p) é a função de densidade marginal para P. Reescrevendo (31), tem-se:
θ(d, p) = θ(d|p)θp (p)
(33)
Na maioria dos casos, considera-se que a probabilidade e P é igual para qualquer
zona da região estudada, ou seja, não se possui qualquer informação adicional sobre P
que o constrinja. Assim, aplicando a definição de estado de informação nula apresentado
anteriormente, tem-se:
θ(d, p) = θ(d|p)µp (p)
(34)
onde µp (p) é a função de densidade de informação nula sobre P no sistema.
Função de densidade a posteriori (PDF)
A conjunção de ρ(x) e θ(x) dá-nos um novo estado de informação sobre o sistema. Assim,
e pela definição expressa por (29), o estado de informação a posteriori é dado pela PDF:
σ(x) =
ρ(x)θ(x)
µ(x)
(35)
Separando a famı́lia de parâmetros X nas subclasses X = (D, P), a equação (35)
pode ser reescrita como
σ(d, p) =
ρ(d, p)θ(d, p)
µ(d, p)
(36)
A partir desta equação, pode-se determinar as PDF’s marginais normalizadas para D
e P:
Z
ρ(d, p)θ(d, p)
σd (d) =
dp
(37)
µ(d, p)
Z
ρ(d, p)θ(d, p)
σp (p) =
dd
(38)
µ(d, p)
Na maior parte dos casos, a informação a priori de D é independente da de P:
ρ(d, p) = ρd (d)ρp (p)
Inserindo, (34) e (39) na equação (38), tem-se:
21
(39)
Z
σp (p) = ρp (p)
ρd (d)θ(d|p)
µp (p)
dd
µ(d, p)
Pela definição de informação independente expressa na equação (39),
a equação anterior simplifica-se da seguinte forma:
Z
ρd (d)θ(d|p)
dd
σp (p) = ρp (p)
µd (d)
(40)
µp (p)
µ(d,p)
=
1
µd (d) ,
logo
(41)
Desta forma, a equação (41) é a solução probabilı́stica completa para uma dada coordenada hipocentral p a determinar.
Aproximação Guassiana
Visto que a distribuição gaussiana é amplamente usada na resolução de inúmeros problemas fı́sicos, apresenta-se de seguida a aproximação gaussiana para as funções (30) e
(31) respectivamente:
(
)
1
T −1
ρ(t|X, Y, Z, T ) = exp − (t − t0 ) Ct (t − t0 )
(42)
2
(
T
)
1
−1
θ(t|X, Y, Z, T ) = exp − t − g(X, Y, Z, T ) CT t − g(X, Y, Z, Y )
(43)
2
onde t0 são os tempos de chegada medidos nas estações; t − g(X, Y, Z, T ) são os tempos de propagação estimados computacionalmente a partir do modelo de velocidades2 ;CT
é a matriz covariante dos erros associados aos tempos de propagação calculados e Ct é
a matriz covariante dos erros associados às medições dos tempos de chegada3 . Note-se
que para uma teoria exacta, isto é, CT é nula (43) toma a probabilidade máxima de 1; para
o nosso caso, isso seria equivalente a termos um modelo de velocidades “perfeito”.
Considerando que para os tempos de propagação estimados (segundo o modelo de
velocidades escolhido), a função de densidade para o estado de informação nula de x é
constante, i.e.:
µ(x) = const.
(44)
então, a PDF para x é determinada por substituição directa das equações (42) e (43) nos
termos da equação (35). Calculando o respectivo integral, obtém-se a solução geral para o
problema espaço-temporal da localização do hipocentro para uma distribuição gaussiana
dos valores medidos nas estações:
(
σ(X, Y, Z, T ) = ρ(X, Y, Z, T )exp
2A
)
1
T
−1
− (t0 − t(x)) (Ct + CT ) (t0 − t(x))
2
expressão refere-se à solução trivial do sistema homogéneo de (21).
que os valores obtidos nas estações seguem uma distribuição gaussiana.
3 Assume-se
22
(45)
Equações de norma L2 implementadas no software NonLinLoc
A equação (45) é implementada no programa NonLinLoc só para as componentes espaciais. A função de densidade de x para as coordenadas espaciais é obtida integrando
(45) no intervalo de [−∞; +∞] o tempo de origem T :
Z ∞
σ(X, Y, Z, T )dT
(46)
σ(x) =
−∞
Assume-se também, que a informação a priori de T na origem não poderá ser independente, contrariamente à informação espacial, que está directamente relacionada com
a configuração tectónica da região ou com a ocorrência de sismos anteriores em zonas
especı́ficas dessa mesma região. Logo,
ρ(X, Y, Z, T ) = ρ(T ).ρ(X, Y, Z) = ρ(X, Y, Z)
(47)
O cálculo do tempo de chegada entre o ponto (X, Y, Z) à estação i é:
gi (X, Y, Z, T ) = hi (X, Y, Z) + T
(48)
onde hi é o tempo de propagação entre o ponto (X, Y, Z) e a estação i.
Inserindo (47) e (48) em (46) e integrando, tem-se
(
)
1
σ(x) = Kρ(x)exp − f (x) , K = const.
2
T
⇒ f (x) = t̂0 − t̂(x) (Ct + CT )−1 t̂0 − t̂(x)
P
j
j wj t o
i
i
t̂0 = t0 − P
j wj
P
j
j wj h
ĥi = hi − P
j wj
X
wi =
wij ;
wij = [(Ct − CT )−1 ]ij
(49)
(50)
j
(51)
onde K é um factor de normalização e f (x) é uma função L2 .
A equação de determinação das coordenadas espaciais x = (x, y, z) implementada
pelo NonLinLoc segue a formulação de Tarantola e Valette descrita na equação (49), recorrendo às diferenças quadráticas entre os tempos observados e os tempos calculados
previamente, dando, portanto, o erro nominal associado a cada observação. Assim, podese reescrever σ(x) da seguinte forma:
2 !
1 X T obsi (x) − T calci (x)
pdf (x) = Kρ(x) −
(52)
2
σi2
obsi
onde σi é o erro associado à obsi . Portanto, a equação (52) é a PDF L2 −RM S (LS-L2) para
as coordenadas espaciais com aproximação Gaussiana implementada pelo NonLinLoc.
A probabilidade máxima para o tempo de origem é determinada pela função Tml (x)
segundo Moser et al. (1992) [12]:
23
XX
Tml (x) =
i
wij (t0 )i − ti (x)
j
XX
i
(53)
wij
j
Função densidade “Equal-Differential-Time” (EDT)
O software NonLinLoc disponibiliza uma PDF com uma norma diferente construı́da pelo
método das hipérboles [14]. Este método é mais robusto que as equações de norma LS-L2
para observações com resı́duos maiores que o seu erro nominal. Tal como nas equações
LS-L2, assume-se que os erros dos tempos de chegada observados e dos tempos de
propagação calculados seguem a distribuição normal.
A PDF EDT implementada é:
"
pdf(EDT ) (x) = Kρ(x)
X
p
obsa obsb
1
σa2 + σb2
2 !#N
[T obsa (x) − T obsb (x)] − [T T calca (x) − T T calcb (x)]
exp −
σa2 + σb2
(54)
onde T obsa e T obsb são os tempos de chegada e T T calca e T T calcb são os tempos de
propagação, respectivamente, para as observações a e b. O somatório é efectuado para
todas as combinações de pares possı́veis. N é o número total de observações. Na exponencial, o primeiro termo em parêntesis rectos expressa a diferença entre os tempos de
chegada observados nos locais a e b; o segundo termo traduz a diferença entre os tempos
de chegada dos tempos de propagação calculados teoricamente para os locais a e b segundo o modelo de velocidades utilizado. Este termo é zero quando ambos os termos são
iguais, e portanto para esses dois pontos a exponencial atingirá o valor máximo de 1. As
coordenadas de x estimadas para o caso em que os termos da diferença são iguais (daı́
a PDF se denominar por “Equal Differential Time” (EDT) são as que melhor satisfazem os
valores temporais obtidos para a e b. Como o somatório sobre as observações é efectuado fora da exponencial, a EDT é maior na região de x que satisfaz o maior número de
pares de observações, não sendo, por isso, sensı́vel a observações com um grande erro
associado [6].
A EDT não determina o tempo de origem T . Ao separar o cálculo da componente
temporal, reduz-se a equação aos parâmetros espaciais 3–D tornando a equação mais
robusta e minimizando erros associados a um possı́vel trade-off entre a profundidade e
tempos de chegada. Este cálculo é efectuado após se ter determinado o valor de x que
satisfaz a solução máxima da PDF.
2.2.2
Estimadores de Gauss [5] [7]
Uma vez que se assume uma estrutura gaussiana para os dados observados, é efectuada
uma estimativa dos parâmetros “tradicionais” da estatı́stica de Gauss: o valor esperado e
os semi-eixos que definem o elipsóide de confiança.
Para uma grelha com xi,j,k , o valor esperado é:
24
E(x) = ∆V
X
xi,j,k σ(xi,j,k )
(55)
i,j,k
∆V é o volume da célula da grelha. Para N amostras distribuı́das segunda a PDF com
localizações xn :
E(x) =
1 X
xn
N n
(56)
onde σ(xi,j,k não são necessários, uma vez que a própria amostra já segue a PDF. Em
ambos os casos, a matriz covariante é dada por,
i
h
(57)
C = E x − E(x) · (x − E(x))T
A elipsóide com 95% de confiança é obtido através da determinação dos valores e
vectores próprios das matrizes quadradas C T C e CC T :
C = U[diag wi ]VT
(58)
T
onde U é a matriz quadrada dos vectores próprios de C C e V é a matriz quadrada dos
vectores próprios de CC T ; [diag wi ] são os valores próprios iguais para ambas as matrizes.
Os semi-diâmetros das elipsóides com 95% de confiança são [7, 815wi ]1/2 , onde 7, 815 é o
valor de 95% confiança para três graus de liberdade.
Estes parâmetros da estatı́stica de Gauss são bons indicadores das incertezas associadas à localização, somente se a PDF não-linear tiver um único máximo e forma elipsoidal.
Os resultados finais irão mostrar isso.
2.2.3
Algoritmo de cálculo da PDF: Grid-Search [3] [5]
Este algoritmo efectua uma procura sistemática da zona onde a PDF assume maiores
valores. Na primeira iteração analisa toda a região, partindo para subsecções desta nas
iterações seguintes onde a PDF é mais elevada. A figura 7 ilustra o método.
As zonas encontradas nas iterações seguintes estão centradas no nodo com o valor hipocentral óptimo segundo a PDF usada. Apesar das zonas determinadas em cada
iteração serem obviamente menores que as das iterações anteriores, as grelhas que as
definem têm mais nodos que as anteriores.
Em cada iteração, o algoritmo avalia o valor da PDF para cada nodo usando ou a
equação (52) ou a (54), conforme a escolha do utilizador. Para isso, ele calcula previamente os tempos de propagação para cada ponto da grelha, assumindo como ponto
inicial o tempo observado em cada estação e guarda o resultado no disco. Para todas as
estações este cálculo é efectuado. Para minimizar o tempo de cálculo, esses tempos são
guardados no disco em ficheiros com grelhas YZ ao invés de o guardar num único ficheiro
com uma malha 3–D. Assim, o programa avalia sistematicamente para um conjunto de “folhas” 2–D com o mesmo x e correspondentes às diferentes observações, os pontos onde
a PDF é maior, construindo a grelha com os novos pontos para a iteração seguinte.
Para a grelha final, os valores de todos os pontos são normalizados assumindo que o
integral da PDF nessa malha é 1.
As vantagens são:
25
Figura 7: Ilustração do algoritmo Grid-Search [5]
1. O método não recorre a derivadas parciais (contrariamente aos métodos lineares),
logo o cálculo não é complexo sendo possı́vel usá-lo em modelos de velocidades
3–D heterogéneos.
2. Fácil implementação: o método é sistemático e linear em cada iteração para a nova
sub-região.
3. O método de leitura em planos 2–D permite que o método seja eficiente na leitura
das grelhas 3–D dos tempos de propagação.
As desvantagens são:
1. O tempo de cálculo é grande.
2. O método não tem boa resolução para grelhas finais muito grandes ou pode originar
truncatura da PDF para grelhas finais pequenas.
3. Requer uma boa escolha do tamanho da grelha e do espaço entre os nodos.
Este método numérico só avalia as coordenadas espaciais do hipocentro.
2.2.4
Algoritmo de cálculo da PDF: Metropolis-Gibbs [3] [5]
O método Metropolis-Gibbs baseia-se no cálculo de um percurso directo no espaço (x, y, z)
que tende na direcção da localização desejada: aquele que maximiza a PDF.
Para cada iteração, o ponto local xactual é perturbado por um vector dx de direcção
arbitrária e com dado comprimento l. A probabilidade σ(xnovo ) é determinada para a
nova localização e comparada com a probabilidade σ(xactual ). Se σ(xnovo ) ≥ σ(xactual ),
26
então a nova localização é aceite com a probabilidade Pnovo = σ(xnovo ); se σ(xnovo ) <
σ(xnovo )
σ(xactual ), então a nova localização é aceite com probabilidade Pnovo = σ(x
. Quando
actual )
uma nova localização é aceite, torna-se na localização xactual e é guardada com o valor
respectivo da PDF no disco.
Para efeitos de localização hipocentral, se a região de estudo for muito grande, a PDF
pode variar bastante dentro desta. Logo, é importante inicializar o algoritmo com um l alto
o suficiente para que o método possa explorar toda a região.
O método Metropolis-Gibbs implementado no NonLinLoc obedece a três fases distintas:
1. Fase de leitura: l é fixo a um valor relativamente alto; o percurso calculado cobre
sensivelmente toda a região de estudo, caminhando para as zonas onde a PDF é
elevada. Os valores da PDF para cada ponto não são guardados.
2. Fase de equilı́brio: l é ajustado considerando a distribuição espacial dos pontos
encontrados na segunda metade da fase de leitura; são determinados os desvios
padrão sx , sy , sz dessa distribuição. Nas iterações seguintes: os desvios padrão são
reavaliados incluindo as coordenadas do local da iteração anterior; l é recalculado
através da expressão fs (sx sy sz /Ns )1/3 , onde Ns é o número de valores da PDF a
serem aceites durante a fase seguinte de registo e fs é um factor de escalonamento.
Os valores de Ns e fs são escolhidos pelo utilizador. Assim, assegura-se que l é
redimensionado para o comprimento correcto das Ns células com lados sx , sy , sz , a
serem avaliadas na fase seguinte. O caminho segue para as zonas com valores de
PDF mais altos.
3. Fase de registo: l é fixo com o valor calculado na última iteração da fase de equilı́brio.
O caminho continua a explorar zonas com valores altos para a PDF. Os valores da
PDF são registados até ao valor máximo Ns . No entanto, nem todos os valores
calculados são guardados; existe um certo intervalo de tempo entre cada registo
para garantir uma independência entre os valores obtidos.
As vantagens deste método são:
1. Não recorre ao calculo de derivadas parciais contrariamente aos métodos lineares,
sendo por isso indicado para modelos de velocidade 3–D complexos, como é o nosso
caso de estudo.
2. Boa precisão para PDF’s de formas irregulares.
3. Cerca de 10x mais lento que os modelos lineares, mas 100x mais rápido que o
método Grid-Search.
As desvantagens são:
1. Para PDF’s irregulares, não-elipsoidais e com vários máximos locais, a cobertura
deste método pode ser inconsistente.
2. A cobertura estocástica deste método pode omitir caracterı́sticas importantes da malha, determinantes para o cálculo do hipocentro. Por exemplo, o método pode seguir
um caminho que o leve ao um máximo local e não ao máximo absoluto da PDF. Por
outro lado, se reconhecer o máximo absoluto, poderá não reconhecer os restantes
máximos locais, não se obtendo uma caracterização completa da PDF.
27
Figura 8: Ilustração do algoritmo MetropolisGibbs [5]
3. Requer uma boa selecção de Ns e fs .
4. Como lê quase por completo a grelha 3–D dos tempos de propagação, corre lentamente para o grande número de observações e para uma malha muito grande.
Este método numérico só avalia as coordenadas espaciais do hipocentro.
2.2.5
Algoritmo de cálculo da PDF: Oct-Tree [5]
O método Oct-Tree usa um processo de divisão octal recursiva da malha 3–D em células
cada vez mais pequenas, onde a densidade da PDF no centro dessas células vai aumentando.
A probabilidade do hipocentro localizar-se em dada célula é:
Pi = Vi [pdf (xi )]
(59)
onde Vi é o volume da célula e xi são as coordenadas no centro desta.
Numa primeira fase, o método avalia o valor da PDF para um dado conjunto inicial
de células com dado volume, de forma a que estas cubram toda a região de estudo. As
probabilidade determinadas são colocadas numa lista LP , que as ordena por magnitude.
A célula com maior P é dividida em oito novas células (ver Figura 9). O valor da PDF
é novamente avaliado para essas novas oito células e os seus valores são inseridos e
ordenados na lista LP . Escolhe-se a célula com P mais elevada e o processo repete-se.
Este método converge rapidamente no espaço 3–D, produzindo uma estrutura octal de
células com valores de PDF especı́ficos. Veja-se a figura 10.
No final, os centros das células construı́das ao longo do processo iterativo são desenhados sobre a estrutura octal, reproduzindo de forma elucidativa a PDF (ver Figura 11).
As vantagens deste método são:
28
Figura 9: Ilustração do algoritmo OctTree: subdivisão octal de uma célula i nas suas
células-filhas [5]
Figura 10: Algoritmo OctTree: ilustração das subdivisões sistemáticas [5]
29
Figura 11: Algoritmo OctTree: representação compacta da PDF [5]
1. Muito mais rápido que o método Grid-Search num factor de 1/100.
2. Mais completo que o modelo Metropolis-Gibbs
3. É um método bastante simples e com poucos parâmetros a definir. Para além disso,
os parâmetros a definir são de fácil determinação.
As desvantagens deste método são:
1. Pode tornar-se lento para um grande número de observações, malhas grandes e/ou
um grande número de células iniciais.
2. Pode eventualmente não localizar máximos locais. É necessário ter atenção aos
parâmetros numéricos que definem a dimensão da grelha.
Este método numérico só avalia as coordenadas espaciais do hipocentro. No entanto,
é possı́vel vir a ser extendido no sentido de efectuar localizações a 4–D que incluam a
coordenada temporal (ver [6], pág.16).
30
3
3.1
Metodologia do trabalho
Modelo de Velocidades
O modelo de velocidades foi construı́do a partir do padrão de curvas de nı́vel extraı́das
das secções XZ do modelo de velocidades de SWIberia proposto por Raphaël Grandin
et.al [1] (ver Figuras 12 e 13). Os limites deste modelo são [34◦ ; 40◦ ]N, [5◦ ; 35◦ ]W e
[0 ; 60]km em profundidade. Como SWIberia não inclui os Açores e a Madeira (ver Figura
14), extendeu-se o modelo até estas regiões. Para extrapolar o modelo à zona atlântica dos
Açores, prolongou-se o padrão das curvas de nı́vel da longitude de 13◦ W até à longitude de
35◦ W. Em relação à região atlântica sul, que inclui a Madeira, assumiu-se que o padrão
das curvas de nı́vel era idêntico à secção XZ de 34◦ N, isto é, à secção mais a sul do
modelo SWIberia; as curvas de nı́vel, das secções a sul de 34◦ N, sofreram apenas um
ligeiro deslocamento para a esquerda de acordo com a geografia da região à superfı́cie.
Para incluir o registo sı́smico da estação de Manteigas (fornecido pelo IM), extendeu-se
o modelo meio grau para norte. Por outro lado, a profundidade de SWIberia toma como
origem o nı́vel médio das águas do mar; atendendo a que há estações a uma altitude
máxima de ∼ 1, 01km, “subiu-se” as coordenadas dos pontos em z = 0km para z =
−1, 01km (a coordenada z aumenta em profundidade). Assim, os limites do modelo de
velocidades implementado neste trabalho são: [32◦ ; 40, 5◦ ]N, [5◦ ; 35◦ ]W e [0 ; ∼ 63]km.
A malha tridimensional, construı́da para o modelo de velocidades, tem um espaço entre nodos de 4km nas três direcções. Este valor foi escolhido por forma a tornar o cálculo
do modelo computacionalmente viável. Para além do tempo de cálculo computacional dispendido, outras condicionantes foram tomadas em consideração. A caracterização dessas
condicionantes é explicada mais à frente no contexto da descrição do software NonLinLoc.
Desta forma, a profundidade máxima do modelo foi avaliada segundo a condição para a
origem de z, explicitada anteriormente, e o passo de 4km da grelha 3–D.
Implementou-se um programa em C para visualização do modelo de velocidades SWIberia. O programa executa a rotina de aplicação de leitura do modelo, disponibilizada pelos
autores, e carrega o resultado para uma outra rotina que permitirá a construção de gráficos
a partir do programa Gnuplot (para consultar código, ver A.3).
Construiu-se um algoritmo tridimensional de sólidos geométricos segundo a sintaxe
do programa NonLinLoc. Definiram-se, então, secções XZ com um dado padrão de
polı́gonos. A forma destes polı́gonos obedece ao padrão das curvas de nı́vel das secções
consideradas. As curvas de nı́vel têm uma variação consecutiva de 2km/s. As secções sucessivas, que variam em Y , contêm somente pequenas variações geométricas por forma
a que as faces de cada secção consecutiva se liguem entre si, construindo o sólido pretendido. X representa a longitude, Y a latitude e Z a profundidade. As pequenas variações
geométricas de secção para secção permitiram “desenhar” as variações do padrão de velocidades ao longo de Y . Portanto, o sólido representa uma parcela da nossa região onde
a velocidade varia constantemente segundo um gradiente definido. Atribui-se também a
cada sólido uma profundidade de referência e, a essa profundidade, faz-se corresponder
uma dada velocidade de referência (para consultar sintaxe do script construı́do, ver A.4).
Este esquema tridimensional foi construı́do dentro de um modelo de velocidades 1–D por
camadas. As velocidades atribuı́das ao modelo 1–D seguem aproximadamente os valores
referidos na tabela 1,pág.1149 de [1].
Este esquema de sólidos geométricos implementado é adaptável a outro modelo de
velocidades, embora não seja totalmente flexı́vel para regiões com um padrão de velocidades muito dı́spar. Para se construir um novo modelo, basta redefinir as coordenadas
31
Figura 12: Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das
ondas P para a secção XZ, com y = 38, 77◦ .
Esta secção corresponde à localização em latitude da estação de Alcochete.
Mapa de curvas de nível de Vp no plano XZ para 38,77° de latitude
Estação de Alcochete (PACT)
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
8
6
4
2
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Figura 13: Modelo SWIberia de Raphaël Grandin et.al [1]: curvas de nı́vel das ondas P
para a secção XZ, com y = 38, 77◦ .
Esta secção corresponde à localização em latitude da estação de Alcochete.
32
dos vértices, colocar a profundidade de referência desejada a cada sólido, atribuir-lhe uma
velocidade e aplicar um gradiente de velocidades ao sólido.
Como se extrapolou o modelo de velocidades do oceano até aos Açores e Madeira, foi
necessário introduzir um factor de correcção nos pontos da malha referentes às coordenadas geográficas das estações respectivas. A tabela de velocidades das ondas sı́smicas
em [1] (pág.1149) indica que, no oceano, a velocidade das ondas P dos 0 − 4 km está no
intervalo 1, 5 6 Vp 6 2 km/s. Por outro lado, consultando a carta geológica de Portugal
e simultâneamente a tabela de velocidades de [1], verifica-se que a velocidade das ondas
P à superfı́cie das ilhas é ≈ 4km/s. Calculando o tempo demorado pela onda a percorrer
4km/s em ambos os casos e subtraindo, obtém-se um factor de atraso τ = 1, 67s. Este
resı́duo é negativo porque o tempo de chegada deveria ser menor. Para efeitos de cálculo
do hipocentro dos eventos de 12/02/2007 e 1/07/2007, este factor é irrelevante pois não
temos fases lidas nessas estações.
O modelo de velocidades de Raphaël Grandin et.al. para SWIberia não contém nenhuma transformação de coordenadas, isto é, 1 ◦ w 110km. O NonLinLoc aplica transformações
de coordenadas aos cálculos efectuados (o tipo de tranformação é escolhido pelo utilizador). Assim, é necessário construir o modelo de velocidades com base nessa transformação
no sentido de georeferenciá-lo com os resultados obtidos através do programa. Uma
vez que o modelo de velocidades é contruı́do com a ligação de secções XZ (contantes
em latitude) e considerando que as faces dos polı́gonos têm de ser planas (i.e., não podem apresentar qualquer tipo de curvatura ou conterem ângulos concâvos), optou-se pela
transformação TransSimple segundo a nomenclatura do NonLinLoc. Essa tranformação
obedece às seguintes expressões:
Figura 14: Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das
ondas P para a secção XY , z = 0km
33
x = long − long(origem) .111, 111. cos[lat(radianos)]
y = (lat − lat(origem) ).111, 111
lat(origem) + y
lat =
111, 111
long(origem) + x
long =
111, 111. cos[lat(radianos)]
(60)
(61)
(62)
(63)
Como se pode verificar pelas equações (60), esta tranformação corrige as distâncias da
longitude em função da latitude. Assim, as latitudes dos pontos de dada secção permanecem constantes entre si, pois somente se aplicou um factor de correcção da latitude comum a todos os pontos desse plano. Evita-se desta forma que surjam curvaturas ao longo
das secções construı́das para o modelo de velocidades lido pelo NonLinLoc. Por exemplo,
a transformação de Lambert, embora muito conveniente para regiões com grande desenvolvimento em longitude face à latitude (como é o caso da nossa região de estudo), provocaria grandes distúrbios na leitura do modelo. Implementou-se um pequeno programa que
efectua a transformação de coordenadas dos pontos do modelo geométrico em coordenadas de Lambert . O objectivo seria averiguar/confirmar o comportamento do programa. O
resultado foi o previsto: o programa, como não alocava devidamente as coordenadas dos
vértices, não os reconhecia, lançando erro como output.
Extraiu-se também do modelo SWIberia o padrão de velocidades para as ondas S. A
Figura 15 exemplifica muito bem a verosimilhança do modelo das ondas S em relação ao
das ondas P na secção XY à superfı́cie. O padrão prolonga-se de forma indêntica em
profundidade. Como se pode verificar, somente a magnitude é que difere. Assim, não foi
necessário construir um modelo de velocidades distinto para as ondas S. Utilizando o valor
de VVPS = 1, 73 (ver [1], pág. 1149), o valor da velocidade das ondas S para cada ponto da
malha é automaticamente determinado.
3.2
Leitura das fases
A leitura das fases passa essencialmente por extrair dos sinais recebidos pelas estações
os tempos de chegada das ondas P e S para um dado evento sı́smco.
Neste trabalho estudaram-se os eventos sı́smicos de 12/02/2007 e de 01/07/2007.
As fases dos tempos de chegada foram analisadas com o programa Seismic Analysis Code (SAC) [26]. O objectivo do tratamento do sinal é basicamente fazer a sua
desconvolução para obter o ground motion. Construiu-se uma macro, segundo a sintaxe
do próprio programa (ver código na seccção A.5), que tratava o sinal da seguinte forma:
• Retirar a média do sinal
• Retirar a inclinação
• Desconvoluir o sinal usando a função de transferência da estação sı́smica receptora
para assim se obter o ground motion do evento sı́smico em análise.
Retirar a média e a inclinação do sinal é remover o valor de offset inerente a este.
No domı́nio da frequência, a média e a inclinação são funções rectas que, ao passarem
para o domı́nio temporal através da transformada de Fourier, provocam grandes picos de
34
Figura 15: Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das
ondas S para a secção XY , z = 0km
amplitude nos pontos onde os sinais estão em fase. Desta forma, torna-se conveniente
removê-las para facilitar a leitura do sinal.
O sinal recebido pela estação é da forma:
Y (ω) = X(ω) ∗ F (ω)
(64)
onde X(ω) é o ground motion e F (ω) a função de resposta da estação ao sinal recebido,
usualmente denominada por função de transferência4 . A função de transferência tem a
seguinte forma:
(ω − ω1 )(ω − ω2 )
F (ω) = CA
(65)
(ω − ω3 )(ω − ω4 )
onde A é uma constante de normalização, C a constante de conversão de bit’s para m/s,
ω1 e ω2 os zeros e ω3 e ω4 os pólos. Sabendo o valor de AC, dos zeros e dos pólos,
define-se F (ω) e desconvolui-se o sinal, obtendo-se o pretendido sinal de ground motion.
Analisou-se os sinais desconvoluı́dos nas três direcções: este-oeste (E), norte-sul (N) e
vertical (Z) (ver Fig. 16). Identificou-se o primeiro sinal relativo às ondas P e S; definiramse as respectivas chegadas como emergente ou impulsiva; caracterizou-se o primeiro movimento como sendo de compressão ou dilatação (ver pág. 5 de [7]); determinaram-se os
tempos de chegada para cada onda através da leitura do eixo xx dos sinais; mediu-se a
coda do sinal, o perı́odo e amplitude pico-a-pico das ondas P . Para medir a amplitude e
o perı́odo, fez-se um corte ao sinal nos primeiros 5 segundos e mediu-se directamente a
amplitude pico-a-pico nesse intervalo de tempo. Este método simples serve para eliminar
4 Por aproximação, considera-se nesta análise o valor de X(ω) como o sinal oriundo da fonte sı́smica. Na
prática, este sinal é a convolução do sinal oriundo da fonte sı́smica com os efeitos de reflecção nas interfaces do
meio que percorre e com os efeitos de atenuação inelásticos.
35
a possibilidade de estarmos a medir amplitudes de ondas P com a sobreposição de ondas
S, de ondas reflectidas e principalmente de ondas de superfı́cie.
Para o evento de 12/02/2007, introduziram-se 8 leituras de ondas P e S e no evento de
01/07/2007, introduziram-se 6 leituras de cada (cada par refere-se ao sinal registado por
uma estação). Atendendo a que a leitura dos tempos de propagação foi feita manualmente,
introduziu-se certamente bastantes erros, principalmente no que respeita às ondas S, pois
estas surgem no final da coda de P e possivelmente em sobreposição com ondas pP ,
SP e ondas de superfı́cie (ver [7], p.2, Fig.1.1-3). Para minimizar este efeito, escolheu-se
para cada evento a direcção que exibia uma leitura mais clara em média para todas as
estações. Assim, para o evento de 12/02/2007 utilizaram-se as leituras dos picos de fase
da componentes vertical; para o evento de 01/07/2007, utilizaram-se as leituras dos picos
de fase da componente Norte-Sul.
3.3
Funcionamento do NonLinLoc
O funcionamento do NonLinLoc baseia-se essencialmente na utilização de comandos
padrão especı́ficos deste software.
De seguida, descreve-se sumariamente o esquema de trabalho com este programa:
1. Vel2Grid: Leitura do modelo de velocidades (construı́do segundo a sintaxe do NonLinLoc). Após o programa ter corrido este comando produz dois ficheiros *.hdr (header ) e *.buf. (binário). O ficheiro binário contém o modelo de velocidades ou de
Slow Len (conforme a escolha do utilizador) e segue a rotina de leitura do NonLinLoc.
A construção do modelo de velocidades somente importa para efeitos de modelação.
O modelo da Slow Len é que é efectivamente usado nos cálculos seguintes. Esta
variável é definida como:
1
Slow Len = . dgrid
(66)
v
onde dgrid = 4km é o valor do passo da malha, igual nas três direcções x, y e z. Uma
vez que a velocidade v é constante dentro de cada célula da malha (ver secção 2.1),
a Slow Len representa o tempo médio de chegada a um dado ponto a partir de um
ponto anterior adjacente.
2. Grid2Time: Determinação dos tempos de propagação numa grelha definida previamente com o modelo da Slow Len para a região de estudo; imediatamente depois,
são determinados os ângulos take-off. Os tempos de propagação são calculados
em relação às coordenadas das estações sı́smicas colocadas como input. Analogamente ao Vel2Grid, o Grid2Time produz 2 pares de ficheiros *.hdr e *.buf: um dos
pares refere-se à malha dos tempos de propagação; o outro refere-se à malha dos
ângulos take-off. Note-se que nesta fase, ainda não estamos a fazer localização hipocentral. Estamos apenas a determinar o padrão dos tempos de propagação, para
todos os pontos da grelha, definida pelo modelo de velocidades (ver secção 2.1.1).
3. NLLoc: Cálculo do(s) hipocentros(s) com base na leitura dos ficheiros com os picos
de fase para as ondas P e S; cada ficheiros refere-se a um evento sı́smico. Este
comando avalia os tempos de chegada observados nas estações e introduz estes
dados no algoritmo de inversão junto com os tempos de propagação, determinados
anteriormente com o comando Gri2Time.
36
Figura 16: Leitura dos tempos de chegada das ondas P e S na componente vertical (Z) do
sinal recebido pela estação de Manteigas (MTE) para o sismo de 12 de Fevereiro de 2007
37
Para cada um dos outputs supra-citados, é possı́vel obter uma visualização dos resultados através do script Grid2GMT. Este comando permite converter os valores obtidos em
ficheiros de leitura interpretáveis pelo GMT, por forma a este gerar os gráficos pretendidos.
Na modelação das velocidades, assegurou-se que o tamanho do ficheiro binário produzido pelo comando Vel2Grid não excedia o aconselhado. Sendo R o tamanho do ficheiro
◦
*.buf, então R . RAM
n◦ est , onde RAM é a memória primária do computador e n est é o
número de estações consideradas em determinada localização.
38
4
Resultados Finais
O modelo de velocidades implementado bem como os sequentes resultados obtidos estão
em coordenadas rectangulares.
4.1
Modelo de Velocidades e Slow Len
De seguida, mostram-se os mapas gerados para o modelo de velocidades implementado
nas secções XY e XZ. Atendendo à dimensão da malha, colocou-se nesta secção somente duas secções XZ. A primeira refere-se à latitude da estação de Alcochete (PACT) e
a segunda representa um corte na zona do Banco de Gorringe. Em anexo, expõe-se detalhadamente as secções XZ explicitamente construı́das para o modelo e, para efeito comparativo, inclui-se também as respectivas curvas de nı́vel do modelo SWIberia (à excepção
das secções no intervalo [32◦ ; 34◦ [ N que foram extrapoladas).
Apresenta-se também o gráfico da secção XY para z = −1.01km da Slow Len determinado a partir do modelo de velocidades.
4.1.1
Secções XY
Figura 17: Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção XY , z = −1.01km
39
Figura 18: Modelo da Slow Len para as ondas P na região SWIberia + região Atlântica
dos Açores e Madeira: secção XY , z = −1.01km
4.1.2
Secções XZ
Compare-se a verosimilhança do padrão representado na Figura 19 com o de SWIberia
para a mesma secção (ver Figura 13). As figuras 20 e 21 confirmam também a expressão
irregular das curvas de nı́vel referente à zona do Banco de Gorringe caracterizada pelo
modelo SWIberia na figura 22.
Observando com atenção a figura 21, verifica-se que há uma descontinuidade no modelo de velocidades entre os 14km e 33km de profundidade: entre os 10km e 14km de
profundidade a velocidade varia entre os 7, 1 e 7, 6km/s; entre os 14km e os 33km, a velocidade diminui para um intervalo entre 6, 6 e 7, 1km/s; para profundidades superiores,
a velocidade volta a subir. Trata-se de uma anomalia na construção do modelo (provavelmente devido a uma definição errónea das cotas dos vértices que poderá provocar a
formação de ângulos concâvos nos polı́gonos ou então a necessidade de construir mais
sólidos para definir as irregularidades desta zona). Não foi possı́vel rectificar em tempo
útil esta anomalia; possivelmente será necessário reconstruir na totalidade esta zona. No
entanto, não é um problema crı́tico, uma vez que a variação de velocidades registada na
camada de 14 − 33km não é elevada. Pela lei se Snell, sin(ic ) = vv21 , onde ic é o ângulo
crı́tico e v1 e v2 são as velocidades nas camadas de incidência e refracção da onda, respectivamente. Assim, neste caso, somente para ângulos de incidência próximos de 90◦ é
que se verifica uma reflexão total. Por isso, é pouco provável esta camada funcionar como
“guia de ondas”, pois somente para um pequeno intervalo de ângulos próximos dos 90◦ é
que o fenómeno da reflexão total se verifica.
40
Figura 19: Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção (parcial) XZ com y = 38, 77◦ .
Esta secção corresponde à localização em latitude da estação de Alcochete. O padrão
nas extremidades do gráfico prolonga-se de igual forma até aos limites do modelo.
Figura 20: Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção (parcial) XZ na linha da costa com y = 36, 5◦ .
O padrão nas extremidades do gráfico prolonga-se de igual forma até aos limites do modelo.
41
Figura 21: Modelo de velocidades das ondas P para a região SWIberia + região Atlântica
dos Açores e Madeira: secção (parcial) XZ com y = 36, 5◦ , evidenciando o padrão de
velocidades irregular na zona do Banco de Gorringe
Mapa de curvas de nível de Vp no plano XZ para 36,5° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Figura 22: Modelo SWIberia de Raphaël Grandin et.al [1]: padrão de velocidades das ondas P para a secção XZ com y = 36, 5◦ , onde se verificam as irregularidades associadas
à zona do Banco de Gorringe.
42
4.2
4.2.1
Tempos de propagação e ângulos take-off teóricos
Tempos de propagação
Em seguida, apresentam-se os tempos de propagação para a estação de Alcochete (PACT)
nas secções XY com z = −1, 01km e XZ com ∼latitude referente à dita estação. Tal como
o esperado, as zonas circundantes, onde se registam menores tempos de propagação,
são as que se encontram mais próximas da estação.
Figura 23: Tempos de propagação calculados para a estação sı́smica de Alcochete:
secção XY , z = −1, 01km.
(PACT: 38◦ 46.40870 N 8◦ 50.11790 W )
Figura 24: Tempos de propagação calculados para a estação sı́smica de Alcochete:
secção XZ, y = 38, 77◦ .
(PACT: 38◦ 46.40870 N 8◦ 50.11790 W )
43
4.2.2 Ângulos take-off
Nesta subsecção, apresenta-se a tı́tulo de exemplo os ângulos take-off calculados para a
estação de Alcochete. Tal como o esperado, a zona onde se observa a verticalização do
padrão de cores corresponde à localização da dita estação. Os resultados referentes às
outras estações encontram-se em anexo (ver secção A.2.2).
Figura 25: Ângulos take-off calculados para a estação sı́smica de Alcochete: secção XZ,
y = 38, 77◦ .
(PACT: 38◦ 46.40870 N 8◦ 50.11790 W )
Os ângulos take-off para os sismos estudados foram determinados juntamente com as
localizações. Para ver os resultados, consultar tabela em anexo A.6.
44
4.3
Localizações Hipocentrais
Neste capı́tulo, analisam-se os melhores resultados para os eventos sı́smicos de 12/02/2007
e 01/07/2007 obtidos com os três diferentes métodos numéricos de avaliação das PDF’s. A
qualidade dos resultados foi comparada com os respectivos dados fornecidos pelo EMSC5
e pelo IM 6 .
De seguida, apresenta-se a informação relevante fornecida por estes dois orgãos:
**********************************Dados do EMSC************************************
EARTHQUAKE on 12/02/2007 at 10:35 (UTC)
AZORES-CAPE ST. VINCENT RIDGE 178 km SW Sagres
MAGNITUDE: Mw 6.1
Data provided by: BRA GFU IMP INGV IRSA LDG LED LJU LVV MAD
MCSM NEIC NEWS OGS RNS SKO SOF THE
Latitude = 35.81 N
Longitude = 10.26 W
Origin Time = 10:35:22.5 (UTC)
Depth = 32 Km
95% confidence ellipse:
- Semi major = 5.1 Km
- Semi minor = 2.1 Km
- Azimuth of major axis = 169 degrees
Number of data used = 487
EARTHQUAKE on 01/07/2007 at 19:03 (UTC) AZORES-CAPE ST. VINCENT RIDGE 302
km W Sagres
MAGNITUDE: mb 4.9
Data provided by: BUC GFZ IMP LJU MAD NEIC ODC ZAMG
Latitude = 36.32 N
Longitude = 12.21 W
Origin Time = 19:03:06.8 (UTC)
5 Os
dados são fornecidos pelo Alert Systems do EMSC; para mais informações consultar [31] e [32]
informação adicional, consultar [33] e [34]
6 para
45
Depth = 2 Km
95% confidence ellipse:
- Semi major = 9.0 Km
- Semi minor = 4.3 Km
- Azimuth of major axis = 171 degrees
Number of data used = 171
***************************************************************************************
**********************************Dados do IM**************************************
12-Fev
Ho 10:35:25.5 ML 5.9
Lat 35.933◦ N Lon 10.495◦ W Prof 37
Região: Mar de Marrocos
Sentido: V Sagres
AxMax 4.2 AxMin 2.4 Ang 119 ErrZ 12 (elipse de confiança (90%), em km)
01-Jul Ho 19:03:15.2 ML 4.9
Lat 36.568◦ N Lon 11.986◦ W Prof 10
Região: Gorringe
Sentido: IV Lisboa
AxMax 5.2 AxMin 3.6 Ang 76 ErrZ 0 (elipse de confiança (90%), em km)
***************************************************************************************
Os resultados estimados pelo EMSC seguem um modelo de velocidades 1–D por camadas7 . Não foi possı́vel apurar o método de cálculo de localizações hipocentrais utilizado
por este organismo. O serviço Alert Systems8 é activado quando no mı́nimo duas estações
lançam o registo de um sismo com magnitude superior a um dado threshold (na Europa o
threshold é v 5.0), dentro da sua região de análise. O sismologista de serviço analisa os
dados fornecidos pelas várias redes que estão ligadas ao EMSC, seleccionando os dados
por forma a eliminar redundâncias (ver [32], p.10); revê os valores de magnitude dados e
calcula a localização com base em todos os tempos dos picos de chegada enviados pelas
várias redes (ver [32], p.17). Uma descrição pormenorizada dos parâmetros utilizados para
efectuar localizações é apresentada no documento [32], p.9. Provavelmente, o método de
cálculo usado é o tı́pico método linear (ver secção 2.2) pois trata-se de um método desenvolvido e automatizado há mais tempo pela maioria dos programas de análise sı́smica.
Relativamente aos dados disponibilizados pelo IM, as localizações foram estimadas
utilizando o programa Hypocent [36]. Este método implementa o esquema iterativo de
Geiger adaptado à técnica dos mı́nimos quadrados. O modelo de velocidades adoptado é
um modelo 1–D baseado no estudo da refracção das ondas em profundidade em Portugal
Continental (ver [35], págs. 576 e 577, cap. S EISMICITY DATA).
7 consultar
8 para
http://www.emsc-csem.org/index.php?page=euromed&sub=velocity
mais informações do EMSC, consultar [31] e [32]
46
4.3.1 Grid-Search
Apresentam-se, de seguida, os resultados obtidos com o método numérico Grid-Search
para os eventos sı́smicos em estudo. As PDF’s vêm representadas com 1000 amostras
(cada amostra corresponde a um ponto vermelho da scatter sample).
O estudo foi efectuado para ambas as PDF’s. No entanto, para o sismo 12/02/2007 não
foi possı́vel a obtenção de um resultado concreto com a PDF LS-L2. O programa produzia
um overflow.
Relativamente ao sismo de 01/07/2007, observa-se uma maior constrição no resultado
da LS-L2 em relação ao da EDT. Este resultado deve-se à diferente natureza topológica
das duas PDF’s.
A análise detalhada destes resultados é descrita na secção 4.3.4.
PDF
TO (hhHmm ss,ss)
PDF HypLoc
LS-L2
-
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Gauss HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Ell95 (km)
-
EDT
10H35 24,38
35,71
-10,75
26,99
35,80
-10,51
28,57
39,90
55,30
185,00
Tabela 1: Localizações hipocentrais determinadas segundo o método Grid-Search para o
evento sı́smico de 12/02/2007.
Parâmetro numérico: n amostras = 1000
Legenda da tabela:
TO –> Tempo de origem; PDF HypLoc –> Coordenadas espaciais do hipocentro obtidas através da
determinação do valor máximo da PDF; Gauss HypLoc –> Coordenadas espaciais do hipocentro obtidas através
dos estimadores de Gauss; Ell95 –> Comprimento dos semi-eixos do elipsóide de erro (95% intervalo de
confiança)
47
Figura 26: Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
Grid-Search e PDF do tipo EDT.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
PDF
TO (hhHmm ss,ss)
PDF HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
Gauss HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
Ell95 (km)
LS-L2
19H03 13,68
35,96
-12,36
14,99
35,96
-12,36
14,99
0,14
0,27
11,10
EDT
19H03 12,83
35,92
-12,19
2,99
32,87
-18,00
29,28
29,30
51,80
394,00
Tabela 2: Localizações hipocentrais determinadas segundo o método Grid-Search para o
evento sı́smico de 01/07/2007.
Parâmetro numérico: n amostras = 1000
Legenda da tabela:
TO –> Tempo de origem; PDF HypLoc –> Coordenadas espaciais do hipocentro obtidas através da
determinação do valor máximo da PDF; Gauss HypLoc –> Coordenadas espaciais do hipocentro obtidas através
dos estimadores de Gauss; Ell95 –> Comprimento dos semi-eixos do elipsóide de erro (95% intervalo de
confiança)
48
Figura 27: Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
Grid-Search e PDF do tipo LS-L2.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha -> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
Figura 28: Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
Grid-Search e PDF do tipo EDT.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
49
4.3.2 Metropolis-Gibbs
Nesta secção, apresentam-se os resultados para o método Metropolis-Gibbs.
Estes resultado são meramente exemplificativos, pois a variação ligeira do step inicial l
provocava alterações bastante acentuadas no método. Como se pode verificar, o valor de
l = 3 para o sismo de 12/02/2007 e l = 65 para o sismo de 01/07/2007. Por isso, não foi
possı́vel encontrar um valor que tornasse este método fiável para a região em estudo.
Utilizou-se o valor de fs = 8.0 porque era o valor recomendado na literatura do programa (ver [5], sec. Control File –> NLLoc –> “LOCSEARCH statement”). Variou-se este
parâmetro, ligeiramente para 9.0 para o sismo de 12/02/2007 e o resultado melhorou ligeiramente (por comparação com os resultados do EMSC).
Também, neste caso, não se obteve resultados com a LS-L2 e, agora, para ambos os
sismos. O programa convergia para valores da PDF nulos.
Para mais detalhes, consultar a secção 4.3.4.
PDF
TO (hhHmm ss,ss)
PDF HypLoc
LS-L2
-
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Gauss HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Ell95 (km)
-
EDT
10H35 24,12
35,71
-10,77
24,44
35,71
-10,75
37,90
6,31
8,08
37,30
Tabela 3: Localizações hipocentrais determinadas segundo o método Metropolis-Gibbs
com fs = 8, 0 e l=3 para o evento sı́smico de 12/02/2007.
Legenda da tabela:
TO –> Tempo de origem; PDF HypLoc –> Coordenadas espaciais do hipocentro obtidas através da
determinação do valor máximo da PDF; Gauss HypLoc –> Coordenadas espaciais do hipocentro obtidas através
dos estimadores de Gauss; Ell95 –> Comprimento dos semi-eixos do elipsóide de erro (95% intervalo de
confiança)
50
Figura 29: Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
Metropolis-Gibbs, fs = 8, 0, l=3 e PDF do tipo EDT.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
PDF
TO (hhHmm ss,ss)
PDF HypLoc
LS-L2
-
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Gauss HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Ell95 (km)
-
EDT
10H35 24,33
35,71
-10,76
30,72
35,70
-10,23
24,08
26,20
37,60
67,00
Tabela 4: Localizações hipocentrais determinadas segundo o método Metropolis-Gibbs
com fs = 9, 0 e l=3 para o evento sı́smico de 12/02/2007.
Legenda da tabela:
TO –> Tempo de origem; PDF HypLoc –> Coordenadas espaciais do hipocentro obtidas através da
determinação do valor máximo da PDF; Gauss HypLoc –> Coordenadas espaciais do hipocentro obtidas através
dos estimadores de Gauss; Ell95 –> Comprimento dos semi-eixos do elipsóide de erro (95% intervalo de
confiança)
51
Figura 30: Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
Metropolis-Gibbs, fs = 9, 0, l=3 e PDF do tipo EDT.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
PDF
TO (hhHmm ss,ss)
PDF HypLoc
LS-L2
-
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Gauss HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
-
Ell95 (km)
-
EDT
19H03 13,61
35,93
-12,20
2,39
35,66
-12,43
33,44
23,20
50,40
100,00
Tabela 5: Localizações hipocentrais determinadas segundo o método Metropolis-Gibbs
com fs = 8, 0 e l=65 para o evento sı́smico de 01/07/2007.
Legenda da tabela:
TO –> Tempo de origem; PDF HypLoc –> Coordenadas espaciais do hipocentro obtidas através da
determinação do valor máximo da PDF; Gauss HypLoc –> Coordenadas espaciais do hipocentro obtidas através
dos estimadores de Gauss; Ell95 –> Comprimento dos semi-eixos do elipsóide de erro (95% intervalo de
confiança)
52
Figura 31: Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
Metropolis-Gibbs, fs = 8, 0, l=65 e PDF do tipo EDT.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
53
4.3.3 Oct-Tree
De seguida, apresentam-se os resultados para o método Oct-Tree.
Os parâmetros numéricos deste método dependem essencialmente das dimensões da
grelha.
Este método mostrou ser mais rápido que o método Grid-Search e ter uma maior estabilidade numérica do que os dois anteriores. Produziu resultados para ambas as PDF’s.
Também neste método se verificou uma maior constrição da LS-L2 em relação à EDT,
mas com um maior trade-off associado.
A análise detalhada destes resultados encontra-se na secção 4.3.4.
PDF
TO (hhHmm ss,ss)
PDF HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
Gauss HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
Ell95 (km)
LS-L2
10H35 27,63
35,62
-10,39
14,96
35,62
-10,39
14,35
2,40
2,50
3,27
EDT
10H35 24,19
35,70
-10,76
29,14
35,70
-10,73
42,56
8,07
23,00
25,40
Tabela 6: Localizações hipocentrais determinadas segundo o método OCT-Tree para o
evento sı́smico de 12/02/2007.
Parâmetros numéricos: no de células iniciais(x,y,z) –> (18,12,4); comprimento mı́nimo da
célula –> 0,1 (depois de atingido este número, o processo numérico pára); no máximo de
células a avaliar –> 2680944 Legenda da tabela:
TO –> Tempo de origem; PDF HypLoc –> Coordenadas espaciais do hipocentro obtidas através da
determinação do valor máximo da PDF; Gauss HypLoc –> Coordenadas espaciais do hipocentro obtidas através
dos estimadores de Gauss; Ell95 –> Comprimento dos semi-eixos do elipsóide de erro (95% intervalo de
confiança).
54
Figura 32: Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
OCT-Tree e PDF do tipo LS-L2.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
Figura 33: Localização hipocentral para o evento sı́smico de 12/02/2007 com o método
OCT-Tree e PDF do tipo EDT.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
55
PDF
TO (hhHmm ss,ss)
PDF HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
Gauss HypLoc
(Lat(◦ ); Long(◦ ); Profundidade (km))
Ell95 (km)
LS-L2
19H03 13,89
35,96
-12,37
14,90
35,96
-12,36
14,73
0,71
1,10
3,30
EDT
19H03 13,16
35,92
-12,17
2,01
32,86
-18,02
29,30
17,50
48,50
370,00
Tabela 7: Localizações hipocentrais determinadas segundo o método OCT-Tree para o
evento sı́smico de 01/07/2007.
Parâmetros numéricos: no de células iniciais(x,y,z) –> (18,12,4); comprimento mı́nimo da
célula –> 0,1 (depois de atingido este número, o processo numérico pára); no máximo de
células a avaliar –> 2680944 Legenda da tabela:
TO –> Tempo de origem; PDF HypLoc –> Coordenadas espaciais do hipocentro obtidas através da
determinação do valor máximo da PDF; Gauss HypLoc –> Coordenadas espaciais do hipocentro obtidas através
dos estimadores de Gauss; Ell95 –> Comprimento dos semi-eixos do elipsóide de erro (95% intervalo de
confiança)
Figura 34: Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
OCT-Tree e PDF do tipo LS-L2.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
56
Figura 35: Localização hipocentral para o evento sı́smico de 01/07/2007 com o método
OCT-Tree e PDF do tipo EDT.
Legenda: Estrela verde –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a azul –> valor esperado determinado pela estatı́stica de Gauss; Elipsóide a azul –> elipsóide de erro determinado pela estatı́stica de Gauss.
4.3.4
Análise dos Resultados
A qualidade dos resultados obtidos e a avaliação do melhor método com a PDF mais
adequada foram efectuadas por comparação com os resultados determinados pelo EMSC
e o IM e analisando os vários exemplos descritos em [6], págs. 30-46.
O método que se revelou mais eficaz foi o OCT-Tree utilizando a PDF do tipo EDT.
Este resultado não foi surpreendente. O método OCT-Tree foi desenvolvido depois dos
outros dois precisamente para colmatar a rapidez computacional e simultaneamente a
estabilidade numérica.9 Investigando numerosos artigos publicados recentemente, cujos
trabalhos envolveram a utilização do NonLinLoc, constata-se rapidamente que a opção é
unânime. O OCT-Tree é o método escolhido pela maior parte dos utilizadores10
A grelha para o modelo de velocidades e, consequentemente, a grelha de localização
utilizada neste trabalho tem dimensões superiores comparativamente à maioria dos casos
abordados em outros trabalhos (para localização regional, utilizando coordenadas rectangulares). O método manteve a rapidez computacional desejada. No entanto, é necessário
atribuir bons valores aos parâmetros numéricos. Estes parâmetros numéricos dependem
essencialmente do tamanho da grelha. Por isso, uma vez efectuada a análise dos me9 O método foi integrado no NonLinLoc em 2001. Em 2000, o artigo [3] que descreve em detalhe o mecanismo
do NonLinLoc não faz referência ao método OCT-Tree. Artigos publicados antes de 2002 referem o MetropolisGibbs como o melhor método ([17],p.314). Em 2003, vem a primeira referência encontrada a OCT-Tree no artigo
[10] pág.14, sendo este o método escolhido.
10 Lista de referências bibliográficas dos trabalhos que utilizaram o software NonLinLoc com o método OCTTree: [20], pág.886, cap. Application to Yellowstone Earthquake Catalog; [21], pág.4, cap.Data analysis; [22],
pág.8, cap.3.1 Method outline; [23], pág.4; [24] pág.1.
57
lhores parâmetros numéricos para eventos sı́smicos com resultados conhecidos, não é
necessário modificar estes parâmetros numéricos em casos futuros.
Para ambos os eventos sı́smicos em análise, o OCT-Tree foi o único método que apresentou resultados para ambas as PDF’s. O valor máximo da PDF não apresentava grandes
discrepâncias nas coordenadas da latitude e longitude, para ambas as PDF’s em ambos
os eventos, mesmo alterando (suavemente) os parâmetros que determinavam a paragem
do cálculo numérico.
É possı́vel que o método OCT-Tree não identifique máximos locais em algumas situações,
mas é um método adaptado a localizar máximos absolutos em PDF’s irregulares e pouco
homogéneas (consultar secção 2.2.5), como foram os casos aqui analisados, e é isso que
essencialmente se pretende. Esta região estará sempre sujeita ao cálculo de PDF’s irregulares devido à distribuição geográfica das estações. Por outro lado, caso hajam poucos
registos para analisar (como foi também o caso do estudos sı́smicos abordados neste
trabalho), a dispersão da PDF tende a acentuar-se.
O método Grid-Search apresentou resultados consistentes com os previstos. Note-se
que é um método numérico bastante simples a nı́vel de modelação, pois só é necessário
definir o no de subdivisões (no de grelhas embutidas) que realiza e, consequentemente,
o número de vezes que determina o máximo da PDF. Optou-se no final, pelo cálculo de
três malhas embutidas por forma a não se perder acuidade e simultaneamente não tornar
método excessivamente moroso. De qualquer forma, é um método que consome bastante
tempo de cálculo relativamente aos restantes. Este método não produziu resultados para a
PDF LS-L2 relativamente ao evento de 12/02/2007. O programa abortava o processo por
não conseguir à trigésima iteração determinar os valores próprios da matriz covariante,
produzindo de seguida um overflow. A existência de muitos erros, possivelmente de medida, associados aos dados registados pelas estações constrói uma matriz covariante (ver
equação (57)) com muitos termos não-nulos, o que dificulta o processo de convergência
numérica na determinação dos valores próprios11 .
O método Metropolis-Gibbs foi desenvolvido para corrigir o problema relacionado com
o excessivo tempo de cálculo associado ao método Grid-Search. Embora bastante rápido,
é um método muito instável numericamente. Essa instabilidade verificou-se no estudo
dos dois eventos sı́smicos. Variando ligeiramente o step inicial l, as discrepâncias das
coordenadas em latitude e longitude do máximo da PDF, bem como as estimativas da
estatı́stica de Gauss, eram bastante elevadas. Este método mostra não ser eficaz para
PDF’s irregulares (consultar secção 2.2.4), sendo sensı́vel a pequenos máximos locais e
convergindo facilmente para eles. A tı́tulo de exemplo, mostram-se na secção anterior os
melhores resultados obtidos com este método. Fez-se variar o factor de escalonamento
fs (o valor recomendado é de 8.0) e obteve-se um bom resultado para fs = 9.0, mas não
foi possı́vel uniformizar os parâmetros deste método para a região em estudo. Tal como
no método Grid-Search, este método não produziu resultados para a PDF LS-L2, neste
caso, para ambos os eventos. No entanto, o problema nesta situação é de origem estritamente numérica. Após a fase de leitura, o valor da PDF determinada era de 0.00e + 00, i.e.,
nulo e inferior ao valor mı́nimo para a PDF definido num parâmetro numérico associado ao
método. A LS-L2 foi testada com este método para alguns valores de step inicial, embora
a procura não tenha sido exaustiva; é possı́vel que para determinados valores de l, este
problema pudesse desaparecer. No entanto, não era uma procura relevante porque provavelmente, entre esses valores, verificar-se-ia a mesma inconsistência numérica observada
na EDT, reforçado pelo facto de a LS-L2 ser mais sensı́vel a erros associados que a EDT
(ver secção 2.2.1, sec. Função densidade “Equal-Differential-Time” (EDT)). Porém, este
11 Para
uma melhor compreensão deste problema, consultar svdcmp numerical recipes para Linux.
58
resultado era previsı́vel. O método annealing de Metropolis [9], mostra historicamente não
ser adequado para regiões muito grandes convergindo facilmente em mı́nimos ou máximos
locais.
Visivelmente a PDF do tipo EDT mostrou ser mais estável que a LS-L2. Considerando os resultados do EMSC como os mais fidedignos, uma vez que utilizam dados
de mais de uma centena de estações, os resultados da EDT foram bastante consistentes nomeadamente em profundidade. Esta coordenada foi a que nos vários métodos
e para parâmetros numéricos variáveis (relativamente aos métodos Metropolis-Gibbs e
OCT-Tree), registou variações bastante significativas. No entanto, esta discrepância foi
bem mais suave para a PDF EDT do que para a LS-L2. Na formulação da EDT, o tı́pico
trade-off que se verifica entre a determinação do tempo de origem e da profundidade foi
minimizado. O facto de se constranger a localização a uma superfı́cie hiperbolóide deformada só para as componentes espaciais, torna a localização independente do tempo
de origem e dos resı́duos associados (consultar [15], págs.440-442 e [11], págs. 657 e
658). Por outro lado, a soma sobre as observações fora da exponencial torna a EDT mais
elevada na região de x que satisfaz o maior número de pares de observações (i.e. onde
[T obsa (x) − T obsb (x)] − [T T calca (x) − T T calcb (x)] ≈ 0), não sendo, por isso, sensı́vel
a observações com um grande erro (e.g. de medida) associado (ver Fig. 36 e Fig.37).
Por outro lado, a LS-L2 procura soluções que satisfaçam todas as observações porque o
somatório das diferenças entre tempos observados e tempos calculados está dentro do
argumento da exponencial. No entanto, como esta equação lida simultaneamente com o
cálculo das coordenadas espaciais e da coordenada temporal, ocorre um maior trade-off
entre o tempo de origem e a profundidade. Uma alteração na profundidade da fonte altera
muito pouco a variância face à alteração do resı́duo ([15], p.440). Assim, são compatı́veis
diferentes localizações em profundidade para o mesmo tempo de origem. Note-se, portanto, que a diferente natureza topológica das duas PDF’s em consideração provoca scatter samples com diferentes formas. Isto porque, o espaço de soluções da EDT é maior
que o da LS-L2: a EDT procura soluções para pares de observaçoes ao passo que a LSL2 procura soluções únicas para todas as observações. Assim, a dispersão da EDT será
maior que a da LS-L2. Devido a esta diferença topológica das PDF’s, também o elipsóide
de erro dos estimadores de Gauss irá apresentar uma maior amplitude na EDT do que na
LS-L2.
O método OCT-Tree com a PDF EDT foi muito bem modelado. Para além dos resultados bastante próximos dos resultados calculados pelo EMSC, a escolha dos parâmetros
numéricos para OCT-Tree e a eficácia da EDT foram assegurados com o resultado especial que se obteve da análise do sismo de 01/07/2007 (ver Fig. 35). Como se verifica, este
sismo ocorreu a uma profundidade muito pequena. A profundidade do oceano atlântico
é bastante variável. No entanto, o mapa de batimetria gerado com os dados batimétricos
da ETOPO (ver A.7), comprovam que a localização não foi erroneamente determinada na
água. Variando ligeiramente o parâmetro numérico referente ao tamanho mı́nimo da célula
a ser avaliada, obtinha-se valores de profundidade visivelmente errados: z = −1, 01km,
isto é, estava-se a obter localizações à superfı́cie do modelo e na água. Modificando ligeiramente este parâmetro, obteve-se para a EDT o valor em profundidade esperado, mas
uma PDF altamente irregular. Note-se que o valor calculado pela estatı́stica de Gauss é
bastante diferente. No entanto, este resultado era o esperado e inclusivamente foi estudado em [6], página 30. Visivelmente a PDF tem dois máximos: um deles é o máximo
absoluto da PDF e o outro é um máximo local da PDF e também o valor determinado pela
estatı́stica de Gauss. Isto acontece porque o hipocentro está muito perto da fronteira da
nossa região de análise (quase à superfı́cie!). Para além disso, localiza-se numa zona de
59
velocidade de elevado contraste: velocidades baixas acima e velocidades altas em baixo.
Como há mais soluções que verificam o máximo local, a estatı́stica de Gauss vai associar
esse máximo ao valor esperado. No entanto, para a PDF o valor máximo atingido é no
ponto mais à superfı́cie na zona de alto contraste de velocidades. Obtém-se assim o valor
correcto da localização.
Como nota final, é importante realçar a inclusão dos tempos de chegada das ondas
S no estudo apresentado: na ausência da sua leitura, produziam-se localizações nos
máximos locais.
Nas tabelas em anexo (ver A.6), encontra-se informação adicional acerca dos eventos
sı́smicos estudados.
60
Figura 36: Visualização 3–D do hipocentro do sismo de 12/02/2007 com o método OCTTree e PDF do tipo EDT.
Legenda: Ponto amarelo –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a ciano –> valor esperado
determinado pela estatı́stica de Gauss; Elipsóide a ciano –> elipsóide de erro determinado pela estatı́stica de
Gauss; linha verde –> resı́duo positivo da onda P ; linha amarela –> resı́duo positivo da onda S; linha magenta
–> resı́duo negativo da onda P ; linha rosa –> resı́duo negativo da onda S.
Para saber a magnitude dos resı́duos das ondas P e S, consultar tabela no anexo A.6.
61
Figura 37: Visualização 3–D do hipocentro do sismo de 01/07/2007 com o método OCTTree e PDF do tipo EDT.
Legenda: Ponto amarelo –> Máximo da PDF; nuvem vermelha –> PDF; Ponto a ciano –> valor esperado
determinado pela estatı́stica de Gauss; Elipsóide a ciano –> elipsóide de erro determinado pela estatı́stica de
Gauss; linha verde –> resı́duo positivo da onda P ; linha amarela –> resı́duo positivo da onda S; linha magenta
–> resı́duo negativo da onda P ; linha rosa –> resı́duo negativo da onda S.
Para saber a magnitude dos resı́duos das ondas P e S, consultar tabela no anexo A.6.
62
5
Considerações finais
Os resultados obtidos mostram que o método probabı́listico não-linear modelado na região
de estudo é bastante fidedigno. Os resultados revelaram que o modelo de velocidades foi
bem implementado e que o método é bastante auspicioso dentro das caracterı́sticas menos favoráveis da nossa região: modelo de velocidades de alto contraste e localização
geográfica das estações pouco homogénea. Para além disso, mostrou-se bastante eficiente no cálculo de eventos com poucos picos de fase.
Possı́veis desenvolvimentos futuros neste trabalho seriam: corrigir a pequena anomalia
do modelo de velocidades no banco de Gorringe (ver secção 4.1.2); extrapolar o método
até ao Norte de Portugal Continental, por forma, a abranger a futura estação do Porto;
testar o método OCT-Tree com a PDF EDT para mais eventos sı́smicos com mais registos
(preferencialmente, com os registos dos Açores e Madeira para ter-se uma distribuição de
estações mais homogénea) para assegurar a escolha dos parâmetros numéricos; fazer a
análise empı́rica do cálculo da magnitude.
63
Referências bibliográficas
[1] G RANDIN, R.; B ORGES, J.; B EZZEGHOUD, M.; C ALDEIRA, B; C ARRILHO, F.: Simulations of strong motion in SWIberia for the 1969 February 28 (Ms = 8.0) and the 1755
November 1 (M ∼ 8.5) earthquakes – I. Velocity Model, Geophys. J. Int. (2007), 171,
1144–1161
[2] P ODVIN, P.; L ECOMTE, I.: Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, Geophys.
J. Int. (1991), 105, 271–284
[3] L OMAX, A.; V IRIEUX, J.; VOLANT, P.; B ERGE -T HIERRY, C.: Probabilistic earthquake
location in 3–D and leyered models, Advances in Seismic Event Location Thurber,
C.H., e N. Rabinowitz (eds), 2000, Kluwer Academic Publishers, Amesterdão, 101–
134
[4] TARANTOLA, A.; VALETTE, B.: Inverse Problems = Quest for Information, Journal of
Geophysics (1982), 50, 159–170
http://www.ipgp.jussieu.fr/~tarantola/Files/Professional/Papers_PDF/IP_
QI_latex.pdf
[5] L OMAX, A.: The NonLinLoc Software Guide, http://alomax.free.fr/nlloc/
[6] L OMAX, A.; M ICHELINI, A.; C URTIS, A.: Earthquake Location, Direct, Global-Search
Methods
http://www.geos.ed.ac.uk/homes/acurtis/Lomax_etal2008.pdf
[7] S TEIN, S.; W YSESSION, M.: An Introduction to Seismology, Earthquakes, and Earth
Structure, Blackwell Publishing, 2003
[8] M OSEGAARD, K.; TARANTOLA, A.: Probabilistic Approach to Inverse Problems,
International Handbook of Earthquake E Engineering Seismology (Part A), Academic
Press, 2002, 237–265
http://www.ipgp.jussieu.fr/~tarantola/Files/Professional/Papers_PDF/
InverseProblemHandbk.pdf
[9] M OSEGAARD, K.; TARANTOLA, A.: Monte Carlo sampling of solutions to inverse problems, Journal of Geophysical Research (1995), 100, 12431–12447
http://www.gfy.ku.dk/~klaus/papers/B4-MT1995-LaTeX.pdf
[10] H USEN, S.; K ISSLING, E.; D EICHMANN, S.; W IEMER, S.; G IARDINI, D.; B AER, M.:
Probabilistic earthquake location in complex three-dimensional velocity models: Application to Switzerland, Journal of Geophysical Research (2003), 108
[11] F ONT, Y.; K AO, H.; L ALLEMMAND, S.; L IU, C.-S.; C HIAO, L.-Y.: Hypocentre determination offshore of eastern Taiwan using the Maximum Intersection method, Geophys. J.
Int. (2004), 158, 655–675
[12] M OSER, T.J.; VAN E CK, T.; N OLET, G.: Hypocenter determination in strongly heterogeneous earth models using the shortest path method., J. Geophys. Res. (1992), 97,
6563–6572
[13] W ITTLINGER, G.; H ERQUEL G.; N AKACHE, T.:Earthquake location in strongly heterogeneous media, Geophys. J. Int. (1993), 115, 759–777.
64
[14] M ILNE, J.: Earthquakes and Other Earth Movements, D. Appelton and Company,
1986, New York, 361pp. http://alomax.free.fr/nlloc/
[15] Z HOU H-w.; Rapid 3-D hypocentral determination using a master station method, J.
Geophys. (1994), Res., 99, 15439-15455
[16] T U P.; Z ISSERMAN A.; M ASON I.; C OX I.: Identification of Events from 3D Volumes of
Seismic Data, Department of Engineering Science Oxford University
http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=0041383712
[17] L OMAX, A.; Z OLLO, A.; C APUANO, P.; V IRIEUX, J.: Precise, absolute earthquake location under Somma-Vesuvius volcano using a new three-dimensional velocity model,
Geophys. J. Int. (2001) 146, 313–331
[18] L EFELDT, M.; G REVEMEYER, I.: Centroide depth and mechanism of trench-outer rise
earthquakes, J. Geophys. Res. (2008) 172, 240–251
[19] K RAFT, T.; WASSERMANN, J.; I GEL, H.: High-precision relocation and focal mechanism
of the 2002 rain-triggered earthquake swarms at Mt Hochstaufen, SE Germany, J.
Geophys. Res. (2006) 2006, 1513-1528
[20] H USEN, S.; S MITH, R.B.: Probabilistic Earthquake Relocation in Three-Dimensional
Velocity Models for the Yellowstone National Park Region, Wyoming, Bulletin of the
Seismological Society of America, Vol. 94, No. 3, pp. 880?896, June 2004
[21] D EICHMANN, N.; B AER, M.: Earthquakes in Switzerland and surrounding regions
1996-2006, Swiss Seismological Service ? ETH Zürich, October 26, 2007
http://histserver.ethz.ch/seismotectonics/reports/intro.pdf
[22] KONSTANTINOU, K.I.; L IN, C.-H.; L IANG: Seismicity characteristics of a potentially
active Quaternary volcano: The Tatun Volcano Group, northern Taiwan, Journal of
Volcanology and Geothermal Research (2006), VOLGEO-03621
http://www.earth.sinica.edu.tw/papers/LinCH/JVGR_2006.pdf
[23] S ATRIANO, C.; L OMAX, A.; Z OLLO, A.: Real-Time Evolutionary Earthquake Location
for Seismic Early Warning, Bulletin of the Seismological Society of America, Vol. 98,
No. 3, June 2008 http://alomax.free.fr/posters/rtloc/2007_IUGG_Satriano.
pdf
[24] T URINO, C.; S CAFIDI, D.: Seismotectonic Analysis of the Complex Fault System named “Garfagnana – Nord” (Northen Tuscany), Dipartimento per lo Studio del Territorio
e delle sue Risorse, Università di Genova http://www.dipteris.unige.it/eventi/
grgweb2008/abstracts/Scafidi_eng.pdf
[25] P EREZ, M.A.; B ANCROFT, J.C.: Finite-difference methods for estimating traveltimes
and raypaths in anisotropic media, CREWES,Department of Geology and Geophysics,
University of Calgary, Calgary, Alberta, T2N IN4
[26] Seismic Analysis Code,
http://www.iris.washington.edu/manuals/sac/SAC_Home_Index.html
12 Consultar
em locais com autorização, e.g. IST
65
[27] B EZZEGHOUD, M.; B ORGES, J.; C ALDEIRA, B: Actividade sı́smica em Portugal
http://www.alentejolitoral.pt/Downloads/Ambiente/Riscos%20S%C3%
ADsmicos/Actividade%20s%C3%ADsmica%20em%20Portugal.pdf
[28] M C K AY, I. C.: Statistical formulae and tables, University of Glasgow, 1997 Edition, p.
7
http://www.discourses.org.uk/statistics/formulae.pdf
[29] M ELLO, F.G. de: Probabilidades e Estatı́stica - conceitos e métodos fundamentais
vol.1, Escolar Editora, 2a edição, 2000
[30] Site do EMSC
http://www.emsc-csem.org
[31] M AZET-R OUX, G: Real Time Seismicity and Alert Systems: the FAQ for improved
services
http://www.emsc-csem.org/Doc/EMSC_RTS-EWS.pdf
[32] M AZET-R OUX, G.; B OSSU, R.; C ARRE ÑO, E.; G UILBERT, J.; E L -B AKI, J.-M.: Report
on EMSC – Real Time Earthquake Information services
http://www.emsc-csem.org/Doc/EMSC_DOCS/EMSC_RTEI_2005_low.pdf
[33] Boletim Sismológico e Preliminar do Continente e da Madeira para Fevereiro de 2007,
IM
[34] Boletim Sismológico e Preliminar do Continente e da Madeira para Julho de 2007, IM
[35] C ARVALHO, J.; R ABETH, T.; C ABRAL, J.; C ARRILHO, F.; M IRANDA, J.M.: Geophysical
characterization of the Ota-Vila Franca de Xira–Lisbon–Sesimbra fault zone, Portugal,
Geophys. J. Int. (2008), 174, 567–584
[36] L IENERT B.R.: Hypocenter 3.2 – A Computer Program for Locating Earthquakes Locally, Regionally and Globally (October, 1994), Hawaii Institute of Geophysics and
Planetology
[37] Optimization – The samped least squares method
http://www.sinopt.com/learning1/optsoft/optimization/optimization.htm
[38] WALDHAUSER & E LLSWORTH hypoDD double-difference earthquake location algorithm
http://www.ldeo.columbia.edu/~felixw/DD.html
[39] K LEIN, Fred W.: User’s Guide to HYPOINVERSE-2000, a Fortran Program to Solve
for Earthquake Locations and Magnitudes, U. S. Geological Survey, 4/2002 version
http://geopubs.wr.usgs.gov/open-file/of02-171/of02-171.pdf
[40] Dados batimétricos da ETOPO
http://www.ngdc.noaa.gov/mgg/global/
66
A
A.1
Anexo
Secções XZ do modelo de velocidades
A.1.1
Intervalo [32◦ ; 34◦ [ N: zona extrapolada
A.1.2
Intervalo [34◦ ; 40◦ ] N
Nota: Para cada par de gráficos, o da esquerda refere-se a uma parte da secção XZ do
modelo implementado no NonLinLoc; o da direita, representa as curvas de nı́vel extraı́das
do modelo SWIberia.
Mapa de curvas de nível de Vp no plano XZ para 34° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
67
Mapa de curvas de nível de Vp no plano XZ para 34,5° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 35,1° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 35,7° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 36,02° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
68
Mapa de curvas de nível de Vp no plano XZ para 36,8° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
69
Mapa de curvas de nível de Vp no plano XZ para 37° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 37,5° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 38° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 38,4° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
70
Mapa de curvas de nível de Vp no plano XZ para 38,6° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 38,8° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para 39,3° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
Mapa de curvas de nível de Vp no plano XZ para (aprox.) 39,65° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
71
Mapa de curvas de nível de Vp no plano XZ para 40° de latitude
8
6
4
2
Profundidade (km)
9
8
7
6
5
4
3
2
1
0
10
20
30
40
50
60
−13 −12 −11 −10 −9 −8 −7 −6 −5
Longitude (°)
A.1.3
Intervalo [40◦ ; 40, 5◦ [ N: zona extrapolada
72
A.2
A.2.1
Tempos de propagação e ângulos take-off teóricos para as estações
sismológicas
Secções XY e XZ dos tempos de propagação
73
74
75
A.2.2
Secções XZ dos ângulos take-off
76
77
A.3
Código C do programa de visualização do modelo G RANDIN et.al.
para a plataforma gráfica Gnuplot
Nota: O código, que se segue, não contém as macros com a sintaxe própria do Gnuplot
para construção de gráficos. Essas macros variam conforme as necessidades do utilizador. Este código apenas transforma os ficheiros binários13 com rotina de aplicação em C,
em ficheiros binários que possam ser lidos pelo programa Gnuplot.
/***Programa de leitura do modelo de Velocidades SWIberia-Grandin
para os sismos de 1755/Novembro, 1969/Fevereiro/28.****/
/*Nota: Para alterar o ficheiros a serem lidos,
basta faze-lo na rotina de abertura do ficheiro binario.
Os nomes dos ficheiros de output tambem podem ser alterados
nas rotinas que geram os mesmos.*/
#include <stdio.h>
#include <stdlib.h>
/*Declaracao da funcao auxiliar writematrix*/
void writematrix(char *filename,
float *x,
float *y,
float **m,
int nx,
int ny);
int main(){
/*Declaracao e inicializacao das variaveis*/
int i, j, k;
float ***p;
int xsize=881;
int ysize=661;
int zsize=61;
FILE *fon;
//eixo ns = NORTE-SUL corresponde a yy
float ns[ysize];
//eixo oe = OESTE-ESTE corresponde a xx
float oe[xsize];
//eixo pr = profundidade corresponde a zz
float pr[zsize];
13 concedido
por J. B ORGES
78
//ponteiro para P com menos uma dimensao para plotar
float **V;
/********************Leitura do modelo SWIberia-Grandin*******************/
/*Alocacao do tensor "p" que guarda os valores da grelha contidos
no ficheiro binario.*/
p=(float ***)malloc(ysize*sizeof(float *));
for(i=0;i<ysize;i++)
{
p[i]=(float **)malloc(zsize*sizeof(float *));
for(j=0;j<zsize;j++)
{
p[i][j]=(float *)malloc(xsize*sizeof(float));
}
}
/*Abre o ficheiro binario do modelo SW: os ficheiros podem ser referentes
a Vp,Vs,pho e Q com a extensao da respectiva variavel.*/
fon=fopen("SWIB_2006.p","rb");
if (fon==NULL)
{
printf("O ficheiro de leitura do modelo n~
ao foi colocado
na directoria corrente.\n");
return -1;
}
/*Le o ficheiro de acordo com a rotina do modelo SW-Grandin.*/
for(i = 0; i < ysize; i++)
{
for(j = 0; j < zsize; j++)
{
for(k = 0; k < xsize; k++)
{
fread(&(p[i][j][k]), 4, 1, fon);
}
}
}
/*Fecha o ficheiro binario dos dados do modelo.*/
fclose(fon);
/*************************Leitura concluida******************************/
79
/*--------Criacao dos vectores para os eixos dos graficos--------*/
/*Criar o "ns" (Latitude), vector do eixo Y*/
for(i=0;i<ysize;i++)
{
ns[i]=34+6.0*((float)i)/((float)ysize-1.0);
}
/*Criar o "oe" (Longitude), vector do eixo X*/
for(k=0;k<xsize;k++)
{
oe[k]=-13.0+8.0*((float)k)/((float)xsize-1.0);
}
/*Criar o "pr" (km), vector do eixo Z*/
for(j=0;j<zsize;j++)
{
pr[j]=60.0*((float)j)/((float)zsize-1.0);
}
/*---------------Criacao dos vectores concluida-------------------*/
/*As seccoes para "plotar" a 2D e 3D, que se seguem,
podem funcionar alternadamente: dao-nos o output pretendido.*/
/*Para alterar o plano a "plotar" no caso 2D,
tem de modificar adequadamente o ciclo for.*/
/*Para alterar os planos a "plotar" no caso 3D,
tem de se modificar adequadamente a rotina de alocacao de V,
associacao de V a P, rotina writematrix e rotina de desalocaç~
ao do V.*/
/********************"Plotar" graficos 2D**************************/
/*Rotina que guarda os valores de p so para uma dimensao.*/
/*O ficheiro gerado esta pronto para ser "plotado" pelo Gnuplot.*/
//fon=fopen("Vp2D_y=0.dat","w");
/*O ficheiro contém 2 colunas:
1a para o eixo dos xx e 2a para o eixo dos yy.*/
/*for(j=0;j<zsize;j++)
{
fprintf(fon,"%g %g\n",pr[j],p[0][j][440]);
}
*/
/* fclose(fon);*/
/*Com k=440, estamos a fazer o perfil do plano XZ no centro deste plano.*/
/*******************Rotina 2D concluida****************************/
80
/********************"Plotar" graficos 3D**************************/
/*Alocacao da matriz das velocidades "V",
que vai corresponder a uma so layer da matriz "p".*/
V=(float **)malloc(xsize*sizeof(float *));
for(k=0;k<xsize;k++)
{
V[k]=(float *)malloc(zsize*sizeof(float));
}
/*Associacao da matriz V aos valores da matriz "p"*/
for(k=0;k<xsize;k++)
{
for(j=0;j<zsize;j++)
{
V[k][j]=p[449][j][k];
}
}
/*Funcao para escrever os dados em formato legivel pelo gnuplot
para serem plotados a 3D.*/
/*Os parametros desta funcao incluem V,
o vector que contem os dados do nosso modelo
e o ficheiro que sera criado para ser lido pelo gnuplot.*/
writematrix("Vp3Dy=449.dat",oe,pr,V,xsize,zsize);
/*Libertar a matriz "V"*/
for(k=0;k<xsize;k++)
{
free(V[k]);
}
free(V);
/*******************Rotina 3D concluida************************/
/*Libertar de memoria a matriz "p"*/
for(i=0;i<ysize;i++)
{
for(j=0;j<zsize;j++)
{
free(p[i][j]);
}
free(p[i]);
}
81
free(p);
return 0;
}
/***********Funç~
ao auxiliar
writematrix**************/
#include <stdio.h>
#include <stdlib.h>
void writematrix(char *filename,
float *x,
float *y,
float **m,
int nx,
int ny)
{
FILE *fp;
int i, r;
float **matriz;
fp=fopen(filename,"wb");
matriz=(float **)malloc((ny+1)*sizeof(float *));
for(i=0;i<ny+1;i++)
matriz[i]=(float *)malloc((nx+1)*sizeof(float));
matriz[0][0]=(float)nx;
for(r=0;r<nx;r++)
matriz[0][r+1]=x[r];
for(i=0;i<ny;i++)
matriz[i+1][0]=y[i];
for(r=0;r<nx;r++)
for(i=0;i<ny;i++)
matriz[i+1][r+1]=m[r][i];
for(i=0;i<ny+1;i++)
fwrite(matriz[i],sizeof(float),nx+1,fp);
for(i=0;i<ny+1;i++)
free(matriz[i]);
free(matriz);
fclose(fp);
}
82
A.4
Sintaxe do script para o modelo de velocidades
VERTEX id_num zloc xloc yloc
EDGE id_num vertex1 vertex2
POLYGON3 id_num n_edges
[NEW_LINE] edge1 edge2 ...
SOLID id_num n_edges depth Vp Vp_grad Vs Vs_grad p_top p_grad
[NEW_LINE] polygon1 polygon2 ...
Legenda: id num - no de identificação do vértice, linha, polı́gono ou sólido
xloc yloc zloc - identificação (x, y, z) em km da localização do vértice
no referencial com origem definida no ficheiro de controlo
n edges - no de linhas que definem o polı́gono
depth - profundidade de referência do sólido
Vp - velocidade da onda P na profundidade de referência do sólido
Vp grad - gradiente da velocidade da onda P no sólido
Vp - velocidade da onda S na profundidade de referência do sólido
Vp grad - gradiente da velocidade da onda S no sólido
p - densidade à profundidade de referência no sólido
p grad - gradiente da densidade no sólido
Nota: Foi considerada como profundidade de referência a cota mais profunda do sólido
83
A.5
Macro para o SAC das leituras dos picos de fase (exemplo)
* Ler o ficheiro
r ../raw/2007.043.10.34.20.4977.ES.EOSO..BHZ.R.SAC
* Remove a media
rmean
* Remove a inclinacao
rtrend
* Desconvolucao
TRANSFER FROM POLEZERO SUBTYPE
/home/lapsis/wfmod/PZ/SAC_PZs_ES_EOSO_BHZ__2000.001.00.00.00.0000
TO NONE FREQ 0.033 0.05 20 30 PREW 2
* Quick and dirty plot
QDP off
* Abre o x
bd x
* Gera grafico
plot
* Grava
w 2007.043.10.34.20.4977.ES.EOSO..BHZ.R.decon.SAC
84
A.6
Informação adicional relativa à localização hipocentral dos sismos de 12/02/2007 e 01/07/2007
Legenda das tabelas:
sta: código da estação
inst: banda de frequência recebida pela estação
cmp: componente do sinal analisada
on: descrição do tipo de chegada da onda (e- emergente; i - impulsiva)
phs: tipo de onda (P ou S)
fm: first motion (C ou U - compressão; D - dilatação)
date: data do evento sı́smico
time: hora, minuto, segundo da recepção da chegada da onda P
error: tipo de erro (GAU - gaussiano)
amp: amplitude
per: perı́odo
pred tt: tempo de chegada previsto
resid: resı́duo da onda P (Tobs -Tprevisto )
weight: peso da onda
sta:
Dist - distância epicentral da estação ao hipocentro (em km)
Azim - ângulo azimutal do epicentro à estação (em graus, no sentido horário a partir do Norte)
ray:
RAz - ângulo take-off azimutal do hipocentro estimado (em graus, no sentido horário a partir do Norte)
RDip - ângulo take-off de profundidade do hipocentro (em graus, na vertical: 0 = para baixo; 180 = para cima)
q - factor de qualidade do ângulo take-off estimado (0 = mau; 10 = excelente)
85
A.7
Batimetria SWIberia + região Atlântica dos Açores e da Madeira
86
Download
Tese 9,7 MB - Técnico Lisboa

ONG AMIGOS DOS BICHOS FEIRA DE ADOÇÃO DE CÃES E

IMPULSO

GESTÃO DE VELOCIDADE NO BRASIL E INICIATIVAS

04 - Física na Veia!

VT COLHEITADEIRA BC6500.indd

Magnetismo e Magnetos Moleculares

Aula1deExercícios

Motores Entre as propostas mecánicas para a gama
