CN114167494A - 一种基于集合经验模态分解的海底地震波降噪方法 - Google Patents

Abstract

本发明提供一种基于集合经验模态分解的海底地震波降噪方法，包括如下步骤：步骤1：建立海底地震波褶积模型；步骤2：根据步骤1的模型设定仿真参数，得到海底地震波仿真数据，采用集合经验模态分解得到海底地震波信号分量，并构建低通滤波算法；步骤3：建立最优光滑降噪模型，根据最优光滑降噪模型设定评价函数对步骤2构建的算法进行评价；步骤4：选取最优算法实现海底地震波降噪。本发明建立褶积模型用来仿真海底地震波信号，在使用集合模态经验分解的方法处理信号后，建立最优光滑降噪模型来选取最优算法，得到高信噪比的海底地震波信号。

Description

一种基于集合经验模态分解的海底地震波降噪方法

技术领域

本发明涉及海底地震波信号降噪领域，具体涉及一种基于集合经验模态分解的海底地震波降噪方法。

背景技术

海底地震波是由舰船航行、海洋风浪、地壳运动等因素激发产生，可以在海底固液分界面远距离传播的一种振动波。通过海底地震仪来探测海底地震波是一种良好的获取海洋信息的方法，但是海底地震波信号受制于背景噪声和传播距离的影响，这就导致海底地震仪探测到的信号往往难以直接分辨，需要滤波降噪处理才能够有效提取海洋信息。

海底地震波信号往往是非线性非平稳信号，传统的线性滤波方法并不能很好的反应信号的非线性关系，滤波得到的信号也往往丢失关键信息。近年来有采用小波变换的方法来对海底地震波滤波，并且取得了不错的效果，如文献《一种海底地震波小波阈值降噪方法》(出处：中国科技信息,2017(17):88-89+14.)。但是由于小波变换需要自行选择基函数和降噪阈值，在自适应性方面就稍显不足，因此往往与其他方法结合使用，较为复杂，如专利《一种基于自适应滤波与小波变换结合的地震波降噪方法》(专利号CN201710334723.9)。上述的文献和专利虽然是有效的海底地震波降噪的方法，但是仍然存在一些问题，这是小波变换方法自身的缺陷导致的，难以完善。

综上所述，需要一种既可以实现海底地震波的有效信息提取，又兼具自适应性及简便性的方法来获取海底地震波信息。集合经验模态分解方法的出现给这种需求提供了一个新的思路，其根据信号自身尺度将信号分解成不同频率的分量，再通过算法可以实现各种水平的降噪效果。其算法简单、自适应性强，在处理非线性非平稳信号方面有巨大优势。采用集合经验模态分解对仿真的海底地震波信号进行滤波处理对于海底地震动数据分析、获取海洋信息有重要意义。

发明内容

本发明的目的是为了建立褶积模型用来仿真海底地震波信号，在使用集合模态经验分解的方法处理信号后，建立最优光滑降噪模型来选取最优算法，得到高信噪比的海底地震波信号的一种基于集合经验模态分解的海底地震波降噪方法。

本发明的目的是这样实现的：

一种基于集合经验模态分解的海底地震波降噪方法，包括如下步骤：

步骤1：建立海底地震波褶积模型；

步骤2：根据步骤1的模型设定仿真参数，得到海底地震波仿真数据，采用集合经验模态分解得到海底地震波信号分量，并构建低通滤波算法；

步骤3：建立最优光滑降噪模型，根据最优光滑降噪模型设定评价函数对步骤2构建的算法进行评价；

步骤4：选取最优算法实现海底地震波降噪。

进一步的，所述步骤1中所述的建立海底地震波褶积模型，其以陆地地震的褶积模型为参考，将其应用于海底情形。

将水下传感器采集到的海底地震波记录看成是由海底不同反射界面的反射系数与地震子波的褶积构成，在同一片海域其反射系数是相同的，所不同的是地震子波的频率。借此仿真海底地震波信号s(t)时将其简化为地震子波与反射系数序列褶积的结果：

s(t)＝ω(t)*r(t)

其中ω(t)称为地震子波，r(t)称为反射系数，*为褶积运算(又叫卷积)。舰船地震波频率主要集中在10-100Hz，地壳运动地震波频率一般不高于50Hz，仿真时采用f0＝50Hz的雷克子波与反射系数褶积构造地震波信号，再根据实际情况加背景噪声将其淹没。

