有限元分析颈椎棘突骨折内固定有效性

沈 彦 王朝阳* 吴 荣 陈 博

1(中国人民解放军第98医院骨科, 浙江 湖州 313000)2(上海交通大学医学院附属瑞金医院伤骨科研究所,上海 200025)

有限元分析颈椎棘突骨折内固定有效性

沈 彦1王朝阳1*吴 荣1陈 博2

1(中国人民解放军第98医院骨科, 浙江 湖州 313000)2(上海交通大学医学院附属瑞金医院伤骨科研究所,上海 200025)

通过对颈椎棘突骨折(累及椎板)内固定治疗有限元模型的建立和分析,明确此种治疗方式对颈椎棘突骨折的有效性。先建立正常全颈椎(C0-T1)的有限元模型并与文献报告进行对比验证,模型验证后,在正常模型基础上建立颈椎棘突骨折(累及椎板)模型,并模拟直型接骨板行内固定治疗,测量并比较颈椎棘突骨折模型及手术内固定模型和原始正常模型在前屈、后伸、左右侧弯、左右旋转6种条件下活动度改变。以及颈椎各结构的应力变化。结果表明，在正常模型上结合临床病例建立的颈椎棘突骨折(累及椎板)外观逼真,生物力学相似度良好。骨折模型部分节段,主要为C7-T1的活动度(前屈+后伸9.20°,左右侧弯5.83°,左右旋转13.12°)较正常模型(前屈+后伸7.11°,左右侧弯4.92°,左右旋转 9.59°)增大,尤其是旋转活动度,模拟植入内固定后稳定性增加(前屈+后伸4.07°,左右侧弯2.21°,左右旋转2.91°),且内固定钢板应力分析提示，承受最大应力值在安全范围。颈椎棘突骨折(累及椎板)及内固定模型可以较好地模拟临床实际病例，通过有限元分析预示，此型骨折存在潜在不稳的可能性,探讨微型棘突钢板在骨折手术治疗中的应用,具有一定的临床参考价值。

颈椎棘突骨折；内固定；有限元分析

引言

颈椎棘突骨折是一类附件骨折,通常也称为Clay shoveler骨折,即“铲土者骨折”,此类骨折对颈椎稳定性的影响较小,所以常常被临床工作者所忽视。虽然颈椎棘突骨折患者保守治疗常常可以获得满意的临床效果,但随着对颈椎棘突骨折的研究逐渐加深,有学者发现，此类骨折的部分类型有潜在颈椎不稳的可能性,特别是累及椎板的棘突骨折,因其损伤暴力较大,常常存在潜在的颈椎不稳及神经损伤[1-2]。临床上一旦发生颈椎不稳和(或)脊髓损伤,需及时给予手术治疗以避免严重并发症的发生。由于此类骨折发生率较低,前人对颈椎棘突骨折的研究大多是个案报道为主[3-5],缺乏大样本数据支持,治疗方案也主要依靠临床医师个人经验决定。故此研究运用有限元方法结合临床实际病例建立颈椎棘突骨折(累及椎板)模型[6-7],探讨此型骨折对颈椎稳定性的影响,并在此基础上模拟自行设计的颈椎微型钢板行内固定手术治疗[8-9],对其进行生物力学分析,探讨此种治疗方式的有效性及安全性,为临床上颈椎棘突骨折的诊断及治疗提供理论依据。

1 材料与方法

1.1 实验对象

青年志愿者(男,22岁,身高175 cm,体重70 kg),既往无颈椎疾病及外伤病史,检查前无明显颈部不适症状,同时给予相关检查排除可能影响本研究的其他相关疾病。本实验经过中国人民解放军第九八医院伦理委员会审核批准,符合医学伦理学要求及规范。

1.2 实验环境

1)硬件环境:GE 宝石 CT (Discovery CT750)及 HP Z800数据处理工作站。

2)软件环境:Mimics17 (医学有限元建模软件)，Geomagic Studio 2012 (逆向重建软件)，Hypermesh12.0 (数据分析前处理软件)，ABAQUS 6.9 (有限元分析处理软件)。

1.3 数据获取及实体模型建立

1.3.1 正常模型建立

