一种匹配追踪地震谱分解方法及装置与流程

文档序号:12611751阅读:467来源:国知局
一种匹配追踪地震谱分解方法及装置与流程

本申请涉及地震资料处理技术领域,尤其是涉及一种匹配追踪地震谱分解方法及装置。



背景技术:

随着油气勘探开发的不断深入,小断层、薄储层和岩性闭圈的识别与解释成为地震资料精细解释的重点。地震谱分解技术不仅可以提高地震资料对薄储层的解释与预测能力,而且可从地震数据中获得更为丰富的地质信息,一经推出便引起业界广泛关注,并得到了快速发展。已在地层或沉积相解释、油气检测等方面取得广泛应用。其基本原理是选取一组在时间域和频率域均有限的基函数,通过观察待分析信号在基函数上的投影,得到信号在时间域和频率域的联合分布。

最初的地震谱分解主要通过短时傅里叶变换(STFT)来实现,后来又逐渐出现了基于小波变换、S变换以及匹配追踪算法(MP)等的地震谱分解方法。其中,基于匹配追踪算法(MP)的地震谱分解方法在信号时频分解时,具有简单、直接和有效的特点,在地震资料解释领域得到了广泛的运用。匹配追踪算法的出现使得利用过完备库来实现信号的稀疏表示在大多数情况下成为可能。但由于匹配追踪方法是一种贪婪算法,且过完备库一般非常庞大,寻找真正的最佳匹配在计算上是非常昂贵和不现实的。目前大多数情况下选择过完备库中相关性最大的当前原子时,采用遍历所有原子选择内积最大的那一个原子的做法,其缺点是不能适应不同特征的地震信号,而且计算效率和分解精度不高。



技术实现要素:

本申请实施例的目的在于提供一种匹配追踪地震谱分解方法及装置,可以提高匹配追踪算法用于地震谱分解的计算效率、分解精度以及适应性。

为达到上述目的,本申请实施例提供了一种匹配追踪地震谱分解方法,所述方法包括:

(1)获取地震信号,并将所述地震信号作为当前信号;

(2)根据所述当前信号,确定待搜索原子的中心频率;

(3)使用匹配追踪算法搜索预设过完备库中满足第一预设条件的原子,选择所述满足第一预设条件的原子中与所述当前信号最相关的原子,获取所述当前信号在所述最相关的原子处的对应投影分量,并获取所述当前信号与所述对应投影分量的信号残差;所述第一预设条件为所述原子的中心频率等于所述待搜索原子的中心频率;

(4)将所述信号残差作为新的当前信号,重复步骤(2)至(4),直至当前得到的所述信号残差小于预设阈值为止;

(5)根据所有所述对应投影分量,得到所述地震数据的地震谱分解结果。

为达上述目的,本申请实施例还提供了一种匹配追踪地震谱分解装置,所述装置包括:

获取模块,用于获取地震信号,并将所述地震信号作为当前信号;

中心频率确定模块,用于根据所述当前信号,确定待搜索原子的中心频率;

搜索模块,用于使用匹配追踪算法搜索预设过完备库中满足第一预设条件的原子,选择所述满足第一预设条件的原子中与所述当前信号最相关的原子,获取所述当前信号在所述最相关的原子处的对应投影分量,并获取所述当前信号与所述对应投影分量的信号残差;所述第一预设条件为所述原子的中心频率等于所述待搜索原子的中心频率;

重复模块,用于将所述信号残差作为新的当前信号,重复执行所述中心频率确定模块至所述重复模块,直至当前得到的所述信号残差小于预设阈值为止;

结果获得模块,用于根据所有所述对应投影分量,得到所述地震数据的地震谱分解结果。

由上述本申请实施例所提供的技术方案可知,本申请实施例在使用匹配追踪算法进行地震谱分解时,通过当前信号确定了过完备库中需要搜索的原子的中心频率,不同于现有技术的搜索所有原子的做法,提高了计算效率。另一方面,本申请实施例只搜索预设过完备库中满足第一预设条件的原子,从中确定用于表示地震信号的原子,最终得到的用于表达地震信号的过完备库与地震信号之间更加匹配,提高了匹配追踪算法地震谱分解的分解精度和对不同地震信号的适应性。

