OpenGeoRock
3.11

Análisis probabilístico

Válido para OGR Slip2D 0.1.235

El factor de seguridad se calcula a partir de parámetros que solo se conocen de forma aproximada. En un análisis probabilístico, algunos de ellos se describen como variables aleatorias y se muestrean muchas veces, con lo que se obtiene cómo varía el factor de seguridad: su media y su dispersión, la probabilidad de que resulte inferior a uno y un índice de fiabilidad. El análisis de sensibilidad, que hace variar los parámetros de uno en uno, se trata al final del capítulo.

3.11.1Variables aleatorias#

Una variable aleatoria identifica un dato de entrada del modelo (qué parámetro de qué objeto) y le asigna una distribución. Pueden hacerse aleatorios los datos de entrada siguientes:

  • cualquier parámetro del modelo de resistencia de un material: c′c' y ϕ′\phi', y los parámetros de todos los modelos no lineales (Sección 3.5);
  • el peso específico y el peso específico saturado de un material, así como su rur_u, su presión intersticial constante, su ϕb\phi^b y su valor de entrada de aire;
  • las propiedades hidráulicas de los materiales que las tienen: permeabilidad saturada, relación y ángulo de anisotropía, parámetros de las funciones no saturadas, contenidos de agua saturado y residual, y almacenamiento específico;
  • la magnitud de las cargas distribuidas y de las lineales;
  • los coeficientes sísmicos khk_h y kvk_v;
  • la posición del nivel freático, definida como un desplazamiento vertical de media nula que desplaza todos los vértices del nivel freático en la misma cantidad.

Cada variable se define con una media μ\mu, que parte del valor actual del parámetro, una desviación típica σ\sigma y un mínimo y un máximo relativos, es decir, distancias por debajo y por encima de la media. La variable se trunca entre los límites

a=μ−Δmin,b=μ+Δmaxa = \mu - \Delta_{min}, \qquad b = \mu + \Delta_{max}

que es la forma habitual de expresar una incertidumbre («10 kPa, más o menos 4»). Si la distancia mínima es menor que la media, además se garantiza que toda muestra sea físicamente posible: un ángulo de rozamiento truncado así nunca resulta negativo. Una variable nueva parte con σ\sigma igual al 10 % de su media y con un margen del 30 % de la media a cada lado, o de 0.1 y 0.3 en valor absoluto si la media es cero. Una variable cuyo margen es nulo, o una variable normal o lognormal sin desviación típica, no es aleatoria y conserva su media; las distribuciones normal, lognormal, beta y gamma requieren todas una desviación típica positiva.

En OGR Slip2D

Las variables aleatorias se definen en Statistics → Random Variables, donde se enumeran todos los parámetros del modelo que pueden muestrearse y se representa la densidad de la variable que se está editando. Cada muestra se escribe en una copia del proyecto, por lo que un cálculo estadístico nunca modifica el modelo. Cuando una definición ya no corresponde al modelo (por ejemplo, porque se ha eliminado después un material), el cálculo señala la variable afectada, en lugar de omitirla sin avisar.

Código: ogr_core/statistics/random_variables.py

3.11.2Distribuciones#

Se dispone de siete distribuciones. La Tabla 3.11.1 da la densidad de cada una y cómo se obtienen sus parámetros a partir de los datos de entrada; aa y bb son los límites de la Ec. (3.11.1).

