Métodos potenciales de prospección, FCAG, 2023.
Para el modelado de anomalías magnéticas, en primer lugar escribimos una subrutina que implemente el método de Talwani-Heirtzler. El resultado que obtendremos de la subrutina de modelado es el campo magnetizante o de inducción magnética $\mathbf{H}=(H_x,0,H_z)$ en cada estación:
\begin{align} H_x(\mathbf{0}) &= 2 \left(m_x\,P+m_z\,Q\right), \\ H_z(\mathbf{0}) &= 2 \left(m_x\,Q-m_z\,P\right); \end{align}considerando, sin pérdida de generalidad, el punto de estación $\mathbf{0}=(x_E=0,z_E=0)$. El vector $\mathbf{m}=(m_x,0,m_z)$ es la polarización magnética o magnetización. Las funciones $P(x_1,x_2,z_1,z_2)$ y $Q(x_1,x_2,z_1,z_2)$ dependen de los vértices que definen el contorno del cuerpo y estan dadas por:
\begin{align} P &\equiv - \cos{\phi}\ \sin{\phi} \ln\frac{r_2}{r_1} + (\theta_2-\theta_1)\ \sin^2{\phi} ,\\ Q &\equiv \phantom{+}\sin^2{\phi}\ \ln{\frac{r_2}{r_1}} + (\theta_2-\theta_1) \cos{\phi}\ \sin{\phi} . \end{align}Donde $r_1 = \sqrt{x_1^2 + z_1^2}$; $r_2 = \sqrt{x_2^2 + z_2^2}$, $\tan{\theta_1} = \dfrac{z_1}{x_1}$; $\tan{\theta_2} = \dfrac{z_2}{x_2}$; $\tan{\phi} = \dfrac{z_2-z_1}{x_1-x_2}$ y
\begin{align} \sin{\phi} &= \frac{z_2-z_1}{\big[(x_1-x_2)^2+(z_2-z_1)^2)\big]^{1/2}};\\ \cos{\phi} &= \frac{x_1-x_2}{\big[(x_1-x_2)^2+(z_2-z_1)^2\big]^{1/2}}. \end{align}Estas expresiones son válidas para cada lado del polígono si se recorre en sentido antihorario. Los subíndices $i=1,2$ en $x_i,z_i$ indican vértices sucesivos. Para obtener el campo magnetizante del polígono, sumamos las contribuciones de cada lado para cada punto de observación.
Las componentes del vector magnetización $\mathbf{m}$ deben ser proyectadas al plano de observación del perfil. Para ello, dados los ángulos de inclinación $I$, declinación $D$, y el azimut $C$ del eje X del perfil, resulta:
\begin{align} m_x &= ||\mathbf{m}|| \cos(I)\cos(C-D),\\ m_z &= ||\mathbf{m}|| \sin(I). \end{align}A raíz de lo anterior, una subrutina para calcular el vector magnetización puede adoptar la siguiente estructura:
def campo_magnetizante(x,m,vertices,I,C,D):
...
mx = m * np.cos(I) * np.cos(C-D)
mz = m * np.sin(I)
H = []
...
# loop sobre estaciones:
for xi in x:
xe, ze = xi # Coordenadas del punto de estación.
Hx, Hz = 0.,0.
...
# loop sobre lados del polígono:
for v in vertices:
x1, x2, z1, z2 = v
x1, x2, z1, z2 = x1-xe, x2-xe, z1-ze, z2-ze # para evaluar en el punto de estación (0,0)
...
Hx += 2.0 *( mx * P + mz * Q )
Hz += 2.0 *( mx * Q - mz * P )
H.append([Hx Hz])
return np.array(H)
Como el campo de inducción magnética $\mathbf{H}$ no es el observable de los métodos potenciales, proyectamos $\mathbf{H}$ en la dirección del campo magnético ambiental $\mathbf{B}$ en el plano del perfil. Obtenemos así la anomalía escalar de intensidad total (total field anomaly, TFA):
$$\Delta T := \mathbf{H}^T\dfrac{\mathbf{B}}{||\mathbf{B}||}= H_x{v_x} + H_z{v_z}.$$Asumimos que la magnitud del campo magnetizante es mucho menor que la magnitud del campo principal del lugar de medición, lo que asemeja la anomalía calculada con las anomalías registradas en el campo.
Las componentes del versor del campo principal $\mathbf{v}=\mathbf{B}/||\mathbf{B}||=(v_x,0,v_z)$ pueden obtenerse de los ángulos de inclinación $I$, declinación $D$ y del azimut $C$ del eje X del perfil por medio de
\begin{align} v_x &= \cos(I)\cos(C-D),\\ v_z &= \sin(I). \end{align}Un programa para obtener la anomalía escalar de intensidad total podría ser:
def anomalia_escalar_de_intensidad_total(H,B):
...
Hx, Hz = H
...
return Hx*vx + Hz*vz
Definimos el modelo y procedemos a la evaluación de la anomalía escalar total de intensidad $\Delta T$. El modelo es parametrizado por el valor de la susceptibilidad $\kappa$, los vértices de la sección transversal $(x_1,x_2,z_1,z_2)$, la orientación del perfil y los vectores involucrados.
En el ejemplo ilustrativo, consideramos magnetización puramente inducida por el campo principal,
$$\mathbf{m}=\kappa\mathbf{B},$$y una susceptibilidad $\kappa=0,012$ (unidades cgs). El campo principal tiene intensidad $||\mathbf{B}||_2=50000$ nT, inclinación $I = 0^\circ$ y declinación $D=10^\circ$. El perfil es paralelo a la dirección N–S (por lo tanto es $C=0^\circ)$.
Sección poligonal: coordenadas en X: [-3000. -3000. 3000. 3000.] coordenadas en Z: [ 3000. 10000. 10000. 3000.] número de vértices: 4
Calculamos el campo de inducción magnética $\mathbf{H}$ producido por el cuerpo en cada estación de medición $x$:
x = ...
H = campo_magnetizante(x,m,vertices,I,D,C)
Finalmente, proyectamos el campo de inducción magnética $\mathbf{H}$ en cada punto de la estación sobre el campo principal, obteniendo la anomalía TFA denotada por $\Delta T$:
TFA = anomalia_escalar_de_intensidad_total(H,B)
Representamos gráficamente los resultados e interpretamos los resultados numéricos. Los vectores no están a escala.
Graficamos de forma interactiva para distintos valores de la inclinación $I$, la susceptibilidad $\kappa$, el techo y la base del cuerpo anómalo rectangular. Los vectores representados no están a escala.
En esta visualización interactiva podemos analizar como la anomalía puede representar localmente un refuerzo del campo principal ($\Delta T>0$), una debilitación del campo principal ($\Delta T<0$) o llegar a ser nula.
Además se hace evidente el efecto del ángulo de inclinación con la forma de la anomalía a lo largo del perfil. Este último detalle se vinculará con la herramienta de procesamiento de reducción al polo (RTP).
interactive(children=(IntSlider(value=90, description='$I$ ($\\circ$)', max=90, min=-90, step=5), FloatSlider(…
Eso es todo por hoy.