附图说明

此处所说明的附图用来提供对本申请实施例的进一步理解,构成本申请实施例的一部分,并不构成对本申请实施例的限定。在附图中:

图1为本申请实施例的一种匹配追踪地震谱分解方法示意图;

图2为本申请实施例的主频为30HZ的Ricker子波示意图;

图3为本申请实施例的小波变换时频分解的结果示意图;

图4为本申请实施例的S变换时频分解的结果示意图;

图5为本申请实施例的常规匹配追踪算法时频分解的结果示意图;

图6为本申请实施例的改进匹配追踪算法时频分解的结果示意图;

图7为本申请实施例的某一道地震道的地震信号示意图;

图8为本申请实施例的对图7所示的地震信号使用CWT分解的结果示意图;

图9为本申请实施例的对图7所示的地震信号使用常规匹配追踪分解的结果示意图;

图10为本申请实施例的对图7所示的地震信号使用改进后的匹配追踪分解的结果示意图;

图11为本申请实施例的改进匹配追踪算法产生的10HZ频率剖面;

图12为本申请实施例的改进匹配追踪算法产生的30HZ频率剖面;

图13为本申请实施例的改进匹配追踪算法产生的50HZ频率剖面;

图14为本申请实施例的改进匹配追踪算法分解获得的10Hz、30Hz和50Hz沿层切片多谱图像合成结果;

图15为本申请实施例的短时快速傅里叶变换(FFT)分解获得的10Hz、30Hz和50Hz沿层切片多谱图像合成结果;

图16为本申请实施例的单CPU情况下改进匹配追踪算法和短时窗FFT算法耗时对比图;

图17为本申请实施例的CPU多线程下改进匹配追踪算法和短时窗FFT算法耗时对比图;

图18为本申请实施例的一种匹配追踪地震谱分解装置示意图。

具体实施方式

为使本申请实施例的目的、技术方案和优点更加清楚明白,下面结合实施例和附图,对本申请实施例做进一步详细说明。在此,本申请实施例的示意性实施例及其说明用于解释本申请实施例,但并不作为对本申请实施例的限定。

为了更清楚的介绍本申请实施例所提供的技术方案,首先在这里对匹配追踪算法做一个介绍。匹配追踪算法是信号稀疏表达的一种方式,实质是将信号在过完备库上进行分解。

给定一个过完备库D∈Rn×k,其中它的每列表示一种原型信号的原子。给定一个信号y,它可以被表示成这些原子的稀疏线性组合。即,信号y可以表示成:y=Dx,或者y≈Dx。过完备库的过完备性指的是原子的个数远远大于信号y的长度。匹配追踪算法的基本思路为:从过完备库D中,选择一个与信号y最匹配的原子构成一个稀疏逼近,并求出信号残差,然后继续选择与信号残差最匹配的原子,反复迭代,信号y可以由这些原子来线性和,再加上最后的残差值来表示。很显然,如果残差值在可以忽略的范围内,则信号y就近似是这些原子的线性组合。

上述过程中,在选择与信号y最匹配的原子时,可以计算信号y与过完备库中每一个原子的内积,选择内积绝对值最大的一个原子,它就是与信号y在本次迭代运算中最匹配的。该过程用公式表达可以为:令信号y∈H(H表示希尔伯特空间),从过完备库中选择一个最匹配的原子,满足其中,r0表示一个过完备库矩阵的列索引。这样信号y就被分解为在最匹配原子的垂直投影分量和残差两个部分,即:

对残差R1f进行同样的分解,并依次对后续每一次得到的残差进行同样的分解,那么第k步可以得到:其中满足可见信号y经过K+1步分解之后为:

上述Rif表示残差,i=1,2,3…

上述就是常规匹配追踪算法的简单介绍,下面结合附图,对本申请实施例的具体实施方式作进一步的详细说明。

参考图1所示,本申请实施例提供的一种匹配追踪地震谱分解方法,可以包括以下步骤。该方法可以用于地震信号的时频分解。

