Problem a in Verso 2

7/21/2019 Problem a in Verso 2

1/26

Prospeco Geofsica I Captulo 8. Inverso de dados emPGEMFernando M. Santos-2006

________________________________________________________________________________________________

1

CAPTULO 8. INVERSO DE DADOS EM PROSPECOGEOELECTROMAGNTICA

8.1. Introduo.8.2. Formulao do problema inverso.

8.3. A informao contida nos dados.

8.4. Resoluo do problema inverso linear.8.4.1. Constrangimentos8.4.2. Covarincia dos parmetros8.4.3. Matrizes da resoluo dos dados e dos parmetros

8.5. Resoluo do problema inverso no linear usando mtodos de minimizaolocal.8.5.1.O mtodo clssico de soluo do PI sobre-determinado (a inverso comoprocesso de optimizao).8.5.2. O PI amortecido (misto). A informao a priori.8.5.3. Os mtodos de regularizao (mtodo de Tikhonov).

8.5.4. O mtodo do gradiente descendente.8.5.5. O mtodo de Newton.8.5.6. O mtodo do gradiente conjugado.

8.6. Resoluo do problema inverso no linear usando mtodos de minimizaoglobal.8.6.1. Simulated annealing (SA).8.6.2. Controlled Random Search(CRS) -Procura aleatria controlada.

8.7. Anlise da soluo do PI.8.8. Ambiguidade (equivalncia) e sensibilidade da soluo.

8.9. Inverso conjunta.

8.9.1. Inverso conjunta de dados de resistividade e MT.8.9.2. Inverso conjunta de dados de GEM e ssmica.

8.10. Concluses___________________________________________________________________

8.1. Introduo

A interpretao dos dados adquiridos numa campanha geofsica constitu um dosmais importantes estdios do trabalho de geofsico. Por interpretao dos dadosentende-se aqui, a complexa (e difcil) tarefa de estabelecer um modelo geofsicoque cumpra duas condies fundamentais: 1) que satisfaa o conjunto de dadosrecolhidos, de acordo com um critrio pr-estabelecido, e 2) seja susceptvel de ter

interpretao geolgica que esteja de acordo com a informao conhecida. Nestecontexto o termo modelodesigna o conjunto de valores das propriedades fsicas dasrochas (designados por parmetros do modelo) que caracterizam o(s) fenmeno(s)fsico em que se baseia a metodologia de prospeco utilizada. Veremos, numpargrafo posterior, que o termo modelo pode ter, no mbito da resoluo doproblema inverso, outro significado.

O estabelecimento do modelo geofsico uma operao complexa que, nos dias quecorrem e cada vez com maior frequncia, faz apelo resoluo do problemainverso, na medida em que se parte dos dados para se obterem estimativas doconjunto de parmetros do modelo, bem como informao relativa aos seus limitesde confiana. O problema inverso (PI) , assim, definido em oposio ao problemadirecto. Neste, conhecem-se os parmetros e calcula(m)-se a(s) resposta(s) do

modelo. O problema directo um problema bem postono sentido em que, dado um


2/26


________________________________________________________________________________________________

2

conjunto de valores para os parmetros e conhecidas as expresses que regem ocomportamento do campo, a soluo nica e pode ser calculada, em princpio,com a preciso que se desejar. O problema inverso , pelo contrrio, um problemamal postopois diferentes modelos podem representar igualmente bem os dados (naprtica o problema inverso no tem soluo nica) e a soluo geralmente

instvel. Estes, e outros problemas, sero abordados em pargrafos futuros.

O problema inverso e a sua resoluo constitui, de facto, um ramo da matemticaaplicada sendo as tcnicas a desenvolvidas usadas na soluo de muitos problemasencontrados na actividade humana (em diferentes ramos de engenharia, em estudoseconmicos, sociais, bio-mdicos, etc.). O presente captulo analisa alguns dasquestes relacionadas com a resoluo do problema inverso em prospeco geo-electromagntica (PGEM), englobando nesta designao, todos os mtodos quefazem uso do campo electromagntico, com especial relevo para os mtodos daresistividade e magneto-telrico. A importncia destes mtodos na geofsica geral ena geofsica aplicada por demais conhecida. Os mtodos GEM so aplicados como objectivo de conhecer a distribuio da resistividade elctrica no subsolo aprofundidades que podem variar de algumas dcimas de metro a dezenas de

quilmetros. Deste modo a sua aplicao abrange reas to diversas como as doreconhecimento geolgico, hidrogeolgico, geotrmico e mineiro. Nos ltimos anostem crescido o interesse na aplicao dos mtodos GEM a problemas ambientais,nomeadamente no reconhecimento e monitorizao de zonas contaminadas.Todavia, os mtodos GEM no se limitam a contribuir para a soluo de algunsproblemas prticos de prospeco: eles permitem, tambm, adquirir informaorespeitante ao interior profundo do nosso planeta.

A importncia dos mtodos GEM facilmente compreendida se se recordar que aresistividade (condutividade) elctrica das formaes rochosas , talvez, apropriedade fsica que maior variaes sofre com as variaes do contedo equalidade da gua existente nos poros e fracturas das rochas, da temperatura e coma presena de alguns materiais geolgicos como sejam as argilas, os metais e a

grafite. Entende-se, ento, que o conhecimento da distribuio da condutividadeelctrica (ou da resistividade elctrica) pode ser um precioso auxiliar na observaoindirecta do interior do nosso planeta.

8.2. Formulao do problema inverso

