Numerical simulation of subduction-induced molten plume: Destruction of overriding plate and its dynamic topographic responses
-
摘要: 洋壳俯冲过程中温度、压力升高和密度差异，可导致俯冲板片熔融柱的快速上涌，并作用在上覆板块岩石圈地幔底部，从而导致岩石圈的破坏减薄以及地表形态的剧烈变化，该过程类似于地幔柱对岩石圈的破坏作用。目前，对于俯冲板片熔融柱的形成及其对岩石圈破坏程度的研究相对较少，特别是地表动力地形变化与深部岩石圈破坏作用之间的响应关系依然不清楚。为此，本文将利用I2VIS有限差分方法，基于质量守恒方程、动量守恒方程以及能量守恒方程，通过给定材料参数和一定边界条件，计算揭示俯冲洋壳在不同时间和不同深度下发生部分熔融并形成俯冲板片熔融柱的过程，从而模拟再现该熔融柱对上覆板块岩石圈的破坏过程，并进一步分析其导致的浅部地表地形变化响应。数值模拟结果显示，在大洋板片俯冲过程中，由俯冲的陆源沉积物以及洋壳形成的混合熔融柱垂向侵蚀岩石圈底部，造成岩石圈减薄。在熔融柱的横向侵蚀过程中，岩石圈地幔熔融范围增加，可达300 km。在地形变化方面，板块俯冲造成大陆前缘受挤压变形，引起构造变形，构造变形范围可达300 km。同时，与俯冲相关形成的熔融柱对岩石圈地幔底部的侵蚀作用逐渐增强，动力地形变化幅度增大，并持续抬升，最终可垂向抬升至4 km。动力地形的变化范围局限在300 km以内，这与岩石圈地幔的破坏范围保持一致。Abstract: In the process of oceanic crust subduction, with the increase in temperature and pressure and the difference in density, the subduction-induced molten plume will rise rapidly and act on the lithospheric mantle bottom of the overriding plate, which may lead to the decrease in the lithospheric damage and the drastic change of surface morphology. This process is similar to the destruction of the lithosphere caused by mantle plume. So far, there have been little studies on the formation of subduction-induced molten plumes and their damage to the lithosphere, especially on the responses of surface dynamic topographic changes to the deep destruction. Based on the conservation equations of matter, momentum and energy, the I2VIS finite difference method is adopted by the authors to calculate and reveal the partial melting of the subducted oceanic crust at different times and depths under given material parameters and boundary conditions. The process of forming a subduction-induced molten plume is obtained, and then the process of molten plume-lithosphere interactions is further simulated, and the response of shallow topographic changes are analyzed. The numerical simulation results show that in the process of oceanic plate subduction, the composite molten plumes formed by subducted terrigenous sediment and oceanic crust eroded the bottom of the lithosphere longitudinally, and resulted in lithospheric thinning. During the transverse erosion of the molten plumes, the melting range of the lithospheric mantle increases up to 300 km. In terms of geomorphic change, plate subduction results in compression deformation of the continental front, which may reach 300 km. At the same time, the erosion of the molten plumes associated with subduction to the bottom of the lithospheric mantle is gradually strengthened, and the dynamic topographic changes increased, while uplifting continued, and ultimately reached a figure of 4 km. The variation range of dynamic topography is limited to 300 km, which is consistent with the damage range of lithospheric mantle.
-
1. 研究背景
S变换[1]是由美国的地球物理学家Stockwell等[2-3]在短时傅里叶变换和小波变换的基础上提出的。与短时傅里叶变换和小波变换相比，S变换既解决了短时傅里叶变换时窗固定的问题，同时具有小波变换的多分辨率分析。由于S变换所具有的独特优势，其在信号分析和处理领域受到广泛关注。由于高斯时窗和频率的函数关系较为单一，降低了S变换在实际应用中的灵活性，因此很多学者从不同角度对S变换进行了推广和改进得到广义S变换[4-8]，使得广义S变换更加适用于地震资料处理、油气检测[9-11]等方面的信号处理。
广义S变换作为时频域的滤波工具，能够压制噪音从而提高地震资料信噪比[12]。赵淑红等[13]首先在频域和时域确定干扰波范围，再利用广义S变换处理VSP数据；李雪英等[14]考虑到地震信号有效频带随时间逐渐变窄，在时频域采用交互式的方法确定高频噪音切除边界；王云专等采用时频谱叠加的方法计算共中心点道集中的滤波器进而压制随机干扰[15]。陈学华等引入λ和p两个参数控制高斯窗函数的变化趋势，在应用过程中通过试验确定最优参数[5]。本文对陈学华[5]的广义S变换方法进行了改进，使其兼顾对高低频端信号的分析能力，能够更加适用于随机噪音压制；其次由于标准S逆变换存在能量泄露，通过修正S逆变换[16-19]，以消除滤波噪音。在压制噪音时，首先计算信号瞬时信噪比，根据瞬时信噪比的高低结合不同处理策略实现噪音压制。通过处理合成数据和实际资料，验证改进广义S正反变换在随机噪音压制方面的有效性。
2. 方法原理
2.1 改进广义S变换
信号u(t)的标准S正变换的表达式为
$$ \int S\left(\tau ,f\right)={\int }_{-\infty }^{\infty }u\left(t\right)w(\tau -t,f){\mathrm{e}}^{-i2{{\text π} }ft}\mathrm{d}t $$ (1) 式中，
$ S\left(\tau ,f\right) $ 是信号$ u\left(t\right) $ 的时频谱，$ w(t,f) $ 是高斯窗函数，其具体表达式为$$ w\left(t,f\right)=\frac{1}{\sigma \sqrt{2\mathrm{{\text π} }}}{\mathrm{e}}^{\frac{-{t}^{2}}{2{\sigma }^{2}}} $$ (2) 式中，
$ \sigma =\dfrac{1}{\left|f\right|} $ 是控制高斯窗宽度的尺度因子。广义S变换就是通过调节尺度因子大小控制对应的高斯窗函数的有效宽度。针对标准S变换高斯窗形态随频率变化趋势单一进而导致S变换时频谱聚焦度低的问题，陈学华等[5]在标准S变换的基础上引入了控制参数λ和p，目的为丰富尺度因子随频率的变化趋势，针对不同的资料提供更多的参数组合。其改进的尺度因子为
$$ \sigma =\frac{1}{\lambda {\left|f\right|}^{p}} $$ (3) 不同的控制参数λ和p影响尺度因子变化的快慢（图1）。不同的λ值对尺度因子的主要影响集中在低频端，即随着λ值的增大，同一频率对应的尺度因子幅值增大。不同的p值对尺度因子的主要影响在高频端，即随着p值的增大，同一频率对应的尺度因子幅值减小。为了进一步说明两个参数对时变信号的聚焦作用，选取理论信号进行效果分析。该信号是由两个不同的宽频带调频信号组成，即使得高斯窗尺度因子随频率的变化范围比较大，可以有效检验方法的适用性。
![]()
图 1 理想时频谱与不同控制参数的广义S变换的时频谱a. 理想时频；b. p=0.7，λ=1对应的时频谱；c. p=0.9, λ=1对应的时频谱；d. p=1，λ=1标准S变换时频谱；e. p=1.2，λ=1对应的时频谱；f. p=1，λ=1.2对应的时频谱；g. p=1，λ=0.9对应的时频谱；h. p=1，λ=0.6对应的时频谱。Figure 1. Ideal time-frequency spectrum and time-frequency spectrum of S transform using different control factorsa. ideal time-frequency spectrum, b. p=0.7, λ=1 corresponding to time-frequency spectrum, c. p=0.9, λ=1 corresponding to time-frequency spectrum, d. p=1, λ=1 corresponding to time-frequency spectrum, e. p=1.2, λ=1 corresponding to time-frequency spectrum, f. p=1, λ=1.2 corresponding to time-frequency spectrum, g. p=1, λ=0.9 corresponding to time-frequency spectrum, h. p=1, λ=0.6 corresponding to time-frequency spectrum.$$\begin{split} y=&{\sin}\left\{3*\left[140{{\text π} }t+40{\sin}\left(3{{\text π} }t\right)\right]\right\}+\\&{\sin}\left\{3*\left[140{{\text π} }t+40{\cos}\left(3{{\text π} }t+120\right)\right]\right\} \end{split}$$ (4) 采用不同控制参数λ和p对上述信号进行广义S变换，分析对应的时频谱聚焦程度。当p值由小变大时，信号时频谱在低频端与高频端的聚焦程度均响应明显。其中低频端的聚焦程度随p值增大而变好，高频端随p值增大而变差，甚至当p＞1时，高频端出现严重失焦（图1a-d）。当p值由小变大时，信号时频谱在低频端的聚焦程度无响应，而高频端的聚焦程度也变化不大（图1c、f、g、h）。广义S变换中参数λ和p只能取固定值，其单一性的变化趋势不能兼顾高低频端的信号聚焦程度。综合以上实验结果可知，p值对时频域聚焦程度影响更大。改进广义S变换主要是更改p值变化趋势，而固定λ=1。
p值由小变大时，低频端与高频端聚焦程度响应不同。高频端信号需要小的p值，而低频端信号需要大的p值。为实现广义S变换在时频域的宽频信号扫描，将参数p改为频率的线性函数，即
$$ p=p\left(f\right)=a+bf $$ (5) 其中，初始值a＞0，变化率b＜0。
在S变换计算过程中，计算的最小频率为0 Hz，最大频率为待分析信号的Nyquist频率，在这个频带范围内规定p值线性递减，同时调整p值的初始值a和变化率b，使p值能够自动匹配信号，增加信号时频分析的聚焦程度。
为了检验该方法的效果，以理论模拟信号为例进行时频分析对比。该信号是由两个频率变化趋势相反的线性调频信号和一个正弦调频信号组成，见公式（6）。图2是标准S变换、广义S变换的优选方案（p=0.8，λ=1）和改进广义S变换的时频谱图（p=0→1，λ=1）。由图2可知，标准S变换在频率方向出现了条带状能量团，模糊了真实时频谱。广义S变换改善了时频谱图，增强了时频谱的聚焦度，而改进的广义S变换其高频的线性调频信号时频谱聚焦度更高，正弦型调频信号能量分布更均匀，分析效果要优于广义S变换。
![]()
图 2 合成信号时频谱图a. 标准S变换，b. 广义S变换（p=0.8，λ=1），c. 改进的广义S变换（p=0→1，λ=1）。Figure 2. Synthetic signal time-frequency spectruma. the standard S transform, b. conventional generalized S transform (p=0.8, λ=1), c. improved generalized S transform (p=0→1, λ=1).$$\begin{split} y=&{\cos}\left(68{{\text π} }t-20{{\text π} }{t}^{2}\right)+{\cos}\left[2{{\text π} }{\sin}\left(5{{\text π} }t\right)+120{{\text π} }t\right]+\\&{\cos}(168{{\text π} }t+28{{\text π} }{t}^{2}) \end{split}$$ (6) 为了进一步检验改进的广义S变换在地震信号分析中的适用性，选取一道添加高斯白噪声的实际地震信号进行分析。在高频噪声存在的情况下，利用上述三种方法得到时频谱的能量分布（图3）。由图3中立体图与平面图对比可知，标准S变换对低频部分的信号分析结果聚焦度比较好，但其高频噪音的时频谱的能量被过分放大，如果用这种方法在时频域采用阈值滤波的方法进行去噪，显然会把部分信号当成噪音同时也无法把高频噪音去掉。广义S变换的时频谱改善了高频噪音部分的聚焦度，但在低频端能量团几乎连在一起，时间分辨率变差。改进的广义S变换和标准S变换在低频端的能量分布相近，在高频端使高频噪音的能量分散聚焦，而且各个时频点的能量比较稳定，没有被过分放大。可以看出改进的广义S变换对地震资料进行分析或者时频域滤波，效果显然要优。
![]()
图 3 含随机噪声的地震信号时频谱图a. 标准S变换，b. 广义S变换（p=0.8, λ=1），c. 改进的广义S变换（p=0→1, λ=1）。Figure 3. Time-frequency spectrum of seismic signal with random noisea. the standard S transform, b. conventional generalized S transform (p=0.8, λ=1), c. improved generalized S transform (p=0→1, λ=1) .2.2 改进S域时频滤波
由于广义S变换具有逆变换，因此，可以利用广义S正反变换实现时频滤波。其常规滤波思路是利用二维滤波函数F(τ，f )与S变换时频谱相乘，对乘积结果沿时间轴积分，再反傅氏变换回时间域。滤波公式为：
$$ {u}_{\mathrm{f}\mathrm{i}\mathrm{l}\mathrm{t}1}\left(t\right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{S}(\tau ,f){F}(\tau ,f){\mathrm{e}}^{i2\mathrm{{\text π} }ft}\mathrm{d}\tau \mathrm{d}f $$ (7) 式中，F(τ, f )是时频域滤波函数。通过计算发现当F(τ, f ) =1时，即不对时频谱作任何处理可以完全重构信号，但当F(τ, f ) ≠1时，实际滤波结果和理想滤波结果存在差异，无法保证时频滤波的准确性[10]，如图4所示，在时频域将合成信号为500～600 ms的时频谱滤掉，反变换回时间域就出现了能量泄露。
![]()
图 4 标准S逆变换时频滤波效果图a. 合成信号，b. 合成信号时频域滤波，c. 时间域滤波效果。Figure 4. Effect of time and frequency filtering of standard S inverse transforma. synthetic signal, b. time and frequency filtering of synthetic signal, c. filtering effect of time field.为了解决这个问题，需要修正变换公式得到精确滤波的改进S域时频滤波公式。
$$\begin{split} S\left(\tau ,f\right)=&{\int }_{-\infty }^{\infty }u\left(t\right)w\left(\tau -t,f\right){\mathrm{e}}^{-i2\mathrm{{\text π} }ft}\mathrm{d}t=\\&{\int }_{-\infty }^{\infty }u\left(t\right)\frac{\lambda {\left|f\right|}^{p}}{\sqrt{2\mathrm{{\text π} }}}{\mathrm{e}}^{\frac{-{\lambda }^{2}{f}^{2p}{(\tau -t)}^{2}}{2}}{\mathrm{e}}^{-i2\mathrm{{\text π} }ft}\mathrm{d}t\end{split}$$ (8) 式中，u(t)是待分析信号，w(τ-t, f)是窗函数。
根据高斯窗函数的表达式，引入一个新的变量x(τ, t)。
$$ x\left(\tau ,t\right)=u\left(t\right){\mathrm{e}}^{\frac{-{\lambda }^{2}{f}^{2p}{(\tau -t)}^{2}}{2}} $$ 对式（8）针对变量t做傅里叶变换：
$$ x\left(\tau ,f\right)={\int }_{-\infty }^{\infty }u\left(t\right){{\mathrm{e}}^{\frac{-{\lambda }^{2}{f}^{2p}{(\tau -t)}^{2}}{2}}\mathrm{e}}^{-i2\mathrm{{\text π} }ft}\mathrm{d}t $$ (9) 对比式（7）和式（9）可得
$$ x\left(\tau ,f\right)=\frac{\sqrt{2\mathrm{{\text π} }}}{\lambda {\left|f\right|}^{p}}S\left(\tau ,f\right) $$ (10) 对式（10）针对变量 f 作反傅里叶变换：
$$ x\left(\tau ,t\right)=\frac{\sqrt{2\mathrm{{\text π} }}}{\lambda }{\int }_{-\infty }^{\infty }\frac{S\left(\tau ,f\right)}{{\left|f\right|}^{p}}{\mathrm{e}}^{i2\mathrm{{\text π} }ft}\mathrm{d}f $$ (11) 由式（8）可知，当τ=t时，u(t)和x(τ, t)之间具有对应关系：
$$ u\left(t\right)=x(t,t) $$ (12) 因此，令τ=t代入式（11）可得：
$$ u\left(t\right)=x\left(t,t\right)=\frac{\sqrt{2\mathrm{{\text π} }}}{\lambda }{\int }_{-\infty }^{\infty }\frac{S\left(t,f\right)}{{\left|f\right|}^{p}}{\mathrm{e}}^{i2\mathrm{{\text π} }ft}\mathrm{d}f $$ (13) 对u(t)进一步展开：
$$ \begin{split}u\left(t\right)=&\frac{\sqrt{2\mathrm{{\text π} }}}{\lambda }{\int }_{-\infty }^{\infty }\frac{S\left(t,f\right)}{{\left|f\right|}^{p}}{\mathrm{e}}^{i2\mathrm{{\text π} }ft}\mathrm{d}f=\\&{\int }_{-\infty }^{\infty }u\left(t\right){\mathrm{e}}^{\frac{-{\lambda }^{2}{f}^{2p}{\left(\tau -t\right)}^{2}}{2}}{\mathrm{e}}^{-i2\mathrm{{\text π} }ft}\mathrm{d}t=\\&{\int }_{-\infty }^{\infty }u\left(\tau \right){\int }_{-\infty }^{\infty }{\mathrm{e}}^{\frac{-{\lambda }^{2}{f}^{2p}{\left(t-\tau \right)}^{2}}{2}}{\mathrm{e}}^{-i2\mathrm{{\text π} }f(t-\tau )}\mathrm{d}f\mathrm{d}\tau =\\&{\int }_{-\infty }^{\infty }u\left(\tau \right)m\left(t-\tau \right)\mathrm{d}\tau =u\left(t\right)\otimes m\left(t\right)\end{split} $$ (14) 式中，
$ \otimes $ 代表褶积，函数m(t)的表达式：$$ m\left(t\right)={\int }_{-\infty }^{\infty }{\mathrm{e}}^{\frac{-{\lambda }^{2}{f}^{2p}{t}^{2}}{2}}{\mathrm{e}}^{i2\mathrm{{\text π} }ft}\mathrm{d}f $$ (15) 由式（14）可知，计算得到的结果是原始信号和信号m(t)的褶积，因此在时频域滤波结果的基础上需要再做一步滤波修正处理，最终的滤波表达式为：
$$ {u}_{filt2}\left(t\right)=\frac{\sqrt{2\mathrm{{\text π} }}}{\lambda }{\int }_{-\infty }^{\infty }\frac{S(t,f)F(t,f)}{{\left|f\right|}^{p}}{\mathrm{e}}^{i2\mathrm{{\text π} }ft}\mathrm{d}f $$ (16) $$ {u}_{\mathrm{f}\mathrm{i}\mathrm{l}\mathrm{t}3}\left(t\right)={\rm{IFT}}\left(\frac{{U}_{\mathrm{f}\mathrm{i}\mathrm{l}\mathrm{t}2}\left(f\right)}{M\left(f\right)}\right) $$ (17) 式中，F(t, f)是时频域滤波函数，Ufilt2(f)是ufilt2(t)的频谱，M(f)是m(t)的频谱，IFT表示反傅里叶变换。式（16）和式（17）是最终的滤波公式。
为了检验对S逆变换的改进效果，分别用标准逆变换、改进逆变换将图4b中的时频谱逆变换，并计算与准确滤波结果的差值，如图5所示。对比结果显示，改进逆变换消除了滤波区信号的畸变，避免产生滤波噪音。
![]()
图 5 准确滤波结果和不同滤波方法计算的差值a. 准确滤波结果，b. 标准逆变换计算差值（红色）和改进逆变换（蓝色）计算差值。Figure 5. Accurate filtering result and difference calculated by different filtering methodsa. accurate filtering result, b. difference calculated by standard S inverse transform (red) and difference calculated by improved S inverse transform (blue).基于上述改进广义S变换和修正的S逆变换，可以将随机噪音压制概括为3步：① 对含噪信号进行改进广义S变换得到信号的S谱；② 在S谱的基础上利用谱比法[20]计算信号的瞬时信噪比；③ 根据瞬时信噪比的不同，采用指数衰减和阈值滤波两种策略压制随机噪音。具体处理步骤如下：
（1）计算信号的S变换时频谱，利用该时频谱计算瞬时信噪比和阈值。阈值的计算思路是：计算初至前两个频率段内噪音的时频谱的平均值，作为初始阈值，具体的频率段可以分成有效频带内和有效频带外两段，计算两个阈值。由于高频随机干扰的时频谱幅值的均值随时间变化不大，而相对浅层数据深层数据的高频成分衰减快，能量弱，高频噪音能量占优，因此在无法利用初至前记录计算阈值时也可利用深层数据较高频段的时频谱计算初始阈值。
$$ {E}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}\left(t,f\right)=\frac{1}{T\times F}\sum _{t={t}_{1}}^{{t}_{2}}\sum _{f={f}_{1}}^{{f}_{2}}\left|S(t,f)\right| $$ (18) 式中，Emean(t, f)是阈值，S(t, f)是时频谱，T是参与计算时间点个数，F是参与计算频率点个数，t1到t2是参与计算的时间段，f1到f2是参与计算的频率段。
（2）实际地震数据需要考虑地层的吸收衰减，且高频衰减更快，因此深部有效频带变窄，有效信号能量弱[21]，此时如果将深部地震记录按照阈值处理，会损失很多有效信号，因此为了保留一部分有效信号以供后续处理和进一步使用，此时对瞬时信噪比低的数据段采用指数衰减的方法处理，即对信噪比小于1的时间点，对该时间点大于有效频带高频端的S谱的幅值进行一个指数衰减的加权，即将S(t, f)乘以一个以f为变量的指数衰减函数H(f)，如式（19）：
$$ H\left(f\right)={\mathrm{e}}^{-\frac{f}{\mu }} $$ (19) 式中，μ是一个衰减系数，f是频率，μ越大对高频衰减越小，反之对高频衰减越大。
（3）当瞬时信噪比大于1时，有效信号能量占优，且随机噪音的时频谱分布在整个频带上，对该时间点所有S谱采用阈值滤波法，当某频率点的幅值小于r倍的Emean(t, f)的时候，定义该点的S谱为噪音，将其充零，当其幅值大于r倍的Emean(t, f)的时候，将其时频谱的幅值减去r倍的Emean(t, f)，并将处理之后的时频谱反变换回时间域。
在实际应用中，应根据实际地震资料的情况，进行反复试验，选取合适的μ和r值，既能压制随机噪音，又能保护高频端的有效信号，为数据拓频留有空间。
3. 合成信号试验
为了检验方法的有效性，首先利用零相位雷克子波模拟48道经地层吸收衰减的地震记录，加入随机白噪声后，利用改进广义S变换对地震记录压制随机噪音（图6）。抽取单道地震模拟数据对比，能够清晰地看出改进广义S变换压制随机噪声的效果，尤其深部低信噪比信号也得到有效的恢复（图7）。分析时频域中滤波效果，可看到有效信号时频谱得以保留，随机噪音时频谱得到有效压制（图8）。
![]()
图 6 含噪地震记录去噪效果图a. 48道模拟地震记录， b. 含白噪声的模拟地震记录， c. 去除噪声的模拟地震记录。Figure 6. Denoising effect for seismic records with noisea. 48 channels simulated seismic record, b. simulated seismic record with white noise, c. simulated seismic record after denoising.![]()
图 7 单道信号对比图a. 加白噪声后与未加噪声单道信号对比， b. 去噪后与未加噪声单道信号对比。Figure 7. Comparison of single channel signalsa. comparison of single channel signals with noise and those without noise, b. comparison of single channel signals without noise and those after noise reduction.![]()
图 8 时频域滤波效果对比a. 不含白噪声模拟地震道的时频谱，b. 含白噪声模拟地震道的时频谱，c. 去除噪声模拟地震道的时频谱。Figure 8. Comparison of filtering effects in time-frequency domaina. time and frequency spectrum of simulated seismic channel without noise, b. time and frequency spectrum of simulated seismic channel with white noise, c. time and frequency spectrum of simulated seismic channel after noise reduction.4. 实际资料处理
为了检验本文改进的广义S变换方法在实际资料中的应用效果，选择高分辨率海洋二维地震叠前剖面进行去噪处理，并将其与传统广义S变换滤波去噪法进行对比。本文首先通过单炮初至前的噪声信号确定去噪阈值，再利用改进的广义S变换进行时频滤波。图9为地震叠前单炮剖面。传统广义S变换参数经过对比后确定为（p=0.6, λ=1），改进后的广义S变换应用参数（p=0.6→1, λ=1）。图10是两种方法压制噪声效果图。从整体上可以看出，两种广义S变换时频域去噪均能达到一定效果，随机噪声得到合理去除。但改进的广义S变换在远偏移距处的弱信号去噪效果更好，深部反射同相轴连续性更优，有效信号的损失较小，对弱信号起到保护作用。
![]()
图 10 两种方法压制噪声效果对比a. 改进广义S变换压噪后剖面，b. 传统广义S变换压噪后剖面，c. 改进广义S变换去除的噪声，d. 传统广义S变换去除的噪声。Figure 10. Comparison of noise suppression effects between two methodsa. denoised profile by improved generalized S transform, b. denoised profile by traditional generalized S transform, c. noise removal by improved generalized S transform, d. noise removal by generalized S transform.5. 结论
本文针对广义S变换高斯窗参数λ和p固定、逆变换存在信号泄露的问题，提出了改进的思路和计算过程，并将其应用于随机噪音压制。通过对合成信号和实际资料的应用，得到以下认识：
（1）改进广义S变换时频谱聚焦度要优于广义S变换，而且参数p与频率f呈线性变化，以适应不同的应用目的；改进后的逆变换可以实现精确滤波，消除了滤波噪音。
（2）利用单炮初至前或深部较高频段时频谱计算阈值，根据每个时刻时频谱的幅值与阈值的差异压制噪音，同时考虑到地震信号的有效频带和瞬时信噪比随时间变化的情况，对瞬时信噪比较低的采用指数衰减的策略压制噪音，可以在保护有效信号的基础上实现去噪目的。
（3）高斯窗参数
$ p $ 的范围选择，应当根据高低频端的聚焦程度而改变。通过试验确定的参数$ p $ ，既能够有效压制噪音又能避免弱有效信号的损失，为后续处理打好基础。 -
![]()
图 1 初始模型设置
a.模型主要区域（4 000 km×670 km）的初始物质场和温度场以及边界条件，图中白色的线条是以400 ℃为梯度绘制的等温线，黄色箭头是大洋俯冲速率和大陆仰冲速率；b.模型中不同颜色代表不同的岩石组成，分别是：1.空气层，2.海水，3.沉积物，4.洋壳，5.上地壳，6.下地壳，7.岩石圈地幔，8.软流圈，9.薄弱带，10、11.熔融的沉积物，12.熔融的洋壳，13.熔融的上地壳，14.熔融的下地壳
Figure 1. Initial model configuration
a.Enlargement (4 000 km×670 km) of the numerical model shows composition field and boundary conditions. The isotherms (white lines) are plotted for each 400 °C increment. Yellow arrow represents the subduction rates of oceanic plate and obduction rates of continental plate，b.The colored grid for different rock types: 1—air; 2—water; 3 —sedimentary cover; 4—oceanic crust; 5—upper continental crust; 6—lower continental crust; 7—lithospheric mantle; 8—athenospheric mantle; 9—weak zone mantle; 10 and 11—partially molten sediment; 12—partially molten oceanic crust;13 and 14—partially molten continental crust(5and 6)
![]()
图 2 与大洋板块俯冲相关形成的熔融柱对岩石圈地幔破坏的物质场演化过程
黑色数字表示实验积累时间，单位是Ma；白色的实线是温度线，数字表示温度，单位 ℃
Figure 2. Material field evolution of mantle plume damage to lithospheric mantle associated with subduction of an oceanic plate
The black number indicates the accumulated time of the experiment in millions of years. The solid white line is the temperature line, the number indicates the temperature, unit ℃
![]()
图 3 不同的时间动力地形瞬时变化（a）及地形随时间演化过程（b）
时间与图2相对应。红色虚线表示在不同时间的地形剖面
Figure 3. a. Different temporal dynamic terrain instantaneous changes (in meters) corresponding to Fig. 2, b. Represents the evolution of terrain over time
The dotted red line represents the topographic profile at different times, which is consistent with a
![]()
图 5 深部俯冲动力学过程与浅部地表变化的响应
（红色虚线表示动力地形扩展范围，黑色虚线表示构造变形范围，虚线重合区带表示地形变化稳定区）
Figure 5. Responses of shallow topographic changes to deep subduction dynamics
(The red dashed line indicates the extension range of dynamic topography, the black dashed line indicates the range of deformation, the dotted line superposition area indicates the stable area of topographic change)
表 1 数值模拟采用的黏滞性流变参数（据文献[35-38]）
Table 1 Parameters of viscous flow in the numerical experiments (after references[35-38])
标号 流变性质 E/(K·J·mol−1) V/(J·MPa−1·mol−1) n AD/(MPa−n·s−1) η0a/(Pa·s) A* 空气/水 0 0 1.0 1.0×10−12 1×1018 B* 湿石英（强） 154 0 2.3 3.2×10−6 1.97×1019 C* 斜长An75（强） 238 0 3.2 3.3×10−6 4.80×1024 D* 斜长石 An75 238 0 3.2 3.3×10−4 4.80×1022 E* 无水橄榄岩 532 8 3.5 2.5×104 3.98×1016 F*b 湿橄榄岩 470 8 4.0 2.0×103 5.01×1020 G*b 长英质熔体 0 0 1.0 2.0×10−9 5.00×1014 H 铁镁质熔体 0 0 1.0 1.0×10−7 1.00×1013 注：:a η0表示为有效黏滞系数，计算公式为：η0=(1/AD)×106n;
b 熔融的长英质熔体表示的是熔融的沉积物和地壳。
下载: 导出CSV
表 2 数值模型中的主要材料参数
Table 2 Parameters of the materials in the numerical models
物质 状态 ρ0
/(kg·m−3)ρe
/(kg·m−3)Cp
/(J·kg−1·K−1)Ka
/(W·m−1·K−1)$T_{\simfont\text{固相}}^{\rm b} $/K $T_{\simfont\text{液相}}^{\rm b} $/K Hr
/(μW·m−3)α
/K−1β
/MPa黏滞性流变参数 塑性性质
Sin (φeff)空气 — 1 — 100 20 — — 0 0 0 A* 0 水 — 1 000 — 3 330 20 — — 0 0 0 A* 0 沉积物
(6 km)固态 2 700 — 1 000 K1 TS1 TL1 2 3×10−5 1×10−5 B* 0.15 熔融 2 500 G* 0.06 上地壳
(14 km)固态 2 700 — 1 000 K1 TS1 TL1 2 3×10−5 1×10−5 B* 0.15 熔融 2 500 G* 0.06 下地壳
(15 km)固态 3 000 — 1 000 K2 TS2 TL2 0.5 3×10−5 1×10−5 C* 0.15 熔融 2 500 G* 0.06 洋壳(8 km) 固态 3 000 3 800 1 000 K2 TS2 TL2 0.25 3×10−5 1×10−5 D* 0.15 熔融 2 900 H* 0.06 岩石圈—软流圈地幔 固态 3 300 — 1 000 K3 — — 0.022 3×10−5 1×10−5 E* 0.6 熔融 2 700 0.06 水化地幔 固态 3 200 — 1 000 K3 — — 0.022 3×10−5 1×10−5 F* 0.6 熔融 2 700 0.06 注：a. K1=[0.64+807/(TK+77)]exp(0.000 04P); K2=[1.18+474/(TK+77)]exp(0.000 04P); K3=[0.73+1 293/(TK+77)]exp(0.000 04P);
b. 当P＜1 200 MPa, TS1=889+17 900/(P+54)+20 200/(P+54)2; 当P＞1 200 MPa, TS1=831+0.06P. TL1=1 262+0.09P; 当P＜1 600 MPa, TS2=973–70 400/(P+354)+778×105/(P+354)2; 当P＞1 600 MPa, TS2=935+0.003 5P+0.000 006 2P2. TL2=1 423+0.105P
下载: 导出CSV
-
[1] Morgan W J. Convection plumes in the lower mantle [J]. Nature, 1971, 230(5288): 42-43. doi: 10.1038/230042a0
[2] Wilson J T. A possible origin of the Hawaiian islands [J]. Canadian Journal of Physics, 1963, 41(6): 863-870. doi: 10.1139/p63-094
[3] Burov E, Guillou-Frottier L, d'Acremont E, et al. Plume head-lithosphere interactions near intra-continental plate boundaries [J]. Tectonophysics, 2007, 434(1-4): 15-38. doi: 10.1016/j.tecto.2007.01.002
[4] Christensen U R, Harder H. 3-D convection with variable viscosity [J]. Geophysical Journal International, 1991, 104(1): 213-220. doi: 10.1111/j.1365-246X.1991.tb02505.x
[5] 李建康, 王登红. 地幔柱数值模拟研究进展[J]. 地质科技情报, 2005, 24(4):13-20. [LI Jiankang, WANG Denghong. Advances in the numerical simulation of the mantle plume [J]. Geological Science and Technology Information, 2005, 24(4): 13-20. doi: 10.3969/j.issn.1000-7849.2005.04.003 [6] 卢记仁. 峨眉地幔柱的动力学特征[J]. 地球学报, 1996, 17(4):424-438. [LU Jiren. Dynamic characteristics of EMEI mantle plume [J]. Acta Geoscientia Sinica, 1996, 17(4): 424-438. [7] 徐义刚. 地幔柱构造、大火成岩省及其地质效应[J]. 地学前缘, 2002, 9(4):341-353. [XU Yigang. Mantle plumes, large igneous provinces and their geologic consequences [J]. Earth Science Frontiers, 2002, 9(4): 341-353. doi: 10.3321/j.issn:1005-2321.2002.04.014 [8] Gorczyk W, Hobbs B, Gessner K, et al. Intracratonic geodynamics [J]. Gondwana Research, 2013, 24(3-4): 838-848. doi: 10.1016/j.gr.2013.01.006
[9] King S D, Ritsema J. African hot spot volcanism: small-scale convection in the upper mantle beneath cratons [J]. Science, 2000, 290(5494): 1137-1140. doi: 10.1126/science.290.5494.1137
[10] Liao J, Gerya T. Influence of lithospheric mantle stratification on craton extension: Insight from two-dimensional thermo-mechanical modeling [J]. Tectonophysics, 2014, 631: 50-64. doi: 10.1016/j.tecto.2014.01.020
[11] Yang S H, Li Z H, Gerya T, et al. Dynamics of terrane accretion during seaward continental drifting and oceanic subduction: Numerical modeling and implications for the Jurassic crustal growth of the Lhasa Terrane, Tibet [J]. Tectonophysics, 2018, 746: 212-228. doi: 10.1016/j.tecto.2017.07.018
[12] van Avendonk H J A, Holbrook W S, Lizarralde D, et al. Structure and serpentinization of the subducting Cocos plate offshore Nicaragua and Costa Rica [J]. Geochemistry, Geophysics, Geosystems, 2011, 12(6): Q06009.
[13] Contreras-Reyes E, Grevemeyer I, Watts A B, et al. Deep seismic structure of the Tonga subduction zone: Implications for mantle hydration, tectonic erosion, and arc magmatism [J]. Journal of Geophysical Research: Solid Earth, 2011, 116(B10): B10103. doi: 10.1029/2011JB008434
[14] Fumagalli P, Stixrude L, Poli S, et al. The 10Å phase: A high-pressure expandable sheet silicate stable during subduction of hydrated lithosphere [J]. Earth & Planetary Science Letters, 2001, 186(2): 125-141.
[15] Irifune T, Kubo N, Isshiki M, et al. Phase transformations in serpentine and transportation of water into the lower mantle [J]. Geophysical Research Letters, 1998, 25(2): 203-206. doi: 10.1029/97GL03572
[16] Ranero C R, Weinrebe W, Grevemeyer I, et al. Tectonic structure of the Middle America Pacific margin and incoming Cocos Plate from Costa Rica to Guatemala[C]//American Geophysical Union, Fall Meeting 2003. AGU, 2003.
[17] Sano A, Ohtani E, Kubo T, et al. Effect of water on garnet-perovskite phase transformation in MORB system[C]//American Geophysical Union, Fall Meeting 2004. AGU, 2004.
[18] Li Z H, Xu Z Q, Gerya T V. Flat versus steep subduction: Contrasting modes for the formation and exhumation of high- to ultrahigh-pressure rocks in continental collision zones [J]. Earth & Planetary Science Letters, 2011, 301(1-2): 65-77.
[19] Campbell I H, Griffiths R W. Implications of mantle plume structure for the evolution of flood basalts [J]. Earth & Planetary Science Letters, 1990, 99(1-2): 79-93.
[20] d’Acremont E, Leroy S, Burov E B. Numerical modelling of a mantle plume: the plume head-lithosphere interaction in the formation of an oceanic large igneous province [J]. Earth and Planetary Science Letters, 2003, 206(3-4): 379-396. doi: 10.1016/S0012-821X(02)01058-0
[21] Pekeris C L. Thermal convection in the interior of the earth [J]. Geophysical Journal, 1935, 3(8): 343-367.
[22] Flament N, Gurnis M, Müller R D. A review of observations and models of dynamic topography [J]. Lithosphere, 2013, 5(2): 189-210. doi: 10.1130/L245.1
[23] Duretz T, Gerya T V. Slab detachment during continental collision: Influence of crustal rheology and interaction with lithospheric delamination [J]. Tectonophysics, 2013, 602: 124-140. doi: 10.1016/j.tecto.2012.12.024
[24] Gerya T V, Yuen D A. Rayleigh-Taylor instabilities from hydration and melting propel ‘cold plumes’ at subduction zones [J]. Earth & Planetary Science Letters, 2003, 212(1-2): 47-62.
[25] Gerya T V, Yuen D A. Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems [J]. Physics of the Earth and Planetary Interiors, 2007, 163(1-4): 83-105. doi: 10.1016/j.pepi.2007.04.015
[26] Vogt K, Gerya T V, Castro A. Crustal growth at active continental margins: Numerical modeling [J]. Physics of the Earth and Planetary Interiors, 2012, 192-193: 1-20. doi: 10.1016/j.pepi.2011.12.003
[27] Gerya T. Future directions in subduction modeling [J]. Journal of Geodynamics, 2011, 52(5): 344-378. doi: 10.1016/j.jog.2011.06.005
[28] Toussaint G, Burov E, Jolivet L. Continental plate collision: Unstable vs. stable slab dynamics [J]. Geology, 2004, 32(1): 33-36. doi: 10.1130/G19883.1
[29] Gorczyk W, Willner A P, Gerya T V, et al. Physical controls of magmatic productivity at Pacific-type convergent margins: Numerical modelling [J]. Physics of the Earth and Planetary Interiors, 2007, 163(1-4): 209-232. doi: 10.1016/j.pepi.2007.05.010
[30] Gerya T V, Meilick F I. Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts [J]. Journal of Metamorphic Geology, 2011, 29(1): 7-31. doi: 10.1111/j.1525-1314.2010.00904.x
[31] Li Z H, Gerya TV. Polyphase formation and exhumation of high- to ultrahigh-pressure rocks in continental subduction zone: Numerical modeling and application to the Sulu ultrahigh-pressure terrane in eastern China [J]. Journal of Geophysical Research: Solid Earth, 2009, 114(B9): B09406.
[32] Liao J, Wang Q, Gerya T, et al. Modeling craton destruction by hydration-induced weakening of the upper mantle [J]. Journal of Geophysical Research: Solid Earth, 2017, 122(9): 7449-7466. doi: 10.1002/2017JB014157
[33] Manglik A, Singh R N. Rheological stratification of the Indian continental lithosphere: Role of diffusion creep [J]. Proceedings of the Indian Academy of Sciences - Earth and Planetary Sciences, 1999, 108(1): 15-21.
[34] Gerya T V, Stöckhert B, Perchuk A L. Exhumation of high-pressure metamorphic rocks in a subduction channel: A numerical simulation [J]. Tectonics, 2002, 21(6): 6-1-6-19.
[35] Kirby S H, Kronenberg A K. Rheology of the lithosphere: selected topics [J]. Reviews of Geophysics, 1987, 25(6): 1219-1244. doi: 10.1029/RG025i006p01219
[36] Kirby S H. Rheology of the lithosphere [J]. Reviews of Geophysics, 1983, 21(6): 1458-1487. doi: 10.1029/RG021i006p01458
[37] Ranalli G, Murphy D C. Rheological stratification of the lithosphere [J]. Tectonophysics, 1987, 132(4): 281-295. doi: 10.1016/0040-1951(87)90348-9
[38] Ranalli G. Rheology of the Earth[M]. 2nd ed. Netherlands: Springer, 1995.
[39] Burg J P, Gerya T V. The role of viscous heating in Barrovian metamorphism of collisional orogens: thermomechanical models and application to the Lepontine Dome in the Central Alps [J]. Journal of Metamorphic Geology, 2005, 23(2): 75-95. doi: 10.1111/j.1525-1314.2005.00563.x
[40] Li Z H, Gerya T V, Burg J P. Influence of tectonic overpressure on P-T paths of HP-UHP rocks in continental collision zones: thermomechanical modelling [J]. Journal of Metamorphic Geology, 2010, 28(3): 227-247. doi: 10.1111/j.1525-1314.2009.00864.x
[41] Li Z H, Xu Z Q, Gerya T, et al. Collision of continental corner from 3-D numerical modeling [J]. Earth and Planetary Science Letters, 2013, 380: 98-111. doi: 10.1016/j.jpgl.2013.08.034
[42] Turcotte B, Schubert J. Geodynamics [J]. Geological curtain: English version, 2002, 450(2): 136-136.
[43] Gerya T V, Yuen D A. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties [J]. Physics of the Earth and Planetary Interiors, 2003, 140(4): 293-318. doi: 10.1016/j.pepi.2003.09.006
[44] Huangfu P P, Wang Y J, Cawood P A, et al. Thermo-mechanical controls of flat subduction: insights from numerical modeling [J]. Gondwana Research, 2016, 40: 170-183. doi: 10.1016/j.gr.2016.08.012
[45] Schmeling H, Babeyko A Y, Ennsa A, et al. A benchmark comparison of spontaneous subduction models-Towards a free surface [J]. Physics of the Earth and Planetary Interiors, 2008, 171(1-4): 198-223. doi: 10.1016/j.pepi.2008.06.028
[46] Larsen T B, Yeun D A. Fast plumeheads: Temperature-dependent versus non-Newtonian rheology [J]. Geophysical Research Letters, 1997, 24(16): 1995-1998. doi: 10.1029/97GL01886
[47] Manga M, Stone H A, O'Connell R J. The interaction of plume heads with compositional discontinuities in the Earth's mantle [J]. Journal of Geophysical Research, 1993, 98(B11): 19979-19990. doi: 10.1029/93JB00441
[48] Van Keken P. Evolution of starting mantle plumes: A comparison between numerical and laboratory models [J]. Earth & Planetary Science Letters, 1997, 148(1-2): 1-11.
[49] Burov E, Guillou-Frottier L. The plume head-continental lithosphere interaction using a tectonically realistic formulation for the lithosphere [J]. Geophysical Journal International, 2005, 161(2): 469-490. doi: 10.1111/j.1365-246X.2005.02588.x
[50] 郭慧丽, 徐佩芬, 张福勤. 华北克拉通及东邻西太平洋活动大陆边缘地区的P波速度结构: 对岩石圈减薄动力学过程的探讨[J]. 地球物理学报, 2014, 57(7):2352-2361. [GUO Huili, XU Peifen, ZHANG Fuqin. P wave velocity structure of the North China Craton and West Pacific active continental margin: exploration for dynamic processes of lithospheric thinning [J]. Chinese Journal of Geophysics, 2014, 57(7): 2352-2361. doi: 10.6038/cjg20140729 [51] 李三忠, 索艳慧, 李玺瑶, 等. 西太平洋中生代板块俯冲过程与东亚洋陆过渡带构造-岩浆响应[J]. 科学通报, 2018, 63(16):1550-1593. [LI Sanzhong, SUO Yanhui, LI Xiyao, et al. Mesozoic Plate Subduction in West Pacific and Tectono-magmatic Response in the East Asian Ocean-Continent Connection Zone [J]. Chinese Science Bulletin, 2018, 63(16): 1550-1593. [52] Yu S Y, Li S Z, Zhang J X, et al. Multistage anatexis during tectonic evolution from oceanic subduction to continental collision: A review of the North Qaidam UHP Belt, NW China [J]. Earth-Science Reviews, 2019, 191: 190-211. doi: 10.1016/j.earscirev.2019.02.016
[53] Yu S Y, Zhang J X, Li S Z, et al. TTG-Adakitic-like (Tonalitic-Trondhjemitic) magmas resulting from partial melting of Metagabbro under high-pressure condition during continental collision in the North Qaidam UHP Terrane, Western China [J]. Tectonics, 2019, 38(3): 791-822. doi: 10.1029/2018TC005259
[54] 郑永飞, 陈仁旭, 徐峥, 等. 俯冲带中的水迁移[J]. 中国科学: 地球科学, 2016, 59(4):651-682. [ZHENG Yongfei, CHEN Renxu, XU Zheng, et al. The transport of water in subduction zones [J]. Science China Earth Sciences, 2016, 59(4): 651-682. [55] Liu S F, Nummedal D, Gurnis M. Dynamic versus flexural controls of Late Cretaceous Western Interior Basin, USA [J]. Earth & Planetary Science Letters, 2014, 389: 221-229.
[56] Liu L J, Spasojević S, Gurnis M. Reconstructing Farallon plate subduction Beneath North America back to the Late Cretaceous [J]. Science, 2008, 322(5903): 934-938. doi: 10.1126/science.1162921
[57] 刘少峰, 王成善. 构造古地理重建与动力地形[J]. 地学前缘, 2016, 23(6):61-79. [LIU Shaofeng, WANG Chengshan. Reconstruction of tectono-paleogeography and dynamic topography [J]. Earth Science Frontiers, 2016, 23(6): 61-79. [58] Burov E, Gerya T. Asymmetric three-dimensional topography over mantle plumes [J]. Nature, 2014, 513(7516): 85-89. doi: 10.1038/nature13703
-
期刊类型引用(3)
1. 刘景良，苏杰龙，戴逸宸，李宇祖，黄永，郑文婷. 基于滑动窗宽优化的LMSSGST识别非平稳信号瞬时频率. 振动与冲击. 2024(19): 183-193 . 
百度学术
2. 李雪英，王鑫. 二阶同步提取变换的沉积旋回界面定位与追踪识别. 黑龙江科技大学学报. 2023(04): 614-621 . 百度学术
3. 吴耀乐，邓敏. 基于残差卷积神经网络的随机噪声压制网络的海洋地震资料去噪技术研究. 中国石油和化工标准与质量. 2022(22): 175-177 . 百度学术
其他类型引用(4)
计量
- 文章访问数: 3299
- HTML全文浏览量: 571
- PDF下载量: 100
- 被引次数: 7
