CN104459795A - 一种深度变密度的地壳伸展系数热校正重力异常反演方法 - Google Patents

Abstract

本发明公开了一种深度变密度的地壳伸展系数热校正重力异常反演方法。本发明在重力反演地壳伸展系数计算中，引入岩石圈热重力异常校正和深度约束策略，根据海底地形、沉积厚度、磁异常、海洋年龄等时线、地震反射和折射等多参数约束的迭代计算方法，从而对由于岩石圈热扰动会造成地壳密度的变化，对重力解释带来歪曲效应进行校正，为认识地壳伸展和破裂过程提供依据。

Description

一种深度变密度的地壳伸展系数热校正重力异常反演方法

技术领域：

本发明涉及地球科学领域，尤其是一种深度变密度的地壳伸展系数热校正重力异常反演方法。

背景技术：

名词解释：

伸展系数：是伸展型盆地的地壳伸展程度的系数。McKenzie(1978)提出了岩石圈均匀拉张模型，又称纯剪切模型，该模型描述了岩石圈对伸展作用的基本响应,定义伸展系数β＝l/b，式中，l为伸展后的岩石圈宽度，b为初始宽度，用以描述地壳的伸展减薄程度。在洋陆边界(COB)附近或者洋脊扩张轴附近，大陆岩石圈伸展在深度上的变化，并且伸展量随深度而增大。

热校正：大陆岩石圈在伸展减薄过程中，产生破裂并形成海洋岩石圈，会导致地温陡升，之后随着岩石圈的冷却，地温将向平衡热梯度释放，这种扰动地温会引起的岩石圈横向密度变化，产生较大的岩石圈热重力异常，因此岩石圈热重力异常校正对重力反演计算正确性影响很大，其热校正非常重要。

重力异常：由地球内部密度分布差异所产生的扰动质量源和大地水准面起伏联合影响而引起的。根据重力异常场可以研究地球的形状、地质构造和重力探矿等科学问题。

岩石圈热均衡时间t：根据采集的海洋磁异常资料，识别出磁条带，推测出海洋年龄等时线，或参考已有的有关海洋年龄等时线资料，计算获得岩石圈热均衡时间t。

Parker-Oldenburg界面反演迭代公式：1973年Parker在地球物理位场界面正演中引入快速傅里叶变换(FFT)(Parker，1973,1974)，Oldenburg根据Parker公式，提出一种频率域的密度界面迭代反演方法(Oldenburg，1974)。

Moho界面：(Mohorovicic Discontinuity)莫霍面，是地壳与地幔间的不连续面，即：地壳与地幔的分界面。它是1909年由前南斯拉夫学者莫霍洛维奇首先发现的。在莫霍面上，地震波的纵波和横波传播速度增加明显，弹性和密度随深度逐渐增加，地幔物质密度、硬度大于地壳。地震波在地球内部的传播速度一般随深度增加，但又不是匀速增加，在某些深度处发生突然变化。地震波在此突然加速或减速(纵波减速时，横波甚至会消失)。这种波速发生突然变化的面叫作不连续面。

在张裂大陆边缘岩石圈研究中，伸展(减薄)系数、Moho深度及地壳厚度的计算对洋陆过渡带(OCT)结构和位置、大陆边缘与海洋岩石圈形成演化起到关键作用，特别是，描述伸展型盆地的地壳伸展程度的伸展系数，其在空间上的变化特征有助于研究盆地岩石圈演化机制，其研究成为国际海洋地质学前沿，为海洋地质学界研究岩石圈伸展和破裂机制提供依据，也为区域油气勘探提供重要帮助。

重力异常作为一种天然地球物理场，在研究大范围、多尺度地壳及岩石圈结构分析中起到重要作用。然而，岩石圈热扰动会造成地壳密度的变化，对重力解释带来歪曲效应，应该进行校正。

目前还没未见对于深度变密度的地壳伸展系数热校正重力反演方法工作的报道。

发明内容：

本发明的目的是提供一种能对由于岩石圈热扰动会造成地壳密度的变化而对重力解释带来歪曲效应进行校正的深度变密度的地壳伸展系数热校正重力异常反演方法