步骤S101,获取地震信号,并将所述地震信号作为当前信号。

其中,所述地震信号可以为地震数据采集得到的原始时域地震信号。为了使得分频的结果更加平滑,在本申请的一个实施例中,所述地震信号可以为经过平滑处理的地震信号,这样就能减轻噪声对分解结果的影响。

步骤S102,根据所述当前信号,确定待搜索原子的中心频率。

由于本实施例中所介绍的改进匹配追踪方法是一个循环执行的过程,每一次执行循环体时都需要一个当前信号。具体的,所述当前信号可以为获取的地震原始信号或者是按照匹配追踪算法执行完每一次循环后,得到的信号残差。所述待搜索原子的中心频率可以为所述当前信号在瞬时能量最高时刻所对应的瞬时频率。在本申请的一个实施例中,可以通过计算得到当前信号的解析信号,并通过解析信号得到瞬时能量最高的时刻,该时刻的瞬时频率可以作为待搜索原子的中心频率。具体的,若当前信号为c(t),则解析信号为c(t)+iH[c(t)],再根据以下公式得到待搜索原子的中心频率。

首先,根据以下公式得到所述解析信号瞬时能量最高的时刻。

tn=argmax||c(t)+iH[c(t)]|| (1)

式中,tn表示瞬时能量最高的时刻,c(t)表示当前信号,H[·]表示希尔伯特变换;

当前信号在tn处的瞬时频率为待搜索原子的中心频率。

式中,ωn表示待搜索原子的中心频率。

步骤S103,使用匹配追踪算法搜索预设过完备库中满足第一预设条件的原子,选择所述满足第一预设条件的原子中与所述当前信号最相关的原子,获取所述当前信号在所述最相关的原子处的对应投影分量,并获取所述当前信号与所述对应投影分量的信号残差;所述第一预设条件为所述原子的中心频率等于所述待搜索原子的中心频率。

所述预设过完备库可以为一种列数大于行数的矩阵,又可称为过完备字典。如果能够确定需要搜索的原子的中心位置,就可以大大减少需要搜索的原子数量,确定待搜索原子的中心频率实质上是将预设过完备库中的原子进行了一个删选,删选得到的原子都是满足第一预设条件的原子,在这些满足第一预设条件的原子中根据匹配追踪算法得到与当前信号最相关的原子,匹配追踪的每一次循环都得到了一个最相关的原子,这些最相关原子可以用来表示地震信号。选择满足第一预设条件的原子中与所述当前信号最相关的原子,可以为找出所有满足第一预设条件的原子中,与当前信号的内积最大的那一个原子。

在本申请的一个实施例中,将过完备库表示为:G={gn(t)},gn(t)表示原子,其中n=0,1,2,3,…

所述当前信号在所述最相关的原子处的对应投影分量可以为当前信号在最相关的原子上的垂直投影an,具体可以表示为:<c(t),gn(t)>gn(t),其中c(t)为当前信号。当前信号与该投影分量的差信号为所述信号残差。

步骤S104,将所述信号残差作为新的当前信号,重复步骤S102至S104,直至当前得到的所述信号残差小于预设阈值为止。

匹配追踪时频分解的目的为将信号表示成原子的稀疏线性组合。每一次迭代过程中,原子的线性组合与原始信号之间都存在一定残差,当残差可以忽略不计时,我们就可以认为地震信号时频分解完成了。所述预设阈值可以为预先设定的一个可以忽略不计的界限,当信号残差小于这个界限时,可以认为地震信号可以通过各个原子的线性组合表示。

步骤S105,根据所有所述对应投影分量,得到所述地震数据的地震谱分解结果。

所述地震谱分解结果可以为地震信号在各个原子上投影分量的和,该过程用公式表达可以为:

式中,s(t)表示原始地震信号,an表示原始地震信号在每个原子上的投影长度。

