灰色關(guān)聯(lián)及非線性規(guī)劃法構(gòu)建傳遞函數(shù)估算黑土水力參數(shù)

王子龍，常廣義，姜秋香，付 強(qiáng)，陳偉杰，林百健，印玉明

灰色關(guān)聯(lián)及非線性規(guī)劃法構(gòu)建傳遞函數(shù)估算黑土水力參數(shù)

王子龍，常廣義，姜秋香※，付 強(qiáng)，陳偉杰，林百健，印玉明

（東北農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，哈爾濱 150030）

土壤水分特征曲線和飽和導(dǎo)水率是重要的水力參數(shù)，為了簡便準(zhǔn)確獲取這些參數(shù)，以松嫩平原黑土區(qū)南部為研究區(qū)域，采集136個(gè)采樣點(diǎn)土樣用于測定不同土層土壤水分特征曲線、飽和導(dǎo)水率以及土壤理化性質(zhì)，并運(yùn)用灰色關(guān)聯(lián)分析確定影響土壤水力參數(shù)的主要土壤理化性質(zhì)，采用非線性規(guī)劃構(gòu)建土壤分形維數(shù)、有機(jī)質(zhì)、干容重、土壤顆粒組成與土壤水分特征曲線、飽和導(dǎo)水率之間的土壤傳遞函數(shù)，并通過與現(xiàn)有土壤傳遞函數(shù)對比分析進(jìn)行精度驗(yàn)證。結(jié)果表明：1）土壤分形維數(shù)是估算土壤水分特征曲線模型參數(shù)和飽和導(dǎo)水率的主要參數(shù)之一，同時(shí)，干容重和有機(jī)質(zhì)含量也在不同土層土壤傳遞函數(shù)中起到重要的作用；2）通過驗(yàn)證分析，不同土層各參數(shù)平均絕對誤差接近于0，均方根誤差值也都較小，其中在不同土層土壤傳遞函數(shù)估算的土壤含水率均方根誤差分別為0.022、0.017 cm3/cm3；3）對比分析其他已存的土壤水分特征曲線和飽和導(dǎo)水率的土壤傳遞函數(shù)，該文構(gòu)建的土壤傳遞函數(shù)均方根誤差值均較小，決定系數(shù)值都在0.66以上，表明估算精度較高，均好于其他方法估算精度，具有良好的區(qū)域適應(yīng)性。綜上，所構(gòu)建的土壤水分特征曲線和飽和導(dǎo)水率土壤傳遞函數(shù)可以用于松嫩平原黑土區(qū)土壤水力參數(shù)估算。

土壤；灰色關(guān)聯(lián)；水分；松嫩平原黑土區(qū)；土壤水分特征曲線；飽和導(dǎo)水率；土壤傳遞函數(shù)

0 引 言

土壤水分特征曲線和飽和導(dǎo)水率是構(gòu)建土壤水分與溶質(zhì)運(yùn)移模型的重要參數(shù)[1-2]。兩者都受土壤質(zhì)地、容重、有機(jī)質(zhì)含量、孔隙度等的影響，不同地區(qū)影響因素也存在明顯的差異，另外，土壤具有空間變異性特點(diǎn)，隨著尺度的變化，各因素的影響強(qiáng)度也隨之發(fā)生變化[3-4]。直接測定土壤水分特征曲線和飽和導(dǎo)水率具有耗費(fèi)人力、花費(fèi)時(shí)間、成本高等特點(diǎn)，對于小尺度可行，但對大尺度而言，往往需要采集大量的樣本，因此，可操作性較低[5]。

研究表明，通過較為容易測定的土壤基本理化性質(zhì)來獲得土壤水力參數(shù)的方法是一種較為簡便有效方法，即土壤傳遞函數(shù)。構(gòu)建土壤傳遞函數(shù)的方法有多元線性和非線性回歸、人工神經(jīng)網(wǎng)絡(luò)、數(shù)據(jù)分類分組等方法[6-8]。近年來，國內(nèi)外學(xué)者針對不同地區(qū)建立了土壤水分特征曲線、飽和導(dǎo)水率的土壤傳遞函數(shù)。Li等[9]建立了封丘地區(qū)土壤水分特征曲線的土壤傳遞函數(shù)，并用于玉米、小麥產(chǎn)量估算，與Vereecken等[10-12]建立的土壤水分特征曲線的土壤傳遞函數(shù)對比分析，結(jié)果表明Li等建立的各參數(shù)傳遞函數(shù)估算精度最好。孫美等[13]建立了北京大興區(qū)土壤飽和導(dǎo)水率傳遞函數(shù)。孫麗等[14]建立了科爾沁沙丘-草甸相間地區(qū)表土飽和導(dǎo)水率的土壤傳遞函數(shù)，并與Wosten等[12,15-17]建立的飽和導(dǎo)水率土壤傳遞函數(shù)對比發(fā)現(xiàn)，建立的土壤傳遞函數(shù)估算精度優(yōu)于Campbell等建立的模型。由此可知，土壤傳遞函數(shù)具有區(qū)域局限性的特點(diǎn)，因此，針對研究區(qū)域建立一套適應(yīng)性和精度高的土壤水力性質(zhì)的傳遞函數(shù)是十分有必要的。

本文以松嫩平原黑土區(qū)南部為研究區(qū)域，把分形維數(shù)引入土壤傳遞函數(shù)當(dāng)中，采用灰色關(guān)聯(lián)度分析和非線性規(guī)劃相結(jié)合的方法，建立適合松嫩平原黑土區(qū)的土壤水分特征曲線和飽和導(dǎo)水率的土壤傳遞函數(shù)，為研究區(qū)土壤水分和溶質(zhì)運(yùn)移提供參數(shù)，同時(shí)為區(qū)域土壤水力性質(zhì)的獲取提供技術(shù)支持。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)為松嫩平原黑土區(qū)南部的哈爾濱市區(qū)、五常市、雙城區(qū)、阿城區(qū)、賓縣、呼蘭區(qū)、巴彥縣，如圖1所示，總面積達(dá)2.8′104km2，約占松嫩平原黑土區(qū)面積的20%。研究區(qū)屬于溫帶大陸性氣候，降水主要集中在 6－9月，主要以水稻、玉米、大豆農(nóng)作物為主，土壤肥沃，適合農(nóng)作物生長。