进一步的，所述步骤2具体为：

步骤2.1：根据步骤1的模型设定仿真参数。

首先选定雷克子波为地震子波，其在时间域的表达式为：

r(t)＝A[1-2(πf₀t)²]exp[-(πf₀t)²]

式中，f₀为峰值频率，t为时间，A为幅值。取仿真参数A＝1，时间域采样间隔为0.001s，采样点数为100点，则总的时间长度为0.1s，采样率1000Hz，频率域采样间隔为10Hz。为了展示方便，将波形向右平移0.5s，即将式改写成：

r(t)＝A{1-2[πf₀(t-50)]²}exp{-[πf₀(t-50)]²}

随机生成500个范围在-1～1之间的数据作为反射系数，根据步骤1将反射系数与雷克子波褶积后再施加5dB的背景噪声。得到600个采样点的数据，取中间500个采样点数据进行步骤2.2。

步骤2.2：采用集合经验模态分解方法将步骤2.1仿真的海底地震波信号分解。

采用集合经验模态分解方法将步骤2.1得到的信号可以分解成8个分量。海底地震波信号可以表示为：

其中F_j(t)表示分解得到的有用信号分量，一共有7个，其频率按照从高到低排列，r₉(t)表示余项，也可表示为频率最低的第8个信号分量，记做F₉(t)。

步骤2.3：根据步骤2.2分解得到的信号分量构建低通滤波算法。

考虑到地震噪音多为高频，则需要滤除掉经过分解得到的高频信号分量，在步骤2.2中得到了8个信号分量，则具体的滤波器算法表达式可以表示为：

进一步的，所述步骤3建立的最优光滑降噪模型。其包括降噪偏差均方根误差(RMSE)和平均光滑度(SMSE)两个指标。

定义降噪偏差均方根误差：对于x1,x2,x3,…,xn一共n个数据，降噪偏差均方根误差可用下式表示：

式中，

为降噪后的信号，x_i为含噪信号。

定义平均光滑度：对于离散函数曲线y(x)，定义其在x0点的光滑程度(SN)为：

式中，h为步长，h＝1。点x0的光滑度越接近零，曲线的点x0的邻域越光滑，数据的波动越小，即滤波曲线越稳定。对于x1,x2,x3,…,xn一共n个数据，定义滤波曲线的所有点的光滑度的绝对值(两端去除2个点)的平均值称为滤波算法的平均光滑度，其满足：

降噪偏差均方根误差反映的是降噪的稳定性，平均光滑度反映的是信号的波动程度。当信号中含有大量噪声时，应当使这两个指标都小。最优光滑降噪函数F(n)表示为：

F(n)＝aRMSE+(1-a)SMSE

其中，a为加权系数，其在本次仿真中取a＝0.8。

进一步的，所述步骤4中选定最优算法为：

经过该算法处理得到的滤波后的地震波数据的信噪比可以达到12.5210dB。

与现有技术相比，本发明的有益效果是：

1.本发明建立了褶积模型用来模拟海底地震波。其以陆地地震的褶积模型为参考，将其应用于海底情形。

2.本发明建立了最优光滑降噪模型作为海底地震波滤波方法的评价函数，综合考虑了滤波的稳定性和信号的波动程度。其也可以作为其他滤波方法的评价参考，具有普适性。

3.本发明采用集合模态经验分解的方法对5dB信噪比的海底地震波仿真信号进行了滤波处理，采用的算法提高了信噪比，能够极大提高海底地震波的信息获取能力。

附图说明

图1为本发明中分析方法流程图；

图2为本发明中海底地震波信号仿真波形图；

图3为本发明中含噪海底地震波信号仿真波形图；

图4为本发明中经过集合模态分解后得到的8个信号分量示意图；

图5为本发明中光滑曲线示意图；

图6为本发明中评价函数图像；

图7为本发明中滤波前后对比效果图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述。

实施例1

一种基于集合经验模态分解的海底地震波降噪方法，所述分析方法包括以下步骤：

步骤1：建立海底地震波褶积模型；

步骤2：根据步骤1的模型设定仿真参数，得到海底地震波仿真数据，采用集合经验模态分解得到海底地震波信号分量，并构建低通滤波算法；