采用GE Discovery 宝石 CT扫描机,在自然状态(取仰卧位,保持颈椎中立)下进行扫描,扫描范围为枕骨至第二胸椎(C0-T2),得到0.625 mm层厚的CT断层扫描影像522张,点阵物理距离为0.799 mm,生成DICOM格式导入数据处理工作站。运用Mimics17把前期临床采集的CT图像数据逆向重建出颈椎的三维点云模型,运用Geomagic Studio 2012软件对逆向重建的数据模型图像进行重建后修饰,转化为外观逼真的三维模型。

1.3.2 分析前处理

颈椎13种韧带主要包括:黄韧带(LF)、齿状突尖韧带(APL)、前纵韧带(ALL)、后纵韧带(PLL)、覆膜(TM)、横韧带(TL)、十字韧带垂直部分(CLV)、寰枕前膜(AAOM)、寰枕后膜(PAOM)、翼状韧带(AL)、关节囊韧带(JC)、棘上韧带(SSL)、棘间韧带(ISL),相关韧带的起止点参考相关论著[10]。各颈椎关节间的相互作用通过非线性通用接触关系来模拟。颈椎小关节以及相应棘突之间的接触,齿状突与寰横韧带、齿状突与枢椎的接触,C3-C7钩椎关节,均设置为默认硬接触。

模型网格划分:皮质骨、松质骨按照实际厚度建立分别采用三维4节点(C3D4) 及三维6节点(C3D6)实体单元建立;同时采用三维减缩积分壳(S4R)单元模拟齿突横韧带,采用双节点SPRINGA弹簧单元模拟其余12种韧带。椎间盘(含髓核)采用辅以沙漏控制的C3D8R单元建模。

颈椎各结构材料属性:椎体(皮质骨、松质骨)以及横韧带，均设置为正交各向异性材料[11-12](见表1)；椎间盘(包括髓核及纤维环)，参考相关文献[13](见表2)；韧带材料，除横韧带之外，其余韧带应用弹塑性材料进行设置(见表3)。

表1 颈椎各椎体及横韧带材料属性

Tab.1 The material properties of conical bone，trabecular bone and transverse ligament

材料弹性模量/MPa泊松比单元数椎体皮质骨15×1030.3112486松质骨0.5×1030.3398023终板0.5×1030.414182韧带横韧带E1=86,E2=60.01

表2 椎间盘属性

设定模型载荷与边界条件:在T1椎下终板约束全部6个自由度作为边界条件来模拟加载。选取

表3 颈椎十二条韧带拟合数据

C0旋转轴上的点设置为加载参考点,将其与椎体上终板节点建立起Distribution Coupling约束,然后根据右手准则在中性点上施加1.5 Nm不同方向的纯扭矩载荷,方向设定上分别为X、Y、Z全局坐标(即:水平面为X-Y平面、冠状面为X-Z平面、矢状面为Y-Z平面)。分别用以模拟颈椎有限元模型受到的6种载荷作用(前屈、后伸、左右侧弯和左右轴向扭转)。

1.4 正常全颈椎有限元模型的对比验证

在载荷及边界条件设定相同的情况下,对正常颈椎模型进行模拟计算,将前屈、后伸、左右侧弯和左右旋转6个活动情况下的相对活动度(range of motion,ROM)与Panjabi等前期所做的实验数据[14]进行比较分析,比较其差异性并对模型进行验证。

1.5 颈椎棘突骨折模型与内固定模型的制备

在正常C0-T1模型的基础上,在C7预期骨折部位,删除对应骨折线的部分皮质骨和松质骨单元,使得棘突骨折块与椎体两部分完全分离,即保证两者之间没有任何节点和单元相连。同时,保留C7与C6及T1之间的棘间韧带和棘上韧带。用SolidWorks建立自行设计的颈椎直型接骨板几何模型,通过在Hypermesh中测量棘突椎板的大致弧度,在SolidWorks中进行预弯;并通过移动命令,移动棘突接骨板至实际手术时的摆放位置,完成接骨板与颈椎的装配。螺钉按照板的弯曲弧度,在植入处垂直于此处的切线方向植入。螺钉与板的相互关系为绑定,即与板接触的部分采用绑定接触关系,分析过程中螺钉与板不会分离。

1.6 模型加载及运算