圖1 研究區(qū)地理位置及采樣點(diǎn)分布

1.2 土壤樣品的采集與測定

土樣的采集時(shí)間為2017年7－9月，在農(nóng)作物未收割之前。采樣方案的設(shè)計(jì)是在綜合考慮有機(jī)質(zhì)含量、土壤質(zhì)地、土壤容重、土地利用、坡度等土壤水力參數(shù)重要影響因素的基礎(chǔ)上，利用相關(guān)分析確定關(guān)鍵影響因子，通過各因子貢獻(xiàn)賦值及加權(quán)疊加獲得空間計(jì)算單元的優(yōu)先級指數(shù)和排序，運(yùn)用合理采樣數(shù)方法計(jì)算采樣數(shù)量，結(jié)合優(yōu)先級排序確定采樣點(diǎn)的空間位置[18]。采樣點(diǎn)的布設(shè)如圖1所示，研究區(qū)總共136個(gè)采樣點(diǎn)，每個(gè)采樣點(diǎn)分別利用環(huán)刀在0～20 cm和>20～40 cm取原狀土、自封袋取散土，每個(gè)采樣點(diǎn)上下2層各取3個(gè)重復(fù)。原狀土用于測定土壤水分特征曲線和土壤干容重。土壤水分特征曲線采用日本HITACHI公司的高速離心機(jī)測定在1、3、5、10、30、50、100、300、400、500、700、900 kPa下的體積含水率；飽和導(dǎo)水率（K）采用定水頭法測定；土壤干容重（bulk density，BD）采用環(huán)刀法測定；散土土樣進(jìn)行風(fēng)干、過2 mm篩之后，采用英國馬爾文儀器公司MS2000激光粒度儀進(jìn)行土壤顆粒分析，按國際制進(jìn)行分類，分為砂粒（sand）（粒徑≥0.02～2 mm）、粉粒（silt）（粒徑≥0.002～0.02 mm）和黏粒（clay）（粒徑<0.002 mm）。分形維數(shù)（）采用王國梁等[19]提出的方法計(jì)算。土壤有機(jī)質(zhì)（soil organic matter，SOM）采用農(nóng)業(yè)行業(yè)標(biāo)準(zhǔn)（NY/T1121.6-2006）[20]方法進(jìn)行測定。

1.3 研究方法

1.3.1 土壤水分特征曲線模型

土壤水分特征曲線的經(jīng)驗(yàn)?zāi)Ｐ洼^多，國內(nèi)外學(xué)者研究表明van Genuchten模型對不同質(zhì)地的土壤適應(yīng)程度較高[21-23]，因此本文選擇van Genuchten模型來擬合實(shí)測土壤水分特征曲線[24]，公式如下：

式中為土壤體積含水率，cm3/cm3；土壤飽和含水率，cm3/cm3；為土壤殘余含水率，cm3/cm3；為土壤吸力，cm；=1-1/為模型參數(shù)，其中數(shù)值上等于土壤進(jìn)氣值的倒數(shù)。

1.3.2 土壤傳遞函數(shù)構(gòu)建方法

目前，為了獲得土壤傳遞函數(shù)的具體表達(dá)式，以便于推廣應(yīng)用，運(yùn)用線性回歸、非線性回歸等統(tǒng)計(jì)方法建立的土壤傳遞函數(shù)較多，但是，這些方法構(gòu)建的土壤傳遞函數(shù)普遍存在極值點(diǎn)處估算值與實(shí)測值嚴(yán)重偏差，導(dǎo)致較大計(jì)算誤差的問題，而非線性規(guī)劃逐漸逼近目標(biāo)函數(shù)的極值點(diǎn)，能夠有效消除極值偏差問題。

本文把灰色關(guān)聯(lián)度分析法（grey relation analysis，GRA）與非線性規(guī)劃（nonlinear programming，NP）[25-26]相結(jié)合構(gòu)建研究區(qū)土壤水力參數(shù)傳遞函數(shù)，首先運(yùn)用灰色關(guān)聯(lián)分析土壤水力參數(shù)與土壤理化性質(zhì)之間的關(guān)系，確定對土壤水力參數(shù)主要的影響因子，然后采用非線性規(guī)劃構(gòu)建van Genuchten模型參數(shù)、、和飽和導(dǎo)水率K與土壤分形維數(shù)、有機(jī)質(zhì)含量、干容重、土壤顆粒組成之間的土壤傳遞函數(shù)。同時(shí)，在0～20 cm和>20～40 cm土層各136個(gè)土樣中各隨機(jī)選取109個(gè)土樣參與土壤傳遞函數(shù)構(gòu)建，剩余27個(gè)土樣用于土壤傳遞函數(shù)的精度 驗(yàn)證。

1.3.3 模型精度評價(jià)指標(biāo)

選擇平均絕對誤差（mean absolute error，MAE）、均方根誤差（root of mean square error，RMSE）、決定系數(shù)2[9]，對構(gòu)建的土壤傳遞函數(shù)精度進(jìn)行評價(jià)。

2 結(jié)果與分析

2.1 土壤理化特性的統(tǒng)計(jì)分析

對研究區(qū)土壤理化特性進(jìn)行統(tǒng)計(jì)分析，結(jié)果如表1所示。由表可知，研究區(qū)0～20 cm土壤中黏粒、粉粒、砂粒體積分?jǐn)?shù)變化范圍分別為11.30%～42.95%、11.65%～48.59%、18.13%～76.18%；>20～40cm土壤中黏粒、粉粒、砂粒的體積分?jǐn)?shù)變化范圍分別為5.49%～48.95%、3.70%～43%、18.61%～90.80%，由國際制土壤質(zhì)地分類標(biāo)準(zhǔn)可知研究區(qū)土壤質(zhì)地以壤黏土和壤土為主，其中0～20和>20～40 cm的黏粒、粉粒、砂粒體積含量均屬于中等變異，隨著土壤深度的增加，土壤黏粒、粉粒、砂粒變異系數(shù)均增加。0～20 cm土壤有機(jī)質(zhì)質(zhì)量分?jǐn)?shù)在4.74%～98.70%范圍內(nèi)，變異系數(shù)為41%；而>20～40 cm的則在6.19%～107%之間，變異系數(shù)為55%，表明隨著土壤深度的增加變異強(qiáng)度增加。0～20 和>20～40 cm土壤干容重均值分別為1.43和1.33 g/cm3，變異程度均為中等變異。另外分形維數(shù)在0～20 cm的均值為2.82，在>20～40 cm的均值為2.80，變異系數(shù)均為2%，為弱變異。

