波致Gibson粉土质海床累积孔压响应的简化分析*

张 琪, 徐继尚

(中国海洋大学海洋地球科学学院，海底科学与探测技术教育部重点实验室，山东 青岛 266100)

研究简报

波致Gibson粉土质海床累积孔压响应的简化分析*

张 琪, 徐继尚**

(中国海洋大学海洋地球科学学院，海底科学与探测技术教育部重点实验室，山东 青岛 266100)

基于累积孔隙水压的控制方程，采用有限差分法求解了波浪作用下Gibson粉土质海床的累积孔隙水压。首先采用有限差分法求解了均匀及双层海床的累积孔隙水压，通过与解析解对比，验证了该方法的准确性。其次针对Gibson粉土质海床累积孔压计算量大的缺点，提出了利用等效替代法获取简化土体性质参数的思路，并设计算例予以验证。计算结果表明：在一定的时间尺度内，可将Gibson土简化为双层土体，简化运算减少了孔压计算量，方便工程应用。

Gibson海床；固结；累积孔隙水压；有限差分法；等效替代法

通过现场观测和试验研究，目前认为波浪作用下海床的动力响应可根据孔隙水压生成机制的不同分为两类：瞬时孔隙水压和累积孔隙水压。瞬时孔隙水压通常呈周期性变化，一般存在振幅衰减和相位滞后[1]。累积孔隙水压指循环剪应力导致松散的海床土体颗粒压缩，造成的孔压随时间增长的现象[2]。与地震作用下的土体液化类似，波致累积孔隙水压通常发生在荷载的初始阶段，如果孔压不能及时消散会引起海床土体的液化破坏。

现有对累积孔隙水压的研究主要是根据不排水循环三轴试验建立孔隙水压增长模式，通过求解修正的Biot固结方程，得到孔压解析解[3]和数值解[4]。 Seed[4]构建了一维有限元模型用于计算海床土的累积孔隙水压，并以1964年日本新潟地震的某一典型砂层为例估算了孔隙水压力的变化过程。McDougal[5]在求解Biot固结方程时对累积孔压比和循环应力比做了线性近似，得到了适用于不同厚度海床的累积孔压解析解。Sumer[6]给出线性推进波作用下均匀海床累积孔隙水压的解析解，并提出了一种预测海床累积孔压的随机游走模型。Cheng[7]采用有限差分法求解了推进波作用下均匀土体的累积孔隙水压，发现海床厚度较大时，孔压对表层土剪切力比较敏感。Jeng[8]利用拉普拉斯变换得到了波浪荷载下无限厚海床累积孔压的简化解，通过变动参数分析发现，相对水深减小或波陡增大时，海床更易发生累积液化。

然而，以上研究针对的都是均匀海床，而实际海床通常是不均匀的。例如，对具有一定沉积历史的海床，固结过程的增密效应和上层土体的压实作用往往导致剪切模量及渗透系数随深度连续变化[9]。因此，工程中用到的Gibson土可能更符合实际情况。Gibson土是指剪切模量及渗透系数随深度连续变化的非均质弹性半空间土体模型[10]。目前关于海床累积液化的研究主要针对的是均匀或双层土体，对诸如Gibson土等非均匀土体的研究较少。刘占阁[11]运用分离变量法和格林函数得到了推进波作用下三层海床累积孔压解析解，不过在分析土况参数对累积孔压的影响时仍将海床视为双层土体。ZHU Jieran[12]采用椭圆余弦波理论计算浅水区域海床面的波压力，引入双曲线模型来描述累积孔压比和循环应力比的关系，进而分析了双层海床的累积孔压响应。

本文采用有限差分法求解了线性推进波作用下Gibson粉土质海床的累积孔隙水压，通过算例分析得到了与Gibson土孔压变化规律相似的简化土体的性质参数，这样在求解Gibson海床累积孔隙水压时可用简化土体近似代替它，方便了工程应用。

1 数值模型

1.1 控制方程

假定土骨架和孔隙流体是可压缩的，土骨架应变足够小，且应力应变关系符合Hooke定律，土中渗流满足Darcy定律。根据应力平衡方程、弹性本构关系及流体连续性条件，Biot[13]建立了三维Biot固结方程。该方程能准确反映孔隙水压消散与土骨架变形的耦合关系，又称为“真三维固结理论”[14]。