Seja ento do o vector contendo o conjunto dos valores observados (valoresmedidos ou calculados a partir de valores medidos de expresses algbricas). Nocaso especfico dos mtodos GEM esses valores so, geralmente, os valores deresistividade aparente e/ou fases calculados para diferentes configuraesgeomtricas do conjunto emissor-receptor e/ou para diferentes frequncias. Cadaelemento de dopode ser descrito como o valor aproximado do verdadeirovalor dagrandeza a medir d, ao qual no se tem acesso, possuindo-se uma estimativa do erro

da medida, tal que

d(do, d ) < d com a condio d

od quando d0

onde representa uma distncia (mtrica) definida no espao dos dados.

Seja mest o vector dos parmetros do modelo adoptado para a interpretao dosdados. No caso dos mtodos GEM o vector mest corresponde a valores deresistividade (ou condutividade) elctrica e, em alguns modelos, a profundidades (localizao) de interfaces que separam meios com propriedades geoelctricasdistintas. No problema inverso o vector mest o vector das incgnitas e a resoluodo problema inverso dever permitir conhecer os elementos deste vector. De acordocom os mtodos adoptados nesta lio, para a resoluo do problema inverso, os


3/26


________________________________________________________________________________________________

3

valores de mestrepresentam uma estimativa dos verdadeirosvalores dos parmetrosm, isto ,

m(mest, m) < m com a condio m

estm quando m0

m representa uma distncia definida no espao dos parmetros do modelo. Nestecontexto o termo modelo usado para designar as relaes matemticas que ligamos parmetros m resposta do modelo d,

d= g (m) .

De acordo com a natureza daquelas relaes usual classificar os problemasinversos em lineares e no lineares. O problema inverso em GEM, atendendo natureza das funes g(m) , claramente, no linear.

O problema inverso pode ainda ser classificado de acordo com a dimensionalidadedo modelo adoptado, podendo ser, ento, designado por problema inverso uni-dimensional (1-D), bi-dimensional (2-D) ou tri-dimensional (3-D). O termodimensionalidade traduz, ento, o modo como a resistividade elctrica variaespacialmente, isto , traduz a dependncia da resistividade elctrica com a posio(x, y, z). Nos modelos 1-D admite-se que a resistividade funo apenas daprofundidadez. Nos modelos 2-D admite-se que a resistividade varia numa secovertical, enquanto nos modelos 3-D se admite que a resistividade funo das trsvariveis espaciais.

Resumindo, o problema inverso consiste na determinao de mest cumprindo duascondies necessrias: 1) os valores da resposta do modelo dc= g (mest), devem sercompatveis com os dados do, isto , d(d

o, dc) < e, 2) deve ser possvel estimaros erros do modelo calculado, isto , deve ser possvel determinar os limites devalidade do modelo calculado (geralmente designados por intervalos de confiana).

Como se fez notar na introduo, as duas condies enumeradas so necessriasmas no so suficientes, pelo menos em geofsica. Uma terceira condio dever serimposta, a saber: o modelo calculado deve ser interpretvel em termos geolgicos,devendo estar de acordo com a informao disponvel. Esta ltima condio, quederiva do facto da soluo do PI no ser nica e de que se trata de dados recolhidosem ambiente geolgico, deve estar sempre presente quando se interpretam dadosgeofsicos.

8.3. A informao contida nos dados

A interpretao de quaisquer dados deve ser feita de acordo com a informao nelescontida. Isto significa que se devem evitar duas situaes: 1) a sobre-parametrizaodos modelos e, 2) o sobre-ajustamentoentre os dados e a respostado modelo. A sobre-parametrizaodo modelo significa que se est a tentar obterinformao sobre estruturas que no esto representadas nos dados. O sobre-ajustamentoentre dados e resposta do modelo significa que o modelo est a incluirestruturas que no correspondem aos dados mas aos erros. Estes problemas soextremamente complexos e esto directamente relacionados com a definio deprofundidade de investigao(ou de volume de investigao) e com o problema daanlise de sensibilidadedo modelo.

De acordo com a teoria geral do comportamento do campo electromagntico emmeios condutores, sabe-se que nos mtodos GEM a profundidade de investigao1

1 A profundidade de investigao ou profundidade efectiva a profundidade para a qual, em meio homogneo e

isotrpico, o mdulo de uma das componentes do campo electromagntico tem um valor igual a 1/e do seu valor superfcie do meio.


4/26


________________________________________________________________________________________________

4

depende da resistividade do meio e da frequncia do campo electromagntico (ouda geometria do dispositivo usado, no caso do mtodo da resistividade).Suponhamos que a investigao se realiza usando sondagens magneto-telricas(MT) ou sondagens Schlumberger. O volume de material investigado depender,portanto, daqueles parmetros e ainda do nmero e disposio espacial das

sondagens. O conceito de volume de material investigadoou de volume efectivoum conceito importante quando se lida com problemas 2-D e 3-D e no temmerecido grande ateno por parte dos investigadores. O tema voltar a serabordado no pargrafo 7.