步骤3：建立最优光滑降噪模型，根据最优光滑降噪模型设定评价函数对步骤2构建的算法进行评价；

步骤4：选取最优算法实现海底地震波降噪。

进一步的，所述步骤1具体为，建立海底地震波褶积模型。将水下传感器采集到的海底地震波记录看成是由海底不同反射界面的反射系数与地震子波的褶积构成，在同一片海域其反射系数是相同的，所不同的是地震子波的频率。借此建立海底地震波褶积模型，仿真海底地震波信号s(t)时将其简化为地震子波与反射系数序列褶积的结果，如下式(1)：

s(t)＝ω(t)*r(t)(1)

其中ω(t)称为地震子波，r(t)称为反射系数，*为褶积运算(又叫卷积)。舰船地震波频率主要集中在10-100Hz，地壳运动地震波频率一般不高于50Hz，仿真时采用f0＝50Hz的雷克子波与反射系数褶积构造地震波信号，再根据实际情况加背景噪声将其淹没。

进一步的，所述步骤2具体包括以下步骤：

步骤2.1：根据步骤1的模型设定仿真参数。

首先选定雷克子波为地震子波，其在时间域的表达式为如下式(2)：

r(t)＝A[1-2(πf₀t)²]exp[-(πf₀t)²] (2)

式中，f₀为峰值频率，t为时间，A为幅值。取仿真参数A＝1，时间域采样间隔为0.001s，采样点数为100点，则总的时间长度为0.1s，采样率1000Hz，频率域采样间隔为10Hz。为了展示方便，将波形向右平移0.5s，即将式改写成如下式(3)：

r(t)＝A{1-2[πf₀(t-50)]²}exp{-[πf₀(t-50)]²} (3)

随机生成500个范围在-1～1之间的数据作为反射系数，根据步骤1将反射系数与雷克子波褶积后再施加5dB的背景噪声。得到600个采样点的数据，取中间500个采样点数据进行步骤2.2。

步骤2.2：采用集合经验模态分解方法将步骤2.1仿真的海底地震波信号分解。

采用集合经验模态分解方法将步骤2.1得到的信号可以分解成8个分量。海底地震波信号可以表示为如下式(4)：

其中F_j(t)表示分解得到的有用信号分量，一共有7个，其频率按照从高到低排列，r₉(t)表示余项，也可表示为频率最低的第8个信号分量，记做F₉(t)。

步骤2.3：根据步骤2.2分解得到的信号分量构建低通滤波算法。

考虑到地震噪音多为高频，则需要滤除掉经过分解得到的高频信号分量，在步骤2.2中得到了8个信号分量，则具体的滤波器算法表达式可以表示为如下式(5)：

进一步的，所述步骤3具体为，建立最优光滑降噪模型。其包括降噪偏差均方根误差(RMSE)和平均光滑度(SMSE)两个指标。

定义降噪偏差均方根误差：对于x1,x2,x3,…,xn一共n个数据，降噪偏差均方根误差可用下式(6)表示：

式中，

为降噪后的信号，x_i为含噪信号。

定义平均光滑度：对于离散函数曲线y(x)，定义其在x0点的光滑程度(SN)为：

式中，h为步长，h＝1。点x0的光滑度越接近零，曲线的点x0的邻域越光滑，数据的波动越小，即滤波曲线越稳定。对于x1,x2,x3,…,xn一共n个数据，定义滤波曲线的所有点的光滑度的绝对值(两端去除2个点)的平均值称为滤波算法的平均光滑度，其满足下式(8)：

降噪偏差均方根误差反映的是降噪的稳定性，平均光滑度反映的是信号的波动程度。当信号中含有大量噪声时，应当使这两个指标都小。最优光滑降噪函数F(n)表示为如下式(9)：

F(n)＝aRMSE+(1-a)SMSE (9)

其中，a为加权系数，其在本次仿真中取a＝0.8，对于给定的权重a，应当使F(n)越小越好。

进一步的，所述步骤4具体为，选定最优算法为下式(10)：

经过该算法处理得到的滤波后的地震波数据的信噪比可以达到12.5210dB。

实施例2

地震的褶积模型指的是在地质勘探领域，地震波通过大地***之后的***响应模型。此时将大地视为一个线性***，震源视为围绕一个固定频率，并且包含一定频率范围的脉冲信号，该脉冲信号会经过地层的调制作用成为探测到的地震信号。