由表1可知，0～20 和>20～40 cm土層值變化幅度較?。煌瑫r(shí)，0～20 和>20～40 cm土層參數(shù)值比較接近，根據(jù)國際制土壤質(zhì)地分類標(biāo)準(zhǔn)可知研究區(qū)2個(gè)深度土層土壤質(zhì)地為黏壤土，土壤質(zhì)地相似。由0～20和>20～40 cm土層參數(shù)值可知，其中<0.1情況占95%以上，表明研究區(qū)土壤平均孔徑較小、質(zhì)地較細(xì)，屬于壤黏土和壤土范疇。

2.2 構(gòu)建模型估算van Genuchten模型參數(shù)及飽和導(dǎo)水率

采用灰色關(guān)聯(lián)法[28]將不同深度土層實(shí)測土壤理化指標(biāo)與不同深度土層土壤水分特征曲線模型參數(shù)、、進(jìn)行定量分析，并對關(guān)聯(lián)度進(jìn)行排序，如表2所示。關(guān)聯(lián)度大小表示比較序列與參考序列密切程度，關(guān)聯(lián)度越大則兩者變化趨勢越接近，反之，則越遠(yuǎn)離。從表2可以看出，隨著土壤深度的增加，土壤水分特征曲線模型參數(shù)與各指標(biāo)之間的關(guān)聯(lián)度大小也在發(fā)生變化。在0～20 cm土層參數(shù)與分形維數(shù)關(guān)聯(lián)度最大，與砂粒含量關(guān)聯(lián)度最小；在>20～40 cm土層參數(shù)與砂粒含量關(guān)聯(lián)度最大，與粉粒含量關(guān)聯(lián)度最小。在0～20 cm土層參數(shù)與分形維數(shù)關(guān)聯(lián)度最大，與有機(jī)質(zhì)含量關(guān)聯(lián)度最??；而在>20～40 cm土層參數(shù)與分形維數(shù)關(guān)系最為密切。另外，在0～20 cm土層參數(shù)與砂粒含量關(guān)聯(lián)度最大；在>20～40 cm與砂粒含量關(guān)聯(lián)度最大，與黏粒含量關(guān)聯(lián)度最小。

表1 土壤理化特性統(tǒng)計(jì)分析

注：樣本數(shù)為136，變異系數(shù)CV>100%時(shí)為強(qiáng)變異，10%≤CV≤100%時(shí)為中等變異，CV<10%時(shí)為弱變異[27]。

Note: Number of samples is 136, Coefficient of variable CV>100% is a strong variation, 10%≤CV≤100% is a medium variation; CV<10% is a weak variation[27].

表2 不同土層van Genuchten模型參數(shù)和飽和導(dǎo)水率與土壤理化特性間的灰色關(guān)聯(lián)度和排序

注：K為飽和導(dǎo)水率。

可知研究區(qū)PTFs估算趨勢最好，Li和WostenNote: Kis saturated hydraulic conductivity.

通過灰色關(guān)聯(lián)度分析結(jié)果可知，在0～20 cm土層參數(shù)與分形維數(shù)、干容重、有機(jī)質(zhì)含量關(guān)聯(lián)度最高，參數(shù)與分形維數(shù)、干容重、黏粒含量關(guān)聯(lián)度最高，參數(shù)與砂粒含量、干容重、分形維數(shù)關(guān)聯(lián)度最高；在>20～40 cm土層參數(shù)與砂粒含量、干容重、有機(jī)質(zhì)含量關(guān)聯(lián)度最高，參數(shù)與分形維數(shù)、有機(jī)質(zhì)含量、粉粒含量關(guān)聯(lián)度最高，參數(shù)與砂粒含量、分形維數(shù)、有機(jī)質(zhì)含量關(guān)聯(lián)度最高。采用非線性規(guī)劃構(gòu)建van Genuchten模型的參數(shù)方程。

在0～20 cm土層模型參數(shù)公式如下：

（5）

ln=3.02-4.2×10-3C+0.93BD-3.39+

ln=-0.69+1.65×10-2S+0.06BD+

0.642(R=0.77 ,<0.01) （7）

在>20～40 cm土層模型參數(shù)公式如下：

=0.08-1.98×10-4S-1.35×10-3BD+

1.11×10-4OM(R=0.75,<0.01) （8）

ln=-0.53-4.10×10-3S-4.70×10-4OM+

1.59+0.24lnS+40.22e-OM

-0.562(R=0.96)<0.01 （9）

式中C、S、S為黏粒、砂粒、粉粒體積百分?jǐn)?shù)，%；BD為容重，g/cm3；為分形維數(shù)；OM為有機(jī)質(zhì)含量，g/kg。

圖2為不同土層所構(gòu)建的土壤傳遞函數(shù)的估算值與實(shí)測值之間的關(guān)系。由圖2可知，土壤傳遞函數(shù)的估算值與實(shí)測值基本落在1:1線附近，表明所構(gòu)建的土壤傳遞函數(shù)估算值與實(shí)測值基本吻合，其中參數(shù)ln對實(shí)測值的估算效果較好。由圖2中0～20 cm和>20～40 cm土層模型參數(shù)、ln、ln的MAE值可知 MAE值接近于0，表明在不同土層所構(gòu)建的土壤傳遞函數(shù)的估算值與實(shí)測值平均絕對偏離程度較小，估算精度較高。從2值來看，在0～20和>20～40 cm土層模型參數(shù)、ln、ln的2值都在0.66以上，可知在0～20 和>20～40 cm土層、ln、ln的土壤傳遞函數(shù)估算值與實(shí)測值相近，估算效果較好。從RMSE值來看，在0～20和>20～40 cm土層模型參數(shù)、ln、ln的RMSE值接近于0，表明在不同土層所構(gòu)建的、ln、ln的土壤傳遞函數(shù)估算值與實(shí)測值吻合程度較高。