Consideremos que os valores a observar (da resistividade aparente, por exemplo),com o uso de mtodos GEM, variam de modo contnuo no espao e em frequncia.Se fosse possvel ter o conhecimento completo do modo como se comportam essesobservveis, seria possvel, em princpio, estimar a distribuio contnua daresistividade elctrica no sub-solo de modo nico (Zhdanov and Keller, 1994). Aaquisio de dados em GEM feita, contudo, de modo discreto, quer espacialmentequer em frequncia. Isto , os elementos do vector do so conhecidos apenas emalguns locais e apenas para um nmero finito de frequncias, do(xi, yi, zi, j), comi=1,., N ej=1,., P. Este facto limita o conhecimento que se pode obter sobre adistribuio dos valores de mest.Os elementos de mestso, geralmente, conhecidosde modo discreto mest(xk,yk,zk), com l=1,.M. Este processo conhecido como aparametrizao do modelo. Tem-se, ento, que o problema inverso, tal comoformulado acima, um caso de problema inverso discreto.

Considere-se um conjunto de N sondagens realizadas para P frequncias em pontoscolineares e igualmente espaados de uma distncia x. De acordo com o teoremada amostragem

2 no se possuir informao sobre a distribuio horizontal daresistividade elctrica que tenha uma dimenso lateral caracterstica inferior a 2x.Nos mtodos GEM o problema de avaliar a dimenso caracterstica pode ser,contudo, mais complexo pois, devido ao facto da profundidade efectiva dependertambm da frequncia do campo electromagntico a aplicao simples do teoremada amostragem pode conduzir a resultados errneos. Por outro lado, como ocomportamento do campo electromagntico obedece a uma lei de difuso (para agama de frequncias utilizadas em GEM), verifica-se uma perda de resoluo com oaumento da profundidade a investigar. Isto deve-se ao facto d a informao obtidareflectir a influncia de um volume de massa rochosa tanto maior quanto menor fora frequncia do campo electromagntico usado (ou quanto maior a distncia entreo emissor e o receptor). O estudo do volume de material que influncia asobservaes, bem como o estudo da importncia relativa dessa influncia geralmente designado por anlise da sensibilidade do modelo e ser abordada nopargrafo 7.

Deve notar-se, no entanto, que a limitao mencionada no ltimo pargrafo podeser usada de modo a rentabilizar as sondagens GEM realizadas. De facto, se osdados adquiridos contm informao referente a um grande volume de materialpoder-se-, ento, aumentar a amostragem espacial (at um certo limite) sem perdasignificativa da informao referente distribuio dos parmetros no subsolo.

8.4. Resoluo do problema inverso linear

Considere-se a resoluo do PI linear em que a relao entre dados e parmetros domodelo seja explcita, isto ,

d =G m (8.4.1)

2

De acordo com este teorema, a frequncia mxima (de Nyquist) que se pode obter com uma amostragem T fN = 1/(2T).


5/26


________________________________________________________________________________________________

5

O mtodo a usar para a resoluo deste PI depende da relao entre N e M. Trscasos podem ocorrer: 1) se N = M o problema completamente determinado, 2) seN > M o problema ser um problema sobredeterminadoe, 3) se N


6/26


________________________________________________________________________________________________

6

componentes sobredeterminada e subdeterminada. Se o valor do parmetro frmuito grande a soluo minimiza a parte subdeterminada da soluo. Se o valor doparmetro fr muito pequeno o erro ser minimizado mas no ser possvelintroduzir informao sobre a parte subdeterminada pelo que ser possvelencontrar vrias solues para o problema. Ter de se procurar, em cada caso, um

valor ptimo para o parmetro . A soluo do problema, minimizando-se (8.4.8),ser,

mest= (GTG + 2I)-1GTdo (8.4.9)

8.4.1. Constrangimentos

H outros tipos de constrangimentos que podem ser usados. Pode-se, por exemplo,pretender que a soluo seja prxima de um valor mdio . Neste caso vir paraL,

L = (m-)T (m-). (8.4.8)

Muitas vezes, contudo, este tipo de constrangimento pode no ser apropriado. Umapossvel alternativa o de considerar-se a soluo comoplana, isto ,

L = [Dm]T[Dm] = mTDTDm (8.4.9)

onde a matriz D dada por:

D=

11

...

...

0....110

0.....011

ou

D=

121

...

...

...1210

...0121

se se pretender que a soluo seja suave. Por vezes conhecem-se algumas relaesque a soluo cumpre, por exemplo, pode conhecer-se o valor mdio da soluo, h 1,isto ,

Fm = h1 (8.4.10)

com F= [1 1 1 ....1] . Se se conhecer o valor de um dos elementos da soluo,mi= h1tem-se, F= [0 0 1 0 ....0].

Pode, ento, definir-se o PI da seguinte forma: calcular mest tal que G mest = do,cumprindo ainda a condio Fm = h.


7/26


________________________________________________________________________________________________

7

Os erros dos dados podem ser usados como constrangimento uma vez que pesama importncia dos dados no processo de inverso. Considere-se que os dados soindependentes obdecendo a uma estatstica Gaussiana e que a covarincia dos dados da forma,

[cov d] = d2

I (8.4.11)

Definindo-se a matriz diagonal W= diag(1/1.... 1/N) o uso dos erros dos dados feito atravs da matriz Gw = W G e do vector dw= W d, vindo a soluo, para oproblema sobredetermindo,

mest= (GwTGw)

-1GwTdw (8.4.12)

8.4.2. Covarincia dos parmetros

Os dados tm erros e esses erros reflectem-se nos valores estimados para osparmetros. Em todos os casos estudados a soluo do tipo,

mest= Md+ v (8.4.11)

A matriz M por vezes designada por inversa generalizada. Seja a covarincia dosdados,

C= [cov d] (8.4.12)

Os valores estimados para os parmetros tero uma covarincia dada por,

[cov m] = M[cov d] MT (8.4.13)

