青藏高原是地球上陆-陆碰撞的典型地区之一,通常认为青藏高原是岗瓦纳大陆与欧亚大陆长期相互作用的结果,其隆升也仅仅是最近几百万年以来的事件. 青藏高原年轻、独特、完整、典型的地质地貌特征已成为国际地球科学的热点研究地区之一. 青藏高原东北缘是印度和欧亚两大板块碰撞作用由近南北方向向东和北东方向转换的重要场所,物质东流的汇聚之处,高原内部块体及边缘块体地壳变形强烈,并在其接触区域形成一系列断裂带(黄全忠,蒲玉太,1989;邹谨敞,吴增益,1989;万天丰,赵维明,2002). 东昆仑—西秦岭活动构造带是青藏高原内一条重要的陆内板块碰撞拼合带,其中位于川北甘南交界的阿尼玛卿缝合带东段是青藏高原东北缘东昆仑南缘缝合带的重要组成部分之一,为甘孜—松潘地块和东昆仑—西秦岭褶皱带两大地质构造单元的分界. 多年来,地学研究人员从不同的研究领域对这些重要构造区进行了大量的研究,特别是采用主动源地震测深方法探测其深部结构和构造,提供地壳运动与地球动力学过程的壳幔深部地震学证据,这对于深入了解青藏高原的隆升变形机制、物质东流等大陆动力学基本过程有着十分重要的意义(许志琴等,1991;边千韬等,1999;张国伟等,2004;杨经绥等,2005;王椿镛等,2003).
深地震测深剖面方法是获得地壳结构信息的重要手段之一. 上地壳结构信息对于了解主要断裂的空间展布性质、深浅构造关系以及地壳变形特征都是十分重要的. 它对于正确了解下地壳乃至上地幔的结构特征也起着不可或缺的作用. 由于在地质结构复杂地区,上地壳通常经受强烈变形,并具有较强的非均匀性,因而要求有更高分辨探测的结果. 地震测深观测技术和资料处理解释方法的快速发展,使我们获得更加精细的上地壳结构模型成为可能(Fuis,1996;Fuis et al, 2001).
2004年在青藏高原东北缘的川北甘南地区穿越阿尼玛卿缝合带东段完成了马尔康—碌曲—古浪深地震测深剖面,并对阿尼玛卿缝合带东段及其两侧的红原—合作段进行了加密炮点和观测点的重点控制. 本文利用这次探测获得的初至折射波资料,采用多种处理解释方法综合分析,获得了阿尼玛卿缝合带东段及其两侧的红原—合作段区域上地壳结构特征.
1 剖面位置和观测系统马尔康—碌曲—古浪深地震测深剖面全长637 km,呈近南北向展布. 其南端起于四川省红原县刷经寺村(102°33.07′E,32°0.0′N),桩号75 km,北端止于甘肃省黄华镇移民村(103°10.08′E,37°40.37′N),桩号712 km. 该剖面由南向北相继跨过松潘—甘孜微块体、阿尼玛卿缝合带东段、东昆仑—西秦岭褶皱带、祁连褶皱带等地质构造单元. 穿越的主要断裂有: 阿坝弧形断裂、阿尼玛卿缝合带内的库赛湖—玛沁断裂、武都—迭部断裂、舟曲—两当断裂以及秦岭地轴北缘断裂等. 其中,阿坝弧形断裂、库赛湖—玛沁断裂、武都—迭部断裂、舟曲—两当断裂分别约在剖面170,280,310 km和330 km桩号附近穿过(图 1).
![]() |
图 1 研究区主要地质构造和深地震测深剖面位置图 ① 阿坝弧形断裂;② 库赛湖—玛沁断裂;③ 武都—迭部断裂;④ 舟曲—两当断裂;⑤ 秦岭地轴北缘断裂;⑥ 日月山—拉脊山断裂;⑦ 毛毛山—南西华山断裂 |
本剖面的观测系统如图 2所示. 沿剖面共进行了13次激发爆破. 由南向北分别在154 km(B1)、205 km(B2)、247 km(B3)、290 km(B4)、377 km(B5)、511 km(B6)、606 km(B7)桩号进行了7次TNT大药量激发,观测点间距2~3 km. 为了获得阿尼玛卿缝合带及其两侧更为精细的上地壳结构,对剖面的红原—合作段(约200~400 km桩号)进行了重点控制,在该段的224 km(S1)、261 km(S2)、277 km(S3)、320 km(S4)、359 km(S5)、380 km(S6)桩号进行了6次TNT中等药量激发,并加密观测,观测点间距约0.8~1 km. 从而形成了以控制阿尼玛卿缝合带为主的较为完整的深地震测深相遇追逐观测系统(图 2).
![]() |
图 2 马尔康—碌曲—古浪宽角反射/折射剖面观测系统 (a)观测系统图;(b)布格重力异常;(c)高程;(d)沿剖面地质构造及地名. ① 阿坝弧形断裂;② 库赛湖—玛沁断裂;③ 武都—迭部断裂;④ 舟曲—两当断裂;⑤ 秦岭地轴北缘断裂;⑥ 日月山—拉脊山断裂 |
本文研究的区域位于阿尼玛卿缝合带及其两侧的红原—合作段内(约200~400 km桩号),其中包括了若尔盖盆地、阿尼玛卿缝合带和西秦岭褶皱带中南段.
2 Pg波走时曲线特征Pg波是地震记录中的初至震相,被认为是来自上地壳顶部覆盖层的回折波或首波. Pg波走时曲线直观地反映了介质速度和界面特征的变化,是构造及其特征判别的震相依据. 当有速度异常和界面起伏时,Pg波走时常表现出显著超前或滞后,视速度则有突变显示. 远炮点的Pg波视速度则反映了与此相应的折射界面的速度.
图 3给出了剖面205,247,258,280 km和356 km桩号炮的Pg波记录截面. 图 3显示,在各炮记录截面上,Pg波大约可在炮检距0~100 km范围内被连续追踪对比. 近炮点Pg波视速度较低,反映其回折波特征. 由近炮点Pg波走时估算的近地表速度,自南向北分别为: 205 km桩号处4.37 km/s,247 km桩号处3.24 km/s,258 km桩号处1.46 km/s,280 km桩号处1.46 km/s以及356 km桩号处3.51 km/s(图 1). 若尔盖盆地内近地表速度明显偏低,反映了表层沉积盖层的特征. 随着炮检距的增大,Pg波视速度逐渐增大,在炮检距约大于60 km之后,视速度趋于稳定,达到5.8~6.0 km/s. 沿剖面Pg波走时曲线起伏较大,考虑到在本文讨论的剖面范围内高程变化最大约为490 m,因此,Pg波走时曲线如此大的起伏变化不会是由于高程差异引起的. 稳定后的Pg波折合到时约在0.8~1.7 s之间变化,反映了与此相应的盖层埋深存在较大的变化. 约160~280 km桩号范围,大致对应于若尔盖盆地,Pg波到时相对于南北两侧表现出明显的超前,反映了若尔盖盆地南侧的阿坝弧形断裂带和北侧的库赛湖—玛沁断裂带区域的基底下陷和若尔盖盆地相对隆升的构造特征. 大约在340~360 km桩号范围内,Pg波到时相对其两侧明显超前,表明基底在阿尼玛卿缝合带北侧隆升之后继续向北下陷. 由图 3可见,约在剖面170,280,315 km和335 km桩号附近,Pg波走时均出现明显突跳,暗示了断裂构造存在的可能性.
![]() |
图 3 205,247,258,280 km和356 km桩号炮的Pg波记录截面图 |
该方法利用接收段所有的初至波走时资料,包括近炮点视速度变化大的折射走时和远炮点稳定的折射走时,重建沿剖面从地表到折射深度层之间的速度结构图象(Vidale, 1988,1990;Hole,1992).
利用炮点位于200~400 km桩号段共10炮的Pg波资料,约900多个走时数据. 反演初始模型的建立主要依据折射地震记录的有关参数. 浅地表速度由近炮点直达波走时估算得到;较深处的初始模型速度则由较远炮点的Pg波视速度给出. 用不同的初始模型进行了反演,以检验解对初始模型的依赖程度. 经过10次反演迭代,得到了沿剖面200~400 km段的上地壳P波速度结构.图 4a给出了P波速度结构图象和一维初始模型,图 4b是相应的射线数分布图. 由图 4b可见,在100~270 km桩号区间8 km深度以上,270~320 km桩号区间10 km深度以上以及350~400 km桩号区间4 km深度以上,射线数大多大于20条,有的区域达80条以上. 由图 4a可见,速度等值线的剧烈变化约发生在2 km深度之下,相对而言,2 km深度之上,速度等值线变化较为平缓. 浅地表速度在阿尼玛卿缝合带南端较低,若尔盖盆地内次之,在阿尼玛卿缝合带以北的西秦岭褶皱带中南段区域最高,反映了沿剖面的地表结构特征.
![]() |
图 4 马尔康—碌曲—古浪剖面Pg波有限差分反演结果 (a)P波速度结构图;(b)射线数分布图. ① 阿坝弧形断裂; ② 库赛湖—玛沁断裂;③ 武都—迭部断裂;④ 舟曲—两当断裂 |
在若尔盖盆地内,P波速度纵向梯度变化不大,横向比较均匀. 阿尼玛卿缝合带则总体表现为明显南倾的低速带,其内部速度结构纵向梯度变化明显,横向上呈现强烈的非均匀性,体现了阿尼玛卿缝合带上地壳受挤压的破碎结构特征,其南倾的产状特征则可能暗示了阿尼玛卿缝合带在深部向其南侧若尔盖盆地扩展的迹象. 在西秦岭褶皱带中南段地区,速度变化较为平缓.
在阿尼玛卿缝合带南端与若尔盖盆地接触区域,约280~300 km桩号范围,速度等值线自地表就表现出强烈的下凹形态,该低速结构随深度加深向南倾斜,逐渐形成了宽度近30 km、深度约达8 km的低速带(如图 4a中A1所示),与其南侧的若尔盖盆地形成明显的速度差异;大约在310 km桩号附近,约3 km深度之下可见一个较小规模的低速带(如图 4a中A2所示),其延伸深度大约6~7 km;大约在330 km桩号附近,约在2.2 km深度之下存在着一个延伸深度不小于6 km的低速带(如图 4a中A3所示).
3.2 Sg波走时的射线追踪反演在获得的地震记录截面上还可以识别Sg波(图 5),它是与Pg相应的横波. 利用Sg震相可以构组上地壳S波速度结构.
![]() |
图 5 380 km桩号炮的Sg波记录截面图 |
Zelt和Smith(1992)的射线追踪反演方法适用于任意类型的体波地震数据,它可以同时获得地壳二维速度结构和界面形态,并能给出反演结果的分辨和误差信息. 利用射线追踪反演方法对炮点位于150~600 km桩号段共12炮,约860个Sg走时数据进行了反演处理.为了可与上节得到的P波速度相比较,初始模型也选择一维速度梯度模型.经过5次反演迭代,得到了沿剖面的S波速度结构图象(图 6a)和相应的射线分布(图 6b). 由图 6b可见,剖面200~400 km桩号范围由于处在加密观测区段,射线数密集,反演结果比较可靠.
![]() |
图 6 马尔康—碌曲—古浪剖面Sg波射线追踪反演结果 (a)S波速度结构;(b)射线分布(图中虚线为模型参数化网格). ① 阿坝弧形断裂; ② 库赛湖—玛沁断裂;③ 武都—迭部断裂;④ 舟曲—两当断裂 |
由图 6a可见,对应于阿尼玛卿缝合带的280~350 km桩号区段内,S波速度整体呈现低速,并且存在较大的横向变化;而在其南北两侧的若尔盖盆地和西秦岭褶皱带中南段区域,S波速度变化则相对平缓. 在阿尼玛卿缝合带南端与若尔盖盆地接触区域,约280~300 km桩号范围,明显可见一个低速条带(如图 6a中B所示),其向下延伸深度达6~7 km.
4 上地壳泊松比分布在获得上地壳P波和S波速度分布的基础上,可进一步求得泊松比分布(刘志等,2005). 本文用于计算泊松比的P波和S波速度结构均由射线追踪反演方法得到,初始模型均选择连续速度模型. 图 7为泊松比分布图. 图中白色虚线为用于计算泊松比的P波速度等值线,它由Pg波走时射线追踪反演得到. 与Pg波走时有限差分反演结果(图 4a)相比,它在横向变化上较为平滑,而与S波速度结构(图 6a)相比,两者的速度变化形态比较接近. 泊松比分布结果主要反映构造差异特征.
![]() |
图 7 马尔康—碌曲—古浪剖面上地壳泊松比分布(图中白色虚线为P波速度等值线) ① 阿坝弧形断裂; ② 库赛湖—玛沁断裂;③ 武都—迭部断裂;④ 舟曲—两当断裂 |
由图 7可见,泊松比分布在8 km深度以上部位呈现明显的横向分区性. 约在280~350 km桩号区域,其泊松比平均值约为0.29,表现为高泊松比分布;而其南北两侧泊松比平均值约为0.25和0.26,表现为相对低值区. 这一中间高两侧偏低的泊松比分布特征,体现了阿尼玛卿缝合带与若尔盖盆地、西秦岭褶皱带中南段上地壳结构和物性的整体差异,也反映了它们之间的构造分界.
在缝合带南端,剖面280 km桩号左右,介质泊松比表现出显著的横向差异,并构成向深部延伸的差异条带(如图 7中C所示),反映了可能的构造单元分界.
5 基底结构沿剖面Pg波在远炮点视速度相对稳定,并接近于5.8~6.0 km/s,表明上地壳上部可能存在速度或速度梯度界面. 为了获得基底界面结构,对Pg波走时进行了射线追踪反演和时间项反演,并将两者的结果做了比较(图 8).
![]() |
图 8 马尔康—碌曲—古浪剖面Pg波射线追踪反演和时间项反演结果 ① 阿坝弧形断裂; ② 库赛湖—玛沁断裂;③ 武都—迭部断裂;④ 舟曲—两当断裂 |
利用炮点位于150~600 km桩号段共12炮的Pg波资料,约1 100个走时数据,采用射线追踪反演方法(Zelt,Smith,1992)进行了处理. 近炮点的Pg波作为回折波对待,远炮点的Pg波被看做是沿基底界面滑行的首波,目的是获得研究区范围内的基底构造形态.
图 8中白线显示的是由射线追踪反演得到的基底界面. 由图 8可见,沿剖面基底起伏较大. 在若尔盖盆地内(约180~280 km桩号),基底界面埋深约在3.0~4.5 km之间,该范围南北两侧基底埋深迅速加深. 这一界面埋深变化反映了若尔盖盆地基底总体上相对于其北侧的阿尼玛卿缝合带和其南侧的阿坝弧形断裂带呈隆起状态,在布格重力图上也反映为局部重力高(图 2b),但在220~240 km桩号附近,基底在总体上隆的背景下存在局部下凹. 在阿尼玛卿缝合带内,基底强烈下凹,最深处接近6 km. 在缝合带北侧,基底经小范围隆升,深度达2 km左右,之后呈向北下陷的界面形态. 为了与射线追踪反演结果比较,对炮点位于150~400 km桩号段共11炮的Pg波远炮点首波走时资料,约800个数据,采用时间项方法进行了反演,上覆介质平均速度由各炮一维速度-深度函数反演结果分段确定,得到了沿剖面基底界面形态(如图 8中黑色点所示). 尽管时间项方法假定了上覆介质均匀且界面倾斜起伏不大,但由其得到的基底界面形态仍与射线追踪反演结果总体一致. 基底下的速度也同时由时间项反演得到,约为5.75 km/s.
6 上地壳主要断裂通过对地震记录Pg波时距曲线特征的分析,可以帮助判定上地壳断裂存在的可能性(Yan et al,2005). 图 9给出了200~400 km桩号段共10炮的Pg波折合时距曲线. 由图可见,在280,315 km和335 km桩号附近,分别相应于库赛湖—玛沁断裂、武都—迭部断裂和舟曲—两当断裂位置,各炮走时曲线均发生明显的转折(或视速度反转),表明这些断裂可能至少向下延伸到基底.
![]() |
图 9 马尔康—碌曲—古浪剖面Pg波折合时距曲线图 ① 阿坝弧形断裂; ② 库赛湖—玛沁断裂;③ 武都—迭部断裂;④ 舟曲—两当断裂 |
Pg波走时有限差分速度反演结果(图 4)显示,在阿尼玛卿缝合带南端与若尔盖盆地接触区域存在明显的低速带结构. 该低速带位置大致与库赛湖—玛沁断裂带位置相当,其断裂特征自近地表就有明显表现,并且至少向下延伸至8 km深度. 310 km桩号附近约3 km深度之下显示的相对较小规模的低速带,可能是武都—迭部断裂的反映,其延伸深度约6~7 km. 330 km桩号附近约2.2 km深度之下存在的低速带结构,则可能与舟曲—两当断裂相对应,该断裂延伸深度不小于6 km. S波速度结构(图 6)和泊松比结构(图 7)对库赛湖—玛沁断裂都有明显的反映,而对武都—迭部断裂和舟曲—两当断裂的反映并不明显,这可能是由于射线追踪方法对横向变化有平滑效应. 因此,它对较弱构造的反映不如有限差分方法得到的结果明显. 从Pg走时射线追踪反演和时间项反演结果中可以看到,库赛湖—玛沁断裂和舟曲—两当断裂位置处基底界面埋深有明显变化,是区隔阿尼玛卿缝合带南北两侧的主要断裂.而武都—迭部断裂发生在基底界面深度变化带上,它与舟曲—两当断裂可能具有同一深部构造背景.
7 讨论和结论 7.1 结论图 10是200~400 km桩号范围的上地壳结构综合图(图 10),图中包括Pg波有限差分速度反演结果、射线追踪反演和时间项反演得到的基底界面结果以及有限差分射线数分布分析结果(徐朝繁等,2006).
![]() |
图 10 上地壳结构综合图 |
由南向北,研究区上地壳结构明显分成3个区域,即: 若尔盖盆地(约200~280 km桩号)、阿尼玛卿缝合带(约280~340 km桩号)和秦岭褶皱带中南段(约340~400 km桩号).
1)速度结构. 若尔盖盆地和西秦岭褶皱带中南段区域,其速度分布较为均匀,横向变化不大,与它们之间的阿尼玛卿缝合带低速断裂结构相比呈现相对稳定的速度结构特征. 阿尼玛卿缝合带主要由库赛湖—玛沁断裂、武都—迭部断裂和舟曲—两当断裂组成,速度与两侧差异强烈,以整体低速背景下呈现非均匀性为特征;与3条断裂对应的区域,上地壳表现为低速条带结构,低速特征最明显的是与若尔盖盆地接触区域的库赛湖—玛沁断裂,舟曲—两当断裂次之,武都—迭部断裂最弱.
2)泊松比结构. 阿尼玛卿缝合带相对其南北两侧的若尔盖盆地和西秦岭褶皱带中南段表现为高泊松比分布,体现了阿尼玛卿缝合带与南北两侧介质的物性差异以及断裂发育的破碎结构特征. 在库赛湖—玛沁断裂两侧,介质泊松比表现出明显的横向差异.
3)基底结构. Pg波走时曲线特征表明,在研究区2~6 km深度范围内存在明显的速度界面,该界面上方P波速度约为5.4~5.6 km/s,下方约为5.75~5.90 km/s. 沿剖面该基底界面横向变化强烈. 射线追踪反演和时间项反演结果显示,若尔盖盆地基底相对其南侧的阿坝弧形断裂带及北侧的阿尼玛卿缝合带呈上隆状态,局部存在起伏,埋深约3.0~4.5 km;西秦岭褶皱带中南段区域基底界面经小范围隆升之后呈向北下陷的形态. 由图 10可见,这两个区域由射线数分布分析方法(徐朝繁等,2006)得到的水平射线密集区埋深也大致与上述两种方法得到的基底深度具有较好地一致性,体现了若尔盖盆地和西秦岭褶皱带中南段区域完整稳定的基底界面构造特征. 而阿尼玛卿缝合带的基底结构复杂,射线追踪反演和时间项反演结果反映该区域基底呈强烈下凹的构造形态,射线数分布分析方法在该区域得到了多个不同深度的水平射线密集区(徐朝繁等,2006),但与射线追踪和时间项反演得到的基底界面并不重合,表明缝合带内基底并非为速度反差明显的界面,这意味着阿尼玛卿缝合带基底经历了强烈的挤压变形,呈现破碎结构特征.
4)断裂. 库赛湖—玛沁断裂是研究区规模较大的主要断裂,在P波和S波速度分布上均显示为向南倾的明显低速条带,延伸深度达8 km左右,断裂两侧基底埋深急剧变化. 武都—迭部断裂和舟曲—两当断裂也以南倾的P波低速分布为特征,但在S波速度分布上却不明显. 它们位于阿尼玛卿缝合带北端同一基底深度变化带上,具有相同的深部构造背景.
7.2 讨论利用深地震折射剖面的Pg和Sg走时资料建立上地壳精细结构图象,不仅对于研究深浅构造关系以及上地壳速度分布有重要意义,也是进一步利用续至震相建立地壳模型的必要条件. 为了使获得的上地壳结构模型尽量可靠,本文对于相同的数据集采用了不同的正反演方法,并且将其结果进行了综合比较. 尽管不同的方法假设前提不同,但获得的结构图象在整体上的一致和差异反映了结果的可接受程度. 基于连续速度模型的有限差分成像可以获得较为详细的速度分布图象,但对于速度界面的判断不够直接. 尽管近水平射线密集程度可以帮助判别可能存在的速度界面(徐超繁等,2006),但因素并非单一. 基于射线追踪的反演方法有利于获得界面的构造形态. 因此它们的结合使用可以得到较为详细的上地壳结构.
[1] |
边千韬,罗小全,陈海泓,等. 1999. 阿尼玛卿蛇绿岩带花岗—英云闪长岩锆石U-Pb同位素定年及大地构造意义[J]. 地质科学,34(4): 420-426.(![]() |
[2] |
黄全忠,蒲玉太. 1989. 四川地震构造[G]//马杏垣主编.中国岩石圈动力学地图集. 北京: 中国地图出版社: 49.(![]() |
[3] |
刘志,张先康,王夫运,等. 2005. 用地震走时反演计算长白山天池火山区二维地壳泊松比[J]. 地震学报,27(3): 324-331.(![]() |
[4] |
万天丰,赵维明. 2002. 论中国大陆的板内变形机制[J]. 地学前缘,9(2): 451-463.(![]() |
[5] |
王椿镛,韩渭宾,吴建平,等. 2003. 松潘—甘孜造山带地壳速度结构[J]. 地震学报,25(3): 229-241.(![]() |
[6] |
许志琴,候立玮,王大可,等. 1991. “西康式”褶皱及其变形机制: 一种新的造山带褶皱类型[J]. 中国区域地质,(1): 1-9.(![]() |
[7] |
徐朝繁,张先康,张建狮,等. 2006. 射线数分布分析法及其在地壳上部复杂结构探测中的应用[J]. 地震学报,28(2): 167-175.(![]() |
[8] |
杨经绥,许志琴,李海兵,等. 2005. 东昆仑阿尼玛卿地区古特提斯火山作用和板块构造体系[J]. 岩石矿物学杂志,24(5): 369-380.(![]() |
[9] |
张国伟,郭安林,姚安平. 2004. 中国大陆构造中的西秦岭—松潘大陆构造结[J]. 地学前缘,11(3): 23-32.(![]() |
[10] |
邹谨敞,吴增益. 1989. 甘肃地震构造[G]//马杏垣主编. 中国岩石圈动力学地图集. 北京: 中国地图出版社: 54.(![]() |
[11] |
Fuis G S. 1996. Images of crust beneath southern California will aid study of earthquakes and their effects[J]. EOS AGU,77: 173-176.(![]() |
[12] |
Fuis G S,Ryberg T,Lutter W J, et al. 2001. Seismic mapping of shallow fault zones in the San Gabriel Mountains from the Los Angeles region seismic experiment, Southern California[J]. J Geophys Res,106(B4): 6 549-6 568.(![]() |
[13] |
Hole J A. 1992. Nonlinear high-resolution three-dimensional seismic travel time tomography[J]. J Geophys Res,97: 6 553-6 562.(![]() |
[14] |
Vidale J E. 1988. Finite-difference calculation of traveltime[J]. Bull Seism Soc Amer,78: 2 062-2 076.(![]() |
[15] |
Vidale J E. 1990. Finite-difference calculations of traveltimes in three dimensions[J]. Geophysics,55: 521-526.(![]() |
[16] |
Yan Zhimei,Robert W C,Jason S. 2005. Seismic refraction evidence for steep faults cutting highly attenuated continental basement in the central Transverse ranges, California[J]. Geophys J Int,160: 651-666.(![]() |
[17] |
Zelt C A,Smith R B. 1992. Seismic traveltime inversion for 2-D crustal velocity structure[J]. Geophys J Int,108: 16-34.(![]() |