Tabla 3.11.1 Las siete distribuciones de OGR Slip2D: densidad f(x)f(x) antes del truncamiento y parametrización de cada una a partir de la media μ\mu, la desviación típica σ\sigma y los límites a, ba,\ b.
DistribuciónDensidadParámetros
Normal1σ2π e−(x−μ)2/2σ2\dfrac{1}{\sigma\sqrt{2\pi}}\,e^{-(x-\mu)^2/2\sigma^2}μ, σ\mu,\ \sigma, los valores introducidos
Uniforme1b−a\dfrac{1}{b - a} en [a, b][a,\ b]solo los límites; no se usa σ\sigma
Triangular2(x−a)(b−a)(c−a)\dfrac{2(x - a)}{(b - a)(c - a)} para x<cx < c, 2(b−x)(b−a)(b−c)\dfrac{2(b - x)}{(b - a)(b - c)} para x≥cx \ge cvértice c=μc = \mu; no se usa σ\sigma
Beta(x−a)p−1(b−x)q−1B(p,q) (b−a)p+q−1\dfrac{(x - a)^{p-1}(b - x)^{q-1}}{B(p,q)\,(b - a)^{p+q-1}} en [a, b][a,\ b]m=μ−ab−am = \frac{\mu - a}{b - a},
v=(σb−a)2≤0.999 m(1−m)v = \left(\frac{\sigma}{b - a}\right)^2 \le 0.999\,m(1 - m),
k=m(1−m)v−1k = \frac{m(1 - m)}{v} - 1,
p=mkp = m k, q=(1−m)kq = (1 - m)k
Exponencialλ e−λx\lambda\,e^{-\lambda x}, x≥0x \ge 0λ=1/μ\lambda = 1/\mu; no se usa σ\sigma
Lognormal1x σln⁡2π e−(ln⁡x−μln⁡)2/2σln⁡2\dfrac{1}{x\,\sigma_{\ln}\sqrt{2\pi}}\,e^{-(\ln x - \mu_{\ln})^2/2\sigma_{\ln}^2}σln⁡2=ln⁡ ⁣(1+σ2/μ2)\sigma_{\ln}^2 = \ln\!\left(1 + \sigma^2/\mu^2\right), μln⁡=ln⁡μ−σln⁡2/2\mu_{\ln} = \ln\mu - \sigma_{\ln}^2/2
Gammaxk−1e−x/θΓ(k) θk\dfrac{x^{k-1}e^{-x/\theta}}{\Gamma(k)\,\theta^{k}}θ=σ2/μ\theta = \sigma^2/\mu, k=μ2/σ2k = \mu^2/\sigma^2

Las distribuciones lognormal y gamma se parametrizan con la media y la desviación típica de la propia variable, que son los datos que maneja el usuario, y no con las de su logaritmo ni con los parámetros de forma y de escala. La distribución beta se ajusta por momentos en [a, b][a,\ b], de modo que conserva μ\mu y σ\sigma siempre que esa desviación típica sea alcanzable en el intervalo (los parámetros de forma no bajan de 10−310^{-3}). En las distribuciones uniforme y triangular, la media solo sirve para situar los límites y, en la triangular, el vértice: si el intervalo es asimétrico, la media de las muestras no coincide con μ\mu. Con μ=10\mu = 10, Δmin=6\Delta_{min} = 6 y Δmax=10\Delta_{max} = 10, las muestras uniformes tienen una media de 12 y las triangulares de 11.33 (Figura 3.11.1).

Densidad de probabilidad de las siete distribuciones de OGR Slip2D para una variable con media 10, desviación típica 3 y mínimo y máximo relativos de 6 y 10, truncada al intervalo de 4 a 20 Normal 4 10 20 Uniforme 4 10 20 Triangular 4 10 20 Beta 4 10 20 Exponencial 4 10 20 Lognormal 4 10 20 Gamma 4 10 20 μ = 10 (a trazos) σ = 3 mín. relativo = 6 máx. relativo = 10 intervalo [4, 20] una escala vertical
Figura 3.11.1 Densidades de las siete distribuciones para una variable con μ=10\mu = 10, σ=3\sigma = 3 y mínimo y máximo relativos de 6 y 10 (es decir, truncada a [4, 20]), tal como las muestrea el programa: cada curva es la densidad del propio programa dividida por F(b)−F(a)F(b) - F(a) (Ec. (3.11.2)). Todos los paneles comparten la misma escala vertical. Cálculo realizado con OGR Slip2D 0.1.235.

3.11.3Truncamiento#

El truncamiento se impone invirtiendo la función de distribución acumulada, y no por rechazo de muestras. Un valor uniforme u∈[0,1)u \in [0, 1) se traslada al tramo de la distribución acumulada comprendido entre los límites y luego se invierte:

x=F−1 ⁣(F(a)+u [F(b)−F(a)])x = F^{-1}\!\Big(F(a) + u\,\big[F(b) - F(a)\big]\Big)

donde

F, F−1F,\ F^{-1}
función de distribución acumulada de la distribución sin truncar, y su inversa