对于褶积模型，如果以脉冲δ(t)的形式激发地震波，理想情况下在传播过程中脉冲δ(t)不会发生变化。但是实际情况中，由于传播过程的复杂性以及噪声的影响，脉冲δ(t)波形会受到调制作用，不再是理想的波形，探测到的地震记录s(t)实际为：

其中，

表示有效波，n(t)表示干扰波，其包含背景噪声和相干噪声两大类，ω(t)称为地震子波，r(t)称为反射系数，

通常将地震记录简化为地震子波与反射系数序列褶积的结果，如下式(2)所示：

s(t)＝ω(t)*r(t) (2)

作为瞬态震源的地震子波包括雷克子波、高斯一阶子波、余弦子波等，震源函数要求在频谱中高于中心频率的部分要迅速衰减，最好选择一个时间脉冲的子波函数，雷克子波在适用性方面优先选用，往往用于建立正演模型、制作合成地震道记录。同样的，在将合成记录与实际处理结果进行对比时，往往默认处理结果的子波是雷克子波。

Ricker子波在时间域的表达式为式(3)：

r(t)＝A[1-2(πf₀t)²]exp[-(πf₀t)²] (3)

式中，f₀为峰值频率，t为时间，A为幅值。

同样的，水下传感器采集到的海底地震波记录也可以看成是由海底不同反射界面的反射系数与地震子波的褶积构成，在同一片海域其反射系数是相同的，所不同的是地震子波的频率。舰船地震波频率主要集中在10-100Hz，地壳运动地震波频率一般不高于50Hz，在此，仿真时先采用取f₀＝50Hz的雷克子波与反射系数褶积构造地震波信号，再加背景噪声将其淹没。对于雷克子波，仿真参数A＝1，时间域采样间隔为0.001s，采样点数为100点，则总的时间长度为0.1s，采样率1000Hz，频率域采样间隔为10Hz。为了展示方便，将波形向右平移0.5s，即将式(3)改写成式(4)的形式：

r(t)＝A{1-2[πf₀(t-50)]²}exp{-[πf₀(t-50)]²} (4)

用随机函数产生的反射系数，取500个采样点。将两者卷积，得到海底地震波仿真数据，其波形图如图2所示。

背景噪声用高斯白噪声代替，其大小用信噪比(SNR)来表示，SNR指的是信号与噪声功率的比值，SNR越大，信号中噪声的含量越小，信号越纯净。满足公式(5)：

SNR＝10log₁₀(p_s/p_z) (5)

式中，p_s为不含噪信号功率，p_z为噪声功率，它们满足：

式(6)和式(7)中，f(n)为不含噪信号，f1(n)为含噪信号。施加5dB高斯白噪声，得到含噪的海底地震波信号仿真波形，如图所示3。

对于集合经验模态分解降噪方法，添加辅助白噪声的幅值标准差为信号标准差的0.2倍，加入新的正态分布白噪声序列的次数为100次。将含噪声的海底地震波信号分解成8个分量，如图4所示。海底地震波信号可以表示为如下式(8)：

其中F_j(t)表示分解得到的有用信号分量，一共有7个，其频率按照从高到低排列，r₉(t)表示余项，也可表示为频率最低的第8个信号分量，记做F₉(t)。

考虑到地震噪音多为高频，则需要滤除掉经过分解得到的高频信号分量，设计如下式(9)低通滤波器：

其中b为滤波器的截断参数，其取值应根据不同信号的滤波要求来优化确定。对于得到的8个信号分量，则具体的滤波器算法表达式可以表示为如下式(10)：

对于滤波器算法，构建最优光滑降噪模型来充当评价函数。其包括降噪偏差均方根误差(RMSE)和平均光滑度(SMSE)两个指标。

定义降噪偏差均方根误差：对于x1,x2,x3,…,xn一共n个数据，满足：

式中，

为降噪后的信号，x_i为含噪信号。降噪偏差均方跟误差越小，降噪后的曲线越接近原始信号，即两个曲线相似度越高，其反映的是降噪的稳定性。但是当去噪很小或者根本没有去噪时，其值仍然很小，仅凭均方跟误差并不能完全反映出去噪的效果。