Se os dados forem independentes (no correlacionados) e com varincia d2, ento3,

[cov d] = d2I (8.4.14)

[cov m] = Md2IMT

No caso do problema sobredeterminado a inversa generalizada ,

M= (GwTGw)

-1GwT (8.4.15)

e tem-se para a covarincia do modelo:

[cov m] = [ (GwTGw)

-1GwT] [ (Gw

TGw)-1Gw

T]T= (GwTGw)

-1 (8.4.16)

ou

[cov m] = d2(GTG)-1 (8.4.17)

Conhecida a covarincia do modelo podem calcular-se os erros associados a cadaparmetro (com 95% de confiana) dado por,

m= 1.96 diag([cov m])1/2 (8.4.18)

3Os dados pesados tm uma covarincia [cov d] = I


8/26


________________________________________________________________________________________________

8

8.4.3. Matrizes da resoluo dos dados e dos parmetros

Quando se resolve o PI pretende-se que a resposta do modelo obtido seja prximados dados observados. Esta condio verifica-se quando a matriz resoluo dosdadosN= G M =I pois,

dc= G mest= G ( M do) = G M do (8.4.15)

Por outro lado, pretende-se que o modelo estimado seja prximo do modeloverdadeiro, isto , que os valores estimados dos parmetros sejam prximos dosverdadeiros,

mverd= Mdo= M(Gmest) = M G mest (8.4.16)

Aquela condio s se verifica quando a matriz da resoluo dos parmetrosR= M G= I.

As matrizes Ne Rpermitem avaliar o resultado obtido na resoluo de um PI. Esta

avaliao, pode ainda ser feita com recurso a um outro parmetro designado porespalhamento(sp) das matrizes Ne Re definidos da seguinte forma:

sp (N) = || N I||2 (8.4.17)

e

sp (R) = || R I||2 (8.4.17)

8.5. Resoluo do problema inverso no linear usando mtodos de minimizaolocal

8.5.1.O mtodo clssico de soluo do PI sobre-determinado (a inverso comoprocesso de optimizao)

Retome-se a formulao do problema inverso em GEM:

conhecido o vector de dados do(d1,,dN), deve determinar-se o vector deparmetros mest(m1,,mM) que, de acordo com o modelo d

c= g (mest), cumpre acondio d(d

o, dc)< .

Atendendo no linearidade do problema inverso em questo, o processo desoluo dever ser um processo iterativo, isto , aceita-se que a resposta do modelo,para pequenas variaes nos parmetros, pode ser desenvolvida em srie de Taylor4,

g(mk) = g(mk-1) + J(mk-1) m+ R(mk-1, m)(8.5.1.1)

onde mk representa o vector de parmetros para a iterao k, m o vector com ascorreces aos parmetros obtidos na iterao k-1 e J a matriz de sensibilidadecom elementos dados por,

Jij= gi/ mj.(8.5.1.2)

4

Com o objectivo de facilitar a escrita o vector dos dados e dos parmetros sero representados, no texto que se segue,por de m, respectivamente.


9/26


________________________________________________________________________________________________

9

As correces aos parmetros podem ser calculadas, no caso de N>M (problemasobre-determinado), usando a tcnica dos mnimos quadrados. Neste caso a funoobjectivo ser, por exemplo,

= d= || d- g(m) ||2.

(8.5.1.3)

Desprezando os termos de ordem superior R(mk-1, m) no desenvolvimento daresposta do modelo em srie de Taylor, as correces aos parmetros do modelo,para a iterao kpodem ser obtidas pela resoluo do sistema de equaes,

(J(mk)TWd-1J(mk)) m= J(mk)TWd

-1(d - g(mk) )(8.5.1.4)

onde a matriz covarincia dos dados dada por, Wd = d2 I. d

2 representa avarincia dos erros dos dados (supostamente no correlacionados) e I a matrizidentidade5.

Estimadas as correces aos parmetros do modelo para a iterao k, os valores dempara a iterao k+1sero dados por,

mk+1= mk+ m(8.5.1.5)

O sistema de equaes (8.5.1.4) pode ser re-escrito da seguinte forma,

m= S+(d - g(mk) )(8.5.1.6)

A matriz S+ a matriz inversa generalizada. Se se considerar a decomposio em

valores singulares da matriz J, tal que J= UV

T

, onde U a matriz ortonormaldos vectores prprios da matriz JJT, Va matriz ortonormal dos vectores prprios damatriz JTJ e a matriz diagonal contendo os valores singulares (raz quadrada dosvalores prprios da matriz JTJ), a inversa generalizada vem,

S+= V-1UT(8.5.1.7)

e

m= V-1UT(d - g(mk) )(8.5.1.8)

Uma anlise cuidada desta expresso mostra que os valores singulares de menormagnitude podero originar instabilidade na soluo do sistema de equaes(8.5.1.8). A estabilizao da soluo pode ser conseguida, atravs da aplicao domtodo de Levenberg-Marquardt, isto , minimizando o efeito dos valoressingulares menos significativos,

m= V UT(d - g(mk) ) ( 2+ 2I)-1(8.5.1.9)

5Aceita-se que os dados tm uma distribuio Gaussiana.


10/26


________________________________________________________________________________________________

10

O parmetro emprico, embora haja alguns procedimentos que geralmente soseguidos quando se programa a expresso (8.5.1.9).

A soluo do problema inverso, tal como acaba de ser descrita pode ser consideradacomo clssica (soluo aos mnimos quadrados), no sentido de que a mais

