柴达木盆地大风山凸起地层压力预测及成因分析

2023-01-18 06:53:46杨韬政刘成林田继先冯德浩李国雄吴育平
岩性油气藏 2023年1期
关键词:干柴大风声波

杨韬政,刘成林,田继先,李 培,冉 钰,冯德浩,李国雄,吴育平

(1.中国石油大学(北京)油气资源育探测国家重点实验室,北京 102249;2.中国石油大学(北京)地球科学学院,北京 102249;3.中国石油勘探开发研究院,北京 100083)

0 引言

孔隙压力预测对于钻井设计、油气成藏等研究具有重要意义,一般通过岩石性质的变化来量化孔隙压力,如声波速度或电阻率的变化[1-2]。孔隙压力预测方法可以分为3 类:①有效应力法。此类方法基于Terzaghi 有效应力原理(上覆地层压力等于岩石有效应力与地层压力之和),通过建立有效应力函数来计算有效应力,进而计算地层压力。代表性方法有等效深度法[3]、Eberhart-Phillips 法[4]、Bowers法[5-6]、Dutta 法[7]及Miller 法[8-9]。②经验公式法。此类方法是根据超压段测井与钻井响应特征的差异来拟合超压段的压力趋势,从而得到地层压力与各个测井参数之间的关系。代表性方法有Eaton法[10-11]、Dc 指数法[12-13]等。③纵波速度法。此类方法也称为地震预测法,通过地震纵波速度与地层压力之间的关系来计算地层压力,代表性方法主要有Fillippone 法[14-15]、刘震法[16]等。目前大多数孔隙压力预测成功的例子基本上都是预测不均衡压实产生的超压[17-20],然而,沉积后流体膨胀机制产生的超压不会造成孔隙度的异常,难以检测和量化孔隙压力。因此,要准确预测孔隙压力,必须要先了解超压产生的机制以及不同超压产生机制影响岩石性质的方式[2,21],然后再根据实际情况选择合适的压力预测模型。

Osborne 等[22]将造成地层超压的成因划分为三大类,即流体体积变化、压应力增加、流体流动及浮力作用。有学者将超压成因分为不均衡压实、流体膨胀、侧向传递及构造加载4 类。进入21 世纪以来,不少学者研究发现成岩作用特别是泥页岩中蒙脱石—伊利石的转化作用对超压的形成也具有重要贡献,且与流体膨胀机制不同,因此提出了蒙脱石向伊利石转化脱水形成的超压。不均衡压实成因超压会导致岩石骨架上的垂直有效应力增加或基本不变,与相同深度的正常压实地层相比,超压层的孔隙度会增加,密度和电阻率均会减小[23-24]。非不均衡压实成因超压,如流体膨胀,由于其形成原因的差异,与流体膨胀相关超压的识别和定量评价则相对困难。

柴达木盆地西部坳陷大风山凸起沉积构造背景复杂,有关地层压力方面的研究尚不充分,严重阻碍了该区井位钻探、油气运移及油气成藏等研究工作。为提高大风山凸起超压预测的准确性,基于已有地质、测井等资料,通过测井曲线组合、声波速度-垂向有效应力交会图和声波速度-密度交会图以及超压综合分析等方法,综合判断大风山凸起超压成因,根据不同层位超压成因的差异,利用平衡深度法预测不均衡压实成因超压,利用伊顿法预测复合成因超压,以期为该区地层超压预测研究提供借鉴。

1 地质概况

柴达木盆地位于青藏高原东北部,面积达12×104km2,是中国典型的内陆高原咸化湖盆[25-26]。柴达木盆地的形成经历了印支运动、燕山运动和喜马拉雅运动共3 个阶段的沉积和构造演化,形成了盆地目前的构造格局[27]。柴达木盆地西部是盆地重点油气勘探区域,具有沉积中心、咸化中心迁移频繁、岩性复杂和多种盐岩并存的特点[28-29]。柴达木盆地西部坳陷(柴西坳陷)可划分为4 个次级构造单元:昆北断阶、茫崖凹陷、大风山凸起和一里坪凹陷(图1a)。柴西坳陷新生界自下而上主要发育路乐河组(E1+2)、下干柴沟组下段(E31)、下干柴沟组上段(E32)、上干柴沟组(N1)、下油砂山组(N21)、上油砂山组(N22)、狮子沟组(N23)和七个泉组(Q1+2)(图1b)。大风山凸起位于柴西坳陷西北部,东邻一里坪坳陷,南至茫崖凹陷,西止月牙山构造,北抵阿尔金山前。大风山凸起内由多个地面构造组成,主要包括红三旱一号、尖顶山、碱山、大风山和长尾梁等地面构造,总面积达4 200 km2。柴西坳陷有很多与超压有关的构造,超压在各个构造带中几乎都有发育,被认为是柴西北油气运移的重要动力来源,且与油气藏的分布密切相关[30-31]。