虽然Biot固结理论与实际情况接近，不过却不便计算[15]。为此，本文仅考虑渗流发生在垂向的情形。当海床厚度大于等于波长的一半，即h≥L/2时，多孔饱和弹性海床瞬时孔隙水压力Posc和剪应力τxz与土况参数无关，表述如下：

Posc=Pb×exp(-kz)×cos(kx+wt)，

(1)

τxz=Pb×kz×exp(-kz)×cos(kx+wt)。

(2)

为了确定波浪作用下海床土孔隙水压发展规律，需要建立不排水条件下孔隙水压力发展模式。现将一维Biot固结理论和不排水条件下的孔压增长模式相结合，得到海床累积孔压Pres的控制方程：

(3)

式中cv(z)为土体深度z处的固结系数，由下式给出：

(4)

其中：G(z)为剪切模量；K(z)为渗透系数；rw为水容重；v为泊松比。公式(3)中，f(z)为累积孔压源项，通常与土体相对密度、初始循环剪应力比等有关，表征了单位时间单位孔隙体积内孔压的生成速率[16]，即：

(5)

(6)

(7)

式中：N为波浪加载次数；Nl为引起土体液化的波浪加载次数，是循环剪应力比的函数：

(8)

α和β是与土体相对密度(Dr)有关的经验参数，本文参考McDougal[5]表1和图3中的土况参数，取α=0.246，β=-0.165，Dr=0.54。当海床厚度不小于半波长时，

f(z)=Azexp(-λz)，

(9)

(10)

(11)

其中：r′为土体浮容重；rs为土容重；K0为土体静压力系数。

1.2 初始及边界条件

假定海水表面存在线性推进波，海床表面水平，为透水边界；海床底部为不透水边界(见图1)。

求解方程 (3)需要给出初始条件和边界条件，由于海床表面为完全透水界面，累积孔压始终为0，即：

(12)

海床基底为不透水边界，垂向孔压梯度为0，即：

(13)

初始时刻，海床内无孔压累积，即：

(14)

在深度z=m处上下分别选取厚度为Δm的土体微元，当Δm很小时，该微元体可视为均质体，渗透系数为K(m)，固结系数为cv(m)，累积孔压为Pres(z,t)；同样深度m至m+Δm处土体微元的渗透系数为K(m+Δm)，固结系数为cv(m+Δm)，累积孔压为Pres′(z,t)。这样，深度z=m-Δm至z=m+Δm的Gibson土可视为固结系数不同的双层土体。由达西定律和流速连续性条件有：

(15)

1.3 构建差分格式

对深度z=m-Δm至z=m的土体，一维Biot固结方程(3)连同边界条件和初始条件(等式(12)-(15))可利用有限差分法求解。本文采用经典显式差分，差分格式构造如下：

(16)

为了提高边界节点的误差阶，对两个一阶中心差商

(17)

(18)

(19)

化简，得：

(20)

再将(17)式和(18)式代入(20)式，有：

(21)

(22)

这样，(15)式变为：

(23)

特别地，当z=h时，垂向孔压梯度为0：

(24)

2 计算结果验证

2.1 均匀海床累积孔隙水压有限差分解与解析解的对比

Clukey[19]根据水槽实验测量了波浪荷载下饱和粉土质海床的累积孔隙水压，本文为了验证均匀海床累积孔压有限差分解的准确性，在设计算例时参考其设置的波浪和土况参数，具体为：线性推进波，周期T=1.76s，波长L=3.473m，水深d=0.5m，波高H=0.22m；均匀海床，泊松比u=0.49，剪切模量G=5.6×105N/m2，K=4×10-8m/s，水容重rw=9806N/m3，土体浮容重r′=8500N/m3，静压力系数K0=0.4，海床厚度h=2m，经验系数α=0.246，经验系数β=-0.165。

考虑到孔压计算模型仅适用于厚度大于半波长的海床，将Clukey水槽实验中的海床厚度h设置为2 m。参考Jeng[8]求取累积孔隙水压使用的解析解方法，将有限差分法得到的累积孔压数值解及解析解对比结果如图2。

2.2 双层海床累积孔隙水压有限差分解与解析解的对比