在此，定义另一个指标：平均光滑度。

对于一条曲线f(t)，若在区间(t0,t1)之间有一阶连续导数，则其图形为一条处处有切线的曲线，且切线随切点的移动而连续转动，这样的曲线称为光滑曲线。此时需要满足条件：该点的左右导数同时存在且相等。对于有一条由曲线P(x)和Q(x)(0≤x≤1)组成的新的曲线R(x)，连接点横坐标为x0，连接点满足P(x0)＝Q(x0)。设左侧为Q(x)，右侧为P(x)，如图5所示。当R′(x)连续，即两条曲线P(x)和Q(x)的连接处的一阶导数相等时，R(x)曲线光滑。

当曲线光滑时，定义曲率κ为：

对于组成的新曲线R(x)，在x0处的曲率

可以表示为

或

它们具有同样的意义。即：

当曲线R(x)光滑即P′(x₀)＝Q′(x₀)时

|P″(x₀)|＝|Q″(x₀)|＝|R″(x₀)| (14)

以上分析同样适用于离散函数。当函数为离散函数时，假设h为采样步长(h＝1)，根据差分方法求离散函数二阶导数定义：

可以得到：

可以用Q″(x₀)-P″(x₀)的分子部分来表示离散曲线在x0点的光滑程度(SN)。

对于同一离散函数曲线y(x)，定义其在x0点的光滑程度SN为：

点x0的光滑度越接近零，曲线的点x0的邻域越光滑，数据的波动越小，即滤波曲线越稳定。对于x1,x2,x3,…,xn一共n个数据，定义滤波曲线的所有点的光滑度的绝对值(两端去除2个点)的平均值称为滤波算法的平均光滑度：

平均光滑度SMSE的值越小说明信号的波动程度越小，但对于一个实际信号并不意味着平均光滑度越小越好，光滑度过小会使信号发生失真。

对于最优光滑降噪模型，应当使降噪偏差均方误差和平均光滑度都小，这样才能说明实现了较好的降噪效果。为使其考虑指标的地位相当，对RMSE和SMSE都标准化，给定加权系数a，构建最优光滑降噪函数F(n)可表示为：

F(n)＝aRMSE+(1-a)SMSE (21)

式中a表示权重。a越大，则表示评价函数更注重***的相似度，a越小，则表示评价函数更注重***的平均光滑度。在信号降噪中，对于给定的权重a，应当使F(n)越小越好。

由于分解得到的信号分量为由高频到低频依次排列的，对于构建的低通滤波器，减去信号分量的频率越高，其SMSE指标越小。当信号中的信号和噪声重叠频带越多，SMSE指标的作用越重要，仿真的信号加入的是高斯白噪声，SMSE指标的作用相对较弱，对于构建的8个滤波算法，更看重算法的RMSE指标，取a＝0.8。对RMSE和SMSE标准化后，得到评价函数如图6所示，其中虚线表示RMSE函数，点线表示SMSE函数，实线表示F(n)函数。由图6可知，从最优降噪光滑函数的角度F(n)出发，构建的算法A₂(t)函数值最小，其含义是指原始信号减去最高频信号分量，最优算法为：

该算法的滤波对比效果如下图7所示，图的左边为不含噪的海底地震波波形图，图的右边为采用最优算法降噪后的海底地震波波形图。对于添加的5dB的噪声的海底地震波信号，降噪后信噪比SNR＝12.5210dB，降噪后的波形和原始信号相似程度得到了提高。

综上所述，本发明通过建立海底地震波的褶积模型对海底地震波信号进行了数据仿真；采用集合经验模态分解方法对信号进行分解后，构建了相应的低通滤波器；通过建立的最优光滑模型来充当评价函数，选出了最优降噪算法，提高了信噪比，增强了海底地震波的信号获取能力。本发明验证了集合模态分解方法在海底地震波滤波方面的应用价值，建立的模型对于其他类型的仿真分析也有参考意义，具有普适性。

Claims (5)

1.一种基于集合经验模态分解的海底地震波降噪方法，其特征是：包括如下步骤：

步骤1：建立海底地震波褶积模型；

步骤2：根据步骤1的模型设定仿真参数，得到海底地震波仿真数据，采用集合经验模态分解得到海底地震波信号分量，并构建低通滤波算法；