本发明在重力反演地壳伸展系数计算中，引入岩石圈热重力异常校正和深度约束策略，根据海底地形、沉积厚度、磁异常、海洋年龄等时线、地震反射和折射等多参数约束的迭代计算方法，从而对由于岩石圈热扰动会造成地壳密度的变化，对重力解释带来歪曲效应进行校正，为认识地壳伸展和破裂过程提供依据。

本发明的深度变密度的地壳伸展系数热校正重力异常反演方法，其特征在于，包括以下步骤：

a、计算岩石圈热重力异常

根据采集的海洋磁异常资料，识别出磁条带，推测出海洋年龄等时线，或参考已有的有关海洋年龄等时线资料，获得岩石圈热均衡时间t：

根据实测海底热流资料，计算模型热扰动温度T_z：

$T_z=\frac{2T_m}{\mathrm{\π}}\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\frac{{(-1)}^{n+1}}{n}\left[\frac{\mathrm{\β}}{\mathrm{n\π}}\sin \left(\frac{\mathrm{n\π}}{\mathrm{\β}}\right)\right]\exp (-\frac{n^{2}t}{\mathrm{\τ}})\sin \left(\frac{\mathrm{n\πz}}{a}\right)$

式中，T_m是岩石圈基准温度＝1300℃，β是岩石圈伸展系数(β＝l/b)，其是将均衡岩石圈厚度l除以初始薄岩石圈厚度b，τ是岩石圈冷却热衰变常数，a是于岩石圈厚度；

根据热扰动引起密度差估算公式：

Δρ＝ραΔT

式中，α是热膨胀系数，ΔT为岩石圈温度异常数值，这里等于T_z的值，ρ岩石圈密度常数＝3300kgm^-3

计算岩石圈热重力异常g_t：

$\begin{array}{c}{g}_t=\frac{8\mathrm{G\α\Δ\ρa}{T}_m}{\mathrm{\π}}\\ \·\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{1}{(2m+1)}\left[\frac{\mathrm{\β}}{(2m+1)}\sin \frac{(2m+1)\mathrm{\π}}{\mathrm{\β}}\right]\\ \·\exp (-\frac{{(2m+1)}^{2}t}{\mathrm{\τ}})\end{array}$

式中，G＝6.67×10^-11m³kg^-1s^-2万有引力常数，a是岩石圈厚度，α是热膨胀系数＝3.28×10^-5℃^-1，ρ岩石圈密度常数＝3300kgm^-3，T_m是岩石圈基准温度＝1300℃，β是岩石圈伸展系数(β＝l/b)，其是将均衡岩石圈厚度l除以初始薄岩石圈厚度b，τ是岩石圈冷却温度衰减常数，t是岩石圈温度均衡时间(Ma)；

b、计算地幔剩余重力异常

根据采集的海水平的重力异常资料，获得自由空间重力g_faa(见图2)，海底地形b(见图3)计算出海底地形重力g_b，沉积厚度s计算出g_s

再根据如下公式计算地幔剩余重力异常

g_mra＝g_faa-g_b-g_s-g_t

式中，g_faa是自由空间重力异常，g_mra是由莫霍面深度所引起地幔剩余重力异常，g_b是由海底地形横向变化所引起重力异常，g_t是岩石圈热重力异常，g_s是由沉积厚度和密度变化所引起重力异常；

c、计算深度变密度的Moho起伏深度

常密度模型的三维重力异常界面反演Parker-Oldenburg迭代公式：