由图1所示的实施例可知,本申请实施例在使用匹配追踪算法进行地震谱分解时,通过当前信号确定了过完备库中需要搜索的原子的中心频率,不同于现有技术的搜索所有原子的做法,提高了计算效率。另一方面,本申请实施例只搜索预设过完备库中满足第一预设条件的原子,从中确定用于表示地震信号的原子,最终得到的用于表达地震信号的过完备库与地震信号之间更加匹配,提高了匹配追踪算法地震谱分解的分解精度和对不同地震信号的适应性。

时频分析时,过完备库中的原子一般都是基于已知形式的小波构建的原子,具体的可以包括雷克子波(Ricker子波)和Morlet子波等。在本申请的一个实施例中,使用雷克子波(Ricker子波)构建过完备库。原子的一般表达式可以为:

式中,Ln表示尺度参数,tn表示小波的位移参数,表示构成原子的母小波。

如果母小波为:

那么,将式(5)代入式(4)可以得到原子表达式为:

公式(6)所示原子的频谱为:

已知Ricker子波的形式为:

式中,f0表示Ricker子波的峰值频率,ω0=2πf0

Ricker子波对应的频谱为:

式中,f表示频率。

根据以上公式可知,利用Ricker子波构建的原子表达式为:

令ωn=ω0/Ln为小波的中心频率,则:

其频谱为:

由式(12)可知,此类原子频谱的相位由小波的中心位置tn确定,其振幅在小波中心频率ωn处最大。本实施例中,在根据以上公式确定原子表达式之后,可以按照图1所示的流程图,进行地震谱分解。

Ricker子波具有形式简单,零相位等优点,在本实施例中,将其作为构建原子的母小波,计算过程更加简单,也更接近爆炸子波。

在本申请的一个实施例中,为了避免用来表示地震信号的各个原子之间的间距过小,引入了一个新的参数——最小原子间距,该变量规定了用于表达地震信号的原子之间可允许的最小间距,这样在确定原子库中与当前信号最相关的原子时,除了要选择内积最大的原子,还需要满足该原子与之前所有被选中用来表示地震信号的原子之间的位置间隔大于或等于预设最小原子间隔。

具体的,在本申请的一个实施例中,S103中所述选择所述满足第一预设条件的原子中与所述当前信号最相关的原子,可以包括以下两个步骤。

(1)确定所述满足第一预设条件的原子中与所述当前信号的内积绝对值最大的原子。

(2)判断该内积绝对值最大的原子与之前确定的最相关的原子之间的间距是否大于等于预设最小原子间隔,若判断结果为是,则该内积绝对值最大的原子为所述最相关的原子,若判断为否,则将该内积绝对值最大的原子从所述满足第一预设条件的原子中剔除,重复以上步骤(1)~(2),直至判断结果为是为止。

在本申请的一个具体实施例中,之前按照图1所示流程,并同时考虑最小原子间隔,得到了用来表示地震信号的原子g1(t)、g2(t)和g3(t),此时,再根据S102中步骤确定了此次循环所要搜索的原子应该满足的第一预设条件,假设预设过完备库中,满足第一预设条件的原子为g4(t)、g5(t)和g6(t)。此时执行S103时,首先确定g4(t)、g5(t)、g6(t)中与当前信号最相关的原子,假设为g4(t),再判断g4(t)与之前确定的每一个原子g1(t)、g2(t)和g3(t)之间的间隔是否大于等于最小原子间隔,若不是,则舍弃g4(t),重新在剩下的原子g5(t)和g6(t)之间确定一个与当前信号最相关的当前原子,重复以上步骤,直至得到满足最小原子间隔的当前原子。

研究表明,原子间距是影响匹配追踪计算效率和分解质量的关键,但遗憾的是常规计算中,并没有考虑这一因素,在本申请上述几个考虑原子间隔的实施例中,将最小原子间隔考虑进去有利于提高分解精度。

在本申请的一个具体实施例中,预设最小原子间隔为10,即按照图1的流程提进行地震谱分解时,只有S103中选择的与当前信号内积最大的原子满足与之前得到的所有最相关原子之间的间隔大于等于10时,才可以留下该内积最大的原子,作为最相关原子,用来表达地震信号。