本研究在三维模型C0上加载 1.5 Nm 的载荷,运行Abaqus6.9 进行有限元计算分析,测量各个方向上的活动度,然后比较相同扭矩下骨折模型及内固定模型与正常模型在屈曲、伸展、左右侧弯和左右旋转工况下颈椎活动度(ROM)。运算并显示模型 Von Mises 应力云图,显示内固定钢板及椎体应力变化。

2 结果

2.1 正常有限元模型的对比验证

建立的正常颈椎有限元模型、骨折和内固定模型见图1～3,正常颈椎有限元模型共包含535 721单元,171 392个节点。正常全颈椎模型与Panjabi等[14]的文献报告结果相比较,本研究建立的正常颈椎有限元模型的生物力学相似度良好,验证结果见表4。

图1 正常全颈椎有限元模型。(a) 前面观；(b) 后面观； (c) 矢状面Fig.1 The three-dimensional finite element model of the normal cervical spine. (a) Norma anterior (b) Norma posterior (c) Sagittal plain

表4 正常全颈椎模型不同工况下各椎节活动度同Panjabi的实验数据对比验证

注：在屈伸、侧弯及旋转工况下，本研究正常模型与Panjabi实验所得的实测结果基本相符，本研究各个节段大多落在Panjabi实验测得的活动度范围内。

Note：Under the conditions of flexion, extension, lateral bending and rotation, the experimental results of this paper are consistent with the experimental results of Panjabi, and most of the segments in this paper fall within the range of activity measured by Panjabi.

图2 颈椎棘突骨折(累及椎板型)有限元模型Fig.2 The finite model of cervical spinous process fracture(extension into the lamina)

图3 加载内固定系统有限元模型。 (a) 侧面观； (b)后面观； (c)装配螺钉示意Fig.3 The finite model of fixation system (a) Lateml view (b) Norma posterior (c) Fixation system extended the screw

2.2 骨折模型、内固定模型同正常模型之间的活动比较

加载1.5 Nm 载荷后正常颈椎模型和棘突骨折、手术内固定3个模型在不同颈椎活动下的活动度 (ROM) 详见表5。

1)骨折模型在C0-C1、C1-C2、C2-C3、C3-C4、C4-C5、C5-C6节段的6个活动度与健康模型相差无几;在C6-C7节段,骨折模型的前屈、侧弯及扭转活动度比正常模型略减小,而后伸活动度略微增大。在C7-T1节段,除后伸活动外,其余活动度均比正常模型要大,且活动度变化范围较为显著(单侧旋转活动度增加36.7%)。

2)在C5-C6节段,内固定模型比正常模型和骨折模型在旋转的活动度略上升,侧弯活动度下降,前屈和后伸仍然保持不变。在C6-C7节段,内固定模型活动度在后伸、侧弯比正常模型略微下降,旋转显著下降。在该节段与骨折模型相比,内固定模型在后伸时降低,在前屈和旋转时变大,在侧弯时基本相同。在C7-T1节段,植入钢板模型活动度比正常和骨折模型都低。整体C0-T1活动度方面,植入钢板模型结果比正常活动度偏小,整体稳定性增加。

2.3 内固定模型颈椎结构的应力变化及内固定装置的安全应力验证

骨折后颈椎应力改变见图4,内固定重建后颈椎及内固定Von Mises应力分布见图5。重建后内固定Von Mises应力主要集中于靠近椎体的两个螺钉孔周围,成同心圆状分布。在屈曲工况下,C7高应力区主要在椎体的椎弓根,以及打钉处的椎板周围,内固定Von Mises应力峰值为452.5 MPa;在后伸工况下,C7受力主要集中在椎体的椎弓根,内固定Von Mises应力峰值为715.2 MPa;在左右侧弯工况下,C7应力集中在同侧椎弓根,内固定Von Mises应力峰值为391.0 MPa;在左右旋转工况下,C7高应力区集中在同侧椎体,内固定Von Mises应力峰值为544.0 MPa。钛合金的屈服应力为880～950 MPa,此分析中最大应力为715.2 MPa,发生在后伸活动下,但仍在钛合金的弹性范围内,因此植入物是安全的。但疲劳强度一般低于屈服强度，持续疲劳应力仍存在断板可能，需进一步生物力学验证。

表5 颈椎棘突骨折模型、内固定模型与正常模型C6-C7、C7-T1节段及全颈椎活动度比较

Tab.5 The ROM of C6-C7、C7-T1 segment and total cervical spine under all conditions in FEM/fracture and FEM/fixation and FEM/intact model