对双层海床，参考刘占阁[11]设置的波浪参数，设计算例，具体为：线性推进波，周期T=7s，波长L=40m，水深d=4m，波高H=2.75m；上层海床，泊松比u=0.49，剪切模量G=52 353 N/m2，K=0.000 1m/s，水容重rw=10 000N/m3，土体浮容重r′=10000N/m3，静压力系数K0=0.4，海床厚度h=10m，经验系数α=0.246，经验系数β=-0.165；下层海床，泊松比u=0.49，剪切模量G=26 177N/m2，K=0.00005m/s，水容重rw=10000N/m3，土体浮容重r′=10000N/m3，静压力系数K0=0.4，海床厚度h=20m，经验系数α=0.246，经验系数β=-0.165。

将累积孔压有限差分数值解与刘占阁的解析解对比结果如图3。

图2和3显示，在不同时刻均匀及双层海床累积孔压的数值解与解析解吻合，说明利用有限差分法求解此类土体的累积孔压是可行的。

3 Gibson粉土质海床累积孔隙水压的简化分析

由于现阶段尚未对波浪作用下Gibson海床的累积孔隙水压响应展开深入研究，缺乏相应的数值解和解析解。因此，工程上在求取该类土体的累积孔压时，只能将其简化为均匀或双层土体，虽然计算简单，不过简化土体的土况与Gibson土差异较大，孔压解可能不太准确。考虑到对实际三维Gibson海床，当海床较厚或土性参数在横向变化剧烈时，采用有限差分法求解累积孔压需要将土体划分为极为精细的网格，这将极大地增加孔压计算量，不便于工程应用。本节基于Gibson海床累积孔压的有限差分解，采用等效替代法得到了简化土体的性质参数，只需获取Gibson海床表面的土性参数即可判断累积孔压的大小及变化规律，既保证了孔压解的准确性，又简化了运算，可为工程上研究此类土体的累积孔压变化规律提供参考。

Jeng[20]基于Byrant[21]和Suzuki[22]取自墨西哥湾的海底沉积物样品所做的固结实验，拟合得到了Gibson海床渗透系数及剪切模量随深度变化的函数，如式(25)和(26)所示：

K(z)=Ks·exp(-4.69z/h)，

(25)

G(z)=Gs·(1+15z/h)。

(26)

式中，Ks和Gs分别表示海床表面的渗透系数和剪切模量，由于Gibson海床的渗透系数及剪切模量均随深度变化，直接确定简化土体的性质参数有些困难。为此，本文首先对渗透系数或剪切模量随深度变化的土体做了简化，进而得到Gibson海床的简化土体的性质参数，现考虑以下3种土体的累积孔压响应：

土体a：渗透系数随深度指数减小，剪切模量不变；

土体b：渗透系数不变，剪切模量随深度线性增大；

土体c：渗透系数随深度指数减小，剪切模量随深度线性增大。

王立忠[23]在研究波浪作用下砂质及粉质海床的孔压响应时发现，砂质海床不会出现孔压累积现象，而粉质海床的孔压累积现象明显。因此，本文在设计土况参数时也只考虑粉质海床的情况，参考王立忠[23]水槽实验选用的粉土海床参数以及王栋[24]在分析波浪作用下海床动力响应时用到的线性推进波参数，设计算例A，具体为：线性推进波，周期T=5.49 s，波长L=40 m，水深d=8 m，波高H=1 m；Gibson海床，泊松比u=0.49，海床面剪切模量Gs=1.1×106N/m2，海床面渗透系数Ks=2.1×10-7m/s，水容重rw=10 000 N/m3，土体浮容重r‵=10 000 N/m3，静压力系数K0=0.4，海床厚度h=24 m，经验系数α=0.246，经验系数β=-0.165。

3.1 土体a累积孔隙水压响应的简化分析

对土体a，假定可以将其等效为均匀土体d，该土体的渗透系数为K，剪切模量与土体a相同，如果它具有与a相似的累积孔压变化规律，那么在分析复杂土体a的孔压响应特点时，可用均匀土体d代替它。图4给出了对土体d取不同的渗透系数，两种土体在波浪作用3 000T时累积孔压随深度的对比曲线，不难看出当K=0.75Ks时，土体d的累积孔压幅值及孔压随深度的变化规律与a最接近。图5显示，如果波浪作用时间不超过3 000T，渗透系数为0.75Ks的土体d的累积孔压幅值与a相差较小，且两种土体孔压随深度的变化趋势相同。