utilizada em GEM, quando N>M (problema sobre-determinado). A metodologiadescrita , por exemplo, frequentemente empregue na resoluo do problemainverso quando se adoptam modelos estratificados (Johansen, 1997).

Exerccio. Considere o conjunto de valores mostrados na tabela seguinte. Calcule arecta (w =a+be) que melhor se aproxima aos pontos, no sentido dos mnimosquadrados.

resistividadedo terreno (e)

resistividadeda gua (w)

421 7.35250 4.69359 6.56

351 4.90285 4.9547 1.5021 0.5347 0.7565 1.8623 1.36

69.64 2.80208 4.27

214.5 5.19

Exerccio. Elabore um algoritmo para a resoluo do PI no caso de sondagenselctricas verticais, assumindo um modelo de terreno estratificado. Sugesto: utilize

a tcnica de SVD e, para o problema directo, utilize o conceito de transformada daresistividade (Captulo 4).

8.5.2. O PI amortecido (misto). A informaoa priori

A maioria dos problemas inversos em geofsica no so nem (completamente) sub-determinados nem sobre-determinados, mas mistos, no sentido em que haverparmetros que sero bem (completamente) determinados e parmetros que sero(completamente) sub-determinados. Os parmetros do modelo podem, ento,dividirem-se em dois conjuntos: os sub- e os sobre-determinados. A soluo do PIpoderia ser obtida determinando-se os parmetros sobre-determinados de acordocom o mtodo dos mnimos quadrados e estimando os parmetros sub-determinados, impondo que estes cumpram uma condio definida a priori. A

forma mais simples de proceder, contudo, a de calcular uma soluo que minimizeuma combinao de de da condio a priori, S (Menke, 1984).

Neste contexto o vector das correces dos parmetros no caso do PI no linearabordado em 4.1., poder ser estimado minimizando-se a funo:

= d+ 2 S = || g(m) d||2+ 2 S.

(8.5.2.1)

O parmetro emprico determina a importncia relativa de de de S. O parmetro dever variar ao longo do processo iterativo, havendo diferentes aproximaesque podem ser seguidas. A condio S depende, geralmente, da informaogeolgica disponvel. No havendo informao pode impor-se, por exemplo, que a


11/26


________________________________________________________________________________________________

11

covarincia dos parmetros seja dada por Wm = m2 I obtendo-se, ento, para a

iterao k o sistema de equaes (Tarantola, 2004),

(J(mk)TWd-1J(mk) + 2Wm

-1) m= J(m

k)TWd-1(d - g(mk) )

(8.5.2.2)

Dependendo da escolha de S (ver Tabela 1), diferentes tipos de informao a prioripodem ser usados com o objectivo de diminuir a importncia dos parmetros sub-determinados na soluo final, contribuindo para a estabilizao desta.

frmula terica representao numrica

[ ] dzzmzm o

0