Toda muestra queda dentro de [a, b][a,\ b] (además se recorta a ese intervalo para evitar errores de redondeo), la densidad dentro del intervalo es f(x)/[F(b)−F(a)]f(x)/[F(b) - F(a)] y, a diferencia del rechazo, la transformación no altera la estratificación del muestreo por hipercubo latino: cada estrato de uu sigue dando exactamente una muestra.

3.11.4Muestreo#

En el muestreo de Monte Carlo (Metropolis y Ulam, 1949) cada valor uniforme se extrae de forma independiente. El muestreo por hipercubo latino (McKay, Beckman y Conover, 1979) divide [0,1)[0, 1) en NN estratos iguales y extrae un valor en cada uno,

ui=i+UiN,i=0,1,…,N−1,Ui∼U[0,1)u_i = \dfrac{i + U_i}{N}, \qquad i = 0, 1, \ldots, N - 1, \qquad U_i \sim \mathcal{U}[0, 1)

y después los baraja de forma independiente para cada variable. De este modo, cada variable recorre toda su distribución de manera uniforme con NN muestras, mientras que el emparejamiento entre variables sigue siendo aleatorio (Figura 3.11.2). Una opción permite, en cambio, compartir una misma estratificación entre todas las variables: la muestra ii cae entonces en el mismo estrato de todas ellas, y todas se mueven juntas. Eso responde a otra pregunta (qué ocurre si todos los parámetros son desfavorables a la vez) y suele aumentar la dispersión del factor de seguridad; no tiene sentido en el muestreo de Monte Carlo, que no tiene estratos.

Los números aleatorios pueden generarse de forma pseudoaleatoria, a partir de una semilla, de modo que el mismo modelo con la misma semilla da las mismas muestras, o de forma aleatoria, con la semilla tomada del reloj, lo que permite comprobar que una conclusión no depende de una semilla concreta.

Diez muestras de dos variables uniformes independientes, extraídas con el generador de muestras de OGR Slip2D (semilla 10116) por muestreo de Monte Carlo y por hipercubo latino Monte Carlo 0 0 0.5 0.5 1 1 u1 u2 Hipercubo latino 0 0 0.5 0.5 1 1 u1 u2
Figura 3.11.2 Diez muestras de dos variables uniformes independientes, extraídas con el generador de muestras de OGR Slip2D (semilla 10116) por muestreo de Monte Carlo (izquierda) y por hipercubo latino (derecha). En el hipercubo latino, cada fila y cada columna de la cuadrícula de 10 × 10 contiene exactamente una muestra.

3.11.5Correlación entre cohesión y ángulo de rozamiento#

Una variable puede declararse correlacionada con otra mediante un coeficiente ρ\rho comprendido entre −1 y 1. El uso previsto es el de la cohesión y el ángulo de rozamiento de un material de Mohr-Coulomb, en el que lo habitual es un coeficiente negativo. OGR Slip2D impone la correlación después del muestreo, reordenando los valores según el procedimiento de Iman y Conover (1982): las dos variables conservan exactamente los valores muestreados y, por tanto, sus distribuciones; solo cambia el emparejamiento. Si rir_i es el rango del ii-ésimo valor de la variable de referencia entre sus NN valores,

zi=Φ−1 ⁣(ri−12N),ti=ρ zi+1−ρ2  εiz_i = \Phi^{-1}\!\left(\dfrac{r_i - \frac12}{N}\right), \qquad t_i = \rho\,z_i + \sqrt{1 - \rho^2}\;\varepsilon_i

donde

Φ−1\Phi^{-1}
inversa de la función de distribución normal estándar
εi\varepsilon_i
números aleatorios normales estándar independientes

y los valores ordenados de la variable correlacionada se reasignan de modo que sus rangos sigan los de tt. Trabajar con las puntuaciones normales de los rangos, en lugar de con los rangos mismos, es lo que acerca el coeficiente obtenido al solicitado, con el que solo coincide de forma aproximada. Las puntuaciones empleadas son Φ−1((r−12)/N)\Phi^{-1}((r - \frac12)/N), una variante de las puntuaciones de van der Waerden Φ−1(r/(N+1))\Phi^{-1}(r/(N + 1)) de Iman y Conover; ambas se diferencian solo en la posición de representación.