不同土層飽和導(dǎo)水率與土壤理化指標(biāo)之間的關(guān)聯(lián)度及排序如表2所示。由表2可知，在0～20 和>20～40 cm土層土壤飽和導(dǎo)水率與土壤有機(jī)質(zhì)、干容重、分形維數(shù)關(guān)聯(lián)度較高，表明土壤飽和導(dǎo)水率在0～20 cm和>20～40 cm土層與土壤有機(jī)質(zhì)、干容重、分形維數(shù)相關(guān)程度最高，采用非線性規(guī)劃法構(gòu)建飽和導(dǎo)水率的土壤傳遞函數(shù)見式(11)和式(12)。

注：樣本數(shù)為136. MAE和RMSE單位同縱坐標(biāo)，下同。

0～20 cm土層：

（2=0.63,<0.01） (11)

>20～40 cm土層：

（2=0.63,<0.01） (12)

圖3為根據(jù)式（11）和（12）計(jì)算的飽和導(dǎo)水率估算值與實(shí)測值的關(guān)系。由圖3可知，除個(gè)別點(diǎn)外，飽和導(dǎo)水率的估算值與實(shí)測值基本分布在1:1線的兩側(cè)。0～20 cm土層，MAE、RMSE、2值分別為0.34 cm/d、0.472 cm/d、0.63；>20～40 cm土層，MAE、RMSE、2值分別為0.47 cm/d、0.691 cm/d、0.61，由此可知，不同土層飽和導(dǎo)水率土壤傳遞函數(shù)估算與實(shí)測值相近，表明飽和導(dǎo)水率土壤傳遞函數(shù)估算精度較高。

2.3 土壤水力參數(shù)傳遞函數(shù)估算模型驗(yàn)證

根據(jù)已構(gòu)建的土壤傳遞函數(shù)，將剩余27個(gè)土樣實(shí)測數(shù)據(jù)代入式（5）～式（12），計(jì)算得出各參數(shù)的估算值，進(jìn)而得到土壤水分特征曲線和飽和導(dǎo)水率。圖4為不同土層根據(jù)van Genuchten模型參數(shù)土壤傳遞函數(shù)估算土壤含水率和飽和導(dǎo)水率土壤傳遞函數(shù)估算值與實(shí)測值關(guān)系。由圖4可知，在0～20 和>20～40 cm土層土壤含水率散點(diǎn)基本都分布在1:1線兩側(cè)；飽和導(dǎo)水率傳遞函數(shù)的估算值在1:1線兩側(cè)分散程度較小，表明構(gòu)建的土壤傳遞函數(shù)估算精度較好。表3為van Genuchten模型參數(shù)和飽和導(dǎo)水率土壤傳遞函數(shù)估算值的MAE、RMSE、2值。由表3知，各土層、ln、ln、K的MAE和RMSE值均較小，表明各參數(shù)土壤傳遞函數(shù)估算值與實(shí)測值相接近。從2值來看，在0～20 cm和>20～40 cm土層、ln、ln、K的2，除0～20 cm土層K的2值為0.68，其他均在0.9以上，表明在不同土層各參數(shù)土壤傳遞函數(shù)估算值與實(shí)測值基本吻合。總體上本文構(gòu)建的研究區(qū)土壤水力參數(shù)傳遞函數(shù)估算精度較高。

圖3 不同土層飽和導(dǎo)水率估算值與實(shí)測值的關(guān)系

圖4 驗(yàn)證階段不同土層土壤含水率θ和飽和導(dǎo)水率Ks估算值與實(shí)測值關(guān)系

2.4 本文模型與其他模型對比

采用多元回歸等方法構(gòu)建的土壤水力參數(shù)傳遞函數(shù)，分別為Li等[9]、Vereecken等[10]、Wosten等[12]土壤水分特征曲線土壤傳遞函數(shù)和Wosten等[12]、孫麗等[14]、Cosby等[16]飽和導(dǎo)水率土壤傳遞函數(shù)與本文采用非線性規(guī)劃構(gòu)建的土壤傳遞函數(shù)（研究區(qū)PTFs）進(jìn)行對比分析。圖5為研究區(qū)27個(gè)點(diǎn)土壤含水率和飽和導(dǎo)水率估算值與實(shí)測值關(guān)系圖。

表3 不同土層van Genuchten模型參數(shù)驗(yàn)證

由圖5可知，3種土壤傳遞函數(shù)對土壤含水率和飽和導(dǎo)水率估算值與實(shí)測值偏差較大，估算效果不理想，研究區(qū)PTFs估算值與實(shí)測值基本分布在1:1線附近。由圖可知，在0～20 cm土層基于Li、Wosten、Vereecken和本文模型4種土壤傳遞函數(shù)估算的體積含水率MAE值分別為0.988、0.373、1.554、0.015 cm3/cm3，可知研究區(qū)PTFs估算與實(shí)測值平均偏離較小，而其他方法平均偏離程度較大；相應(yīng)的RMSE值分別為0.076、0.054、0.108、0.022 cm3/cm3，表明Vereecken估算效果最差；從2值來看，研究區(qū)PTFs2為0.92，而Li、Vereecken和Wosten的2均小于0.7，表明構(gòu)建的PTFs估算值與實(shí)測值吻合程度較高。在>20～40 cm土層Li、Wosten、Vereecken和本文模型4種土壤傳遞函數(shù)估算的體積含水率MAE值分別為1.879、2.712、1.284、0.021 cm3/cm3，表明Li、Vereecken和Wosten 3種土壤傳遞函數(shù)對實(shí)測值估算誤差較大；相應(yīng)的RMSE值分別為0.087、0.133、0.066、0.017 cm3/cm3，研究區(qū)PTFs偏差最小，Wosten偏差最大；其他區(qū)域3種土壤傳遞函數(shù)2值均小于0.73，而研究區(qū)2值為0.94，可知研究區(qū)PTFs估算趨勢最好，Li和Wosten估算趨勢較差。由誤差分析可知，本文構(gòu)建的土壤水分特征曲線的土壤傳遞函數(shù)估算效果均為最優(yōu)。

