Upgrade to Pro — share decks privately, control downloads, hide ads and more …

PhD defense

PhD defense

Slides from my PhD thesis defense "Forward modeling and inversion of gravitational fields in spherical coordinates".

More information at http://www.leouieda.com

Leonardo Uieda

April 29, 2016
Tweet

More Decks by Leonardo Uieda

Other Decks in Science

Transcript

  1. modelagem direta
    e
    inversão
    de campos gravitacionais em
    coordenadas esféricas
    Observatório Nacional – 29 de Abril de 2016

    View Slide

  2. Leonardo Uieda
    Valéria C F Barbosa (orientadora)

    View Slide

  3. A COORDENADORA DE PÓS-GRADUAÇÃO ADVERTE:
    Efeitos colaterais dessa defesa incluem:
    discussão
    discórdia
    crítica
    choro
    aprovação
    reprovação

    View Slide

  4. introdução

    View Slide

  5. View Slide

  6. View Slide

  7. View Slide

  8. conhecido

    View Slide

  9. conhecido
    topografia
    batimetria

    View Slide

  10. conhecido
    topografia
    batimetria
    crosta média
    manto médio

    View Slide

  11. View Slide

  12. bacias

    View Slide

  13. bacias
    intrusões

    View Slide

  14. bacias
    intrusões
    crosta
    oceânica

    View Slide

  15. bacias
    intrusões
    crosta
    oceânica
    Moho

    View Slide

  16. bacias
    intrusões
    crosta
    oceânica
    Moho
    heterogeneidades

    View Slide

  17. View Slide

  18. inferir

    View Slide

  19. GRACE – EUA/Alemanha
    2002 - 2017

    View Slide

  20. GOCE – ESA
    2009 - 2013

    View Slide

  21. distúrbio da
    gravidade
    (GGM05G)

    View Slide

  22. dados ✔

    View Slide

  23. dados ✔
    subsuperfície
    aproximação esférica

    View Slide

  24. problema inverso

    View Slide

  25. problema inverso

    View Slide

  26. problema inverso
    modelagem
    direta

    View Slide

  27. problema inverso
    modelagem
    direta
    otimização

    View Slide

  28. problema inverso
    modelagem
    direta
    otimização regularização

    View Slide

  29. programação

    View Slide

  30. 1 2
    3

    View Slide

  31. programa A programa B
    Novo método

    View Slide

  32. programa B
    modelagem direta
    aprox. esférica
    Novo método

    View Slide

  33. fatiando a terra
    modelagem direta
    otimização
    regularização, etc
    modelagem direta
    aprox. esférica
    Novo método

    View Slide

  34. fatiando a terra
    modelagem direta
    otimização
    regularização, etc
    método: inversão não-linear rápida
    aplicação: Moho da América do Sul
    modelagem direta
    aprox. esférica

    View Slide

  35. Introdução
    Tesseroids
    Fatiando a Terra
    Inversão Moho
    Conclusão

    View Slide

  36. Introdução
    Tesseroids
    Fatiando a Terra
    Inversão Moho
    Conclusão

    View Slide

  37. modelagem direta em coordenadas esféricas

    View Slide

  38. GEOPHYSICS
    Submetido: Mar 2015
    R1: Nov 2015
    Aceito: Abr 2016
    Geophysical Software and Algorithms

    View Slide

  39. tesseroide = prisma esférico

    View Slide

  40. g
    z
    =Gρ∭
    λ
    1
    ϕ
    1
    r
    1
    λ2
    ϕ2
    r
    2
    K
    z
    dr d ϕd λ

    View Slide

  41. g
    z
    =Gρ∭
    λ
    1
    ϕ
    1
    r
    1
    λ2
    ϕ2
    r
    2
    K
    z
    dr d ϕd λ
    cte.

    View Slide

  42. g
    z
    =Gρ∭
    λ
    1
    ϕ
    1
    r
    1
    λ2
    ϕ2
    r
    2
    K
    z
    dr d ϕd λ
    cte.
    densidade

    View Slide

  43. g
    z
    =Gρ∭
    λ
    1
    ϕ
    1
    r
    1
    λ2
    ϕ2
    r
    2
    K
    z
    dr d ϕd λ
    cte.
    densidade kernel
    (Grombein et al., 2013)

    View Slide

  44. Quadratura
    Gauss-Legendre

    View Slide


  45. a
    b
    f (x)dx≈
    b−a
    2

    i=1
    N
    W
    i
    f (x
    i
    )

    View Slide


  46. a
    b
    f (x)dx≈
    b−a
    2

    i=1
    N
    W
    i
    f (x
    i
    )

    View Slide


  47. a
    b
    f (x)dx≈
    b−a
    2

    i=1
    N
    W
    i
    f (x
    i
    )
    ordem

    View Slide


  48. a
    b
    f (x)dx≈
    b−a
    2

    i=1
    N
    W
    i
    f (x
    i
    )
    pesos
    ordem

    View Slide


  49. a
    b
    f (x)dx≈
    b−a
    2

    i=1
    N
    W
    i
    f (x
    i
    )
    pesos
    raízes de P
    N
    (x)
    ordem

    View Slide


  50. Ω
    K
    z
    (r ,ϕ,λ)d Ω

    View Slide


  51. Ω
    K
    z
    (r ,ϕ,λ)d Ω
    A ∑
    i=1
    N
    r

    j=1
    N
    ϕ

    k=1
    N
    λ
    W
    i
    W
    j
    W
    k
    K
    z
    (r
    i
    ,ϕj
    ,λk
    )

    View Slide

  52. A ∑
    i=1
    N
    r

    j=1
    N
    ϕ

    k=1
    N
    λ
    W
    i
    W
    j
    W
    k
    K
    z
    (r
    i
    ,ϕj
    ,λk
    )

    Ω
    K
    z
    (r ,ϕ,λ)d Ω
    massa pontual

    View Slide

  53. View Slide

  54. View Slide

  55. View Slide

  56. Ku (1977)
    dist. até
    observação
    dist. entre
    massas
    <

    View Slide

  57. View Slide

  58. + massas

    View Slide

  59. + massas + tesseroides

    View Slide

  60. + massas + tesseroides

    View Slide

  61. View Slide

  62. if d/L
    λ
    < D:

    View Slide

  63. if d/L
    λ
    < D:
    Divide em λ

    View Slide

  64. if d/LФ < D:

    View Slide

  65. if d/LФ < D:
    Divide em Ф

    View Slide

  66. if d/Lr < D:

    View Slide

  67. if d/Lr < D:
    Divide em r

    View Slide

  68. View Slide

  69. View Slide

  70. View Slide

  71. algoritmo

    View Slide

  72. View Slide

  73. pilha

    View Slide

  74. pilha
    adicionar no
    topo da pilha

    View Slide

  75. pilha
    “pop”

    View Slide

  76. pilha
    if d/Lx < D:
    Divide em x

    View Slide

  77. pilha
    if d/Lx < D:
    Divide em x

    View Slide

  78. pilha

    View Slide

  79. pilha
    “pop”

    View Slide

  80. pilha
    if d/Lx < D:
    Divide em x

    View Slide

  81. pilha
    if d/Lx < D:
    Divide em x

    View Slide

  82. pilha

    View Slide

  83. pilha
    “pop”

    View Slide

  84. pilha
    if d/Lx < D:
    Divide em x

    View Slide

  85. pilha
    if d/Lx < D:
    Divide em x

    View Slide

  86. pilha
    g total += g tesseroide

    View Slide

  87. pilha

    View Slide

  88. pilha
    “pop”

    View Slide

  89. pilha
    if d/Lx < D:
    Divide em x

    View Slide

  90. pilha

    View Slide

  91. pilha
    vazia
    fim
    g total

    View Slide

  92. D
    d/Lx < D

    View Slide

  93. D=1 4 tesseroides

    View Slide

  94. D=2 38 tesseroides

    View Slide

  95. D=6 936 tesseroides

    View Slide

  96. View Slide

  97. View Slide

  98. Melhor D?

    View Slide

  99. View Slide

  100. 1. pólo h=2km 1°x 1°x 1km
    2. equador h=2km 1°x 1°x 1km
    3. pólo h=260km 1°x 1°x 1km
    4. pólo h=2km 30°x 30°x 1km
    Casca x Quadratura

    View Slide

  101. View Slide

  102. 0.1% erro
    0.1% erro
    0.1% erro

    View Slide

  103. D
    1.5
    1
    8

    View Slide

  104. gzz

    View Slide

  105. software

    View Slide

  106. 29 programas
    linguagem C
    open-source: BSD license

    View Slide

  107. linha de comando

    View Slide

  108. tesseroids.leouieda.com

    View Slide

  109. 34 citações
    (google scholar)
    tesseroids.leouieda.com

    View Slide

  110. conclusão

    View Slide

  111. Discretização adaptativa:
    algoritmo melhor definido pilha x recursivo

    View Slide

  112. Discretização adaptativa:
    algoritmo melhor definido pilha x recursivo
    Quantificação do erro:
    D ótimo Ds diferentes: V, gz, gzz

    View Slide

  113. Discretização adaptativa:
    algoritmo melhor definido pilha x recursivo
    Quantificação do erro:
    D ótimo Ds diferentes: V, gz, gzz
    Implementação:
    open-source usuários desenvolvedores

    View Slide

  114. programa B
    modelagem direta
    aprox. esférica
    Novo método

    View Slide

  115. Introdução
    Tesseroids
    Fatiando a Terra
    Inversão Moho
    Conclusão

    View Slide

  116. modelagem direta e inversa para geofísica
    fatiando a terra

    View Slide

  117. SciPy 2013
    “Modeling the Earth with Fatiando a Terra”

    View Slide

  118. SciPy 2013
    “Modeling the Earth with Fatiando a Terra”
    v0.1

    View Slide

  119. SciPy 2013
    “Modeling the Earth with Fatiando a Terra”
    v0.1
    YouTube

    View Slide

  120. fatiando.org

    View Slide

  121. ~950 downloads/mês
    7 citações
    fatiando.org

    View Slide

  122. View Slide

  123. código
    github.com

    View Slide

  124. View Slide

  125. View Slide

  126. View Slide

  127. View Slide

  128. histórico

    View Slide

  129. View Slide

  130. View Slide

  131. View Slide

  132. View Slide

  133. juntar código
    disciplinas, teses
    dissertações, etc

    View Slide

  134. biblioteca
    funções, classes, etc

    View Slide

  135. fácil de aprender
    rápido de implementar
    open-source
    “bateries included”

    View Slide

  136. 1 código
    ++ usuários

    View Slide

  137. 1 código
    ++ programadores

    View Slide

  138. a biblioteca

    View Slide

  139. fatiando. gridder
    mesher
    utils
    vis
    gravmag
    seismic
    inversion
    Módulos

    View Slide

  140. fatiando. gridder
    mesher
    utils
    vis
    gravmag
    seismic
    inversion
    Módulos

    View Slide

  141. from fatiando.mesher import Prism
    modelo = [
    Prism(0, 1000, 0, 2000, 1500, 2500,
    props={'magnetization': [4,-1,3]})]
    from fatiando.vis import myv
    myv.figure()
    myv.prisms(modelo)
    myv.show()
    In:
    Out:

    View Slide

  142. from fatiando import gridder
    area = (-5000, 5000, -5000, 5000)
    shape = (50, 50)
    x, y, z = gridder.regular(area, shape, z=0)
    print(x)
    print(x.size)
    In:
    Out: [-5000. -5000. -5000. ... 5000. 5000. 5000.]
    2500

    View Slide

  143. from fatiando.gravmag import prism
    mag = prism.tf(x, y, z, modelo,
    inc=-60, dec=20)
    import matplotlib.pyplot as plt
    plt.tricontourf(y, x, mag, 50, cmap="RdBu_r")
    plt.colorbar().set_label('nT')
    In:
    Out:

    View Slide

  144. from fatiando.gravmag import transform
    up = transform.upcontinue(x, y, mag, shape,
    1500)
    plt.tricontourf(y, x, up, 20, cmap="RdBu_r")
    plt.colorbar().set_label('nT')
    In:
    Out:

    View Slide

  145. from fatiando.mesher import Tesseroid
    mod = [Tesseroid(-60, -50, -25, -20, 2000, 0,
    props={'density': 2670})]
    myv.figure()
    myv.tesseroids(mod, scale=(1, 1, 500))
    myv.earth()
    myv.continents()
    In:
    Out:

    View Slide

  146. lon, lat, h = gridder.scatter(
    (-70, -40, -35, -10), 500, z=250e3)
    from fatiando.gravmag import tesseroid
    g = tesseroid.gz(lon, lat, h, mod)
    plt.tricontourf(lon, lat, g, 50, cmap="Reds")
    plt.colorbar()
    In:
    Out:

    View Slide

  147. fatiando. gridder
    mesher
    utils
    vis
    gravmag
    seismic
    inversion
    Módulos

    View Slide

  148. fatiando. gridder
    mesher
    utils
    vis
    gravmag
    seismic
    inversion
    Módulos

    View Slide

  149. fatiando.inversion

    View Slide

  150. y
    i
    =a x
    i
    +b

    View Slide

  151. y
    1
    =a x
    1
    +b
    y
    2
    =a x
    2
    +b

    y
    N
    =a x
    N
    +b

    View Slide

  152. [y
    1
    y
    2

    y
    N
    ]=
    [x
    1
    1
    x
    2
    1
    ⋮ ⋮
    x
    N
    1
    ][a
    b
    ]

    View Slide

  153. [y
    1
    y
    2

    y
    N
    ]=
    [x
    1
    1
    x
    2
    1
    ⋮ ⋮
    x
    N
    1
    ][a
    b
    ]
    ¯
    d

    View Slide

  154. [y
    1
    y
    2

    y
    N
    ]=
    [x
    1
    1
    x
    2
    1
    ⋮ ⋮
    x
    N
    1
    ][a
    b
    ]
    ¯
    d ¯
    ¯
    A

    View Slide

  155. [y
    1
    y
    2

    y
    N
    ]=
    [x
    1
    1
    x
    2
    1
    ⋮ ⋮
    x
    N
    1
    ][a
    b
    ]
    ¯
    d ¯
    p
    ¯
    ¯
    A

    View Slide

  156. [y
    1
    y
    2

    y
    N
    ]=
    [x
    1
    1
    x
    2
    1
    ⋮ ⋮
    x
    N
    1
    ][a
    b
    ]
    ¯
    d ¯
    p
    ¯
    ¯
    A
    dados
    preditos
    matriz
    Jacobiana
    vetor de
    parâmetros

    View Slide

  157. ϕ=‖¯
    do−¯
    d‖2

    View Slide

  158. ^
    ¯
    p
    ϕ=‖¯
    do−¯
    d‖2

    View Slide

  159. ^
    ¯
    p
    ϕ=‖¯
    do−¯
    d‖2

    View Slide

  160. ^
    ¯
    p
    ϕ=‖¯
    do−¯
    d‖2
    minimizar

    View Slide

  161. ^
    ¯
    p
    ϕ=‖¯
    do−¯
    d‖2
    minimizar
    função
    do ajuste
    (misfit)

    View Slide

  162. específico genérico

    View Slide

  163. específico genérico
    y
    i
    =a x
    i
    +b
    ¯
    ¯
    A

    View Slide

  164. específico genérico
    y
    i
    =a x
    i
    +b
    ¯
    ¯
    A
    ϕ=‖¯
    do−¯
    d‖2
    minimizar

    View Slide

  165. específico genérico
    y
    i
    =a x
    i
    +b
    ¯
    ¯
    A
    ϕ=‖¯
    do−¯
    d‖2
    minimizar
    usuário
    implementa

    View Slide

  166. específico genérico
    y
    i
    =a x
    i
    +b
    ¯
    ¯
    A
    ϕ=‖¯
    do−¯
    d‖2
    minimizar
    fatiando.inversion

    View Slide

  167. from fatiando.inversion import Misfit
    class Regressao(Misfit):
    def __init__(self, x, y):
    Misfit.__init__(self, data=y, nparams=2,
    islinear=True)
    self.x = x
    def predicted(self, p):
    a, b = p
    return a*self.x + b
    def jacobian(self, p):
    A = np.empty((self.ndata,
    self.nparams))
    A[:, 0] = self.x
    A[:, 1] = 1
    return A

    View Slide

  168. from fatiando.inversion import Misfit
    class Regressao(Misfit):
    def __init__(self, x, y):
    Misfit.__init__(self, data=y, nparams=2,
    islinear=True)
    self.x = x
    def predicted(self, p):
    a, b = p
    return a*self.x + b
    def jacobian(self, p):
    A = np.empty((self.ndata,
    self.nparams))
    A[:, 0] = self.x
    A[:, 1] = 1
    return A

    View Slide

  169. from fatiando.inversion import Misfit
    class Regressao(Misfit):
    def __init__(self, x, y):
    Misfit.__init__(self, data=y, nparams=2,
    islinear=True)
    self.x = x
    def predicted(self, p):
    a, b = p
    return a*self.x + b
    def jacobian(self, p):
    A = np.empty((self.ndata,
    self.nparams))
    A[:, 0] = self.x
    A[:, 1] = 1
    return A

    View Slide

  170. from fatiando.inversion import Misfit
    class Regressao(Misfit):
    def __init__(self, x, y):
    Misfit.__init__(self, data=y, nparams=2,
    islinear=True)
    self.x = x
    def predicted(self, p):
    a, b = p
    return a*self.x + b
    def jacobian(self, p):
    A = np.empty((self.ndata,
    self.nparams))
    A[:, 0] = self.x
    A[:, 1] = 1
    return A
    y
    i
    =a x
    i
    +b

    View Slide

  171. from fatiando.inversion import Misfit
    class Regressao(Misfit):
    def __init__(self, x, y):
    Misfit.__init__(self, data=y, nparams=2,
    islinear=True)
    self.x = x
    def predicted(self, p):
    a, b = p
    return a*self.x + b
    def jacobian(self, p):
    A = np.empty((self.ndata,
    self.nparams))
    A[:, 0] = self.x
    A[:, 1] = 1
    return A
    ¯
    ¯
    A
    y
    i
    =a x
    i
    +b

    View Slide

  172. específico genérico
    y
    i
    =a x
    i
    +b
    ¯
    ¯
    A
    ϕ=‖¯
    do−¯
    d‖2
    minimizar
    usuário
    implementa

    View Slide

  173. específico genérico
    y
    i
    =a x
    i
    +b
    ¯
    ¯
    A
    ϕ=‖¯
    do−¯
    d‖2
    minimizar
    fatiando.inversion

    View Slide

  174. reg = Regressao(x, yo)
    reg.fit()
    print(reg.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]

    View Slide

  175. reg = Regressao(x, yo)
    reg.fit()
    print(reg.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]
    a b

    View Slide

  176. reg = Regressao(x, yo)
    reg.fit()
    print(reg.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]
    plt.plot(x, yo, 'ok')
    plt.plot(x, reg.predicted(), '-r')
    In:
    Out:
    a b

    View Slide

  177. reg.config('newton', initial=[1, 1]).fit()
    print(reg.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]

    View Slide

  178. reg.config('newton', initial=[1, 1]).fit()
    print(reg.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]
    reg.config('levmarq', initial=[0, 0]).fit()
    print(reg.estimate_)
    In:
    Out: [ 20.39562568 4980.21718193]

    View Slide

  179. reg.config('newton', initial=[1, 1]).fit()
    print(reg.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]
    reg.config('levmarq', initial=[0, 0]).fit()
    print(reg.estimate_)
    In:
    Out: [ 20.39562568 4980.21718193]
    reg.config('acor',
    bounds=[0, 1000, 0, 10e5]).fit()
    print(reg.estimate_)
    In:
    Out: [ 20.39543101 4980.22844959]

    View Slide

  180. ϕ(¯
    p)=‖¯
    do−¯
    d‖2
    minimizar

    View Slide

  181. ϕ(¯
    p)=‖¯
    do−¯
    d‖2
    minimizar instável

    View Slide

  182. Γ(¯
    p)=ϕ(¯
    p)+μθ(¯
    p)
    ϕ(¯
    p)=‖¯
    do−¯
    d‖2
    minimizar instável
    minimizar

    View Slide

  183. Γ(¯
    p)=ϕ(¯
    p)+μθ(¯
    p)
    ϕ(¯
    p)=‖¯
    do−¯
    d‖2
    minimizar instável
    escalar
    minimizar
    função
    regularizadora

    View Slide

  184. from fatiando.inversion import Damping
    phi = Regressao(x, yo)
    gamma = phi + 10e-7*Damping(nparams=2)
    In:

    View Slide

  185. from fatiando.inversion import Damping
    phi = Regressao(x, yo)
    gamma = phi + 10e-7*Damping(nparams=2)
    In:
    Γ(¯
    p)=ϕ(¯
    p)+μθ(¯
    p)

    View Slide

  186. from fatiando.inversion import Damping
    phi = Regressao(x, yo)
    gamma = phi + 10e-7*Damping(nparams=2)
    In:
    Γ(¯
    p)=ϕ(¯
    p)+μθ(¯
    p)
    gamma.fit()
    print(gamma.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]

    View Slide

  187. from fatiando.inversion import Damping
    phi = Regressao(x, yo)
    gamma = phi + 10e-7*Damping(nparams=2)
    In:
    Γ(¯
    p)=ϕ(¯
    p)+μθ(¯
    p)
    gamma.fit()
    print(gamma.estimate_)
    In:
    Out: [ 20.395431 4980.22844991]
    gamma.config('acor',
    bounds=[0, 1000, 0, 10e5]).fit()
    print(gamma.estimate_)
    In:
    Out: [ 20.14654754 4992.15040076]

    View Slide

  188. conclusão

    View Slide

  189. biblioteca: funções e classes
    fatiando a terra

    View Slide

  190. biblioteca: funções e classes
    modelagem, processamento, visualização
    fatiando a terra

    View Slide

  191. biblioteca: funções e classes
    modelagem, processamento, visualização
    inversão: fatiando.inversion
    fatiando a terra

    View Slide

  192. biblioteca: funções e classes
    modelagem, processamento, visualização
    inversão: fatiando.inversion
    Simples Complexos (reais)
    fatiando a terra

    View Slide

  193. biblioteca: funções e classes
    modelagem, processamento, visualização
    inversão: fatiando.inversion
    Simples Complexos (reais)
    Caching
    Matrizes esparsas
    BLAS, LAPACK, MKL
    fatiando a terra

    View Slide

  194. Usuários: downloads, citações

    View Slide

  195. Usuários: downloads, citações
    Colaboradores: locais e externos

    View Slide

  196. Open-source: BSD license
    Usuários: downloads, citações
    Colaboradores: locais e externos

    View Slide

  197. Open-source: BSD license
    Usuários: downloads, citações
    Colaboradores: locais e externos
    fatiando.org

    View Slide

  198. Open-source: BSD license
    Usuários: downloads, citações
    Colaboradores: locais e externos
    fatiando.org
    github.com/fatiando

    View Slide

  199. fatiando.org

    View Slide

  200. fatiando a terra
    modelagem direta
    otimização
    regularização, etc
    modelagem direta
    aprox. esférica
    Novo método

    View Slide

  201. Introdução
    Tesseroids
    Fatiando a Terra
    Inversão Moho
    Conclusão

    View Slide

  202. Inversão não-linear
    em
    coordenadas esféricas
    com aplicação na Moho da América do Sul
    com aplicação na Moho da América do Sul
    rápida

    View Slide

  203. inferir

    View Slide

  204. inferir

    View Slide

  205. Grid de observações

    View Slide

  206. ¯
    do

    View Slide

  207. 1 tesseroide para cada

    View Slide

  208. parâmetros = h
    h

    View Slide

  209. ¯
    p

    View Slide

  210. Estimar
    a partir de ¯
    do
    ¯
    p

    View Slide

  211. inversão
    não-linear

    View Slide

  212. d
    i
    =f
    i

    p)

    View Slide

  213. d
    i
    =f
    i

    p)
    modelagem direta
    (tesseroides)

    View Slide

  214. ¯
    r= ¯
    do−¯
    d
    Resíduos

    View Slide

  215. φ(¯
    p)=‖¯
    r‖2
    Minimizar

    View Slide

  216. Gauss-Newton

    View Slide

  217. ¯
    p0

    View Slide

  218. ¯
    p0
    ¯
    p1=¯
    p0+ ¯
    Δ p0

    View Slide

  219. ( ¯
    ¯
    AkT ¯
    ¯
    Ak) ¯
    Δ pk= ¯
    ¯
    AkT [ ¯
    do−¯
    d(¯
    pk)]
    ¯
    p0
    ¯
    p1=¯
    p0+ ¯
    Δ p0

    View Slide

  220. ( ¯
    ¯
    AkT ¯
    ¯
    Ak) ¯
    Δ pk= ¯
    ¯
    AkT [ ¯
    do−¯
    d(¯
    pk)]
    ¯
    p0
    ¯
    p1=¯
    p0+ ¯
    Δ p0
    A
    i j
    =
    ∂ f
    i
    ∂ p
    j
    Jacobiana

    View Slide

  221. ( ¯
    ¯
    AkT ¯
    ¯
    Ak) ¯
    Δ pk= ¯
    ¯
    AkT [ ¯
    do−¯
    d(¯
    pk)]
    ¯
    p0
    ¯
    p1=¯
    p0+ ¯
    Δ p0
    A
    i j
    =
    ∂ f
    i
    ∂ p
    j
    Jacobiana
    mal posto

    View Slide

  222. Regularização

    View Slide

  223. θ(¯
    p)=‖¯
    ¯
    R ¯
    p‖2
    Suavidade

    View Slide

  224. Γ(¯
    p) = φ + μθ
    Função objetivo

    View Slide

  225. Γ(¯
    p) = φ + μθ
    Função objetivo
    ajuste

    View Slide

  226. Γ(¯
    p) = φ + μθ
    Função objetivo
    ajuste
    regularização

    View Slide

  227. Γ(¯
    p) = φ + μθ
    Função objetivo
    ajuste
    regularização
    balanço

    View Slide

  228. Gauss-Newton
    ¯
    Δ p=( ¯
    ¯
    AT ¯
    ¯
    A + μ ¯
    ¯
    RT ¯
    ¯
    R)−1[
    ¯
    ¯
    AT
    ¯
    rk−μ ¯
    ¯
    RT ¯
    ¯
    R ¯
    pk ]

    View Slide

  229. Gauss-Newton
    ¯
    Δ p=( ¯
    ¯
    AT ¯
    ¯
    A + μ ¯
    ¯
    RT ¯
    ¯
    R)−1[
    ¯
    ¯
    AT
    ¯
    rk−μ ¯
    ¯
    RT ¯
    ¯
    R ¯
    pk ]
    custoso
    (computacionalmente)

    View Slide

  230. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r

    View Slide

  231. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r

    View Slide

  232. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r

    View Slide

  233. Bott
    (1960)

    View Slide

  234. ¯
    Δ p= ¯
    r
    2πG Δρ

    View Slide

  235. ¯
    Δ p= ¯
    r
    2πG Δρ
    ∂ platôde Bouguer
    ∂h

    View Slide

  236. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r

    View Slide

  237. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r

    View Slide

  238. rápido
    pouca memória
    converge

    View Slide

  239. instável
    regularização

    View Slide

  240. Silva et al.
    (2014)

    View Slide

  241. Bott
    Gauss-Newton
    caso
    particular

    View Slide

  242. Bott
    Gauss-Newton
    caso
    particular

    View Slide

  243. caso
    particular
    ¯
    Δ p= ¯
    r
    2πG Δρ
    ¯
    Δ p=( ¯
    ¯
    AT ¯
    ¯
    A)−1 ¯
    ¯
    AT
    ¯
    r

    View Slide

  244. caso
    particular
    ¯
    Δ p= ¯
    r
    2πG Δρ
    ¯
    Δ p=( ¯
    ¯
    AT ¯
    ¯
    A)−1 ¯
    ¯
    AT
    ¯
    r

    View Slide

  245. caso
    particular
    ¯
    Δ p= ¯
    r
    2πG Δρ
    ¯
    Δ p=( ¯
    ¯
    AT ¯
    ¯
    A)−1 ¯
    ¯
    AT
    ¯
    r
    A
    ii
    =2πGΔρ
    A
    ij
    =0 para i≠ j

    View Slide

  246. (2πG Δρ 0 ⋯ 0
    0 2πGΔρ ⋯ 0
    ⋮ ⋮ ⋱ ⋮
    0 0 ⋯ 2πGΔρ
    )
    ~
    A

    View Slide

  247. nesse
    trabalho

    View Slide

  248. Gauss-Newton
    ¯
    Δ p=( ¯
    ¯
    AT ¯
    ¯
    A)−1 ¯
    ¯
    AT [¯
    do−¯
    d(¯
    pk)]

    View Slide

  249. ¯
    Δ p=( ¯
    ¯
    AT ¯
    ¯
    A)−1 ¯
    ¯
    AT [¯
    do−¯
    d(¯
    pk)]
    Gauss-Newton
    tesseroides

    View Slide

  250. ¯
    Δ p=(
    ~
    AT ~
    A)−1 ~
    AT [¯
    do−¯
    d(¯
    pk)]
    Gauss-Newton
    tesseroides
    ~
    A
    ii
    =2πGΔρ
    i

    View Slide

  251. regularização
    suavidade

    View Slide

  252. ¯
    Δ p=(
    ~
    AT ~
    A+μ ¯
    ¯
    RT ¯
    ¯
    R)−1[
    ~
    AT
    ¯
    rk−μ ¯
    ¯
    RT ¯
    ¯
    R ¯
    pk ]

    View Slide

  253. esparsas
    ¯
    Δ p=(
    ~
    AT ~
    A+μ ¯
    ¯
    RT ¯
    ¯
    R)−1[
    ~
    AT
    ¯
    rk−μ ¯
    ¯
    RT ¯
    ¯
    R ¯
    pk ]

    View Slide

  254. esparsas
    ~99.8% tempo de computação
    ~
    AT
    ¯
    rk−μ ¯
    ¯
    RT ¯
    ¯
    R ¯
    pk ]
    ¯
    Δ p=(
    ~
    AT ~
    A+μ ¯
    ¯
    RT ¯
    ¯
    R)−1[

    View Slide

  255. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r

    View Slide

  256. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r
    Bott

    View Slide

  257. 1. Construir
    2. Sistema linear
    3. Calcular
    ¯
    ¯
    A
    ¯
    r
    matrizes esparsas
    Bott

    View Slide

  258. rápido
    pouca memória
    converge

    View Slide

  259. instável
    regularização

    View Slide

  260. instável
    regularização

    View Slide

  261. implementação

    View Slide

  262. scipy (matrizes esparsas)

    View Slide

  263. estimar
    hiperparâmetros

    View Slide

  264. −Δρ +Δρ

    View Slide

  265. z
    ref

    View Slide

  266. Γ(¯
    p) = φ + μθ
    Função objetivo
    ajuste
    regularização
    balanço

    View Slide

  267. z
    ref
    μ
    Δρ

    View Slide

  268. validação cruzada

    View Slide

  269. validação cruzada
    μ

    View Slide

  270. validação cruzada
    μ
    z
    ref
    Δρ

    View Slide

  271. validação cruzada
    μ
    h
    ref
    Δρ

    View Slide

  272. Separar os dados
    ¯
    do

    View Slide

  273. Separar os dados
    ¯
    d
    inv
    o
    ¯
    do

    View Slide

  274. Separar os dados
    ¯
    d
    inv
    o ¯
    d
    test
    o
    ¯
    do

    View Slide

  275. para em :
    μn
    [μ1,
    …,μN
    ]

    View Slide

  276. para em :
    μn
    inversão: ¯
    d
    inv
    o ^
    ¯
    pn
    [μ1,
    …,μN
    ]

    View Slide

  277. para em :
    μn
    inversão: ¯
    d
    inv
    o
    ¯
    d
    test
    n
    ^
    ¯
    pn
    prever
    ^
    ¯
    pn
    [μ1,
    …,μN
    ]

    View Slide

  278. para em :
    μn
    inversão: ¯
    d
    inv
    o
    ¯
    d
    test
    n
    ^
    ¯
    pn
    prever
    ^
    ¯
    pn
    MSE
    n
    =
    ‖¯
    d
    test
    o −¯
    d
    test
    n‖2
    N
    test
    [μ1,
    …,μN
    ]

    View Slide

  279. MSE
    μ

    View Slide

  280. MSE
    μ

    View Slide

  281. MSE
    μ
    melhor μ

    View Slide

  282. Separar os dados
    ¯
    d
    inv
    o ¯
    d
    test
    o
    ¯
    do

    View Slide

  283. dado completo

    View Slide

  284. dado inversão

    View Slide

  285. dado teste

    View Slide

  286. dado teste dado inversão

    View Slide

  287. validação cruzada
    μ
    z
    ref
    Δρ

    View Slide

  288. validação cruzada
    μ
    z
    ref
    Δρ

    View Slide

  289. estimativas
    pontuais
    ¯
    z
    s
    o

    View Slide

  290. Δρ
    z
    ref
    l, m

    View Slide

  291. para cada e :
    Δρ
    m
    z
    ref ,l

    View Slide

  292. para cada e :
    Δρ
    m
    inversão: ¯
    d
    inv
    o ^
    ¯
    pl ,m
    z
    ref ,l

    View Slide

  293. para cada e :
    Δρ
    m
    inversão: ¯
    d
    inv
    o
    ¯
    z
    s
    l ,m
    ^
    ¯
    pl ,m
    interpolar ^
    ¯
    pl ,m
    z
    ref ,l

    View Slide

  294. para cada e :
    Δρ
    m
    inversão: ¯
    d
    inv
    o
    ¯
    z
    s
    l ,m
    ^
    ¯
    pl ,m
    interpolar ^
    ¯
    pl ,m
    MSE
    l,m
    =
    ‖¯
    z
    s
    o−¯
    z
    s
    l,m‖2
    N
    s
    z
    ref ,l

    View Slide

  295. Δρ
    z
    ref

    View Slide

  296. Δρ
    z
    ref

    View Slide

  297. Δρ
    Melhor e
    Δρ z
    ref
    z
    ref

    View Slide

  298. ^
    ¯
    p

    View Slide

  299. sintético

    View Slide

  300. Moho
    CRUST1.0
    Δρ=350 kg.m−3
    z
    ref
    =30 km

    View Slide

  301. View Slide

  302. View Slide

  303. View Slide

  304. View Slide

  305. Δρ=350 kg.m−3
    z
    ref
    =30 km
    verdadeiro

    View Slide

  306. View Slide

  307. View Slide

  308. View Slide

  309. View Slide

  310. América
    do Sul

    View Slide

  311. View Slide

  312. View Slide

  313. Assumpção et al. (2013)

    View Slide

  314. View Slide

  315. View Slide

  316. View Slide

  317. 35 km

    View Slide

  318. 35 km

    View Slide

  319. 35 km

    View Slide

  320. 35 km

    View Slide

  321. 35 km

    View Slide

  322. View Slide

  323. View Slide

  324. View Slide

  325. View Slide

  326. View Slide

  327. View Slide

  328. View Slide

  329. View Slide

  330. View Slide

  331. View Slide

  332. View Slide

  333. View Slide

  334. View Slide

  335. View Slide

  336. View Slide

  337. conclusão

    View Slide

  338. Baseado em Bott (1960) e Silva et al. (2014)
    Tesseroides
    Gauss-Newton + Regularização
    Matrizes esparsas
    Validação cruzada
    Novo método
    μ Δρ z
    ref

    View Slide

  339. Compatível com soluções anteriores (grav e sismo)
    Correções (topo e sedimentos) apropriadas
    ~6 km desvio padrão com sísmica
    Diferença grande concentrada nos Andes
    Resolução pode ser falsa
    Depende de correções corretas
    Difícil estimar o erro
    Moho América do Sul

    View Slide

  340. Compatível com soluções anteriores (grav e sismo)
    Correções (topo e sedimentos) apropriadas
    ~6 km desvio padrão com sísmica
    Diferença grande concentrada nos Andes
    Resolução pode ser falsa
    Depende de correções corretas
    Difícil estimar o erro
    Moho América do Sul

    View Slide

  341. Introdução
    Tesseroids
    Fatiando a Terra
    Inversão Moho
    Conclusão

    View Slide

  342. conclusões
    finais

    View Slide

  343. programa A programa B
    Novo método

    View Slide

  344. fatiando a terra
    modelagem direta
    otimização
    regularização, etc
    método: inversão não-linear rápida
    aplicação: Moho da América do Sul
    modelagem direta
    aprox. esférica

    View Slide

  345. Algoritmo aprimorado
    Quantificação do erro
    Usuários
    Desenvolvedores
    Desenvolvimento estagnado

    View Slide

  346. Algoritmo aprimorado
    Quantificação do erro
    Usuários
    Desenvolvedores
    Desenvolvimento estagnado

    View Slide

  347. fatiando a terra
    Biblioteca em Python
    Vários módulos
    Pesquisa + educação
    Usuários + desenvolvedores
    Além dessa tese

    View Slide

  348. fatiando a terra
    Facilitar a participação:
    Guia do desenvolvedor
    Guia do usuário
    Workshop
    Videos?

    View Slide

  349. Inversão rápida + Moho Am. do Sul
    Método novo + eficiente
    Approx. esférica
    Python + Fatiando a Terra
    Depende das correções
    Incerteza

    View Slide

  350. Inversão rápida + Moho Am. do Sul
    Método novo + eficiente
    Approx. esférica
    Python + Fatiando a Terra
    Depende das correções
    Incerteza

    View Slide

  351. código + dados + figuras
    pinga-lab.org
    leouieda.com

    View Slide

  352. https://xkcd.com/1403

    View Slide