En OGR Slip2D

El análisis se configura en Analysis → Project Settings → Statistics, donde se fijan Probabilistic analysis, Analysis type (Global Minimum por defecto), Sampling method (Monte Carlo por defecto) y Number of samples (1000), y en Random Numbers, donde se fijan Generation (pseudoaleatoria por defecto), Seed (10116) y la estratificación compartida del hipercubo latino (desactivada). La misma semilla rige también las búsquedas que generan superficies al azar. La correlación de una variable se define junto a su distribución, en Statistics → Random Variables.

Código: ogr_core/statistics/distributions.py · ogr_core/project/settings.py

3.11.6Mínimo global y talud completo#

El tipo mínimo global, que es el que se aplica por defecto, ejecuta primero el análisis determinista para localizar la superficie crítica de cada método y después vuelve a evaluar esa superficie para cada muestra. El círculo (o la polilínea) es fijo, pero cada muestra determina sobre él su propia masa deslizante, con los ajustes de búsqueda del análisis determinista, de modo que cada muestra requiere un cálculo de equilibrio límite por método. Cada método conserva su propia superficie crítica.

El tipo talud completo repite la búsqueda entera para cada muestra, por lo que no se presupone la posición de la superficie crítica: el factor de seguridad de una muestra es el mínimo que encuentra su búsqueda. Esto supone una búsqueda completa por muestra. Además, el cálculo acumula los factores de todas las superficies que evalúa y da la superficie probabilística crítica, que es la de mayor probabilidad de rotura entre las evaluadas al menos cinco veces (si hay empate, la de menor índice de fiabilidad) y no tiene por qué coincidir con la superficie crítica determinista. Las superficies se identifican de una muestra a otra por su geometría, redondeada a 0.5 m.

En ambos tipos, una muestra cuya superficie no puede resolverse, o es inadmisible (Sección 3.3), carece de factor de seguridad: se excluye de los estadísticos y se contabiliza según la causa, y el cálculo avisa de la pérdida (en el mínimo global, cuando se pierde más del 20 % de las muestras).

En OGR Slip2D

Si hay una norma de diseño activa, cada muestra se prepara igual que en el análisis determinista: los valores muestreados se escriben en una copia del modelo y a esa copia se le aplican después los coeficientes parciales, de modo que los estadísticos describen el factor de sobredimensionamiento y su probabilidad de ser inferior a uno (Sección 3.12).

Código: ogr_core/statistics/probabilistic.py · ogr_slip2d/analysis_runner.py

3.11.7Probabilidad de rotura e índice de fiabilidad#

La probabilidad de rotura es la fracción de muestras válidas cuyo factor de seguridad es inferior a uno:

PF=n(F<1)NvP_F = \dfrac{n(F < 1)}{N_v}

donde

NvN_v
número de muestras con un factor de seguridad válido y admisible
n(F<1)n(F < 1)
número de ellas con FF estrictamente inferior a uno

El índice de fiabilidad mide la distancia entre el factor de seguridad medio y la rotura, expresada en desviaciones típicas (Cornell, 1969). OGR Slip2D lo da en dos formas:

β=μF−1σF\beta = \dfrac{\mu_F - 1}{\sigma_F}
βln⁡=μln⁡Fσln⁡F\beta_{\ln} = \dfrac{\mu_{\ln F}}{\sigma_{\ln F}}

donde

μF, σF\mu_F,\ \sigma_F
media y desviación típica de los factores de seguridad muestreados (la desviación típica, con Nv−1N_v - 1 en el denominador)
μln⁡F, σln⁡F\mu_{\ln F},\ \sigma_{\ln F}
los mismos estadísticos de ln⁡F\ln F, calculados con los factores de seguridad positivos

La forma lognormal aplica la misma definición a ln⁡F\ln F, para el que la rotura corresponde a ln⁡F=0\ln F = 0; suele describir mejor el factor de seguridad, ya que este no puede ser negativo. Un muestreo sin dispersión da β=±∞\beta = \pm\infty. La probabilidad de rotura, por sí misma, no presupone ninguna forma para la distribución de FF.