2)()( 2

, )( ioi i mm

[ ] dzzm

0

2)( 2

1 )( ii i mm

dzzm

0

)( 0,)( 221 + i ii mm

[ ]

[ ]0,

)()(

)()(0 22

2

+

dzzmzm

zmzm

o

o 0,

)(

)(22

,

2,

+

i ioi

ioi

mm

mm

[ ]

[ ]0,

)(

)(0 22

2

+

dzzm

zm 0,

)(

)(22

1

21 +

i ii

ii

mm

mm

Tabela 1. Resumo dos funcionais estabilizadores definidos para o caso de modelos 1-D. 1(Tikhonov andArsenin, 1977) minimiza a norma L2 (mnimos quadrados) do desvio do modelo relativamente a ummodelo a priori. 2(Constable et al., 1987) minimiza a rugosidade do modelo no sentido dos mnimosquadrados. 3(Rudin et al., 1992) minimiza a variao total dos parmetros. indiferente magnitude dasvariaes dos parmetros podendo preservar variaes bruscas na soluo. Como no derivvel em

pontos onde o gradiente nulo, introduz-se um valor estabilizador .4

(Portniaguine and Zhdanov, 1999)Para suficientemente pequeno, este estabilizador penaliza do mesmo modo qualquer desvio dosparmetros, relativamente ao modelo a priori, independentemente do seu valor. 5 (Portniaguine andZhdanov, 1999) Para suficientemente pequeno, penaliza qualquer variao sbita nos parmetros,independentemente do seu valor. Este estabilizador tende a produzir modelos 1-D com um nmeromnimo de camadas.

Exerccio. Considere o perfil de anomalia gravimtrica dado na tabela seguinte.Supondo que a anomalia devida ao preenchimento de uma bacia sedimentar(material de massa volmica de =300 kg/m3) calcule a topografia da base da bacia,impondo a condio de esta ter uma variao suave. Sugesto: considere, para cadaobservao, a aproximao gi=2Gihipara o clculo da anomalia gravimtrica noponto de observao ionde hirepresenta a espessura da placa de massa volmica Isob o ponto de observao.

coordenadas g (mGal)-300 0.15-250 0.12-200 -0.56-150 -0.90-100 -0.95-50 -1.100 -1.2350 -1.35

100 -1.46150 -0.96200 0.02250 0.17

300 0.29


12/26


________________________________________________________________________________________________

12

8.5.3. Os mtodos de regularizao (mtodo de Tikhonov)

O PI tal como analisado no pargrafo anterior, pode ser considerado um caso

particular dos mtodos no lineares de regularizao desenvolvidos por Tikhonov(Tikhonov and Arsenin, 1977). A ideia bsica desta metodologia a de calcular,entre as possveis solues do PI, aquela que mais robusta. A escolha destasoluo feita usando-se um funcional S(m) designado porfuncional estabilizador(a Tabela 1 resume alguns dos funcionais mais utilizados em problemas de GEM).A soluo obtida impondo-se que o funcional S (na expresso 4.2.1) seja mnimo,

S(m) = mnimo(8.5.3.1)

e que se cumpra a condio de que a resposta do modelo esteja de acordo com asobservaes (dados). Esta ltima condio , geralmente, traduzida por,

d= || g(m) - d||2= .(8.5.3.2)

O problema inverso fica ento formulado do seguinte modo: conhecido o vector dedados d(d1,,dN), deve determinar-se o vector de parmetros m(m1,,mM) que, deacordo com o modelo dc= g (m), minimiza a funo S sujeito condio d= .Trata-se, ento, de um problema de optimizao que pode ser resolvido definindo-se uma funo objectivo mais geral (global) e recorrendo-se ao mtodo dosmultiplicadores de Lagrange. O problema inicial transforma-se, ento, no problema:

estimar mque minimiza a funo objectivo = d + S .

O parmetro de regularizao , determinado impondo-se a condio

d= || g(m) d||2=

(8.5.3.3)

onde m a soluo do PI que minimiza a funo objectiva global para um dado .

Definindo,

d= || Wd( g(m) d) ||2

(8.5.3.4)

S = || Wm(m mo) ||

2

(8.5.3.5)

onde Wm uma matriz diagonal que controla a importncia do afastamento decada parmetro ao modelo a priorimo. Aqui, Wd= (1/d) I. Desprezando os termosde ordem superior R(mk-1, m) no desenvolvimento em srie de Taylor da respostado modelo, as correces aos parmetros do modelo, para a iterao k podem serobtidas pela resoluo do sistema de equaes,

(J(mk)TWdTWdJ(m

k) + WmT Wm) m= J(m

k)TWdTWd(d- g(m

k) ) - WmT

Wm(mk mo) . (8.5.3.6)

Exerccio. Suponha-se um modelo que relaciona as observaes yi(xi) com a

distribuio da propriedade dj(xj) atravs de:


13/26


________________________________________________________________________________________________

13

[ ]=

= +=

7

1 22 .)500()(.628M

j ji

ji

xx

dy

Obtenha-se a distribuio de dconsiderando-se os valores deyda tabela (unidadesarbitrrias), impondo uma variao suave a d.

xI yi-600 2.51-400 2.99-200 3.24

0 3.31200 3.24400 2.99600 2.51

Exerccio. Elabore um algoritmo para a resoluo do PI no caso de sondagenselctricas verticais, assumindo que a resistividade varia suavemente emprofundidade.

8.5.4. O mtodo do gradiente descendente

O mtodo do gradiente descendente um dos mtodos mais usados na resoluo deproblemas de minimizao de funcionais. Calcule-se a variao de primeira ordemdo funcional ,

= 2 (Wdg(m))T(Wd( g(m) d)) + 2 (Wmm)

T(Wm(m mo))(8.5.4.1)

Como,g(m) = Jm

(8.5.4.2)

obtm-se, ento

g(m) = - k l(m)(8.5.4.3)

onde k um real positivo a determinar e l(m) representa o gradiente do funcional,

l(m) = JTWd

2( g(m) d)) + Wm2(m mo)

(8.5.4.4)

De acordo com o processo iterativo, os parmetros para a iterao k+1 seroestimados por,

mk+1= mk k l(mk)

(8.5.4.5)

O coeficiente k calculado a partir da condio de que a funo

(mk+1) = (mk k l(mk)) deve ser mnima relativamente a k.


14/26


________________________________________________________________________________________________

14

8.5.5. O mtodo de Newton

Considere-se a primeira iterao do processo iterativo, em que

m1= m0+ m

(8.5.5.1)

Para o clculo do vector que contm as correces m aos parmetros queconstituem o modelo inicial m0define-se a funo objectivo global,

= || Wd( g(m0+ m) d) ||2+ || Wm(m

0+ m mo) ||2 = || Wd( g(m

0) + Jm d) ||2+

+ || Wm(m0+ m mo) ||

2(8.5.5.2)

A condio necessria para que a funo objectivo seja minimizada que a variaode primeira ordem do funcional seja nula. Esta condio cumpre-se quando,

JTWd2( g(m0) + Jm d) + Wm

2(m0+ m mo) = 0(8.5.5.3)

Daqui se obtm,

m= - H1-1l(m

0)(8.5.5.4)

A matriz H1 designada por matriz Hessiana regularizada e definida por,

H1= JTWd

2J + Wm2

(8.5.5.5)

De acordo com o algoritmo de Newton o vector dos parmetros, na iterao k+1ser calculado por,

mk+1= mk Hk-1l(m

k)(8.5.5.6)

Este mtodo converge, em geral, mais rapidamente que o mtodo do gradientedescendente. Contudo, o clculo da matriz Hessiana inversa pode no ser muitofcil. Com o objectivo de resolver esta dificuldade frequente usar-se umaaproximao, Gk verdadeira inversa da matriz Hessiana calculando-se ascorreces aos parmetros de acordo com,

mk+1= mk kkGkl(m

k)(8.5.5.7)

onde o parmetro kk calculado como se exps em pargrafo anterior.

8.5.6. O mtodo do gradiente conjugado