正常模型骨折模型内固定模型C6-C7C7-T1totalC6-C7C7-T1totalC6-C7C7-T1total前屈+后伸/(°)6.017.1199.495.969.20101.196.034.0796.03左右侧弯/(°)5.284.9266.425.045.8367.255.152.2162.81左右旋转/(°)6.409.59142.474.8013.12144.405.712.91133.75

3 讨论

铲土者骨折是一类应力所致棘突撕脱骨折,最早见于铲土者工人。当工人用力挥动铁铲时,常常由于背部肌肉瞬间用力,作用力较大,加上棘突骨块相对薄弱,故可引起下颈椎以及上胸椎的棘突撕脱性骨折[15]。随着工业化社会的到来,这类由体力劳动引起的单纯棘突撕脱骨折越来越少见,但铲土者骨折这个概念却被沿用至今,人们仍然习惯性地将发生在颈椎棘突的骨折都称为铲土者骨折,然而颈椎棘突骨折和铲土者骨折并不是完全相等的两个概念,铲土者骨折这个概念是被误用及泛用的,颈椎棘突骨折除了体力劳动所致单纯棘突撕脱骨折外,也包括其他致伤因素(如直接打击等)引起的颈椎棘突骨折,但通常因为致伤暴力较大,并不能单一当作临床意义较小的铲土者骨折来诊断及治疗。

颈椎棘突骨折目前公认的治疗方法是保守治疗,通过颈围固定4～6周疼痛会显著减轻[16],3个月左右可以进行正常的日常活动,1年左右可以获得一个稳定的临床康复。以前临床医生对此类骨折患者采取手术治疗通常是因为残留疼痛和骨折不愈合等因素[17-18],但骨折后并发症(颈椎不稳及脊髓损伤)是目前棘突骨折尽早行手术治疗的主要原因。Matar LD[1]和泰德安[2]等认为累及椎板的颈椎棘突骨折因致伤暴力相对较大,引起骨质较致密的棘突椎板交界处断裂,破坏相应节段的后方韧带复合体(项韧带、棘间韧带及黄韧带)和双侧关节突关节损伤移位,破坏局部稳定性,存在潜在不稳的可能。同时,Mirvis等[19]认为累及椎板骨折存在脊髓损伤的可能性,故特将此型骨折称之为非典型颈椎棘突骨折。一旦发现颈椎不稳和(或)脊髓损伤,尽早手术治疗可以有效避免严重并发症发生,增加颈椎整体稳定性。

本研究通过模拟其他学者的建模方法[6-7],在计算机上运用多个建模软件进行全颈椎的模型重建,建模完成后并对模型的测试数据与文献报告数据进行有效性验证,通过验证并进行下一步有限元分析。本研究建立的颈椎C7棘突骨折(累及椎板)以及微型接骨板内固定有限元模型较好地模拟了临床骨折及实际手术病例,并进行活动度比较和应力分析。骨折模型在一定程度上反映了颈椎棘突骨折(累及椎板)潜在不稳的可能性,临床上需临床医生结合体征及影像学资料做出准确具体的判断,针对损伤较为严重的骨折类型(累及椎板的棘突骨折),一旦发现颈椎不稳和(或)脊髓损伤,需及时给予手术固定等治疗方法避免出现严重并发症。

笔者自行设计的颈椎棘突微型接骨板在有限元分析中体现了较好的骨折部位适应性,且应力分析在安全范围内,内固定方式安全有效,得到了与临床一致的结果。本研究并未纳入对内固定疲劳强度的安全性验证，疲劳强度一般低于屈服强度,从本研究来看，后伸和旋转载荷下内固定仍存在断板可能,故临床上应颈托固定4～6周,限制病人旋转及后伸活动。下阶段笔者将联合医疗器械公司对自行设计的颈椎棘突微型接骨板进行进一步安全性验证，包括体外强度测试，并期待进入临床应用。

本研究的不足之处：在研究中使用的是有限元分析方法，即使目前的技术能非常逼真地反映颈椎的生物力学特性，但并不能像真正的人体试验一样具有完全一致的生物特性，在建模和分析过程中都存在一定的实验误差。本研究主要是通过有限元