步骤3：建立最优光滑降噪模型，根据最优光滑降噪模型设定评价函数对步骤2构建的算法进行评价；

步骤4：选取最优算法实现海底地震波降噪。

2.根据权利要求1所述的一种基于集合经验模态分解的海底地震波降噪方法，其特征是：所述步骤1中所述的建立海底地震波褶积模型，其以陆地地震的褶积模型为参考，将其应用于海底情形。

将水下传感器采集到的海底地震波记录看成是由海底不同反射界面的反射系数与地震子波的褶积构成，在同一片海域其反射系数是相同的，所不同的是地震子波的频率。借此仿真海底地震波信号s(t)时将其简化为地震子波与反射系数序列褶积的结果：

s(t)＝ω(t)*r(t)

其中ω(t)称为地震子波，r(t)称为反射系数，*为褶积运算(又叫卷积)。舰船地震波频率主要集中在10-100Hz，地壳运动地震波频率一般不高于50Hz，仿真时采用f0＝50Hz的雷克子波与反射系数褶积构造地震波信号，再根据实际情况加背景噪声将其淹没。

3.根据权利要求1所述的一种基于集合经验模态分解的海底地震波降噪方法，其特征是：所述步骤2具体为：

步骤2.1：根据步骤1的模型设定仿真参数。

首先选定雷克子波为地震子波，其在时间域的表达式为：

r(t)＝A[1-2(πf₀t)²]exp[-(πf₀t)²]

式中，f₀为峰值频率，t为时间，A为幅值。取仿真参数A＝1，时间域采样间隔为0.001s，采样点数为100点，则总的时间长度为0.1s，采样率1000Hz，频率域采样间隔为10Hz。为了展示方便，将波形向右平移0.5s，即将式改写成：

r(t)＝A{1-2[πf₀(t-50)]²}exp{-[πf₀(t-50)]²}

随机生成500个范围在-1～1之间的数据作为反射系数，根据步骤1将反射系数与雷克子波褶积后再施加5dB的背景噪声。得到600个采样点的数据，取中间500个采样点数据进行步骤2.2。

步骤2.2：采用集合经验模态分解方法将步骤2.1仿真的海底地震波信号分解。

采用集合经验模态分解方法将步骤2.1得到的信号可以分解成8个分量。海底地震波信号可以表示为：

其中F_j(t)表示分解得到的有用信号分量，一共有7个，其频率按照从高到低排列，r₉(t)表示余项，也可表示为频率最低的第8个信号分量，记做F₉(t)。

步骤2.3：根据步骤2.2分解得到的信号分量构建低通滤波算法。

考虑到地震噪音多为高频，则需要滤除掉经过分解得到的高频信号分量，在步骤2.2中得到了8个信号分量，则具体的滤波器算法表达式可以表示为：

。

4.根据权利要求1所述的一种基于集合经验模态分解的海底地震波降噪方法，其特征是：所述步骤3建立的最优光滑降噪模型。其包括降噪偏差均方根误差(RMSE)和平均光滑度(SMSE)两个指标。

定义降噪偏差均方根误差：对于x1,x2,x3,…,xn一共n个数据，降噪偏差均方根误差可用下式表示：

式中，

为降噪后的信号，x_i为含噪信号。

定义平均光滑度：对于离散函数曲线y(x)，定义其在x0点的光滑程度(SN)为：

式中，h为步长，h＝1。点x0的光滑度越接近零，曲线的点x0的邻域越光滑，数据的波动越小，即滤波曲线越稳定。对于x1,x2,x3,…,xn一共n个数据，定义滤波曲线的所有点的光滑度的绝对值(两端去除2个点)的平均值称为滤波算法的平均光滑度，其满足：

降噪偏差均方根误差反映的是降噪的稳定性，平均光滑度反映的是信号的波动程度。当信号中含有大量噪声时，应当使这两个指标都小。最优光滑降噪函数F(n)表示为：

F(n)＝aRMSE+(1-a)SMSE

其中，a为加权系数，其在本次仿真中取a＝0.8。

5.根据权利要求1所述的一种基于集合经验模态分解的海底地震波降噪方法，其特征在是：所述步骤4中选定最优算法为：

经过该算法处理得到的滤波后的地震波数据的信噪比可以达到12.5210dB。

Images