Histograma de 1000 factores de seguridad muestreados: las barras por debajo de F = 1 forman la probabilidad de rotura; se marcan la media y el índice de fiabilidad 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 0 40 80 120 muestras F F = 1 PF μ = 1.30 β·σ
Figura 3.11.3 Histograma esquemático de 1000 factores de seguridad, muestreados con el programa a partir de una distribución lognormal de media 1.30 y desviación típica 0.20 (semilla 10116). Las barras sombreadas (F<1F < 1) suman PF=0.047P_F = 0.047; con μF=1.30\mu_F = 1.30 y σF=0.198\sigma_F = 0.198 resultan β=1.51\beta = 1.51 y βln⁡=1.65\beta_{\ln} = 1.65, tal como los calcula el programa a partir de estas muestras.
Ventana Statistics de OGR Slip2D con el histograma de 1 000 factores de seguridad muestreados del talud de dos capas, el umbral de rotura, la media y la línea de resumen
Figura 3.11.4 La ventana Statistics tras una ejecución sobre el talud de dos capas del apartado de inicio rápido (ejemplo 4), con la cohesión de la arcilla (normal, media 12 kPa, desviación típica 3 kPa), su ángulo de rozamiento (24°, 2.4°) y el ángulo de rozamiento de la arena limosa (32°, 2°) como variables aleatorias, cada una truncada a tres desviaciones típicas, y 1 000 muestras de Monte Carlo sobre la superficie crítica determinista. Con Bishop: μF=1.305\mu_F = 1.305, σF=0.113\sigma_F = 0.113, PF=0.5 %P_F = 0.5\,\% (cinco muestras por debajo de 1), β=2.69\beta = 2.69 y βln⁡=2.99\beta_{\ln} = 2.99.

3.11.8Convergencia#

Las estimaciones se estabilizan a medida que aumenta el número de muestras. OGR Slip2D registra el factor de seguridad medio acumulado y la probabilidad de rotura acumulada a lo largo de la secuencia de muestras, en unos 50 puntos (cada N/50N/50 muestras y en la última). Con ellos se dibuja el gráfico de convergencia de un cálculo: si la media o la probabilidad siguen cambiando en la última muestra, hacen falta más muestras. Una probabilidad de rotura pequeña exige muchas muestras para poder estimarse, pues se trata de un recuento de sucesos raros.

Gráfico de convergencia de la ventana Statistics: factor de seguridad medio acumulado y probabilidad de rotura acumulada en función del número de muestras
Figura 3.11.5 Convergencia del cálculo de la Figura 3.11.4. La media acumulada se mantiene a menos de 0.003 de su valor final a partir de unas 300 muestras; la probabilidad de rotura, con cinco muestras que rompen de 1 000, sigue dando un salto con cada nueva rotura. Para estimar una probabilidad de este orden hacen falta muchas más muestras.

3.11.9Análisis de sensibilidad#

El análisis de sensibilidad no es probabilístico: la forma de las distribuciones no interviene y solo se usan el mínimo y el máximo de cada variable. Cada variable se hace variar, de una en una, a lo largo de su intervalo, dividido en nn subintervalos iguales,

xj=a+(b−a) jn,j=0,1,…,nx_j = a + (b - a)\,\dfrac{j}{n}, \qquad j = 0, 1, \ldots, n

y en cada valor se recalcula el factor de seguridad sobre la superficie crítica de cada método, con todos los demás parámetros en los valores del modelo determinista. Con el valor por defecto n=50n = 50 resultan 51 evaluaciones por variable y método. Las variables se ordenan por la amplitud de variación del factor de seguridad en su intervalo; para cada una, el programa da además el punto en que la curva corta un factor de seguridad objetivo, por interpolación lineal entre los dos valores que lo acotan, y el eje horizontal puede expresarse como porcentaje de cada intervalo, de modo que variables con unidades distintas compartan un mismo gráfico. Los valores para los que la superficie no tiene un factor de seguridad admisible no se representan en la curva.

En OGR Slip2D

Sensitivity analysis y el número de Sensitivity intervals (50 por defecto) se encuentran en la página Analysis → Project Settings → Statistics. El análisis recorre las mismas variables aleatorias que el probabilístico y, si hay una norma de diseño activa, a cada modelo se le aplican los coeficientes parciales igual que al determinista.

Código: ogr_core/statistics/sensitivity.py

3.11.10Referencias#