Métodos potenciales de prospección, FCAG, 2024.
Macros $\LaTeX$ $\newcommand{\B}[1]{\mathbf{#1}}$ $\def\transpose{\top}$
La anomalía escalar de intensidad total (TFA) de un prisma recto rectangular semi infinito en la dirección vertical y con magnetización uniforme, se modela utilizando la expresión de la ecuación 10 del artículo científico de Bhattacharyya (1964):
$$ \begin{align*} \dfrac{F(x,y,z=0)}{I_p} = \left[ \dfrac{\alpha_{23}}{2} \log{(\dfrac{r_0-\alpha_1}{r_0+\alpha_1})} +\dfrac{\alpha_{13}}{2} \log{(\dfrac{r_0-\beta_1}{r_0+\beta_1})} -\alpha_{12} \log{(r_0+h)} \\ -lL\tan^{-1}{(\dfrac{\alpha_1\beta_1}{\alpha_1^2+r_0h+h^2})} -mM\tan^{-1}{(\dfrac{\alpha_1\beta_1}{r_0^2+r_0h-\alpha_1^2})} +Nn \tan^{-1}{(\dfrac{\alpha_1\beta_1}{r_0h})} \right] \Big|^{\alpha_u} _{\alpha_l}\Big|^{\beta_u} _{\beta_l}\: . \end{align*} $$$F(x,y,0)$ es la TFA medida en $z=0$. $I_p$ es la magnitud de la magnetización, $h$ la profundidad del cuerpo debajo del nivel de observación. Los valores $l,m,n$ son los cosenos directores del campo principal y $L,M,N$ los cosenos directores del vector polarización. Luego $\alpha_{12}=Lm+Ml$, $\alpha_{13}=Ln+Nl$ y $\alpha_{23}=Mn+Nm$. Por último $\alpha_1=\alpha-x$, $\beta_1=\beta-y$ y $r^2_0=\alpha_1^2+\beta_1^2+(h-z)^2$, donde $\alpha_u,\alpha_u$ y $\beta_l,\beta_l$ son los límites superiores e inferiores de $\alpha_1$ y $\beta_1$, respectivamente. Las coordenadas del punto de observación son, en la notación original, $x,y,z$. Los vértices del cuerpo vienen dados por $\alpha_l,\alpha_u$, ($\alpha_l<\alpha_u$) en la dirección $X$, $\beta_l,\beta_u$ ($\beta_l<\beta_u$) en la dirección $Y$ y $h$ en la dirección $Z$. El sistema de coordenadas considera $X\to \textrm{Norte}$, $Y\to\textrm{Este}$, $Z\to\textrm{Profundidad}.$ Entonces, $z>0$ hacia abajo.
La expresión puede ser utilizada para obtener la TFA de un cuerpo con sección vertical finita $z_2-z_1$. Para ello, consideramos un prisma semi infinito con la cara superior a una profundidad $z_1$ del plano de referencia y otro prisma semi infinito con cara superior a una profundidad $z_2>z_1$. Obtenemos la anomalía de restar al resultado obtenido para $z_1$ el resultado para $z_2$.
Podemos cambiar la notación para modernizar la expresión original del artículo científico utilizando $\hat{m}_x,\hat{m}_y,\hat{m}_z=n,m,l$ y $\hat{F}_x,\hat{F}_y,\hat{F}_z=N,M,L$ para los cosenos directores involucrados. También podríamos escribir $x_p,y_p,x_p$ para las coordenadas del punto de observación y utilizar $x_1,x_2,y_1,y_2,z_1,z_2$ para las coordenadas de los vértices del prisma. Con ello, podemos escribir $\alpha_1=x-x_p$, $\beta_1=y-y_p$, $\alpha_{12}=\hat{m}_x\hat{F}_y+\hat{m}_y\hat{F}_x$,$\alpha_{13}=\hat{m}_x\hat{F}_z+\hat{m}_z\hat{F}_x$, $\alpha_{23}=\hat{m}_y\hat{F}_z+\hat{m}_z\hat{F}_y$ y $r^2_0=(x-x_p)^2+(y-y_p)^2+(z-z_p)^2$. Con esta notación $x \in \{x_1,x_2\}$, $y \in \{y_1,y_2\}$ y $z=z_1$ para la cara superior y $z=z_2$ para la inferior. La anomalía pasaría a ser $\Delta T(x,y,z=0)$.
Los cosenos directores de un vector $\mathbf{v}$ pueden ser descritos por: $\hat{v}_x,\hat{v}_y,\hat{v}_z=\cos{I}\cos{D},\cos{I}\cos{D},\sin{D}$.
Para obtener el resultado en unidades SI multiplicamos la expresión por $C_m =\mu_0 / 4 \pi = 10^{-7}$
Nm $^{-2}$ para obtener la anomalía en Tesla si utilizamos magnetización en unidades de Am $^{-1}$.
$F(\alpha,\beta) \Big|^{\alpha_u}_{\alpha_l}\Big|^{\beta_u}_{\beta_l} = +F(\alpha_u,\beta_u)+F(\alpha_l,\beta_l) - [ F(\alpha_u,\beta_l) + F(\alpha_l,\beta_u) ]$. Al momento de programar la expresión, esta suma podría hacerse con cuatro llamados a la expresión de $F()$ con el signo apropiado.
Generamos el campo de anomalía total (TFA) $\B{\Delta T}$ en unidades de nT para un prisma recto rectangular. El prisma tiene vértices en las coordenadas $x_1,x_2=-1.5,+1.5$ km, $y_1,y_2=-0.5,+0.5$ km y $z_1,z_2=0, 2$ km. Consideramos un intensidad de magnetización $I_p=1$ Am$^{-1}$ y $I_m=-30^\circ$ y $D_m=-20^\circ$ para el vector de magnetización, y $I_f=-60^\circ$ y $D_f=+23^\circ$, para el campo principal. En este ejemplo el vector magnetización es diferente al campo principal, lo que indica la existencia de una componente remanente.
El área de relevamiento va de -7 km, a +7 km en ambas direcciones. Evaluamos la anomalía de intensidad total en $100\times 100$ puntos de observación con igual espaciamiento en cada dirección y elegimos $z_p$ a 300 m sobre el plano horizontal. Notemos que dado que el sistema de coordenadas es con z>0 hacia abajo, en el código debemos utilizar $z_p=-300$ m para ser consistentes.
Generamos los puntos de observación:
n = 100 # [--]:número de puntos de observación en cada dirección.
h = -300 # [m]: altura de las observaciones (Z>0 hacia abajo).
X0,Y0, = 7e3,+7e3 # [m]: límites del área de relevamiento.
x = np.linspace(-X0,+X0,n) # [m]: coordenada en X del punto de observación.
y = np.linspace(-Y0,+Y0,n) # [m]: coordenada en Y del punto de observación.
z = h * np.ones_like(x) # [m]: coordenada en Z del punto de observación.
dx [m]: 141.41414141414143 , dy [m]: 141.41414141414143 , nx: 100 , ny: 100 , N: 10000 Altura de observación [m]: 300
Asignamos los parámetros físicos y geométricos del modelo:
I_f, D_f = -30, -20 # [deg]: dirección del campo geomagnético principal
I_m, D_m = -60, +23 # [deg]: dirección del vector magnetización (caso de magnetización remanente).
m = 1. # [A/m]: magnitud del vector magnetización.
x1,x2 = -1500, +1500 # [m]: coordenadas de los vértices del prisma.
y1,y2 = -500, +500 # [m]: coordenadas de los vértices del prisma.
z1,z2 = 0, +2000 # [m]: coordenadas de los vértices del prisma.
Dirección del campo principal (I,D) [deg]: -30 , -20 Dirección del vector magnetizante (I,D) [deg]: -60 , 23 Vector de intensidad [A/m] [ 0.46025243 0.19536556 -0.8660254 ] magnitud [A/m] 1.0 Coordenadas de los vértices [m]: x1, x2: -1500 1500 y1, y2: -500 500 z1, z2: 0 2000 . Espesor [m]: 2000
Procedemos a calcular la TFA:
def tfa_prisma(...):
CM = 10. ** (-7) # henry/m (SI)
T2NT = 10. ** (9) # Tesla a nanotesla
...
tfa = ...
...
return tfa * CM * T2NT # [nT]
tfa = tfa_prisma(x, y, z, x1, x2, y1, y2, z1, z2, I_m, D_m, I_f, D_f, m) # [nT] y utilizando subrutina tfa_prisma()
Finalmente, visualizamos:
Visualizamos la situación de magnetización puramente inducida:
Por último, el caso de la reducción al polo. Aquí las inclinaciones de los campos involucrados son $90^\circ$ y las declinaciones $0^\circ$.
Consideramos la situación donde el prisma no tiene sus caras alineadas con los ejes de coordenadas horizontales. Realizamos los pasos indicados en el artículo científico. Consideramos el prisma rotado $55^\circ$ en sentido antihorario y el caso de magnetización puramente inducida.
Generamos un código para realizar la rotación de coordenadas, luego calculamos la TFA y graficamos.
Eso es todo por hoy.
Para vectorizar coordenadas aisladas de puntos sobre una cuadrícula. El resultado son dos listas con las coordenadas en x
y en y en la cuadrícula.
x = [1,2,3]
y = [1,2,3]
X = np.repeat(x,len(x))
Y = np.tile(x,len(x))
print(X)
print(Y)
[1 1 1 2 2 2 3 3 3]
[1 2 3 1 2 3 1 2 3]
arbitrary polarization, Geophysics, 29(4), 517, doi: 10.1190/1.1439386.