图1 柴达木盆地大风山凸起构造单元划分(a)及地层岩性综合柱状图(b)Fig.1 Distribution of tectonic units(a)and stratigraphic column(b)of Dafengshan uplift in Qaidam Basin

2 地层压力预测

Hottmann 等[3]在研究墨西哥湾盆地页岩超压时,发现了超压段与常压段之间有效应力的差异。根据孔隙度相等(声波时差相等)的两点,其岩石骨架有效应力相等这一基本假设预测异常高压段的压力,即平衡深度法。图2 为平衡深度法原理示意图,图中B 点位于常压段,A 点位于欠压实段。其关系式为

图2 平衡深度法原理示意图(据文献[7]修改)Fig.2 Schematic diagram showing the principle of equilibrium depth method

式中:σA和σB分别为A 点和B 点处的有效应力,MPa;SA和SB分别为A 点和B 点处的上覆地层压力,MPa;PfA和PfB分别为A 点和B 点处的孔隙流体压力,MPa;ρA和ρB分别为A 点和B 点处的上覆地层平均密度,g/cm3;hA为A 点的深度,m;hB为A点对应的等效深度(B 点的深度),m。

使用平衡深度法计算地层压力的关键在于正常压实趋势线和等效深度的求取,主要包括3 个步骤:

(1)在测井资料中筛选厚度大于2 m 的泥岩,并拟合出该井的泥岩密度趋势线,进而求取任一点上覆地层的平均密度:

(2)作泥岩声波时差-深度交会图,并据此建立正常压实曲线,求A 点的等效深度hB:

(3)联立式(4)、式(5)和式(7)可以得到A 点的PfA的计算公式为

式中:hi为深度为i时的岩石密度,g/cm3;Δt为A 点处的声波时差,μs/m;Δt0为地表的声波时差;μs/m;c为压实系数(常数),由正常压实下的声波速度-深度函数关系得到。

使用平衡深度法计算大风山凸起F2 井地层孔隙压力,建立的正常压实趋势线方程为

根据F2 井的密度测井资料拟合的密度与深度的关系式为

式中:Y1为声波时差,μs/m;Y2为密度,g/cm3;X为深度,m。

平衡深度法的计算结果(表1)显示,对研究区下油砂山组的预测效果较好,总体预测误差小于6.00%,平均误差仅为4.33%;对上干柴沟组和下干柴沟组上段预测效果较差,误差均大于10.00%,最大误差可达18.56%。平衡深度法预测的是不均衡压实作用对超压的贡献,而对非不均衡压实成因和复合成因的超压预测效果较差[32]。不同成因的超压与油气藏形成和分布的关系不同,压力预测所使用的方法也有所不同。平衡深度法的计算结果表明,研究区下油砂山组可能为不均衡压实,而上干柴沟组和下干柴沟组上段则为复合成因或非不均衡压实成因。因此,为解决研究区不同层位压力预测误差大的问题,需要先分析不同层位间超压成因的差异,最后选择合适的方法进行压力预测。

表1 柴达木盆地大风山凸起F2 井平衡深度法计算结果与重复地层测试压力误差对比Table 1 Comparison of error between calculation results by equilibrium depth method and measured formation pressure of well F2 in Dafengshan uplift,Qaidam Basin

3 地层超压成因分析

3.1 测井曲线组合法