图4 颈椎棘突骨折后各椎体及C7单椎体应力云图。 (a) 后伸位椎体应力云图； (b) 后伸位C7单椎应力云图； (c) 前屈位椎体应力云图； (d) 前屈位C7单椎应力云图； (e) 侧屈位椎体应力云图； (f) 侧屈位C7单椎应力云图； (g) 旋转位椎体应力云图； (h) 旋转位C7单椎应力云图Fig.4 The Von Mises stress distribution of the cervical spine and C7 vertebral body under cervical spinous process fracture. (a) The Von Mises stress distribution of the cervical spine under flexion； (b) The Von Mises stress distribution of C7 vertebral body under flexion； (c) The Von Mises stress distribution of the cervical spine under extension； (d) The Von Mises stress distribution of C7 vertebral body under extension； (e) The Von Mises stress distribution of the cervical spine under lateral bending； (f) The Von Mises stress distribution of C7 vertebral body under lateral bending； (g) The Von Mises stress distribution of the cervical spine under rotation； (h) The Von Mises stress distribution of C7 vertebral body under rotation

图5 颈椎棘突骨折内固定术后各椎体及内固定装置应力云图。 (a) 后伸位椎体应力云图； (b) 后伸位内固定应力云图； (c) 前屈位椎体应力云图； (d) 前屈位内固定应力云图； (e) 侧屈位椎体应力云图； (f) 侧屈位内固定应力云图； (g) 旋转位椎体应力云图； (h) 旋转位内固定应力云图Fig.5 The Von Mises stress distribution of the vertebral body and internal fixation devices of the model of cervical spinous process fracture condition. (a) The Von Mises stress distribution of the vertebral body under flexion； (b) The Von Mises stress distribution of internal fixation devices under flexion； (c) The Von Mises stress distribution of the vertebral body under extension； (d) The Von Mises stress distribution of internal fixation devices under extension； (e) The Von Mises stress distribution of the vertebral body under lateral bending； (f) The Von Mises stress distribution of internal fixation devices under lateral bending； (g) The Von Mises stress distribution of the vertebral body under rotation； (h) The Von Mises stress distribution of internal fixation devices under rotation

方法比较同一个体在骨折前后的活动度变化，属于配对设计，活动度验证本模型同文献数据之间的差异并没有对结果造成明显影响，因此可以忽略活动度验证对结果的偏差。

4 结论

保守治疗并不是所有颈椎棘突骨折的通用治疗方式,临床工作者需结合临床实际和患者要求决定棘突骨折治疗方案,手术内固定治疗颈椎棘突骨折(累及椎板)为一种安全有效的治疗方式,可应用于实际骨折病例。

(致谢：感谢李强教授及解放军98医院骨五科、放射科的医生们提供的有力帮助，感谢上海瑞金医院骨研所提供有限元分析平台)

[1] Matar LD,Helms CA, Richardson WJ. “Spinolaminar breach”:an important sign in cervical spinous process fractures[J].Skeletal Radiol,2000,29(2):75-80.

[2] 秦德安,张佐伦,崔新刚，等.颈椎棘突骨折与颈椎稳定性[J].中国骨与关节损伤杂志,2006,9(21): 687-689.

[3] Han SR, Sohn MJ. Twelve contiguous spinous process fracture of cervico-thoracic spine[J].Korean J Spine,2014, 11(3):212-213.

[4] Kazanci A,Gurcan O,Gurcay AG, et al. Six-level isolated spinous process fracture of the thoracic vertebrae (clay-shoveler′s fracture) and a review of the literature[J]. Neurol India, 2015,63(2):223-224.

[5] Olivier EC,Muller E,Janse Van Rensburg DC.Clay-shoveler fracture in a paddler: A case report[J].Clin J Sport Med, 2016,26(3):69-70.

[6] 陈金水,倪斌,陈博,等.寰枢椎脱位三维非线性有限元模型的建立和分析[J].中国脊柱脊髓杂志,2010,20(9):749-753.

[7] 郭群峰,陈方经,倪斌,等.带有颅底的全颈椎三维有限元模型的建立及分析[J].中国脊柱脊髓杂志,2014,24(6):550-554.

[8] 蔡贤华,王威,王志华,等.不同前路内固定方式治疗枢椎椎体横行骨折合并Hangman骨折稳定性的有限元分析[J].中国脊柱脊髓杂志,2014,24(3):257-265.