为了探究在不同的波浪条件下，当土体a的性质参数发生变化时，渗透系数为0.75Ks的土体d与a的累积孔压是否接近。参考刘红军[25]研究黄河三角洲土体强度时所使用的波浪及粉土质海床参数，设计算例B，具体为：线性推进波，周期T=7 s，波长L=46.16 m，水深d=5 m，波高H=1 m；Gibson海床，泊松比u=0.33，海床面剪切模量Gs=107N/m2，海床面渗透系数Ks=10-8m/s，水容重rw=10 000 N/m3，土体浮容重r‵=10 000 N/m3，静压力系数K0=0.4，海床厚度h=30 m，经验系数α=0.246，经验系数β=-0.165。

现作出当K=0.75Ks时，土体a和d在算例B所示的波浪及土况条件下累积孔压随深度的对比曲线：

可见，对不同的波浪和土况条件，当K=0.75Ks时，如果波浪作用时间不超过3 000T，土体d的累积孔压幅值与a相差较小，且孔压随深度的变化趋势与a相同，均表现出随深度先增大后减小的特征。因此，可用渗透系数为0.75Ks的均匀土体来代替它。

3.2 土体b累积孔隙水压响应的简化分析

通过对比均匀土体与土体b累积孔压随深度的分布曲线发现，两种土体累积孔压相差较大。因此，均匀土体可能不适合作为土体b的等效土体。本文假定可以将土体b等效为双层土体e，土体e的渗透系数与b相同，其上层土的剪切模量为G1，下层土的剪切模量为G2。大量的数值实验表明，对算例A所示的波浪及土况条件，当G1=1.4Gs，G2=6.5Gs，上层海床厚度为5 m时，如果波浪作用时间不超过3 000T，土体e和b在不同时刻的累积孔压变化规律相同，且孔压值相差较小(见图7、8)，因此，可将土体b简化为双层土体。

为了探究在不同的波浪条件下，当土体b的性质参数发生变化时，是否仍可将其简化为上层渗透系数为1.4Gs，下层渗透系数为6.5Gs，上下层分界深度为5 m的双层土体。参考算例B的波浪及土况参数，绘制两种土体累积孔压随深度的对比曲线(见图9)。

图9显示，对算例B所示的波浪及土况条件，当上层海床厚度为5m，G1=1.4Gs，G2=6.5Gs时，如果波浪作用时间不超过3 000T，土体b与e的累积孔压值接近，且两种土体具有相同的孔压变化规律。试验证明，可对土体b做这种简化分析。

3.3 土体c累积孔隙水压响应的简化分析

基于得到的土体a和b的简化土体的性质参数，在讨论土体c累积孔压的变化规律时，假定能够选用某一双层土体f来替代它，且f满足以下3个条件：

(1)上下层的渗透系数相同，其值为土体c海床面渗透系数的0.75倍，即K=0.75Ks。

(2)上下层剪切模量分别为土体c海床面剪切模量的1.4倍和6.5倍，即G1=1.4Gs，G2=6.5Gs。

(3)上下层分界深度为5 m。

为了验证上述假设是否合理，参考算例A和B中的波浪及土况参数，绘制两种土体累积孔压随深度的对比曲线，如下：

可见，当双层土体f满足上述假设时，对不同的波浪及土况条件，如果波浪作用时间小于3 000T，土体f与c的累积孔隙水压值相差较小，且孔压随深度的变化规律与c相同。因此，可以将土体c简化为双层土体，且该双层土体的性质参数满足上述3个条件。当波浪作用时间大于3 000T，虽然土体f与c累积孔压随深度的变化规律相同，不过两种土体孔压相差较大。图11显示当波浪作用时间为8 000T，对算例A所示的波浪及土况条件，6 m以深处土体f累积孔压随深度的分布曲线明显偏离土体c。因此，不宜选用土体f作为c的等效土体。

4 结语

本文提出并验证了利用等效替代法计算线性推进波作用下Gibson粉土质海床累积孔隙水压的思路，验证结果表明：在3 000个波周期内，对渗透系数和剪切模量分别满足(25)式和(26)式的Gibson海床，可将其简化为双层土体，该土体上下层的分界深度为5 m；双层土体上下层的渗透系数相同，均为海床面渗透系数的0.75倍；上下层土体的剪切模量分别为海床面剪切模量的1.4和6.5倍。等效替代法的优势在于只需获取Gibson海床表面处的土性参数即可预测不同深度位置的累积孔压，因而能极大地减少孔压运算量，方便工程应用。