在本申请的一个实施例中,为了验证以上本申请实施例所提的改进匹配追踪算法的有效性,进行了相应的数值验证试验。设计一个主频为30HZ,延续长度为200ms,子波峰值在70ms处的Ricker子波,然后分别用小波变换CWT,S变换,常规匹配追踪算法以及图1所示实施例所提供的改进匹配追踪算法对该Ricker子波进行时频分析,来测试各个分解算法的时频局部性是否良好以及是否具有较高的时间频率域分辨率。图2是前面描述的主频为30HZ的Ricker子波,图3至图6分别为小波变换CWT,S变换,常规匹配追踪算法以及本申请实施例所提的改进匹配追踪算法进行时频分解的结果。图3至图6中纵坐标均为时间,单位是ms。从图3至图6可以看出,分解后的图像中,两种匹配追踪算法得到的结果在时间和频率域的分辨率比小波变换CWT和S变换都要高,时频局部性较好。

虽然由于本实施例中只输入了一个Ricker子波,稀疏表示的原子数目较少,而且间距都比较大,改进的匹配追踪算法和常规匹配追踪算法结果类似,但是图3至图6依然可以表明,改进的匹配追踪算法时频局部性良好。

地震数据经过高精度时频分解后,结合图形学、三维可视化等技术可精确的对储层、油气等特殊地质现象开展精细解释。将本申请实施例所提出的改进后的匹配追踪算法运用到来自于某探区碳酸盐岩储层中,地震信号采样率1ms,主频约35HZ,频带宽度约7-85HZ。图7为某一道地震道的地震信号,横坐标为时间(单位ms),纵坐标为振幅。图8为对图7所示的地震信号使用CWT分解的结果,图9为对图7所示的地震信号使用常规匹配追踪分解的结果,图10为对图7所示的地震信号使用改进后的匹配追踪分解的结果,图8至图10中显示出不同的反射率,在时频域可以精确的识别单个反射目标。图8至10中,横坐标为时间,单位ms,纵坐标为频率。对比图8-10看可以得出改进后的匹配追踪算法较CWT方法、常规匹配追踪具有较高的分解精度。

在本申请的一个实施例中,图11-13是用如图1所示的改进匹配追踪算法产生的频率剖面,图11-13中横坐标表示地震剖面中的记录号,纵坐标表示时间,单位s。图11-13分别对应10Hz频率剖面、30Hz频率剖面和50Hz频率剖面。在图11-13中1.85s、2.25s处反射振幅从10HZ到50HZ逐渐增强并随后减弱,30HZ左右最强。这种现象的出现要归结到储层含气特征,储层含气后,使共振频率朝能量较高的一侧迁移,当照亮这些共振频率时,就很容易识别。本实施例说明,如图1所示的改进的匹配追踪方法可以帮助实现在频率域分辨单个目标反射。

在本申请的另一个实施例中,利用如图1所示的改进的匹配追踪方法识别一些特殊的地质目标,如河道、岩性异常体、低频阴影等,并可进一步利用时频分解结果开展基于图像学的多谱图像合成,离散型的频谱分解数据也可以通过动画显示帮助解释人员了解潜在的地质目标,并对这些潜在目标进行分析解释。这可以在部分程度上解决利用单一谱分量不能展示全部信息的缺陷,提高解释精度,极大地挖掘有效信息。图14和图15为10Hz、30Hz和50Hz沿层切片多谱图像合成结果。其中,图14为使用本申请实施例所提的改进的匹配追踪方法分解获得的结果,图15为使用短时快速傅里叶变换(FFT)分解获得的结果。合成后图中储层得到进一步精细刻画,可以很好的对其局部特征结合多种数据进行精细分析。对比图14和图15可知,图14中右上区域内地层的非均质性刻画更加清晰,可以精细识别多个溶洞,另外可以看到北西南东方向的早期断裂被后期发育的北东向断裂所切割,断裂系统更加清晰,特征明显,更加有利于精细解释。