由圖5可知，Wosten和Sun估算的K相對于1:1線偏離程度較大，估算效果較差。表4為研究區(qū)飽和導(dǎo)水率土壤傳遞函數(shù)與Wosten、Sun和Cosby 3種飽和導(dǎo)水率土壤傳遞函數(shù)的誤差分析。在0～20 cm土層由4種土壤傳遞函數(shù)估算的K的MAE值可知，4種土壤傳遞函數(shù)估算效果較好；從RMSE值來看，Wosten、Sun、Cosby 3種傳遞函數(shù)估算研究區(qū)飽和導(dǎo)水率RMSE都大于0.3 m/d，研究區(qū)PTFs估算精度最高，而Wosten估算精度最差；由2值可知，Wosten、Sun和Cosby的2都小于0.45，可知估算值與實(shí)測值變化趨勢偏不一致，而研究區(qū)PTFs估算值與實(shí)測值相吻合。在>20～40 cm土層，從4種土壤傳遞函數(shù)MAE值來看，研究區(qū)PTFs估算值與實(shí)測值最為接近；就RMSE值而言，Wosten、Sun和Cosby估算精度較差，而研究區(qū)PTFs估算精度較高；由2值可知，研究區(qū)PTFs的2值最大，最能反映實(shí)測值變化趨勢。通過對比分析可知，本文構(gòu)建的飽和導(dǎo)水率土壤傳遞函數(shù)估算精度最高、效果最好。

圖5 不同土層土壤含水率和飽和導(dǎo)水率估算值與實(shí)測值對比

表4 4種方法估算的土壤飽和導(dǎo)水率誤差分析

3 討 論

由于受到土壤容重、土壤質(zhì)地、土地利用、有機(jī)質(zhì)含量等多種因素的影響，土壤水分特征曲線和飽和導(dǎo)水率具有較強(qiáng)的空間變異性、尺度效應(yīng)和不確定性，因此，導(dǎo)致了土壤傳遞函數(shù)的區(qū)域性特點(diǎn)。

劉繼龍等[4]研究表明在楊凌地區(qū)影響van Genuchten模型參數(shù)的因素隨著不同土層而發(fā)生變化；朱安寧等[29]研究表明在封丘地區(qū)van Genuchten模型參數(shù)與土壤質(zhì)地、有機(jī)質(zhì)含量、容重等存在相關(guān)性；Vereecken等[10]對比利時(shí)地區(qū)研究結(jié)果表明van Genuchten模型參數(shù)與砂粒、黏粒、容重、有機(jī)碳這幾種影響因素關(guān)系密切，這與本文van Genuchten模型參數(shù)的主要影響因素存在一定差異，這可能是由于研究區(qū)不同而造成的。分形維數(shù)是衡量土壤質(zhì)地粗細(xì)程度的一個(gè)指標(biāo)，本文引入土壤分形維數(shù)這一影響因素，研究表明在0～20 和>20～40 cm土層土壤殘余含水率與干容重、有機(jī)質(zhì)含量、砂粒含量密切相關(guān)，主要是由于隨著土壤容重的變化孔隙度也發(fā)生變化，土壤持水能力也隨之發(fā)生變化；而土壤有機(jī)質(zhì)含量增加會使土壤結(jié)構(gòu)發(fā)生變化，小孔隙增多，比表面積增大，對水的吸附力加大，從而使土壤的持水能力增強(qiáng)；砂粒含量增多會使土壤大孔隙增多，從而使土壤持水能力變?nèi)鮗30]。參數(shù)是與土壤進(jìn)氣值有關(guān)的參數(shù)，在0～20 cm土層參數(shù)與砂粒含量、干容重、分形維數(shù)密切相關(guān)，而在>20～40 cm土層參數(shù)與砂粒含量、有機(jī)質(zhì)、分形維數(shù)關(guān)系密切。參數(shù)是與土壤孔徑分布有關(guān)的參數(shù)，其值與分形維數(shù)呈負(fù)相關(guān)[31]，本文研究表明不同土層參數(shù)與分形維數(shù)關(guān)系最為密切。同時(shí)把分形維數(shù)作為其中1個(gè)自變量與其他重要影響因素共同構(gòu)建土壤傳遞函數(shù)，獲得了較高的預(yù)測精度。另外，把分形維數(shù)作為唯一自變量構(gòu)建van Genuchten模型參數(shù)的土壤傳遞函數(shù)的可行性有待進(jìn)一步研究。

由飽和導(dǎo)水率與土壤理化特性灰色關(guān)聯(lián)度分析可知，影響飽和導(dǎo)水率的主要因素為分形維數(shù)、容重、有機(jī)質(zhì)。這與趙云鵬等[32-34]研究結(jié)果不一致，說明不同地區(qū)飽和導(dǎo)水率影響因素也存在差異。本文情況可能是研究區(qū)黑土有機(jī)質(zhì)含量高，可以促進(jìn)團(tuán)聚體形成，從而改善土壤結(jié)構(gòu)，導(dǎo)致飽和導(dǎo)水率發(fā)生改變[35]。

本文采用非線性規(guī)劃建立研究區(qū)土壤水分特征曲線和飽和導(dǎo)水率土壤傳遞函數(shù)，有效解決了極值點(diǎn)處出現(xiàn)較大偏差的問題，使預(yù)測值估算逼近實(shí)測值，提高了預(yù)測精度。但是，非線性規(guī)劃與線性回歸、非線性回歸等方法相似，都存在不能準(zhǔn)確描述土壤水力參數(shù)與土壤理化性質(zhì)之間的物理機(jī)制，建立的模型不具有物理意義。

由于土壤水分特征曲線采用高速離心機(jī)測定，測定過程中會導(dǎo)致土壤容重發(fā)生變化，從而產(chǎn)生一定的誤差。另外，一些采樣點(diǎn)的預(yù)測精度偏低，可能是由于采樣過程中采集了較為松散或者板結(jié)的土壤，導(dǎo)致土壤水力測參數(shù)的測定出現(xiàn)偏差。

4 結(jié) 論