需要说明的是，在3 000个波周期内，Gibson土与双层土体的累积孔压随深度的变化规律相同，且孔压幅值相近，而当波浪作用时间超过3 000个波周期时，两种土体孔压幅值相差较大，因此，只有当波浪作用时间小于3 000个波周期时才能对Gibson土做上述等效替代。另外，上述等效替代对渗透系数和剪切模量满足(25)式和(26)式的Gibson土体是有效的，尚不清楚是否可以将其他土体简化为双层土体，这也是论文下一步要做的工作。

[1] Madsen O S. Wave induced pore pressures and effective stresses in a porous bed [J]. Geotechnique,1978, 28(4): 377-393.

[2] Nago H, Maeno S, Matsumoto T, et al. Liquefaction and densification of loosely deposited sand bed under water Pressure variation[C]. Singapore: Proceeding of the Third International offshore and Polar Engineering Conference， 1993: 578-584.

[3] Sumer B M, Fredsøe J. The mechanics of scour in the marine environment [M]. Singapore: World Scientific, 2002.

[4] Seed H B, Rahman M S. Wave induced pore pressure in relation to ocean floor stability of cohesionless soils [J]. Marine Geotechnology, 1978, 3(2): 123-150.

[5] Mcdougal W G, Tsai Y T, Liu P L F, et al. Wave induced pore water pressure accumulation in marine soils [J]. Journal of Offshore Mechanics and Arctic Engineering, 1989, 111: 1-11.

[6] Sumer B M, Cheng N S. A random-walk model for pore pressure accumulation in marine soils[C]. Brest: The 9th International Offshore and Polar Engineering Conference， 1999: 521-526.

[7] Cheng L, Sumer B M, Fredsøe J. Solution of pore pressure build up due to progressive waves [J].International Journal for Numerical and Analytical Methods in Geomechanics, 2001, 25(9): 885-907.

[8] Jeng D S, Seymour B R. Simplified analytical approximation for pore-water pressure buildup in marine sediments [J]. Journal of Waterway Port, Coastal, and Ocean Engineering, 2007, 133(4): 309-321.

[9] Gibson R E. Some results concerning displacements and stresses in a non-homogeneous elastic half-space [J]. Geotechnique, 1967, 17(1): 58-67.

[10] 宋岩新． 砂质海床上海底管线稳定性的数值分析[D] ． 大连： 大连理工大学， 2008: 41-48. Song Yanxin. Numerical Analysis of Stability of Submarine Pipelines on a Sandy Seabed [D]. Dalian: Dalian University of Technology, 2008: 41-48.

[11] 刘占阁， 栾茂田， 王忠涛． 波浪作用下成层海床孔隙水压力响应的简化分析 [J]． 地震工程与工程振动， 2006,26(2): 156-160. Liu Zhan-ge, Luan Maotian, Wang Zhongtao. Simplified procedure for evaluation of residual pore pressure response in layered seabed under wave loading [J]. Earthquake Engineering and Engineering Dynamics, 2006, 26(2): 156-160.

[12] Zhu J R. Cnoidal water wave induced pore pressure accumulation in the Seabed [J]. Advances in Petroleum Exploration and Development, 2013, 6(1): 58-64.

[13] Biot M A. General theory of three-dimensional consolidation [J]. Journal of Applied Physics, 1941, 12(2): 155-164.

[14] 钱佳欢， 殷宗泽． 土工原理与计算 [M]． 北京： 中国水利水电出版社，1996. Qian Jiahuan, Yin Zongze. Principle and Calculation of Soil Engineering [M]. Beijing: China Water & Power Press, 1996.

[15] 黄杰卿． 非线性大应变固结理论与试验研究[D]． 杭州： 浙江大学， 2015: 20-23. Huang Jieqing. Theoretical and Experimental Study on Nonlinear Finite Strain Consolidation [D]. Hangzhou: Zhejiang University, 2015: 20-23.

[16] Sumer B M, Ozgur KIRCA V S, Fredsoe J. Experimental validation of a mathematical model for seabed liquefaction under waves [J]. International Journal of Offshore and Polar Engineering, 2012, 22(2): 133-141.