[9] 柳超,王前,张杰峰,等.不同内固定手术方式治疗寰枢椎复合骨折稳定性的有限元分析[J].中国脊柱脊髓杂志,2015,25(10):904-911.

[10] 贾连顺, 袁文. 现代脊柱外科学[M].北京:人民军医出版社,2007,21-24.

[11] Brolin K, Halldin P.Development of a finite element model of the upper cervical spine and a parameter study of ligament characteristics[J].Spine (Phila Pa 1976),2004,29(4):376-385.

[12] Shirazi-Adl A,Ahmed AM, Shrivastava SC. Mechanical response of a lumbar motion segment in axial torque alone and combined with compression[J].Spine (Phila Pa 1976), 1986,11(9):914-927.

[13] El-Rich M,Arnoux PJ,Wagnac E, et al.Finite element investigation of the loading rate effect on the spinal load-sharing changes under impact conditions[J].J Biomech,2009,42(9):1252-1262.

[14] Panjabi MM,Crisco JJ,Vasavada A,et al.Mechanical properties of the human cervical spine as shown by three-dimensional load-displacement curves[J]. Spine (Phila Pa 1976), 2001,26(24):2692-2700.

[15] 李连华, 孙天胜.什么是Clay shoveler骨折?[J].中国脊柱脊髓杂志,2010,20(6):528.

[16] Umredkar A,Sura S, Mohindra S. Multiple contiguous isolated spinous process fracture (Clay-Shoveler′s fracture) of the cervicodorsal spine[J]. Neurol India,2011,595:788-789.

[17] Gershon-Cohen J, Glauser F.Whiplash fractures of cervicodorsal spinous processes: resemblance to shoveler′s fracture[J]. J Am Med Assoc,1954,1556:560-561.

[18] Hirsh LF,Duarte LE,Wolfson EH, et al. Isolated symptomatic cervical spinous process fracture requiring surgery. Case report[J].J Neurosurg,1991,75(1):131-133.

[19] Mirvis SE, Shanmuganathan K. Trauma radiology: Part V. Imaging of acute cervical spine trauma[J]. J Intensive Care Med,1995,10(1):15-33.

Finite Element Analysis of Internal Fixation in Cervical Spinous Process Fracture

Shen Yan1Wang Chaoyang1*Wu Rong1Chen Bo2

1(DepartmentofOrthopedics,ThePLA98thHospital,Huzhou313000,Zhejiang,China)2(DepartmentofOrthopedicsInstitute,RuijinHospital,ShanghaiJiaotongUniversitySchoolofMedicine,Shanghai200025,China)

The aim of this work is to establish a three-dimensional finite element model of cervical spinous process fracture (extension into the lamina) after internal fixation and assure its availability to surgical treatment. Based on a finite element model of a normal cervical spine, a finite element model of cervical spinous process fracture (extension into the lamina) was developed according to the clinical case, with which an internal fixation therapy on fracture model was simulated. The range of motion (ROM) under flexion, extension, lateral-bending and axial rotation were measured and analyzed in the normal and fracture model and fixation model. The force loading was applied on both vertebra and internal fixation to elucidate the safety of this surgical therapy. It was shown that the finite element model of cervical spinous process fracture (extension into the lamina) had a high similarity and profile to the clinical case. The range of motion (ROM) on C7-T1 segment under each movement in fracture model (flexion+extension 9.20°, lateral-bending 5.83°, axial rotation 13.12°) was larger than that in the normal model(flexion+extension 7.11°, lateral-bending 4.92°, axial rotation 9.59°), especially in the rotation movement. The whole cervical vertebra was more stable under internal fixation (flexion+extension 4.07°, lateral-bending 2.21°, axial rotation 2.91°) with the safety of internal fixation system. In conclusion, the finite element model could be used to simulate the biomechanics of cervical spinous process fracture (extension into the lamina) to indicate the potential for delayed instability.

cervical spinous process fracture; internal fixation; finite element analysis

10.3969/j.issn.0258-8021. 2017. 04.010

2016-05-25， 录用日期:2017-02-27

南京军区医学科技创新项目(MS009)；湖州市自然科学基金(2016YZB03)

R318

A

0258-8021(2017) 04-0456-08

*通信作者(Corresponding author)，E-mail: wcyang1228@163.com