1）本文以松嫩平原黑土區(qū)南部為研究區(qū)域，采用灰色關(guān)聯(lián)分析法對土壤理化性質(zhì)指標(biāo)與van Genuchten模型參數(shù)和飽和導(dǎo)水率進(jìn)行分析，土壤分形維數(shù)、有機(jī)質(zhì)含量、干容重及土壤顆粒組成都在不同土層各參數(shù)中起到至關(guān)重要作用，其中分形維數(shù)同各參數(shù)的關(guān)聯(lián)度較大。

2）本文把分形維數(shù)作為自變量，并采用非線性規(guī)劃建立研究區(qū)不同土層van Genuchten模型參數(shù)和飽和導(dǎo)水率的土壤傳遞函數(shù)，通過對各參數(shù)土壤傳遞函數(shù)進(jìn)行驗(yàn)證，在0～20 cm和>20～40 cm土層殘余含水率土壤進(jìn)氣值倒數(shù)ln土壤孔徑分布指數(shù)ln飽和導(dǎo)水率K的平均相對誤差（mean relative error，MAE）值接近于0；而lnln、K的均方根誤差（root mean square error，RMSE）值均較??；決定系數(shù)2值基本都在0.66以上。表明所構(gòu)建的傳遞函數(shù)估算精度高，可用于研究區(qū)土壤水力參數(shù)的研究。

3）對比Li、Vereecken、Wosten土壤水分特征曲線土壤傳遞函數(shù)以及Wosten、Sun、Cosby飽和導(dǎo)水率土壤傳遞函數(shù)估算研究區(qū)土壤水力參數(shù)，在0～20和>20～40cm土層Li、Vereecken、Wosten土壤水分特征曲線土壤傳遞函數(shù)對研究區(qū)土壤含水率估算MAE和RMSE值均大于研究區(qū)傳遞函數(shù)誤差值，且2值都小于0.73。在0～20和>20～40 cm土層Wosten、Sun、Cosby飽和導(dǎo)水率土壤傳遞函數(shù)估算研究區(qū)飽和導(dǎo)水率RMSE都大于0.3 m/d，估算較差，而2都小于0.45，都大于研究區(qū)構(gòu)建傳遞函數(shù)誤差，因此，表明本文構(gòu)建的土壤水力參數(shù)傳遞函數(shù)估算效果明顯優(yōu)于其他方法，且區(qū)域適應(yīng)性較強(qiáng)。

[1] Merdun H，Cinar O，Meral R，et al. Comparison of artificial neural network and regression pedotransfer functions for prediction of soil water retention and saturated hydraulic conductivity[J].Soil and Tillage Research, 2006, 90(1): 108－116.

[2] 廖凱華，徐紹輝，吳吉春，等. 不同土壤轉(zhuǎn)換函數(shù)預(yù)測砂土非飽和導(dǎo)水率的對比分析[J]. 水科學(xué)進(jìn)展，2013，24(4)：560－567. Liao Kaihua, Xu Shaohui, Wu Jichun, et al.Estimating unsaturated hydraulic conductivity of sandy soils using different pedotransfer functions[J]. Advances in Water Science, 2013, 24(4): 560－567. (in Chinese with English abstract)

[3] Zeleke T B, Si B C. Characterizing scale-dependent spatial relationships between soil properties using multifractal techniques[J]. Geoderma, 2006, 134(3): 440－452.