相对于其它分解算法,匹配追踪算法极其耗时,如何提高计算效率,分解效率和数据量存在密切关系,为了进一步提高计算效率,在本申请的一个实施例中,在使用如图1所示改进的匹配追踪算法的基础上,针对匹配追踪算法的特点,采用CPU多线程策略,极大地提高了计算效率。以Xeon系列E5-2670处理器为例,其主频2.30Hz,采用CPU多线程策略,极大地提高了计算效率。图16为单CPU情况下对比图1所示改进的匹配追踪算法与短时窗FFT算法,二者同时利用单线程对不同的样点数进行谱分解,在CPU主频一定的情况下,随着参与计算样点数增加FFT与本方法的耗时都在近似线性增加,但本申请实施例所提的改进的匹配追踪算法的计算耗时的速度增大更快,相比较FFT在较多样点数时更加耗时,这是由于本方法计算时使用大量冗余算子所影响。

在CPU主频一定的情况下固定参与计算的样点数,通过CPU多线程进行并行改造,对两种算法进行耗时对比,如图17所示,从图中可以看到随着参与运算的线程的增多,两种算法的计算效率得到大幅度提升,但随着线程的增多耗时性能提升速率相对降低。本申请实施例所提的改进的匹配追踪算法与FFT算法都具有较好的并行化特征,通过多线程并行加速能够达5-6倍的加速效果,且运算的样点数越多,并行加速效果越明显。

本申请实施例中还提供了一种匹配追踪地震谱分解装置,如下面的实施例所述。由于该装置解决问题的原理与一种匹配追踪地震谱分解方法相似,因此该装置的实施可以参见一种匹配追踪地震谱分解方法实施,重复之处不再赘述。

如图18所示,本申请实施例所提供的一种匹配追踪地震谱分解装置,可以包括以下几个模块。

获取模块1801,用于获取地震信号,并将所述地震信号作为当前信号。

中心频率确定模块1802,用于根据所述当前信号,确定待搜索原子的中心频率。

搜索模块1803,用于使用匹配追踪算法搜索预设过完备库中满足第一预设条件的原子,选择所述满足第一预设条件的原子中与所述当前信号最相关的原子,获取所述当前信号在所述最相关的原子处的对应投影分量,并获取所述当前信号与所述对应投影分量的信号残差;所述第一预设条件为所述原子的中心频率等于所述待搜索原子的中心频率。

重复模块1804,用于将所述信号残差作为新的当前信号,重复执行所述中心频率确定模块至所述重复模块,直至当前得到的所述信号残差小于预设阈值为止。

结果获得模块1805,用于根据所有所述对应投影分量,得到所述地震数据的地震谱分解结果。

由上述本申请装置的实施例可知,本申请实施例在使用匹配追踪算法进行地震谱分解时,通过当前信号确定了过完备库中需要搜索的原子的中心位置和中心频率,不同于现有技术的遍历所有原子的做法,提高了计算效率。另一方面,本申请实施例只使用了事先建立的过完备库中满足第一预设条件的原子去表示地震信号,相当于根据地震信号重新得到了与该信号相匹配的新的过完备库,提高了匹配追踪算法的适应性和分解精度。

本申请实施例中所描述的方法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动磁盘、CD-ROM或本领域中其它任意形式的存储媒介中。示例性地,存储媒介可以与处理器连接,以使得处理器可以从存储媒介中读取信息,并可以向存储媒介存写信息。可选地,存储媒介还可以集成到处理器中。处理器和存储媒介可以设置于ASIC中,ASIC可以设置于用户终端中。可选地,处理器和存储媒介也可以设置于用户终端中的不同的部件中。

在一个或多个示例性的设计中,本申请实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现,这些功能可以存储与电脑可读的媒介上,或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如,这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置,或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外,任何连接都可以被适当地定义为电脑可读媒介,例如,如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电缆、双绞线、数字用户线(DSL)或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片(disk)和磁盘(disc)包括压缩磁盘、镭射盘、光盘、DVD、软盘和蓝光光盘,磁盘通常以磁性复制数据,而碟片通常以激光进行光学复制数据。上述的组合也可以包含在电脑可读媒介中。

以上所述的具体实施例,对本申请的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本申请实施例的具体实施例而已,并不用于限定本申请的保护范围,凡在本申请的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本申请的保护范围之内。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1