[17] Seed H B, Martin P O, Lysmer J. The Generation and Dissipation of Pore Water Pressure During Soil Liquefaction [M]. University of California: College of Engineering, 1975: 166-172.

[18] 李荣华， 刘播．微分方程数值解法 [M]． 北京：高等教育出版社，1998. Li Ronghua, Liu Bo. Numerical Methods of Differential Equations [M]. Beijing: Higher Education Press, 1998.

[19] Clukey E C. Laboratory and field investigation of wave-sediment interaction [R]. New York: Cornell University, 1983.

[20] Jeng D S, Cheng L. Wave-induced seabed response around a pipe laid on a pore-elastic seabed [J]. Journal of Offshore Mechanics and Arctic Engineering, 1999, 121(4): 227-236.

[21] Bryant W R, Hottman W, Trabant P. Permeability of unconsolidated and consolidated marine sediments [J]. Marine Geotechnology, 1975, 1(1): 1-14.

[22] Suzuki H, Ando K, Kitahara M, et al. Shear modulus profile measurement of the shallow-water seabed in Japan [C]∥ Yokahama: Process of International Conference on Geotechnical Engineering for Coastal Development-Theory and Practice on Soft Ground， 1991:3-6.

[23] 王立忠， 潘冬子， 潘存鸿， 等． 波浪对海床作用的试验研究[J]． 土木工程学报， 2007, 40(9): 101-109. Wang Lizhong, Pan Dongzi, Pan Cunhong， et al. Experimental investigation on wave-induced response of seabed [J]. China Civil Engineering Journal, 2007, 40(9):101-109.

[24] 王栋． 波浪作用下海床动力响应与液化的数值分析[D]． 大连： 大连理工大学， 2002: 54-58. Wang Dong. Numerical Analysis for Dynamic Response and Liquefaction Potential of Seabed under Wave Loading [D]. Dalian: Dalian University of Technology, 2002: 54-58.

[25] 刘红军， 张民生， 许国辉， 等． 波浪作用下海床土体强度非均匀化数值分析[J]．岩土工程学报， 2008, 30(6): 924-929. Liu Hongjun, Zhang Minsheng, Xu Guohui， et al. Numerical analysis of wave-induced non-uniformity of soil strength [J]. Chinese Journal of Geotechnical Engineering, 2008, 30(6): 924-929.

责任编辑 徐 环

Simplified Analysis of Wave-Induced Residual Pore Pressure Response in Gibson Silty Seabed

ZHANG Qi, XU Ji-Shang

(College of Marine Geosciences, Key Lab of Submarine Sciences & Prospecting Techniques, Ministry of Education, Ocean University of China, Qingdao 266100, China)

Finite difference method is employed to calculate the wave-induced residual pore pressure in Gibson silty seabed based on the governing equation of residual pore pressure. Finite difference method is employed to calculate the wave-induced residual pore pressure in the uniform and two-layer seabed, and analytic solutions are used to validate the method. Considering the large computation effort of the residual pore pressure in the Gibson silty seabed, equivalent substitution method is proposed to obtain the properties of simplified soil, and numerical examples are designed to prove the effectiveness of this method. Numerical results indicate that Gibson soil could be reduced to two-layer soil in a definite time scale, and equivalent substitution method has the advantage of decreasing the computation complexity, which is convenient for engineering application.

Gibson seabed; consolidation; residual pore pressure; finite-difference method; equivalent substitution method

国家自然科学基金项目(51479182；41006024) 资助

2015-09-13；

2015-10-19

张 琪(1990-)，男，硕士生。E-mail:noweiciji@126.com

** 通讯作者： E-mail: jishangxu@ouc.edu.cn

TU431

A

1672-5174(2017)04-073-08

10.16441/j.cnki.hdxb.20150318

张琪, 徐继尚. 波致Gibson粉土质海床累积孔压响应的简化分析[J]. 中国海洋大学学报(自然科学版), 2017, 47(4): 73-80.

ZHANG Qi, XU Ji-Shang. Simplified analysis of wave-induced residual pore pressure response in Gibson silty seabed[J]. Periodical of Ocean University of China, 2017, 47(4): 73-80.

Supported by the National Natural Science Foundation of China (51479182，41006024)