声波时差、密度、电阻率等测井曲线的组合特征研究是分析沉积盆地超压成因的一种方法,其判断依据是3 条测井曲线反转的同步性[33]。声波时差、密度、电阻率3 条测井曲线的反转包括以下3种情况[33]:①3 条测井曲线同步反转,可以判断为不均衡压实超压;②3 条测井曲线反转不同步或者密度不变或略有减小,指示为生烃增压等流体膨胀成因或压力传导成因;③3 条测井曲线均不发生反转,则超压可能为构造挤压成因。就大风山凸起F2 井而言,超压带上下界面处密度测井位于正常压实趋势线上,而声波时差、电阻率在超压界面下则出现异常,因此可以初步判断其超压的顶界面约在1 600 m。在通过声波时差、电阻率和密度测井的反转变化来分析超压成因时,实际上只考虑了单一成因造成这3 条测井曲线变化的结果,而对于复合成因形成的超压在测井曲线组合上的变化并未考虑。F2 井声波时差和电阻率的反转深度为1 600 m,密度曲线的反转深度则为2 400 m(图3)。一方面超压发育深度较小,地层难以被充分地机械压实;另一方面声波时差随着埋深的增加依然在缓慢偏离正常趋势。由此可以判断埋深为2 400 m 及以上地层的超压成因为不均衡压实,而埋深为2 400 m 及以下地层同样存在不均衡压实作用对地层超压有贡献的可能性。密度测井反转的界面可以作为流体膨胀开始对地层超压起作用的标志,这说明埋深为2 400 m及以下地层的超压为复合成因。

图3 柴达木盆地大风山凸起F2 井声波时差(a)、密度(b)和电阻率(c)测井的变化趋势Fig.3 Variation trend of acoustic time difference(a),density(b)and resistivity(c)logging of well F2 in Dafengshan uplift,Qaidam Basin

3.2 交会图版法

在声波速度-垂向有效应力交会图中,正常压实和不均衡压实成因超压会落在加载曲线上,而流体膨胀、成岩作用(蒙脱石—伊利石转化作用)、压力传递和构造挤压形成的超压均会落在卸载曲线上[34]。Bowers[6]指出声波速度与电阻率测井反映的是岩石的传导属性,而中子和密度测井反映的是岩石的体积属性,并将超压的成因划分为加载成因和卸载成因。实际上,如果地层同时保留了加载成因和卸载成因形成的超压,则在声波速度-垂向有效应力交会图中,超压点依然会落在卸载曲线上。实测超压计算的有效应力投影到交会图上的结果显示,研究区下油砂山组的散点落在加载曲线上,再次证实了其超压成因为不均衡压实,而对上干柴沟组和下干柴沟组上段的散点较为分散,无法准确地落在卸载曲线上。由此判断不均衡压实对上干柴沟组和下干柴沟组上段超压存在一定的贡献,同时也存在卸载成因的贡献(图4)。

图4 柴达木盆地大风山凸起F2 井声波速度-垂向有效应力交会图Fig.4 Cross plot of acoustic velocity and vertical effective stress of well F2 in Dafengshan uplift,Qaidam Basin

声波速度-密度交会图中不均衡压实(图5a 中BC 段)和构造挤压(图5a 中DC 段)形成的超压落在加载曲线上,两者在加载曲线上的延伸方向相反,而其他成因的超压如成岩作用(图5a中CF 段)、流体膨胀作用(图5a 中CH 段)和复合成因(图5a中CE 或CG 段)则会落在加载曲线之外。大风山凸起F2 井下油砂山组的点全部位于加载曲线上,而上干柴沟组和下干柴沟组上段的点则分布较为零散,小部分点位于加载曲线上,大部分点则位于加载曲线之外,且分布没有明显的规律(图5b)。如果将上干柴沟组和下干柴沟组上段的超压分析为单一成因,则很难解释其在声波速度-密度交会图上散点变化异常的原因。出现这一散点趋势的原因很可能是同一层位不同深度段各超压成因的贡献率存在差异,当不均衡压实作用和构造挤压作用的贡献率更大时,曲线会更加偏向加载曲线,而当流体膨胀等卸载成因超压的贡献率更大时,散点则会愈加偏离加载曲线。因此,研究区下油砂山组的超压成因可以判断为不均衡压实,而上干柴沟组和下干柴沟组上段的超压则为包括不均衡压实作用的复合成因,至于卸载成因是成岩作用、压力传递还是生烃膨胀作用则需要进一步分析。

图5 不同超压成因的交会图版(a)和柴达木盆地大风山凸起F2 井声波速度-密度交会图(b)Fig.5 Cross plot of different overpressure causes(a)and cross plot of acoustic velocity and density of well F2 in Dafengshan uplift,Qaidam Basin(b)