$\left\{\begin{array}{c}F\left[\mathrm{\Δh}(x,y)\right]=-\frac{F\left[g_{\mathrm{mra}}(x,y)\right]e^{\left|k\right|d_{\mathrm{ref}}}}{2\mathrm{\πG\Δ}{\mathrm{\ρ}}_0}-\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\frac{{\left|k\right|}^{n-1}}{n!}F\left[\mathrm{\Δh}{(x,y)}^n\right]\\ \mathrm{\Δ}{\mathrm{\ρ}}_0={\mathrm{\ρ}}_m-{\mathrm{\ρ}}_c\end{array}\right.$

式中，G＝6.67×10^-11m³kg^-1s^-2万有引力常数，g_mra是实测地幔剩余重力异常，F[ ]是傅里叶变换，k是傅里叶域波数，Δh是Moho起伏深度，d_ref是Moho参考平均深度，ρ_m是地幔密度，ρ_c是地壳密度，x和y分别是观测点的x向和y向的坐标；

假设壳幔界面的密度差是随深度指数变化

Δρ(z)＝Δρ₀e^-μz

深度变密度的Moho起伏迭代计算公式为：

$F\left[\mathrm{\Δh}(x,y)\right]=-\frac{F\left[g_{\mathrm{mra}}(x,y)\right]e^{\left|k\right|d_{\mathrm{rfe}}}}{2\mathrm{\πG\Δ}{\mathrm{\ρ}}_0}-\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\frac{{|k-\mathrm{\μ}|}^{n-1}}{n!}F\left[\mathrm{\Δh}{(x,y)}^n\right]$

式中，Δρ是随深度成指数变化的双层界面密度差，μ是衰减系数，z是计算点之下的界面深度；

d、计算伸展系数

根据实测海底地形资料、前面获得的Moho起伏深度以及、根据有关资料的推测地壳中火山增厚部分等，综合计算获得伸展系数β：

$\left\{\begin{array}{c}\mathrm{\β}=\frac{c{t}_{\mathrm{ref}}}{c{t}_{\mathrm{now}}-c{t}_{\mathrm{mag}}}\\ d=d_{\mathrm{ref}.}+\mathrm{\Δh}\\ \mathrm{ct}=d-b\end{array}\right.$

式中，d_ref.是Moho参考平均深度，Δh是Moho起伏深度，b海底地形，d是Moho绝对深度，ct_ref参考地壳厚度，ct_mag是火山增厚部分，ct_now是计算现今地壳厚度，等于ct地壳厚度的值。

由于岩石圈热扰动会造成地壳密度的变化，对重力解释带来歪曲效应，应该进行热重力异常校正。同时，考虑到壳幔界面在部分垂向变化比较大，应该使用深度变化的界面密度差。为此，本发明提出一种在重力反演地壳伸展系数计算中，引入岩石圈热重力异常校正和深度变化密度的计算策略，根据海底地形、沉积厚度、磁异常、海洋年龄等时线、地震反射和折射等资料的迭代计算方法。

本发明的方法具有以下特点：

(1)本发明根据海底地形、沉积厚度、重力异常、磁异常、海洋年龄等时线、地震反射和折射等资料，使用重力异常迭代反演方法计算伸展系数。

(2)在重力异常迭代反演中，引入岩石圈热重力异常校正，参考海洋年龄等时线等资料，减少由于岩石圈热扰动会造成地壳密度的变化对重力解释带来歪曲效应影响。

(3)Moho计算中考虑壳幔界面在部分垂向变化比较大情况时，使用深度变化的界面密度差，减少模型计算与实际情况误差。

(4)在地壳厚度计算中，考虑火山增厚的部分影响的计算。

使用迭代计算可以充分利用所有涉及的参数变化所带来的影响，因此可以大限度使用更多约束条件，提高对反演参数的地质合理性，从而对由于岩石圈热扰动会造成地壳密度的变化，对重力解释带来歪曲效应进行校正，为认识地壳伸展和破裂过程提供依据。

附图说明：

图1是本发明的深度变密度的地壳伸展系数热校正重力异常反演方法的流程示意图；

图2是海洋重力异常图；

图3是海底地形图；

图4是Moho起伏深度；

图5是按照常规方法没有热校正和变密度计算出的地壳厚度；

图6是按照本发明的方法经过热校正和变密度计算出的地壳厚度；

图7是按照本发明的方法经过热校正和变密度计算出的伸展系数。

具体实施方式：

以下实施例是对本发明的进一步说明，而不是对本发明的限制。

实施例1：

参照图1，本实施例的深度变密度的地壳伸展系数热校正重力异常反演方法，包括以下步骤：

a、计算岩石圈热重力异常

根据采集的海洋磁异常资料，识别出磁条带，推测出海洋年龄等时线，或参考已有的有关海洋年龄等时线资料，获得岩石圈热均衡时间t：

根据实测海底热流资料，计算模型热扰动温度T_z：

$T_z=\frac{2T_m}{\mathrm{\π}}\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\frac{{(-1)}^{n+1}}{n}\left[\frac{\mathrm{\β}}{\mathrm{n\π}}\sin \left(\frac{\mathrm{n\π}}{\mathrm{\β}}\right)\right]\exp (-\frac{n^{2}t}{\mathrm{\τ}})\sin \left(\frac{\mathrm{n\πz}}{a}\right)$

式中，T_m是岩石圈基准温度＝1300℃，β是岩石圈伸展系数(β＝l/b)，其是将均衡岩石圈厚度l除以初始薄岩石圈厚度b，τ是岩石圈冷却热衰变常数，a是于岩石圈厚度；

根据热扰动引起密度差估算公式：

Δρ＝ραΔT

式中，α是热膨胀系数，ΔT为岩石圈温度异常数值，这里等于T_z的值，ρ岩石圈密度常数＝3300kgm^-3

计算岩石圈热重力异常g_t：

$\begin{array}{c}{g}_t=\frac{8\mathrm{G\α\Δ\ρa}{T}_m}{\mathrm{\π}}\\ \·\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{1}{(2m+1)}\left[\frac{\mathrm{\β}}{(2m+1)}\sin \frac{(2m+1)\mathrm{\π}}{\mathrm{\β}}\right]\\ \·\exp (-\frac{{(2m+1)}^{2}t}{\mathrm{\τ}})\end{array}$

式中，G＝6.67×10^-11m³kg^-1s^-2万有引力常数，a是岩石圈厚度，α是热膨胀系数＝3.28×10^-5℃^-1，ρ岩石圈密度常数＝3300kgm^-3，T_m是岩石圈基准温度＝1300℃，β是岩石圈伸展系数(β＝l/b)，其是将均衡岩石圈厚度l除以初始薄岩石圈厚度b，τ是岩石圈冷却温度衰减常数，t是岩石圈温度均衡时间(Ma)；

b、计算地幔剩余重力异常

根据采集的海水平的重力异常资料，获得自由空间重力g_faa(见图2)，海底地形b(见图3)计算出海底地形重力g_b，沉积厚度s计算出g_s

再根据如下公式计算地幔剩余重力异常

g_mra＝g_faa-g_b-g_s-g_t

式中，g_faa是自由空间重力异常，g_mra是由莫霍面深度所引起地幔剩余重力异常，g_b是由海底地形横向变化所引起重力异常，g_t是岩石圈热重力异常，g_s是由沉积厚度和密度变化所引起重力异常；

c、计算深度变密度的Moho起伏深度

常密度模型的三维重力异常界面反演Parker-Oldenburg迭代公式：

$\left\{\begin{array}{c}F\left[\mathrm{\Δh}(x,y)\right]=-\frac{F\left[g_{\mathrm{mra}}(x,y)\right]e^{\left|k\right|d_{\mathrm{ref}}}}{2\mathrm{\πG\Δ}{\mathrm{\ρ}}_0}-\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\frac{{\left|k\right|}^{n-1}}{n!}F\left[\mathrm{\Δh}{(x,y)}^n\right]\\ \mathrm{\Δ}{\mathrm{\ρ}}_0={\mathrm{\ρ}}_m-{\mathrm{\ρ}}_c\end{array}\right.$

式中，G＝6.67×10^-11m³kg^-1s^-2万有引力常数，g_mra是实测地幔剩余重力异常，F[ ]是傅里叶变换，k是傅里叶域波数，Δh是Moho起伏深度，d_ref是Moho参考平均深度，ρ_m是地幔密度，ρ_c是地壳密度，x和y分别是观测点的x向和y向的坐标；

假设壳幔界面的密度差是随深度指数变化

Δρ(z)＝Δρ₀e^-μz

深度变密度的Moho起伏迭代计算公式为：

$F\left[\mathrm{\Δh}(x,y)\right]=-\frac{F\left[g_{\mathrm{mra}}(x,y)\right]e^{\left|k\right|d_{\mathrm{rfe}}}}{2\mathrm{\πG\Δ}{\mathrm{\ρ}}_0}-\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\frac{{|k-\mathrm{\μ}|}^{n-1}}{n!}F\left[\mathrm{\Δh}{(x,y)}^n\right]$

式中，Δρ是随深度成指数变化的双层界面密度差，μ是衰减系数，z是计算点之下的界面深度；

本步骤是根据三维重力异常界面反演Parker-Oldenburg迭代公式，假设壳幔界面的密度差是随深度指数变化，获得壳幔界面的密度差是随深度指数变化，使用B计算的地幔剩余重力异常g_mra，作为重力异常反演数据进行计算获得Moho起伏深度Δh。

d、计算伸展系数

根据实测海底地形资料、前面获得的Moho起伏深度Δh以及、根据有关资料的推测地壳中火山增厚部分等，综合计算获得伸展系数β：

$\left\{\begin{array}{c}\mathrm{\β}=\frac{c{t}_{\mathrm{ref}}}{c{t}_{\mathrm{now}}-c{t}_{\mathrm{mag}}}\\ d=d_{\mathrm{ref}.}+\mathrm{\Δh}\\ \mathrm{ct}=d-b\end{array}\right.$

式中，d_ref.是Moho参考平均深度，Δh是Moho起伏深度，b海底地形，d是Moho绝对深度，ct_ref参考地壳厚度，ct_mag是火山增厚部分，ct_now是计算现今地壳厚度，等于ct地壳厚度的值。

将图2和图3中的数据代入上述方法中的公式中进行计算，得到相应的结果，如图4、5、6和7所示，式中图4是Moho起伏深度，图5是按照常规方法没有热校正和变密度计算出的地壳厚度，图6是按照本发明的方法经过热校正和变密度计算出的地壳厚度，图7是按照本发明的方法经过热校正和变密度计算出的伸展系数。

大陆岩石圈在伸展减薄过程中，产生破裂并形成海洋岩石圈，会导致地温陡升，之后随着岩石圈的冷却，地温将向平衡热梯度释放，这种扰动地温会引起的岩石圈横向密度变化，产生较大的岩石圈热重力异常。考虑热校正计算Moho面对比与附近地震资料估算结果更接近，未经深度变密度和热校正的计算地壳厚度结果会比地震估算的厚，厚算了近几公里，从而引入干扰。由此可见，在计算地壳伸展系数中，可以考虑使用深度变密度的岩石圈热重力异常校正的重力反演计算比较真实反映实际情况。

Claims (1)

1.一种深度变密度的地壳伸展系数热校正重力异常反演方法，其特征在于，包括以下步骤：

a、计算岩石圈热重力异常

根据实测海底热流资料，计算模型热扰动温度T_z：

$T_z=\frac{2T_m}{\mathrm{\π}}\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\frac{{(-1)}^{n+1}}{n}\left[\frac{\mathrm{\β}}{\mathrm{n\π}}\sin \left(\frac{\mathrm{n\π}}{\mathrm{\β}}\right)\right]\exp (-\frac{n^{2}t}{\mathrm{\τ}})\sin \left(\frac{\mathrm{n\πz}}{a}\right)$

式中，T_m是岩石圈基准温度＝1300℃，β是岩石圈伸展系数(β＝l/b)，其是将均衡岩石圈厚度l除以初始薄岩石圈厚度b，τ是岩石圈冷却热衰变常数，a是于岩石圈厚度，t是岩石圈热均衡时间；

根据热扰动引起密度差估算公式：

Δρ＝ραΔT

式中，α是热膨胀系数，ΔT为岩石圈温度异常数值，这里等于T_z的值，ρ岩石圈密度常数＝3300kgm^-3

计算岩石圈热重力异常g_t：

$\begin{array}{c}{g}_t=\frac{8\mathrm{G\α\Δ\ρa}{T}_m}{\mathrm{\π}}\\ \·\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{1}{(2m+1)}[\frac{\mathrm{\β}}{(2m+1)\mathrm{\π}}\sin \frac{(2m+1)\mathrm{\π}}{\mathrm{\β}}\\ \·\exp (-\frac{{(2m+1)}^{2}t}{\mathrm{\τ}})\end{array},]$

式中，G＝6.67×10^-11m³kg^-1s^-2万有引力常数，a是岩石圈厚度，α是热膨胀系数＝3.28×10^-5℃^-1，ρ岩石圈密度常数＝3300kgm^-3，T_m是岩石圈基准温度＝1300℃，β是岩石圈伸展系数(β＝l/b)，其是将均衡岩石圈厚度l除以初始薄岩石圈厚度b，τ是岩石圈冷却温度衰减常数，t是岩石圈温度均衡时间(Ma)；

b、计算地幔剩余重力异常

根据采集的海水平的重力异常资料，获得自由空间重力g_faa，海底地形b计算出海底地形重力g_b，沉积厚度s计算出g_s

再根据如下公式计算地幔剩余重力异常

g_mra＝g_faa-g_b-g_s-g_t

式中，g_faa是自由空间重力异常，g_mra是由莫霍面深度所引起地幔剩余重力异常，g_b是由海底地形横向变化所引起重力异常，g_t是岩石圈热重力异常，g_s是由沉积厚度和密度变化所引起重力异常；

c、计算深度变密度的Moho起伏深度

常密度模型的三维重力异常界面反演Parker-Oldenburg迭代公式：

$\left\{\begin{array}{c}F\left[\mathrm{\Δh}(x,y)\right]=-\frac{F\left[g_{\mathrm{mra}}(x,y)\right]e^{\left|k\right|d_{\mathrm{ref}}}}{2\mathrm{\πG\Δ}{\mathrm{\ρ}}_0}-\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\frac{{\left|k\right|}^{n-1}}{n!}F\left[\mathrm{\Δh}{(x,y)}^n\right]\\ \mathrm{\Δ}{\mathrm{\ρ}}_0={\mathrm{\ρ}}_m-{\mathrm{\ρ}}_c\end{array}\right.$

式中，G＝6.67×10^-11m³kg^-1s^-2万有引力常数，g_mra是实测地幔剩余重力异常，F[]是傅里叶变换，k是傅里叶域波数，Δh是Moho起伏深度，d_ref是Moho参考平均深度，ρ_m是地幔密度，ρ_c是地壳密度，x和y分别是观测点的x向和y向的坐标；

假设壳幔界面的密度差是随深度指数变化

Δρ(z)＝Δρ₀e^-μz

深度变密度的Moho起伏迭代计算公式为：

$F\left[\mathrm{\Δh}(x,y)\right]=-\frac{F\left[g_{\mathrm{mra}}(x,y)\right]e^{\left|k\right|d_{\mathrm{ref}}}}{2\mathrm{\πG\Δ}{\mathrm{\ρ}}_0}-\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\frac{{|k-\mathrm{\μ}|}^{n-1}}{n!}F\left[\mathrm{\Δh}{(x,y)}^n\right]$

式中，Δρ是随深度成指数变化的双层界面密度差，μ是衰减系数，z是计算点之下的界面深度；

d、计算伸展系数

根据实测海底地形资料、前面获得的Moho起伏深度以及、根据有关资料的推测地壳中火山增厚部分等，综合计算获得伸展系数β：

$\left\{\begin{array}{c}\mathrm{\β}=\frac{{\mathrm{ct}}_{\mathrm{ref}}}{{\mathrm{ct}}_{\mathrm{now}}-{\mathrm{ct}}_{\mathrm{mag}}}\\ d=d_{\mathrm{ref}.}+\mathrm{\Δh}\\ \mathrm{ct}=d-b\end{array}\right.$

式中，d_ref.是Moho参考平均深度，Δh是Moho起伏深度，b海底地形，d是Moho绝对深度，ct_ref参考地壳厚度，ct_mag是火山增厚部分，ct_now是计算现今地壳厚度，等于ct地壳厚度的值。