No mtodo do gradiente conjugado, tal como mo mtodo do gradiente descendente,calculam-se os parmetros na iterao k+1usando a expresso,

mk+1= mk kk l(mk)(8.5.6.1)


15/26


________________________________________________________________________________________________

15

A diferena reside no modo de calcular os gradientes l(mk). Neste mtodo ogradiente na iterao k+1 calculado como sendo uma combinao linear dogradiente descendente calculado para esta iterao e o gradiente da iteraoanterior:

l(mk+1) = l(mk+1) + kl(mk)(8.5.6.2)

O valor do parmetro kk calculado de acordo com o exposto anteriormente. Oscoeficientes k so calculados por (Zhdanov and Keller, 1994),

k= (lT(mk) l(mk) )/( lT(mk-1) l(mk-1) ) .(8.5.6.3)

A vantagem do mtodo do gradiente conjugado que, em princpio, pode evitar-seque o processo de convergncia fique preso num mnimo local, convergindo maisrapidamente para o mnimo global.

8.6. Resoluo do problema inverso no linear usando mtodos de minimizaoglobal

Os mtodos de resoluo do PI expostos no pargrafo 4 so adequados quando atopologia da funo objectivo tem apenas um mnimo. Esta situao no , contudo,a mais frequente, conduzindo a que a soluo obtida dependa do modelo inicial e doparmetro de estabilizao. Com o objectivo de evitar, ou minimizar, estasdependncias recorre-se por vezes ao uso de mtodos de resoluo do PI que podemser considerados como mtodos de minimizao global. Analisam-se, em seguida,dois destes mtodos: o simulated annealing e o mtodo controlled randomsearch (Price, 1977) que traduzimos aqui por procura aleatria controlada. Estesmtodos podem ser considerados como variantes do mtodo de Monte Carlo. Aideia bsica a de investigar (aleatoriamente) o espao de parmetros (toexaustivamente quanto possvel) levando lentamente a investigao a centrar-se,mas no exclusivamente, em torno do mnimo absoluto. Para cada conjunto deparmetros (modelo) determina-se o valor da funo objectivo retendo-se osmodelos que apresentem valores inferiores a um valor estabelecido. Uma anliseestatstica deste conjunto restrito de modelos permite a estimativa das funes dedistribuio a posterioridos parmetros.

Os mtodos de minimizao globaltm, relativamente aos mtodos de minimizaolocalexpostos no pargrafo 4, algumas vantagens. Entre estas assinala-se o facto deno ser necessrio o conhecimento da matriz das sensibilidades, que difcil decalcular para alguns modelos. A maior desvantagem reside no facto de consumiremmuito tempo de clculo. Este facto faz com que a sua aplicao prtica, em PGEM,se restrinja a modelos 1-D e, em alguns casos, a modelos 2-D.

8.6.1. Simulated annealing (SA)

O nome deste algoritmo traduz a semelhana existente entre o algoritmo de SA e oprocesso fsico de arrefecimento de metais. O mtodo inicia a procura da soluo doPI explorando aleatoriamente o espao de parmetros. Lentamente o espao deparmetros explorado restringido ao sub-espao onde se encontra o mnimoabsoluto. Este processo semelhante ao arrefecimento do metal e o seu controlo feito atravs de uma varivel designada, nos algoritmos SA, por temperatura.

O algoritmo de Metropolis et al. (1953) , tambm, um algoritmo iterativo. Aiterao k inicia-se pela gerao aleatria de uma soluo (com base na soluo da

iterao k-1) m

k

e do clculo da respectiva funo objectiva

k

. Se

k

<

k-1

asoluo aceite e o processo iterativo continua pela gerao de uma nova soluo.


16/26


________________________________________________________________________________________________

16

Se k> k-1a aceitao, ou no, da soluo mkobriga ao clculo da funo P() =exp((k-1- k)/T) e gerao aleatria do nmero real r, tal que r[0,1]. A soluomk aceite se r < P() prosseguindo o processo iterativo (mas mantendo Tconstante). T a temperatura, que no incio do processo iterativo deve ser elevada,diminuindo (de acordo com uma lei estabelecida) ao longo do processo.

Nas fases iniciais, e porque T elevado, so aceites solues muito diferentes. medida que T diminui a aceitao de solues torna-se mais exigente e apenas asmelhores solues so aceites levando o algoritmo a convergir. O processo iterativotermina quando (por exemplo) a variao da funo objectivo entre iteraes inferior a um valor preestabelecido.

Para o mesmo problema, o algoritmo SA pode ser corrido vrias vezes, comcondies iniciais diferentes, de modo a obter-se um conjunto (estatisticamente)significativo de solues. A anlise estatstica deste conjunto de solues permitecalcular o modelo mdio e a matriz de covarincia a posterioridos parmetrosdo modelo usando, respectivamente (Mrinal et al., 1993):

= m P(m)(8.6.1.1)

C= (m- ) (m- )TP(m)(8.6.1.2)

onde as somas so sobre os modelos estimados e P(m) a densidade deprobabilidade a posterior,

P(m) exp ((m))(8.6.1.3)

Exerccio. Resolva os exerccios do pargrafo 8.4.1 usando o mtodo SA.

Exerccio. Elabore um algoritmo para a resoluo do PI no caso de sondagensmagneto-telricas, assumindo um modelo de terreno estratificado e usando SA.

8.6.2. Controlled Random Search(CRS) -Procura aleatria controlada