3.3 综合分析法

3.3.1 不均衡压实作用分析

声波时差的异常增大被认为是不均衡压实成因超压的重要标志。近年来,越来越多的学者提出声波时差的偏离并不能指示超压为不均衡压实成因[23-24,33,35]。平衡深度法的计算结果和交会图的分析显示,研究区下油砂山组的超压是不均衡压实作用产生的,上干柴沟组和下干柴沟组上段存在声波时差异常增大的现象,但使用平衡深度法估算的地层压力小于实测值,这说明不均衡压实是成因之一但不是唯一成因,即还存在其他增压作用。根据沉积速率分析,研究区F2 井存在2 个快速沉积时期,即下干柴沟组上段沉积时期和狮子沟组—七个泉组沉积时期(图6)。其中,下干柴沟组上段沉积时期虽然沉积速率很大,达到了600 m/Ma,但由于在沉积时埋藏较浅,地层孔隙的排水速率足够大,使得地层保持静水压力状态,因此这一时期的超压不是地层不均衡压实所形成的。狮子沟组—七个泉组沉积时期,沉积速率约为280 m/Ma,地层的快速沉积加上长时间的沉积使得研究区下油砂山组、上干柴沟组和下干柴沟组上段新生界达到一定的深度,完全具备形成不均衡压实作用超压的条件。

图6 柴达木盆地大风山凸起F2 井地层沉积速率Fig.6 Sedimentation rate of well F2 in Dafengshan uplift,Qaidam Basin

3.3.2 构造挤压作用分析

近年来,构造挤压作用对地层超压存在影响的观点被广泛接受[36-37]。构造挤压作用的数值模拟结果显示,构造增压作用与不均衡压实增压作用的机制类似,相当于横向上的压实[38]。相关研究表明在地层具有良好封闭性的条件下,强烈的构造挤压作用增加的剩余压力占总剩余压力的比例可达50%[39]。研究区新生界始终处于周缘山系的挤压构造背景下,受喜山运动早、中、晚3 期构造运动的影响,形成了目前的构造格局,尤其是喜山运动晚期的挤压作用最为强烈。喜马拉雅晚期,受印度板块持续向北剧烈挤压碰撞的影响,柴达木盆地在北南向挤压下抬升,柴达木盆地西部地区发生了以地层缩短和构造旋转为主的构造变形。因此,构造挤压作用是大风山新生界凸起的重要成因之一。

综上所述,构造挤压作用对整个研究区的地层超压均有贡献。因此,研究区下油砂山组地层超压成因为不均衡压实和构造挤压作用,但在使用平衡深度法计算时预测的效果较好。造成这一现象的原因可能是:在地层埋藏较浅时,由于地层的压实程度不高,较强的构造应力在横向上的压实不仅可以增加地层压力,也可以改变地层的孔隙;在地层埋藏较深时,由于地层的压实已经达到了一定的程度,构造应力无法造成孔隙的弹性回弹,只会导致地层压力的异常增大。使用平衡深度法计算研究区下油砂山组地层压力时,实际上计算的是不均衡压实和构造挤压作用共同产生的超压,因为这2 种增压作用的结果均导致了孔隙度的异常增大。

3.3.3 压力传递作用分析

地下剩余孔隙压力的重新分布称为压力传递,压力传递包括侧向传递和垂向传递[40]。通常认为渗透性砂岩层中的异常压力来自于邻近异常高压泥质岩层的压力传递。砂岩储层中的超压也可以来自更远或更深部超压源的压力传递,正如油气的运移聚集可以远离烃源岩一样。国外大部分学者认为超压不是一成不变的,而是一个可以通过断层或裂缝传递的具有瞬时性特征的水力学现象[2]。压力传递本身并不是产生超压的原生机制,但是它同样会影响地下孔隙压力的分布,导致超压的成因机制难以识别。