[4] 劉繼龍，馬孝義，張振華，等. 基于聯(lián)合多重分形的土壤水分特征曲線土壤傳遞函數(shù)[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2012，34(3)：51－55. Liu Jilong, Ma Xiaoyi, Zhang Zhenhua, et al.Pedotransfer functions of soil water retention curve based on joint multifractal[J]. Transactions of the Chinese Society for Agricultural Machinery, 2012, 34(3): 51－55. (in Chinese with English abstract)

[5] 黃元仿，李韻珠. 土壤水力性質(zhì)的估算—土壤轉(zhuǎn)換函數(shù)[J]. 土壤學(xué)報(bào)，2002，39(4)：517－523. Huang Yuanfang, Li Yunzhu. Estimation of soil hydraulic properties—pedo-transfer functions[J]. Acta Pedologica Sinica, 2002, 39(4): 517－523. (in Chinese with English abstract)

[6] Liu Z X, Shu Q S, Wang Z Y. Apply pedo-transfer functions to simulate spatial heterogeneity of cinnamon soil water retention characteristics in western liaoning province[J]. Water Resources Management, 2007, 21(10): 1751－1762.

[7] Schaap M G, Leij F J, van Genuchten M T. Neural network analysis for hierarchical prediction of soil hydraulic properties[J]. Soil Science Society of America Journal, 1998, 62(4): 847－855.

[8] Minasny B, Mcbratney A. The neuro-m method for fitting neural network parametric pedotransfer functions[J]. Soil Science Society of America Journal, 2002, 66(2): 352－361.

[9] Li Y，Chen D, White R E, et al. Estimating soil hydraulic properties of Fengqiu County soils in the North China Plain using pedo-transfer functions[J]. Geoderma, 2007, 138(3): 261－271.

[10] Vereecken H, Maes J, Feyen J, et al. Estimating the soil moisture retention characte ristic from texture, bulk density, and carbon content[J]. Soil Science, 1989, 148(6): 389－403.

[11] Schaap M G, Leij F J, Genuchten M Th. ROSET TA: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions[J]. Journal of Hydrology, 2001, 251(3): 163－176.

[12] Wosten J H M, Lilly A, Nemes A, et al. Development and use of a database of hydraulic properties of Europ pean soil[J]. Geoderma, 1999, 90(3): 169－185.

[13] 孫美，張曉琳，馮紹元，等. 基于交叉驗(yàn)證的農(nóng)田土壤飽和導(dǎo)水率傳遞函數(shù)研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2014，45(10)：147－152. Sun Mei, Zhang Xiaolin, Feng Shaoyuan, et al. Pedo-transfer function for saturated hydraulic conductivity of agricultural soil based on cross-validation[J]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(10): 147－152. (in Chinese with English abstract)

[14] 孫麗，劉廷璽，段利民，等. 科爾沁沙丘-草甸相間地區(qū)表土飽和導(dǎo)水率的土壤傳遞函數(shù)[J]. 土壤學(xué)報(bào)，2015，52(1)：68－76. Sun Li, Liu Tingxi, Duan Limin, et al. Prediction of saturated hydralic conductivity of surface soil in sand-dune-and- meadow inte rlaced region of Horqin with pedo-transfer functions method.[J]. Acta Pedologica Sinica, 2015, 52(1): 68－76. (in Chinese with English abstract)

[15] Campbell G S. Soil physics with basic: Transport models for soil-plant system[J]. Journal of Hydrology, 1987, 90(3): 359－360.

[16] Cosby B J, Hornberger G M, Clapp R B, et al. A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils[J].Water Resources Research, 1984, 20(6): 682－690.

[17] Saxton K E, Rawls W J, Romberger J S, et al. Estimating generalized soil-water characteri- stics from texture[J].Soil Science Society of America Journal,1986, 50(4): 1031－1036.

[18] 王子龍，陳偉杰，付強(qiáng)，等. 基于優(yōu)先級指數(shù)的土壤采樣設(shè)計(jì)方法研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2018，49(7)：244－251. Wang Zilong, Chen Weijie, Fu Qiang, et al. Method of soil sampling design based on priority index[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(7): 244－251. (in Chinese with English abstract)

[19] 王國梁，周生路，趙其國. 土壤顆粒的體積分形維數(shù)及其在土地利用中的應(yīng)用[J]. 土壤學(xué)報(bào)，2005，42(4)：545－550. Wang Guoliang, Zhou Shenglu, Zhao Qiguo. Volume fractal dimension of particles and it’s applications to land use[J]. Acta Pedologica Sinica, 2005, 42(4): 545－550. (in Chinese with English abstract)

[20] 中華人民共和國農(nóng)業(yè)部. 土壤檢測第6部分：土壤有機(jī)質(zhì)測定：NY/T 1121.6-2006[S]. 北京：中國農(nóng)業(yè)出版社，2006.

[21] Wim M Cornelis, Jan Ronsyn, Marc Van Meirvenne. Evaluation of pedotransfer functions for predicting the soil moisture retention curve[J]. Soil Science, 2001(65): 638－648.

[22] 陳俊英，柴紅陽，Leionid G，等. 再生水水質(zhì)對斥水和親水土壤水分特征曲線的影響[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2018，34(11)：121－127. Chen Junying, Cai Hongyang, Leionid G, et al.Impact of treated waste water quality on repellent and wettable soil water characteristic curve[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(11): 121－127. (in Chinese with English abstract)

[23] 姚姣轉(zhuǎn)，劉廷璽，王天帥，等. 科爾沁沙地土壤水分特征曲線傳遞函數(shù)的構(gòu)建與評估[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2014，30(20)：98－108. Yao Jiaozhuan, Liu Tingxi, Wang Tianshuai, et al.Development and evaluation of pedo-transfer functions of soil water characteristic curves in Horqin sandy land[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(20): 98－108. (in Chinese with English abstract)

[24] van Genuchten M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]Soil Science Society of America Journal, 1980, 44(44): 892－898.

[25] 郭向紅，孫西歡，馬娟娟. 基于混合遺傳算法估計(jì) van Genuchten 方程參數(shù)[J]. 水科學(xué)進(jìn)展，2009，20(5)： 677－682. Guo Xianghong, Sun Xihuan, Ma Juanjuan. Parametric estimation of the van Genuchten's equation based on hybrid genetic algorithm[J]. Advances in Water Science, 2009, 20(5): 677－682. (in Chinese with English abstract)

[26] 付強(qiáng)，蔣睿奇，王子龍，等. 基于改進(jìn)螢火蟲算法的土壤水分特征曲線參數(shù)優(yōu)化[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2015，31(11)：117－122. Fu Qiang, Jiang Ruiqi, Wang Zilong, et al. Optimization of soil water characteristic curves parameters by modified firefly algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(11): 117－122. (in Chinese with English abstract)

[27] 付強(qiáng)，顏培儒，李天霄，等. 凍融期不同覆蓋和氣象因子對土壤導(dǎo)熱率和熱通量的影響[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2017，33(20)：98－105. Fu Qiang, Yan Peiru, Li Tianxiao, et al. Influence of different coverage and meteorological factors on soil thermal conductivity and heat flux during freezing and thawing period[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(20): 98－105. (in Chinese with English abstract)

[28] Deng J L.Introduction to Grey-system theory[J]. The Journal of Grey System, 1989, 1(1): 1－24.

[29] 朱安寧，張佳寶，陳效民，等. 封丘地區(qū)土壤傳遞函數(shù)的研究[J]. 土壤學(xué)報(bào)，2003，40(1): 53－58. Zhu Anning, Zhang Jiabao, Chen Xiaomin, et al. Study on pedo-transfer function in Fengqiu[J]. Acta Pedologica Sinica, 2003, 40(1): 53－58. (in Chinese with English abstract)

[30] 鄭子成，李衛(wèi)，李廷軒，等. 基于分形理論的土壤水分特征曲線研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2012，43(5)：49－54. Zheng Zicheng, Li Wei, Li Tingxuan, et al. Soil water retention curve based on fractal theory in greenhouse soil[J]. Transactions of the Chinese Society for Agricultural Machinery, 2012, 43(5): 49－54. (in Chinese with English abstract)

[31] 徐紹輝，劉建立. 估計(jì)不同質(zhì)地土壤水分特征曲線的分形方法[J]. 水利學(xué)報(bào)，2003，34(1)：78－82. Xu Shaohui, Liu Jianli. Fractal approach for estimating soil water retention curves of various textures[J]. Journal of Hydraulic Engineering, 2003,34(1): 78－82. (in Chinese with English abstract)

[32] 趙云鵬，白一茹，王幼奇，等. 寧夏引黃灌區(qū)土壤飽和導(dǎo)水率空間分異特征[J]. 北方園藝，2017(8)：166－171. Zhao Yunpeng, Bai Yiru, Wang Youqi, et al. Spatial variability of soil saturated hydraulic conductivity in the yellow river irrigated area in ningxia[J]. Northern Horticulture, 2017(8): 166－171. (in Chinese with English abstract)

[33] 李孝良，陳效民，周煉川，等. 西南喀斯特地區(qū)土壤飽和導(dǎo)水率及其影響因素研究[J]. 灌溉排水學(xué)報(bào)，2008，27(5)：74－76. Li Xiaoliang, Chen Xiaomin, Zhou Lianchuan, et al.Soil saturated hydraulic conductivity and its influential factors in southwest karst region of china[J]. Journal of Irrigation and Drainage, 2008, 27(5): 74－76. (in Chinese with English abstract)

[34] 趙春雷，邵明安，賈小旭. 黃土高原北部坡面尺度土壤飽和導(dǎo)水率分布與模擬[J]. 水科學(xué)進(jìn)展，2014，25(6)： 806－815. Zhao Chunlei, Shao Ming'an, Jia Xiaoxu.Distribution and simulation of saturated soil hydraulic conductivity at a slope of northern loess plateau[J]. Advances in Water Science, 2014, 25(6): 806－815. (in Chinese with English abstract)

[35] 張耀方，趙世偉，王子龍，等. 黃土高原土壤團(tuán)聚體膠結(jié)物質(zhì)的分布及作用綜述[J]. 中國水土保持科學(xué)，2015，13(5)：145－150. Zhang Yaofang, Zhao Shiwei, Wang Zilong, et al.Distribution and function of cementing materials of soil aggregates on the Loess plateau，western China[J]. Science of Soil and Water Conservation, 2015, 13(5): 145－150. (in Chinese with English abstract)

Constructing pedo-transfer functions based on grey relational and nonlinear programming to estimate hydraulic parameters in black soil

Wang Zilong, Chang Guangyi, Jiang Qiuxiang※, Fu Qiang, Chen Weijie, Lin Baijian, Yin Yuming

(150030,)

Soil water characteristic curve and saturated hydraulic conductivity are important parameters for constructing soil moisture and solute transport models, both of which are affected by soil texture, bulk density, organic matter content, porosity, etc. In order to simply and accurately obtain the 2 parameters above, a total of 136 soil sampling points in the southern part of the black soil area in the Songnen Plain were designed, and the soil samples in 2 layers of these sampling points were collected to measure the soil water characteristic curve, saturated hydraulic conductivity and soil physical and chemical properties, the grey correlation analysis was used to determine the main soil physical and chemical properties affecting soil hydraulic parameters.The results showed: 1) in 0-20 cm soil layer, the residual moisture contentof van Genuchten model had the higher correlation with fractal dimension, dry bulk density and organic matter content, and the correlation between parameter pore distribution index and fractal dimension, dry bulk density and clay content was the highest. The parameter of reciprocal soil air influx had the higher correlation with sand content, dry bulk density and fractal dimension. In soil layer from 20 to 40 cm, the parameter of residual moisture contentof van Genuchten model had the higher correlation with sand content, dry bulk density and organic matter content; the parameter pore distribution index and fractal dimensions had the highest correlation; The correlation between fractal dimension, organic matter content and silt content and parameter pore distribution index was the highest, and the correlation between parameter reciprocal soil air influx and sand content, fractal dimension and organic matter content was the highest. 2) The correlation between soil saturated hydraulic conductivity and soil organic matter, dry bulk density and fractal dimension in 0-40 cm soil layer was higher. 3) The pedo-transfer functions of the soil water characteristic curve and the saturated hydraulic conductivity based on soil fractal dimension, organic matter, dry bulk density and soil particle composition were established by using the nonlinear optimal programming method. The accuracy of the pedo-transfer functions established in the study was verified by comparing with the existing common pedo-transfer functions. The results showed that soil fractal dimension was one of the primary parameters in constructing soil pedo-transfer functions to predict soil water characteristic curve and saturated hydraulic conductivity. In addition, soil dry bulk density and organic matter content also played an important role in pedo-transfer functions of different soil layers. Based on the verification analysis, the mean absolute errors of the estimated soil hydraulic parameters in different soil layers were all close to 0, and the root mean square error values of the predicted soil hydraulic parameters in different soil layers were also small. Specially, the root mean square errors of soil moisture content estimated by pedo-transfer functions established in the study were 0.022 and 0.017 cm3/cm3in different soil layers, respectively. On the basis of comparing with the other existing common pedo-transfer functions of soil water characteristic curve and saturated hydraulic conductivity, the pedo-transfer functions of soil hydraulic parameters established in the study had high estimation accuracy, the root mean square error values were small; the determination coefficient values were all above 0.66, which was better than the other existing common pedo-transfer functions. In summary, the pedo-transfer functions of the soil water characteristic curve and saturated water conductivity established in the study can be used to estimate the soil hydraulic parameters in the black soil area of the Songnen Plain, and also provide technical support for acquisition of regional soil hydraulic properties.

soils; gray correlation; water content; black soil area of Songnen Plain; soil water characteristic curve; saturated hydraulic conductivity; pedo-transfer functions

10.11975/j.issn.1002-6819.2019.10.008

S152.7

A

1002-6819(2019)-10-0060-09

2018-10-23

2019-03-10

國家自然科學(xué)基金（51579045）；東北農(nóng)業(yè)大學(xué)“學(xué)術(shù)骨干”基金資助項(xiàng)目（16XG10）

王子龍，教授，博士，主要從事寒區(qū)農(nóng)業(yè)水土資源高效利用研究。Email：wangzilong@neau.edu.cn

姜秋香，副教授，博士，主要從事水土資源高效利用和管理研究。Email：jiangqiuxiang2017@163.com

王子龍，常廣義，姜秋香，付 強(qiáng)，陳偉杰，林百健，印玉明. 灰色關(guān)聯(lián)及非線性規(guī)劃法構(gòu)建傳遞函數(shù)估算黑土水力參數(shù)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(10)：60－68. doi：10.11975/j.issn.1002-6819.2019.10.008 http://www.tcsae.org

Wang Zilong, Chang Guangyi, Jiang Qiuxiang, Fu Qiang, Chen Weijie, Lin Baijian, Yin Yuming. Constructing pedo-transfer functions based on grey relational and nonlinear programming to estimate hydraulic parameters in black soil[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(10): 60－68. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2019.10.008 http://www.tcsae.org