Na aplicao do algoritmo (iterativo) de Price (1977) principia-se pela geraoaleatria de L modelos com a condio de que os parmetros pertenam a umsubdomnio do domnio (permitido) dos parmetros. Isto significa que osparmetros tm um limite inferior e superior,

mI< mest

< ms. (8.6.2.1)

Percebe-se facilmente que quanto maior for L mais exaustiva se torna a procura domnimo absoluto. Calcula-se a funo objectivo para cada modelo identificando-se o modelo mmaxque apresentar o valor mximo para = max. Na iterao k, umnovo modelo mk estimado tendo-se em ateno os parmetros de (M+1) modelosescolhidos aleatoriamente entre os L modelos. O vector dos parmetros do novomodelo dado por,

mk= 2 c mM+1)(8.6.2.2)


17/26


________________________________________________________________________________________________

17

onde c o centride dos M primeiros modelos escolhidos, mM+1) o vector dosparmetros correspondente ao (M+1)-simo modelo escolhido sendo M o nmerode parmetros a determinar. Se os elementos de mk obedecerem aosconstrangimentos impostos e se o valor da correspondente funo objectivo forinferior a max o anterior modelo mmax retirado dos L modelos iniciais e

substitudo pelo novo modelo. Se o novo modelo no satisfizer os constrangimentosou se a correspondente funo objectivo for maior que max, o novo modelo no considerado procedendo-se a nova seleco aleatria de (M+1) modelos. Asubstituio de um modelo de L por um novo modelo constitu uma iterao. Oprocesso iterativo termina quando as funes objectivo de todos os L modelostiverem um valor inferior a um valor definido .

O conjunto de L possveis solues estimadas pelo mtodo CRS formam umcluster de modelos na vizinhana do mnimo absoluto. O modelo representadopelo centride deste cluster, m*, dever ser mais significativo que qualquer dosmodelos do conjunto e considerado como a melhor estimativa dos parmetros domodelo. A matriz de covarincia dos parmetros pode ser estimada por (Silva andHohmann, 1983),

Cij= =

L

kL 11

1(mik mI

*)(mjk mj*)

(8.6.2.3)

onde mikrepresenta a k-sima estimativa do parmetro i e m I* a melhor estimativa

do parmetro I, dada por

mI*=

=

L

kL 1

1mIk

*.

(8.6.2.4)

A matriz de correlao R obtida normalizando-se C,

Rij= Cij/(CiiCjj)1/2

(8.6.2.5)

Exerccio. Resolva o exerccio do pargrafo 8.4.3 usando o mtodo CRS.

8.7. Anlise da soluo do PI

A anlise da soluo do PI deve permitir uma avaliao da qualidade da soluoobtida. Assim deve ser possvel:

-

determinar a covarincia e/ou correlao dos parmetros do modelo,-

determinar o intervalo de confiana dos parmetros.

A anlise destas grandezas no oferece dificuldade de maior, quando se trata demodelos obtidos usando o mtodo dos mnimos quadrados (geralmente problemas1-D). Nos casos em que o PI resolvido usando-se algoritmos de regularizao (2-D e 3-D), o problema bastante mais difcil, devido densa parametrizaoutilizada, no havendo nenhuma metodologia definida. usual, contudo, proceder-se a uma anlise da sensibilidade do modelo com o objectivo de estudar-se aestabilidade da soluo e definirem-se os parmetros que so realmenteimportantes. Este tema ser abordado no pargrafo 7.Considere-se, ento, a soluo clssica do PI, obtida usando-se o mtodo dosmnimos quadrados. Assume-se que os dados obedecem a uma distribuio


18/26


________________________________________________________________________________________________

18

Gaussiana. Aceita-se, ainda, que os parmetros do modelo a priori, usado comoinformao a priori, obedece tambm a uma distribuio Gaussiana.

Matriz de covarincia

Assumindo que os dados os parmetros do modelo a priori no estocorrelacionados, a matriz de covarincia do modelo estimado pode ser calculada por(Menke, 1984),

C= Jk-gCd(Jk

-g)T+ (I R) Cm(I R)T

(8.7.1)

onde, Cd e Cm so as matrizes de covarincia dos dados e dos parmetros domodelo a priorie,

Jk-g= CmJk

T(Cd+ JkCmJkT)-1

(8.7.2)

R= Jk-gJk(8.7.3)

sendo a matriz sensibilidade Jkcalculada na iterao final, k.

Intervalo de confiana e valores extremos dos parmetros

Seja moa soluo ptima para um dado PI. A esta soluo corresponde um valor(mnimo) da funo objectivo, o. O valor da diferena entre os dados e a respostadaquele modelo ser bo= d

o g(mo). Para um modelo m, na vizinhana de motem-se,

m= mo+

b= bo+ ce= o+

Seja a funo objectivo definida por = || b Jm||2. Facilmente se mostra que

= cTc= TV 2VT(8.7.4)

Esta a equao de um hiper-elipside com centro em mo e com semi-eixos decomprimento (1/2-1) paralelos aos vectores prprios v. A forma quadrtica detem uma distribuio 2com um nmero de graus de liberdade igual ao nmero

M de parmetros. A superfcie de confiana (1-) % ser, ento, limitada no espaode parmetros pelo elipside6

TV 2VT= M, 2.

(8.7.5)

Com o objectivo de estudar os valores extremos dos parmetros, defina-se (noespao dos parmetros) o funcional linear:

6A superfcie de confiana (1-) obtida considerando os parmetros mtais que

Documents

Problem a in Verso 2