剩余压力是指地层压力与静水压力的差值,在沉积盆地中,当深、浅部地层存在剩余压力差时,深、浅部的超压体通常会通过开启的断裂连通,使得压力趋于新的平衡,最终达到平衡所需要的时间与断裂的输导性能、纵向剩余压力梯度及断裂垂向持续开启的时间有关[41]。大风山凸起是一个典型的随着埋深增大剩余压力逐渐增大的区域,从F2 井剩余压力的计算结果来看,以埋深2 500 m 为界可以分为2 个剩余压力系统:上部剩余压力系统和下部剩余压力系统(图7a)。上部剩余压力系统的剩余压力较小,且变化幅度较小;下部剩余压力系统的剩余压力与深度具有较好的线性关系,且剩余压力较大。从剩余压力的分布来看,压力传递作用主要发生在下部剩余压力系统中。大风山凸起地震解释剖面显示,断层可以作为压力传递的传导条件,当断层处于开启状态时,过剩压力的差异使得压力向上传递,导致浅部压力增加,深部超压减小(图7b)。因此,超压传递沿断层的垂向传递可能是研究区上干柴沟组和下干柴沟组上段超压的成因之一。

图7 柴达木盆地大风山凸起F2 井剩余压力特征(a)与地质剖面(b)Fig.7 Residual pressure characteristics(a)and geological profile(b)of well F2 in Dafengshan uplift,Qaidam Basin

3.3.4 生烃作用分析

有机质在热演化过程中生成烃类的体积会大于原本的干酪根体积,导致生成的烃类挤压原本存在的流体,占据了一定的孔隙空间,从而增加了孔隙压力。生烃作用将相对密度较大的固态干酪根转化为液态烃、气体烃或者油气裂解为分子量较小的轻烃,产生了密度较小的孔隙流体,使原有的孔隙流体体积膨胀总量可达25%[42]。大风山凸起的实测镜质体反射率为0.4%~1.3%,开始生烃(Ro>0.5%)的起始深度约为2 300 m(图8a),反映烃源岩达到成熟阶段。大风山凸起TOC数据分析结果显示,上干柴沟组TOC最大值为2.10%,平均值仅为0.46%,绝大部分样品点的TOC值低于烃源岩的生烃下限;下干柴沟组上段TOC最大值为1.54%,平均值仅为0.38%,大多数数据点的TOC小于0.5%(图8b)。因此,虽然上干柴沟组和下干柴沟组上段已经达到了生烃阶段,但由于有机质丰度没有达到认定为烃源岩的标准,故生烃作用在大风山凸起的影响不明显。

图8 柴达木盆地大风山凸起实测镜质体反射率(a)与总有机碳含量(b)Fig.8 Measured vitrinite reflectance(a)and total organic carbon content(b)in Dafengshan uplift,Qaidam Basin

3.3.5 成岩作用分析(蒙脱石—伊利石转化脱水)

黏土岩中蒙脱石向伊利石转化脱水的反应被认为是超压形成的重要机制[43]。这是由于蒙脱石晶体结构中含有丰富的层间水,而在转化为伊利石的过程中伴随着层间水的释放,孔隙流体的体积增加,从而形成超压。对柴达木盆地西部20 多口井各深度段不同层位黏土矿物X 射线的分析结果显示:柴西地区新生界中黏土矿物的主要成分为伊利石,其次为高岭石和绿泥石,而蒙脱石的含量极低[44],表明柴西地区黏土岩基本上不存在蒙脱石向伊利石的转化。柴达木盆地西部伊利石、高岭石和绿泥石含量随深度的增加并没有明显的变化(图9),由此表明黏土矿物成分的转化对柴西地区异常高压的形成贡献不大。

图9 柴达木盆地西部黏土矿物含量随深度变化特征(据文献[44]修改)Fig.9 Variation characteristics of clay mineral content with depth in western Qaidam Basin

综上所述,通过测井曲线组合、声波速度-垂向有效应力交会图、声波速度-密度交会图的分析,结合研究区实际地质条件的综合分析,认为大风山凸起下油砂山组地层超压成因为不均衡压实和构造挤压作用;上干柴沟组和下干柴沟组上段地层超压成因为不均衡压实、构造挤压作用和压力传递。

4 压力预测的改进

柴达木盆地大风山凸起上干柴沟组和下干柴沟组上段超压成因为复合成因,因此利用平衡深度法预测该段超压时会存在较大的误差。Eaton 公式是建立在泥页岩压实理论基础上的一种半定量经验关系式,其计算是基于压力异常区域声波时差等参数偏离正常压实趋势的程度,用实际参数值与正常趋势值的差值或比值作为直观反映孔隙压力变化的参数。Eaton 法中的伊顿指数N 具有成因意义,对于复合成因的超压仍能进行压力预测。Eaton总结以往的研究成果,提出了利用泥页岩声波时差、电阻率、电导率和密度的压力计算公式,后人称之为Eaton公式[10-11],计算公式为

式中:Pf为孔隙流体压力,MPa;S为上覆地层压力,MPa;Pn为静水压力,MPa;Δtn为正常压实趋势线上的声波时差,μs/m;Δt为实际值测井的声波时差,μs/m;N为伊顿指数。

根据大风山凸起的实际资料,也可以选择电阻率、电导率和密度代替声波时差进行压力计算,计算的形式与式(11)一致。在已建立的研究区F2 井正常压实趋势线的基础上,通过不断调整伊顿指数的值来减小预测结果与实测结果的误差,最终得出在伊顿指数N=4.2 时,预测结果达到最佳。

利用平衡深度法计算大风山凸起F2 井下油砂山组地层压力的效果较好(图10a),而利用伊顿法计算上干柴沟组和下干柴沟组上段地层压力的效果较好(图10b)。因此,对不同超压成因的层位使用不同的计算方法来预测单井地层压力值。分段预测地层压力的误差分析结果显示:相较于利用平衡深度法计算单井的压力,分段预测使得预测效果有了明显的改善,所有预测值与实测值的误差均小于7.00%,平均误差为4.30%(表2)。

表2 柴达木盆地大风山凸起F2 井分段计算地层压力误差对比Table 2 Comparison of formation pressure error in section calculation of well F2 in Dafengshan uplift,Qaidam Basin

图10 柴达木盆地大风山凸起F2 井平衡深度法计(a)和伊顿法(b)计算结果Fig.10 Calculation results by equilibrium depth method(a)and Eaton method(b)of well F2 in Dafengshan uplift,Qaidam Basin

5 超压对油气成藏的意义

超压因其在含油气盆地中分布的普遍性、与油气藏形成关系的密切性以及对钻井安全的重要性,一直是油气地质与勘探研究的热点问题[45]。超压的预测可以为油藏描述和储量估算提供必要的数据,也是井位设计的基础,对于确保安全钻井作业至关重要[46-47]。在泥岩中,油气的运移会受到细小孔径中毛管阻力的束缚,而只有当泥岩与邻近储集层和输导层孔隙流体间的压力差超过了油气运移的阻力时,油气才能从母岩中排出。因此,对于油气成藏来说,超压是油气运移的动力。柴达木盆地大风山凸起上干柴沟组和下干柴沟组上段为主力烃源岩。平面上连续的砂体和纵向上的断层是油气运移良好的疏导体系,是大风山凸起良好的油气运移通道。地层压力预测是计算剩余压力的基础,而剩余压力则是研究区油气运移的动力,可以指示油气运移的方向、估算油气运移的距离,因此对于油气成藏的研究具有重要指导意义。

6 结论

(1)柴达木盆地大风山凸起下油砂山组地层超压成因为不均衡压实和构造挤压作用,而上干柴沟组和下干柴沟组上段超压成因为不均衡压实、构造挤压作用和超压传递。

(2)大风山凸起不同层位超压成因不同,利用平衡深度法对下油砂山组地层压力的预测效果较好,而伊顿法则可以很好地预测上干柴沟组和下干柴沟组上段的地层压力。大风山凸起压力预测的改进对于地层压力分布特征、油气成藏等研究具有重要作用,也为其他地区地层压力预测的研究具有一定的借鉴意义。

(3)超压是大风山凸起油气运移的主要动力,而平面上的连续砂体和纵向上的断层是油气运移良好的疏导体系,过剩压力的分布特征对于大风山凸起的油气成藏有重要指导意义。

猜你喜欢
干柴大风声波
一捆干柴(外一首)
诗歌月刊(2021年6期)2021-07-01 13:59:45
母亲这把“干柴”
文苑(2020年6期)2020-06-22 08:41:46
爱的声波 将爱留在她身边
中国宝玉石(2018年3期)2018-07-09 03:13:58
大风吹(二)
幼儿100(2017年31期)2017-11-27 02:37:46
大风吹(一)
幼儿100(2017年28期)2017-10-27 01:45:49
声波杀手
自适应BPSK在井下钻柱声波传输中的应用
人小鬼大狄仁杰
“声波驱蚊”靠谱吗
一堆干柴