基于貝葉斯最大熵和多源數(shù)據(jù)的作物需水量空間預(yù)測

王景雷，康紹忠，孫景生，陳智芳，宋 妮

?

基于貝葉斯最大熵和多源數(shù)據(jù)的作物需水量空間預(yù)測

王景雷1,2，康紹忠1,3※，孫景生2，陳智芳2，宋 妮2

（1. 西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院，楊凌 712100； 2. 中國農(nóng)業(yè)科學(xué)院農(nóng)田灌溉研究所，新鄉(xiāng) 453002； 3. 中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，北京 100083）

作物需水量是灌溉工程規(guī)劃、設(shè)計和管理的重要基礎(chǔ)數(shù)據(jù)，充分利用多源數(shù)據(jù)和先驗知識，快速經(jīng)濟(jì)地獲取精度較高的區(qū)域作物需水量對于區(qū)域水資源的優(yōu)化配置具有重要意義。為精確預(yù)測作物需水量，該文以長系列實際監(jiān)測和校核作物系數(shù)后計算得到的作物需水量為硬數(shù)據(jù)，利用硬數(shù)據(jù)確定獲得最大熵的約束條件，根據(jù)軟數(shù)據(jù)獲取渠道的不同（部分年份缺失的站點數(shù)據(jù)、文獻(xiàn)中獲得的數(shù)據(jù)、利用灌溉試驗數(shù)據(jù)庫中的作物需水量資料，采用協(xié)同克立格方法獲得的數(shù)據(jù)、考慮主要地形因子和主要氣象要素的影響，采用主成分分析和地理加權(quán)回歸（geographically weighted regression，GWR）方法獲得作物需水量數(shù)據(jù)以及遙感數(shù)據(jù)），提出不同來源軟數(shù)據(jù)的概率密度函數(shù)表達(dá)方法，采用貝葉斯最大熵（Bayesian maximum entropy，BME）方法對不同來源的作物需水量信息進(jìn)行有機(jī)整合。結(jié)果表明：除硬數(shù)據(jù)+文獻(xiàn)軟數(shù)據(jù)外，其他數(shù)據(jù)整合呈現(xiàn)一致結(jié)果。華北地區(qū)冬小麥作物需水量在豫南地區(qū)較小，中部地區(qū)黃河北岸有連片的相對高值區(qū)，山東需水量相對較高，冀東北的樂亭、唐山附近有相對低值區(qū)。除硬數(shù)據(jù)+文獻(xiàn)軟數(shù)據(jù)比不整合的精度低9.41%外，其他軟數(shù)據(jù)源均可不同程度地提高整合效果，硬數(shù)據(jù)+克立格軟數(shù)據(jù)、硬數(shù)據(jù)+GWR軟數(shù)據(jù)和硬數(shù)據(jù)+除文獻(xiàn)數(shù)據(jù)外的其他軟數(shù)據(jù)分別比不整合的精度提高85.33%、85.75%和91.69%。對考慮地形、氣象等要素的多源數(shù)據(jù)進(jìn)行整合可更好地反映冬小麥作物需水量空間分布的細(xì)節(jié)，顯著提高估算精度，為稀疏監(jiān)測站點地區(qū)水土資源的精準(zhǔn)管理和優(yōu)化配置提供數(shù)據(jù)支撐。

數(shù)據(jù)處理；回歸；整合；作物需水量；貝葉斯理論；硬數(shù)據(jù)；軟數(shù)據(jù)；先驗信息

0 引 言

作物需水量是水土資源平衡計算、灌溉工程規(guī)劃設(shè)計與運行管理中不可或缺的基礎(chǔ)數(shù)據(jù)[1]。國內(nèi)外提出了很多測定和估算作物需水量的方法[2-6]。這些方法在其各自的應(yīng)用范圍內(nèi)雖然都有一定的精度，但也有各自的特點：通過布設(shè)監(jiān)測站點測定的作物需水量精度較高，但受人力、物力、技術(shù)手段及外部大環(huán)境影響，數(shù)據(jù)存在缺失現(xiàn)象[7]，如20世紀(jì)80年代建立的灌溉試驗數(shù)據(jù)庫，雖具有一定的系統(tǒng)性，但不少站點在個別時間仍有缺失現(xiàn)象。利用常規(guī)氣象觀測資料采用作物系數(shù)法計算的作物需水量，盡管也具有一定精度[8]，但同樣存在數(shù)據(jù)缺失問題。利用遙感獲取的作物需水信息多是不同空間粒度的瞬時面數(shù)據(jù)，要想獲得多年平均的整個生育期需水量必須進(jìn)行時空尺度轉(zhuǎn)換，在時空尺度轉(zhuǎn)換過程中常常產(chǎn)生一定程度的不確定性[9]。另外，研究論文和課題報告也是作物需水信息來源之一，但其通常是某一具體年份的歷史數(shù)據(jù)，缺乏系統(tǒng)性。綜上所述，作物需水量數(shù)據(jù)來源多，存在時空錯位，系列長度短等問題，導(dǎo)致不少數(shù)據(jù)難以滿足灌溉工程規(guī)劃、設(shè)計及精細(xì)管理的實際需求。

Christakos等[10-13]提出的貝葉斯最大熵（Bayesian maximum entropy，BME）方法由于具有堅實的理論基礎(chǔ)，既可以利用實際監(jiān)測的硬數(shù)據(jù)，又可以充分利用研究對象的內(nèi)在結(jié)構(gòu)和機(jī)理（先驗知識），在土壤科學(xué)、環(huán)境等領(lǐng)域得到了快速應(yīng)用：Douaik等[14]將實驗室分析的土壤電導(dǎo)率與由野外探針獲得的土壤表面電導(dǎo)率，采用線性模型獲取的具有一定不確定性的數(shù)據(jù)進(jìn)行融合，結(jié)果顯示，采用 BME預(yù)測精度比只使用室內(nèi)監(jiān)測數(shù)據(jù)的普通克里格方法精度高；徐英等[15]利用農(nóng)田土壤水分的采樣數(shù)據(jù)與貝葉斯神經(jīng)網(wǎng)絡(luò)法生成的數(shù)據(jù)，采用BME方法對農(nóng)田土壤水分和養(yǎng)分進(jìn)行空間預(yù)測，取得了較好的估算效果；Lee等[16]利用歷史監(jiān)測的最小溫度數(shù)據(jù)，將完全沒有缺失的數(shù)據(jù)和部分缺失的數(shù)據(jù)進(jìn)行融合，評估鳳凰城市的熱島效應(yīng)，結(jié)果表明，與未充分利用缺失數(shù)據(jù)的傳統(tǒng)線性克里格方法相比，精度提高了35.28%。Yu等[17]將實際測定的臺灣省臺北市城區(qū)的PM10和PM2.5數(shù)據(jù)作為硬數(shù)據(jù)，根據(jù)不同土地利用類型面積與PM10和PM2.5含量之間的關(guān)系構(gòu)建回歸模型，計算獲得的PM10和PM2.5作為軟數(shù)據(jù)，對臺北市的空氣質(zhì)量進(jìn)行了時空分析；Zhang等[18]利用華北地區(qū)2003－2012年40個站點監(jiān)測得到的年平均降雨量數(shù)據(jù)和8個考慮距海岸線距離、高程等因素得到的插補(bǔ)數(shù)據(jù)進(jìn)行融合，得到了年降雨量的空間分布。

灌溉工程的規(guī)劃設(shè)計要求數(shù)據(jù)系列不宜少于30 a[19]，但目前嚴(yán)格滿足該要求的站點卻很有限，與此同時，國內(nèi)相關(guān)部門進(jìn)行了大量研究試驗，采用不同的方法或技術(shù)手段獲取了大量數(shù)據(jù)，如何充分利用這些寶貴的數(shù)據(jù)一直是讓人頭疼的難題。由于BME既可以利用嚴(yán)格滿足規(guī)范設(shè)計要求，精度較高的硬數(shù)據(jù)，又可以根據(jù)監(jiān)測屬性所蘊含的先驗知識對不同來源的軟數(shù)據(jù)進(jìn)行硬化處理，因此BME理論的提出及在其他領(lǐng)域的應(yīng)用為不同來源作物需水量的有機(jī)整合提供了可行的技術(shù)手段和可資借鑒的數(shù)據(jù)整合框架。本文研究目的是在充分借鑒與作物需水量密切相關(guān)的土壤、氣象等領(lǐng)域軟數(shù)據(jù)獲取及軟硬數(shù)據(jù)整合經(jīng)驗的基礎(chǔ)上，充分考慮數(shù)據(jù)獲取來源及影響作物需水量的主要因素，將不同來源的數(shù)據(jù)采用BME方法進(jìn)行有機(jī)整合，以期提高大區(qū)域作物需水量的估算精度。

1 貝葉斯最大熵方法概述

依據(jù)BME理論，目前能夠獲取的作物需水信息通常分為2類[10,20-21]：1）廣義知識，主要包括作物需水信息及其影響因子的變化規(guī)律和統(tǒng)計數(shù)據(jù)等；2）特定知識，是針對某一特定時間、特定空間而言，例如為了研究某一尺度作物需水信息的空間分布，在研究空間內(nèi)獲取的監(jiān)測點的作物需水量數(shù)據(jù)就是針對這個點的特定知識。特定知識既包含硬數(shù)據(jù)也包括軟數(shù)據(jù)[10]。硬數(shù)據(jù)理論上要求是真實的沒有誤差的數(shù)據(jù)，但由于目前監(jiān)測技術(shù)的限制，監(jiān)測誤差較小可以忽略不計的數(shù)據(jù)也可作為硬數(shù)據(jù)使用[21]。而軟數(shù)據(jù)則指包含了一定不確定性的數(shù)據(jù)，如根據(jù)專家經(jīng)驗估算的數(shù)據(jù)，或是根據(jù)相關(guān)理論派生的數(shù)據(jù)、使用精度不高的儀器設(shè)備監(jiān)測的數(shù)據(jù)等。廣義知識與特定知識的集合=∪，構(gòu)成了作物需水的全部信息，本文的研究目的就是利用BME對這些不同來源的作物需水量信息進(jìn)行有效整合。BME分析包括3個階段：先驗階段（知識獲?。?、中間階段（知識表達(dá)及整合）和后驗階段（信息處理）[22-24]。

1.1 先驗階段

該階段是基于最大熵原理從豐富的廣義知識中找出包含信息量最大的先驗分布，同時在先驗階段將定性的知識轉(zhuǎn)化為定量的數(shù)學(xué)約束。其主要任務(wù)是利用廣義知識求得最大信息熵，保證將最豐富的信息融入到以后的估計過程中[20-21]。基于Shannon定義的信息熵概念[25]

式中=[hard,soft,χ]表示區(qū)域上每個點的作物需水量數(shù)據(jù)集，hard表示硬數(shù)據(jù)點的需水量數(shù)據(jù)集，soft表示所有軟數(shù)據(jù)點的需水量數(shù)據(jù)集，χ表示待估點的作物需水量數(shù)據(jù)集，f表示作物需水量x的先驗概率密度函數(shù)。其數(shù)學(xué)期望可表示為

（2）

要從先驗知識中獲取最大的信息量，就必須熵值的數(shù)學(xué)期望最大，也即在一定的約束條件下，使式（2）需達(dá)到最大值。因此這一階段的主要目的是在這些約束條件下，確定式（2）取得最大值時的先驗概率密度函數(shù)

式中g取決于先驗知識；c為約束條件的總個數(shù)，與χ總個數(shù)有關(guān)；為配分函數(shù)，其作用是正則化約束，可由式（4）計算

（4）

μ是第1個拉格朗日乘子，可由式（5）求得

1.2 中間階段

該階段的主要任務(wù)是以適當(dāng)?shù)男问奖硎驹谔囟〞r間、特定空間獲取的作物需水信息，這些信息包括硬數(shù)據(jù)和軟數(shù)據(jù)，兩者構(gòu)成了特定知識。硬數(shù)據(jù)通常用hard=[1,…,]T表示，軟數(shù)據(jù)用soft=[χ+1,…,]T表示，特定知識的數(shù)據(jù)用data=[hard,soft]T表示[20]。

1.2.1 硬數(shù)據(jù)

硬數(shù)據(jù)嚴(yán)格說來是沒有任何誤差的實際測定數(shù)據(jù)，但對于作物需水量來說很難做到。因此，此處的硬數(shù)據(jù)是指相對準(zhǔn)確，且不確定性很小的數(shù)據(jù)。華北地區(qū)是中國重要的冬小麥生產(chǎn)基地，但受季風(fēng)氣候影響，小麥生育期降雨少，需要補(bǔ)充灌溉，為了指導(dǎo)冬小麥的科學(xué)灌溉，各省區(qū)布設(shè)了數(shù)量不等的試驗站進(jìn)行需水量和需水規(guī)律試驗，建成了灌溉需水量數(shù)據(jù)庫。由于全生育期的作物需水量多用于灌排工程的設(shè)計及管理，根據(jù)國家標(biāo)準(zhǔn)《灌溉與排水工程設(shè)計規(guī)范》要求[19]，在進(jìn)行灌溉設(shè)計保證率計算時，計算系列不宜少于30 a，但中國不少灌溉試驗站點試驗系列并沒有達(dá)到規(guī)范要求，由于本文的目的是進(jìn)行不同來源數(shù)據(jù)的整合，并且FAO-56推薦的利用常規(guī)氣象要素計算作物需水量的作物系數(shù)法在中國華北地區(qū)也具有一定的精度，因此本文的硬數(shù)據(jù)也包括利用常規(guī)氣象要素（最長系列：1954－2013年）計算得到的計算系列超過30 a的冬小麥需水量數(shù)據(jù)。華北冬小麥適宜種植區(qū)內(nèi)共有52個站點滿足硬數(shù)據(jù)的條件，具體分布如圖1所示，由于這里的硬數(shù)據(jù)是實際測定或計算精度高的數(shù)據(jù)，因此52個點上作物需水量取這些硬數(shù)據(jù)值的概率為1，也即Prob[hard=hard]=1。

1.2.2軟數(shù)據(jù)

1）文獻(xiàn)軟數(shù)據(jù)

國內(nèi)不少專家進(jìn)行了黃淮海冬麥區(qū)的作物需水量試驗，以關(guān)鍵詞“華北地區(qū)”和“冬小麥需水量”從中國知網(wǎng)（http://www.cnki.net）數(shù)據(jù)庫中全文檢索，查到相關(guān)文獻(xiàn)211條，剔除不相關(guān)的文獻(xiàn)后，剩余28篇論文，涉及到河南19個站點、河北5個站點、山東9個站點，山西1個站點，北京2個站點。這些數(shù)據(jù)盡管有些是試驗得到的實測數(shù)據(jù)，但由于這些數(shù)據(jù)與試驗數(shù)據(jù)庫、氣象數(shù)據(jù)估算得到的數(shù)據(jù)時間或空間上不一致，并且這些數(shù)據(jù)系列均小于30 a，因此對這些數(shù)據(jù)進(jìn)行處理時，采用下列原則：如果同一站點的數(shù)據(jù)較多，則按區(qū)間數(shù)據(jù)處理，其最大最小值按所有試驗或估算年份中的最大值、最小值，設(shè)第個站點的最大值為u，最小值為l，則Prob(lu)=1。若站點與硬數(shù)據(jù)點重復(fù)，在數(shù)據(jù)整合時予以刪除。

2）部分年份數(shù)據(jù)缺失的軟數(shù)據(jù)

在華北冬小麥適宜種植區(qū)內(nèi)，有5個站點（定陶、墾利、平度、費縣、陵縣）的試驗數(shù)據(jù)缺失較多，計算數(shù)據(jù)系列沒有達(dá)到30 a，這些數(shù)據(jù)盡管蘊含了較多作物需水量信息，但卻存在一定的不確定性，因此將這些數(shù)據(jù)作為包含不確定信息的軟數(shù)據(jù)（站點分布詳見圖1），這些數(shù)據(jù)簡稱部分年份缺失軟數(shù)據(jù)。

為了構(gòu)建缺失數(shù)據(jù)的概率密度函數(shù)f，根據(jù)統(tǒng)計學(xué)原理，設(shè)第年的作物需水量為x，則從總體樣本（年份）中選出任意（<）個樣本，則個樣本的數(shù)學(xué)期望可表示為

var()=(2)(1 ?/) （7）

式中2為總體方差，在實際應(yīng)用中，樣本數(shù)多為有限樣本，通常用它的無偏估計量s來代替，由于相關(guān)規(guī)范規(guī)定計算系列不少于30 a，因此取30。

由于=(x?)/var()滿足分布[13]，則軟數(shù)據(jù)的概率密度函數(shù)f可表示為：

（9）

式中f為分布的函數(shù)。這里以定陶站為例說明方差的修正過程，定陶站共有19 a的數(shù)據(jù)，根據(jù)式（6）可計算出這19 a的數(shù)學(xué)期望（平均值）為420.02 mm，根據(jù)式（8）可得到有限樣本的2=638.35，利用式（9）可得到定陶站的概率密度函數(shù)為f(soft)=0.285f(0.285x?119.67)。

3）基于灌溉試驗數(shù)據(jù)庫的克立格軟數(shù)據(jù)

由于作物需水量近似認(rèn)為區(qū)域化變量，大多服從正態(tài)分布[7]，其概率密度函數(shù)表示為

以“灌溉試驗數(shù)據(jù)庫”中67個站點的（1961－1990年）作物需水量及其生育期內(nèi)主要氣象要素為數(shù)據(jù)源，利用ArcGIS地統(tǒng)計分析模塊的協(xié)同克立格工具，以日照時數(shù)、風(fēng)速和高程作為協(xié)同因子對作物需水量進(jìn)行空間化處理，可得到整個華北冬小麥種植區(qū)的冬小麥需水量和標(biāo)準(zhǔn)差的空間分布柵格圖。為了避免選擇的點過于集聚，利用ArcGIS數(shù)據(jù)管理工具中的隨機(jī)取樣命令在整個華北冬麥區(qū)設(shè)置如下取樣條件：選取數(shù)目500個，相鄰站點距離不小于20 km，獲得站點數(shù)目274個，然后利用ArcGIS的空間分析模塊的采樣工具獲得這274個站點的平均值（數(shù)學(xué)期望）和標(biāo)準(zhǔn)差，得到平均值和標(biāo)準(zhǔn)差后，根據(jù)式（10）得到軟數(shù)據(jù)的概率密度函數(shù)。這些軟數(shù)據(jù)簡稱克立格軟數(shù)據(jù)。

4）基于主成分分析和地理加權(quán)回歸的作物需水量軟數(shù)據(jù)

空間軟數(shù)據(jù)雖然包含有不確定性的數(shù)據(jù)，但是，它一般與估計的目標(biāo)變量具有不同程度的相關(guān)性，蘊藏著豐富的信息。例如作物需水量與氣象因子密切相關(guān)，則作物需水量的空間分布很有可能受到氣象要素主導(dǎo)，根據(jù)影響作物需水的主導(dǎo)因子與作物需水量的關(guān)系獲得的數(shù)據(jù)就可以作為軟數(shù)據(jù)。由于回歸模型可以考慮地形和氣象因素對作物需水量空間分布的影響，因此本文考慮研究區(qū)內(nèi)67個站點（1961－2004年）的作物需水量為因變量，站點的宏觀地形要素（經(jīng)度、緯度和高程）、微觀地形要素（坡度、坡向和遮蔽度）和對應(yīng)生育期的主要氣象要素（最高氣溫、最低氣溫、平均氣溫、風(fēng)速、相對濕度、降雨、蒸發(fā)、日照時數(shù)、日較差）為自變量構(gòu)建回歸模型獲取軟數(shù)據(jù)。由于這些影響因子之間存在相關(guān)性，而相關(guān)性的存在違背了傳統(tǒng)的全域回歸分析中影響因子相互獨立的基本假定，因此首先采用主成分分析（principal component analysis，PCA）方法消除這些影響因素之間的共線性，選取其中的5個主成分，然后采用地理加權(quán)回歸方法（geographically weighted regression，GWR）得到不同空間的局域回歸方程，以與協(xié)同克立格同樣的限制條件獲得的274個樣點作為驗證點，獲得274個樣點的平均值和標(biāo)準(zhǔn)差，具體過程和步驟詳見文獻(xiàn)[27]。在獲得樣點的平均值和標(biāo)準(zhǔn)差后，根據(jù)式（10）獲得這些樣點的概率密布函數(shù)，由于GWR方法采用的是局域回歸分析，與全局回歸分析相比，可以更好地反映空間變異性[27]，利用此方法獲得的軟數(shù)據(jù)簡稱GWR軟數(shù)據(jù)。

5）遙感軟數(shù)據(jù)

隨著遙感、遙測和信息技術(shù)的快速發(fā)展，大量傳感器應(yīng)用于作物需水量的監(jiān)測，軟數(shù)據(jù)的獲取變得更為經(jīng)濟(jì)、高效。吳炳方等[28]利用MODIS遙感數(shù)據(jù)，采用ET-Watch模型對海河流域2002－2005年的蒸散進(jìn)行了估算。茌偉偉等[29]采用ET-Watch模型和MODIS數(shù)據(jù)對大興地區(qū)2002－2008年冬小麥蒸散進(jìn)行了遙感反演并對SWAT模型得到的結(jié)果進(jìn)行了驗證。陳鶴等[30]采用2006－2008年的MODIS數(shù)據(jù)，利用SEBS模型對位山灌區(qū)2 a的冬小麥蒸散進(jìn)行了估算，由于難以獲取冬小麥整個生育期的數(shù)字化空間分布圖，因此只選取有具體估算值及標(biāo)準(zhǔn)誤差的站點作為遙感軟數(shù)據(jù)。雖然遙感估算的需水量也具有一定的精度，但往往僅僅是1 a或幾年的數(shù)據(jù)，很難構(gòu)成系列，所以本文按區(qū)間軟數(shù)據(jù)進(jìn)行處理。例如河北館陶縣董固監(jiān)測的冬小麥需水量分別為338.1、320、252.5和320 mm，則認(rèn)為館陶董固附近冬小麥的耗水量位于區(qū)間[252.5,338.1]的概率為1，也即Prob(252.5

1.3 后驗階段

該階段的任務(wù)是綜合中間階段獲取的特定知識，根據(jù)所有知識來修正根據(jù)廣義知識得到的先驗概率密度函數(shù)，從而得到后驗條件概率密度函數(shù)，最終實現(xiàn)廣義知識和特定知識的的整合，一般通過貝葉斯條件化來實現(xiàn)

后驗概率密度函數(shù)通常用下式表示

1.4 整合結(jié)果及不確定性評價

整合結(jié)果用后驗概率的數(shù)學(xué)期望來表示，其數(shù)學(xué)期望為

當(dāng)軟數(shù)據(jù)為概率性分布形式時，其數(shù)學(xué)期望可用式（14）[14,20]表示

（14）

整合結(jié)果的不確定性通常用方差來衡量

由于研究區(qū)內(nèi)滿足長系列的站點數(shù)目較少，不宜采用一組數(shù)據(jù)參與建模、另一組數(shù)據(jù)進(jìn)行驗證的方式，因此這里采用交叉驗證（Leave-One-Out）方法[22,31-32]進(jìn)行整合結(jié)果的精度評價。該方法的具體步驟是每次建模時留下1個樣本作為驗證集，其他-1個樣本作為訓(xùn)練集，得到驗證點的預(yù)測誤差，經(jīng)過次重復(fù)，即可得到以下5個評判指標(biāo)的具體值：均方根誤差（root mean square error，RMSE）、平均方差（average variance，AVAR）、標(biāo)準(zhǔn)平均誤差（standardized mean error，SME）、誤差絕對值的平均值（mean of absolute error，ABSMPE）和標(biāo)準(zhǔn)誤差的均方根（root mean square of standardized error，RMSSE）。各個指標(biāo)的具體計算公式詳見文獻(xiàn)[22]。通常來說，RMSE和AVAR指標(biāo)越小，SME和ABSMPE越接近0，RMSSE越接近1，表明整合效果越好，RMSSE小于1表明整合值大于實際值，大于1則表明整合值小于實際值。

另外，為了比較不同數(shù)據(jù)源的整合效果，借鑒Lee[24]于2005年提出的方法，該方法建議采用不同數(shù)據(jù)源得到的方差減少的百分比表示整合精度的變化，計算公式見式（16）。

式中σ為方差，1,…,4分別代表不同數(shù)據(jù)源。

2 結(jié)果與分析

為了說明不同數(shù)據(jù)的整合效果，分別利用硬數(shù)據(jù)、硬數(shù)據(jù)+文獻(xiàn)軟數(shù)據(jù)、硬數(shù)據(jù)+克立格軟數(shù)據(jù)、硬數(shù)據(jù)+GWR軟數(shù)據(jù)和硬數(shù)據(jù)+除文獻(xiàn)數(shù)據(jù)外的其他軟數(shù)據(jù)進(jìn)行了華北地區(qū)冬小麥需水量的空間化處理并加以比較分析。

2.1 華北地區(qū)冬小麥需水量的空間分布

根據(jù)不同的數(shù)據(jù)源，采用BMELIB軟件，得到了作物需水信息整合結(jié)果，不同數(shù)據(jù)源的華北冬麥區(qū)作物需水量空間分布如圖2所示。

從圖2可知，5種數(shù)據(jù)源得到的空間分布結(jié)果，整體上呈現(xiàn)豫南地區(qū)作物需水量小，中部地區(qū)黃河北岸有連片的相對高值區(qū)，山東需水量相對較高，冀東北的樂亭、唐山附近有相對低值區(qū)。圖2b與圖2a相比，整體趨勢相似，僅在豫西、豫西南、豫北的安陽和山東濟(jì)南、惠民等地出現(xiàn)了高值區(qū)，榆社、邢臺、長治等附近低值區(qū)連成一片。出現(xiàn)這種狀況的主要原因是豫西、豫西南的數(shù)據(jù)主要來自文獻(xiàn)[33]，該文獻(xiàn)是通過作物模型計算得到的2003－2012年的冬小麥需水量，豫西的盧氏和欒川的冬小麥需水量分別達(dá)到了515和539 mm，豫西南的西峽也達(dá)到了481 mm，初步懷疑該地區(qū)出現(xiàn)這種現(xiàn)象是2003－2012年冬小麥平均需水量較1954－2013年的平均需水量偏高造成的。為了佐證這個觀點，筆者分析了具有長系列數(shù)據(jù)的豫西的孟津、豫南的信陽和固始這2個時間尺度的年平均冬小麥需水量，發(fā)現(xiàn)這3個站點2003－2012年平均冬小麥需水量分別達(dá)到了500.2、418.4和402.5 mm，而1954－2013年均需水量則分別只有487.4、391.7和390.1 mm，這些地方一方面硬數(shù)據(jù)缺乏，另一方面軟數(shù)據(jù)計算時段的需水量值又較硬數(shù)據(jù)計算時段高，這正是圖2b與圖2a相比，在豫西、豫西南出現(xiàn)大面積高值區(qū)的原因。由于文獻(xiàn)數(shù)據(jù)超過30 a序列的并不多，因此安陽附近很小面積的高值區(qū)（是滑縣年份較短的數(shù)據(jù)造成[34]）以及邢臺、榆社、長治的低值區(qū)[35-36]均可能是由于文獻(xiàn)數(shù)據(jù)與研究系列的不一致造成的。圖2c、圖2d與圖2a相比，豫南地區(qū)整體上變化不大，山東的中部及西北部高值區(qū)連成一片，樂亭、唐山附近的低值區(qū)變小，北京以北、保定西北等研究區(qū)邊緣部分出現(xiàn)了部分高值區(qū)，發(fā)生這些變化的主要原因是研究區(qū)邊緣滿足規(guī)范要求的硬數(shù)據(jù)站點稀疏，同時附近的硬數(shù)據(jù)站點作物需水量相對較低，導(dǎo)致圖2a的研究區(qū)北部出現(xiàn)較大的低值區(qū)，加之圖2c、圖2d所利用的克立格軟數(shù)據(jù)和GWR軟數(shù)據(jù)是根據(jù)“灌溉試驗數(shù)據(jù)庫”提供的不同系列長度的歷史數(shù)據(jù)通過空間化處理獲得的，而通過查閱“灌溉試驗數(shù)據(jù)庫”發(fā)現(xiàn)，天津附近、山西北部的冬小麥需水量均超過了500 mm，而這些數(shù)據(jù)在目前仍具有一定的參考價值，因此導(dǎo)致圖2c、圖2d與圖2a相比在研究區(qū)北部變化較大。同時圖2c、圖2d與圖2a相比，研究區(qū)北部變得相對破碎，特別是北京-保定-石家莊-長治一線破碎程度增加較多，保定-饒陽低值區(qū)不再連成一片，被分割開，這些地方由于軟數(shù)據(jù)的利用，使得數(shù)據(jù)源有所增加，減弱了單純利用有限的硬數(shù)據(jù)進(jìn)行空間化處理造成的平滑效應(yīng)，更能反映作物需水量空間分布的變異性。由于圖2e是將除了文獻(xiàn)的其他所有的軟數(shù)據(jù)參與整合，因此冬小麥需水量空間分布圖細(xì)碎化更嚴(yán)重，表明對冬小麥需水量的空間分布解釋能力更強(qiáng)。

2.2 不同數(shù)據(jù)源整合結(jié)果比較

分別將文獻(xiàn)軟數(shù)據(jù)、克立格軟數(shù)據(jù)、GWR軟數(shù)據(jù)、除文獻(xiàn)數(shù)據(jù)外的其他軟數(shù)據(jù)與硬數(shù)據(jù)進(jìn)行整合，為了比較整合效果，沒有軟數(shù)據(jù)參與整合的相關(guān)指標(biāo)也列于表1。從表1可知，文獻(xiàn)數(shù)據(jù)參與整合后，6個指標(biāo)中，其中SME、ABSMPE、RMSSE指標(biāo)值均大于僅有硬數(shù)據(jù)的處理，其整合效果比沒有軟數(shù)據(jù)參與整合的估算精度（值）低9.41%，除此之外，其他不同來源數(shù)據(jù)的6個指標(biāo)均要比單純利用硬數(shù)據(jù)的空間處理效果要好。在整合精度方面，克立格軟數(shù)據(jù)、GWR軟數(shù)據(jù)、除文獻(xiàn)數(shù)據(jù)外的其他軟數(shù)據(jù)參與整合比整合前的精度分別提高85.33%、85.75%和91.69%均高于85%。利用文獻(xiàn)軟數(shù)據(jù)與硬數(shù)據(jù)整合的標(biāo)準(zhǔn)差指標(biāo)均偏大，其原因可能是文獻(xiàn)數(shù)據(jù)空間點相對比較集中，并且計算或通過試驗獲得的監(jiān)測值與文中所用的硬數(shù)據(jù)監(jiān)測和分析值時間尺度不一致所造成，這也從側(cè)面說明盡管軟數(shù)據(jù)允許具有一定的誤差，但在利用軟數(shù)據(jù)進(jìn)行整合時，必須考慮軟數(shù)據(jù)的獲取及處理方式，在利用區(qū)間數(shù)據(jù)處理較少時間序列的軟數(shù)據(jù)時，必須考慮所選時段在整體時段中的代表性（其最大、最小區(qū)間與整個分析時間序列的區(qū)間盡可能相近）。盡管其他數(shù)據(jù)源5個指標(biāo)相差不大，但相對來說，利用除文獻(xiàn)以外的軟數(shù)據(jù)要比單純利用一種或其他某幾種的整合效果要好些，硬數(shù)據(jù)+GWR軟數(shù)據(jù)要比硬數(shù)據(jù)+克立格軟數(shù)據(jù)和硬數(shù)據(jù)+文獻(xiàn)軟數(shù)據(jù)的整合效果好，這可能是GWR方法可以綜合考慮微地形和主要氣象因素的影響，能夠提供微地形或氣象因子發(fā)生突變地方的冬小麥需水量的細(xì)節(jié)信息有關(guān)[27]。

表1 不同來源數(shù)據(jù)整合效果比較

3 結(jié)論與討論

本文將實際測定和計算精度較高的長系列作物需水量監(jiān)測數(shù)據(jù)作為硬數(shù)據(jù)，查閱相關(guān)文獻(xiàn)得到的站點數(shù)據(jù)、灌溉試驗數(shù)據(jù)庫中，缺測數(shù)值較多的站點數(shù)據(jù)和利用不同的空間化處理結(jié)果獲得的站點數(shù)據(jù)作為軟數(shù)據(jù)，根據(jù)硬數(shù)據(jù)得到的均值或方差作為最大熵的約束條件，求取先驗概率密度函數(shù)的具體表達(dá)方式，探討了各種不同來源軟數(shù)據(jù)及其概率密度的表達(dá)方法，根據(jù)貝葉斯條件概率公式對先驗概率密度函數(shù)進(jìn)行修正，得到后驗概率分配函數(shù)，實現(xiàn)了多源作物需水信息的整合。研究結(jié)果表明，除了文獻(xiàn)軟數(shù)據(jù)外，軟硬數(shù)據(jù)的整合比整合前精度提高均超過85%，具有較明顯的整合效果。盡管不少專家和研究機(jī)構(gòu)對華北地區(qū)冬小麥進(jìn)行了大量研究，積累了相關(guān)文獻(xiàn)，但這些數(shù)據(jù)多是根據(jù)具體任務(wù)針對某一特定條件得到的數(shù)據(jù)，并且缺乏長系列的連續(xù)觀測，對于本文多年平均需水量的整合具有一定影響，整合效果反而不如整合前，這說明在進(jìn)行數(shù)據(jù)整合前，必須加以甄別及預(yù)處理，否則可能適得其反。由于本文分析的是華北地區(qū)多年平均需水量，而遙感數(shù)據(jù)獲取的多是瞬時的作物需水量，盡管一些專家也曾對華北地區(qū)的冬小麥整個生育期的需水量進(jìn)行過相關(guān)研究，但由于研究序列大都較短，并且只有很少站點具有明確的具體值，因此本文沒有單獨將其作為軟數(shù)據(jù)參與數(shù)據(jù)整合效果的比較分析，但這不影響遙感數(shù)據(jù)是快速獲取軟數(shù)據(jù)的性質(zhì)，特別是對于時間尺度短、空間尺度大的區(qū)域，遙感數(shù)據(jù)仍是重要的可以加以整合利用的軟數(shù)據(jù)。

[1] 陳玉民，郭國雙，王廣興，等. 中國主要作物需水量與灌溉[M]. 北京：水利電力出版社，1995.

[2] 孫景生，熊運章，康紹忠. 農(nóng)田蒸發(fā)蒸騰的研究方法與進(jìn)展[J]. 灌溉排水，1993，13(4)：36－38. Sun Jingsheng, Xiong Yunzhang, Kang Shaozhong. Progress and research method in evapotransipiration in cropland[J]. Irrigation and Drainage, 1993, 13(4): 36－38. (in Chinese with English abstract)

[3] 許迪，劉鈺. 測定和估算田間作物騰發(fā)量方法研究綜述[J]. 灌溉排水，1997，16(4)：54－59.Xu Di, Liu Yu. Review of study on calculating and estimating methods on crop evapotanspiration[J]. Irrigation and Drainage, 1997, 16(4): 54－59. (in Chinese with English abstract)

[4] Rana G, Katerji N. Measurement and estimation of actual evapotranspiration in the field under Mediterranean climate:A review[J]. European Journal of Agronomy, 2000, 13(23): 125－153.

[5] 王景雷，孫景生，付明軍，等. 區(qū)域作物需水量估算存在的主要問題及可能的解決途徑[J]. 節(jié)水灌溉，2005(3)：4－7.

Wang Jinglei, Sun Jingsheng, Fu Mingjun, et al. Existing problems and solving methods of regional crop water requirements evaluation[J]. Water Saving Irrigation, 2005(3): 4－7. (in Chinese with English abstract)

[6] Lin Wenjing, Su Zhongbo, Dong Hua, et al. Regional evapotranspiration estimation based on MODIS products and surface energy balance system (SEBS) in Hebei plain Northeastern China[J]. Journal of Remote Sensing, 2008, 12(4): 663－672.

[7] 王景雷，孫景生，劉祖貴，等. 作物需水量觀測站點的優(yōu)化設(shè)計[J]. 水利學(xué)報，2005，36(2)：225－231. Wang Jinglei, Sun Jingsheng, Liu Zugui, et al. Optimal arrangement of observation stations for crop water requirement[J]. Journal of Hydraulic Engineering, 2005, 36(2): 225－231. (in Chinese with English abstract)

[8] 雷志棟，羅毅，楊詩秀，等. 利用常規(guī)氣象資料模擬計算作物系數(shù)的探討[J]. 農(nóng)業(yè)工程學(xué)報，1999，15(3)：119－122. Lei Zhidong, Luo Yi, Yang Shixiu, et al. Calculation of crop coefficient with meteorological data[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 1999, 15(3): 119－122. (in Chinese with English abstract)

[9] 劉元波. 遙感反演區(qū)域蒸散的不確定性研究[C]//2009第五屆蘇皖兩省大氣探測、環(huán)境遙感與電子技術(shù)學(xué)術(shù)研討會專輯. 江蘇省氣象學(xué)會、安徽省氣象學(xué)會，2009：2.

[10] Christakos G. A Bayesian/maximum-entropy view to the spatial estimation problem[J]. Mathematical Geology, 1990, 22(7): 763－777.

[11] Christakos G, Li X. Bayesian maximum entropy analysis and mapping: A farewell to kriging estimators?[J]. Mathematical Geosciences, 1998, 30(4): 435－462.

[12] Serre M L, Christakos G. Modern geostatistics: Computational BME analysis in the light of uncertain physical knowledge–the Equus Beds study[J]. Stochastic Environmental Research and Risk Assessment, 1999, 13(1): 1－26.

[13] Christakos G, Serre M L. BME analysis of spatiotemporal particulate matter distributions in North Carolina[J]. Atmospheric Environment, 2000, 34(20): 3393－3406.

[14] Douaik A, Van Meirvenne M, Tóth T. Soil salinity mapping using spatio-temporal kriging and Bayesian maximum entropy with interval soft data[J].Geoderma, 2005, 128(23): 234－248.

[15] 徐英，夏冰. 綜合BME和BNN法的農(nóng)田土壤水分與養(yǎng)分分布空間插值[J]. 農(nóng)業(yè)工程學(xué)報，2015，31(16)：119－127. Xu Ying, Xia Bing. Spatial interpolation of soil moisture and nutrients using BME combined with BNN[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(16): 119－127. (in Chinese with English abstract)

[16] Lee S J, Balling R, Gober P. Bayesian maximum entropy mapping and the soft data problem in urban climate research[J]. Annals of the Association of American Geographers, 2008, 98(2): 309－322.

[17] Yu H L, Wang C H, Liu M C, et al. Estimation of fine particulate matter in Taipei using land use regression and Bayesian maximum entropy methods[J]. International Journal of Environmental Research and Public Health, 2011, 8(6): 2153－2169.

[18] Zhang Fushen, Yang Zhitao, Zhong Shaobo, et al. Exploring mean annual precipitation values (2003—2012) in a specific area (36°N-43°N, 113°E -120°E) using meteorological, elevational, and the nearest distance to coastline variables[J]. Advances in Meteorology, 2015, 2016(4): 1－13.

[19] 灌溉與排水工程設(shè)計規(guī)范：GB 50288-99[S]. 北京：中國計劃出版社，1999.

[20] Douaik A. Evaluation of the Space-time Variability of Soil Salinity by Statistical, Geostatistical and Bayesian Maximum Entropy Methods[D]. Ghent: Ghent University, 2005.

[21] 張楚天. 貝葉斯最大熵方法時空預(yù)測關(guān)鍵問題研究與應(yīng)用[D]. 武漢：華中農(nóng)業(yè)大學(xué)，2016. Zhang Chutian. Research on Key Issues of Bayesian Maximum Entropy Spatiotemporal Prediction and Its Application[D]. Wuhan: Huazhong Agricultural University, 2016. (in Chinese with English abstract)

[22] Tesar M L L. A Comparison of Spatial Prediction Techniques Using Both Hard and Soft Data[D]. Lincoln: University of Nebraska-Lincoln, 2011.

[23] Serre M L. Environmental Spatiotemporal Mapping and Groundwater Flow Modeling Using the BME and ST Methods[D]. Chapel Hill: University of North Carolina, 1999.

[24] Lee S J. Models of Soft Data in Geostatistics and Their Application in Environmental and Health Mapping[D]. Chapel Hill: University of North Carolina, 2005.

[25] Shannon C E. A mathematical theory of communication[J]. The Bell System Technical Journal, 1948, 27: 379－423.

[26] Manly B F J. Statistics for Environmental Science and Management[M]. Boca Raton: CRC Press, 2008.

[27] Wang J L, Kang S Z, Sun J S, et al. Estimation of crop water requirement based on principal component analysis and geographically weighted regression[J]. Science Bulletin, 2013, 58(27): 3371－3379.

[28] 吳炳方，熊雋，閆娜娜，等. 基于遙感的區(qū)域蒸散量監(jiān)測方法：ETWatch[J]. 水科學(xué)進(jìn)展，2008，19(5)：671－678. Wu Bingfang, Xiong Jun, Yan Nana, et al. ET-Watch for monitoring regional evapotranspiration with remote sensing[J]. Advances in Water Sciences, 2008, 19(5): 671－678. (in Chinese with English Abstract)

[29] 茌偉偉，劉鈺，劉玉龍，等. 基于SWAT模型的區(qū)域蒸散發(fā)模擬及遙感驗證[J]. 中國水利水電科學(xué)研究院學(xué)報，2013，11(3)：167－175. Chi Weiwei, Liu Yu, Liu Yulong, et al. SWAT-based simulation for regional evapotranspiration and verification by means of remote sensing[J]. Journal of China Institute of Water Resources and Hydropower Research, 2013, 11(3): 167－175. (in Chinese with English abstract)

[30] 陳鶴，楊大文，劉鈺，等. 基于遙感模型的華北平原農(nóng)田區(qū)蒸散發(fā)量估算[J]. 灌溉排水學(xué)報，2014，33(4)：21－25. Chen He, Yang Dawen, Liu Yu, et al. Cropland evapotranspiration estimation based on remote sensing model in the north china plain[J]. Journal of Irrigation and Drainage, 2014, 33(4): 21－25. (in Chinese with English abstract)

[31] 于雷，洪永勝，耿雷，等. 基于偏最小二乘回歸的土壤有機(jī)質(zhì)含量高光譜估算[J]. 農(nóng)業(yè)工程學(xué)報，2015，31(14)：103－109. Yu Lei, Hong Yongsheng, Geng Lei, et al. Hyperspectral estimation of soil organic matter content based on partial least squares regression[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(14): 103－109. (in Chinese with English abstract)

[32] 陳鵬飛，楊飛，杜佳. 基于環(huán)境減災(zāi)衛(wèi)星時序歸一化植被指數(shù)的冬小麥產(chǎn)量估測[J]. 農(nóng)業(yè)工程學(xué)報，2013，29(11):124－131.

Chen Pengfei, Yang Fei, Du Jia. Yield forecasting for winter wheat using time series NDVI from HJ satellite[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(11): 124－131. (in Chinese with English abstract)

[33] 陳丹丹，李國強(qiáng)，張杰，等. 基于DSSAT 模型的河南省冬小麥需水量分析[J]. 河南農(nóng)業(yè)科學(xué)，2017，46(2)：153－160. Chen Dandan, Li Guoqiang, Zhang Jie, et al. Analysis on water demand of winter wheat in henan province based on DSSAT model[J]. Journal of henan agricultural sciences, 2017, 46(2): 153－160. (in Chinese with English abstract)

[34] 孫現(xiàn)領(lǐng). 基于冬小麥產(chǎn)量與蒸發(fā)量相關(guān)性研究的安陽節(jié)水農(nóng)業(yè)探討[D]. 重慶：西南大學(xué)，2012. Sun Xianling. Studies Based on Winter Wheat Yield and the Evaporation Correlation Water-saving Agriculture in Anyang[D]. Chongqing: Southwest University, 2012. (in Chinese with English abstract)

[35] 梁薇. 冬小麥經(jīng)濟(jì)節(jié)水灌溉制度的研究[D]. 邯鄲：河北工程大學(xué)，2007. Liang Wei. Study on Economical Water-saving Irrigation Schedule for Winter Wheat[D]. Handan: Hebei University of Engineering, 2007. (in Chinese with English abstract)

[36] 楊計剛. 河北省臨漳縣冬小麥科學(xué)灌溉制度研究[D]. 邯鄲：河北工程大學(xué)，2010.Yang Jigang. The Study of Winter Wheat Science Irrigation in Linzhang City Heibei Province[D]. Handan: Hebei University of Engineering, 2010. (in Chinese with English abstract)

Spatial prediction of crop water requirement based on Bayesian maximum entropy and multi-source data

Wang Jinglei1,2, Kang Shaozhong1,3※, Sun Jingsheng2, Chen Zhifang2, Song Ni2

(1.,,712100,; 2.,,453002,; 3.,,100083,)

Crop water requirement is an important basic data for planning, design and management of irrigation engineering. Obtaining high-precision regional crop water requirement prediction using multi-source data and the priori knowledge has great significance for optimal allocation of regional water resources. In the paper, multisource data was integrated using the Bayesian maximum entropy (BME) method for crop water requirement prediction. A long series of crop water requirement measured and calculated by using the crop coefficients adjusted for actual measurement, were taken as the hard data. The soft data included the missing data in partial years for some stations, literature data, the Kriging interpolation data considering the main influence elements of crop water requirement, the crop water requirement data based on the principal component analysis (PCA) and geographically weighted regression (GWR) method, and the remote sensing data. For the soft data from different sources, the expression method of probability density function was put forward and the crop water requirement information from different sources was well integrated using the BME method. Hard data for the period of 1954-2013were collected from measurements from the irrigation stations in North China. Soft data for winter wheat in North China were also collected by searching literatures and the others. The results showed that spatial distribution of crop water requirement in North China was almost consistent for the hard data, hard data + Kriging soft data, hard data + GWR soft data and hard data + the soft data except for literature data. In the southern Henan had smaller crop water requirement, but the middle part (northern part of the Yellow River) of the North China was relatively high. The crop water requirement was relatively high in Shandong but low in the northeast of Hebei such as Leting, Tangshan. The results from hard data + literature soft data were slightly different from the others and the difference might be because the time periods used were different. In general, the integration accuracy of hard data + literature soft data was 9.41% lower than that based on hard data only. Hard data integrated with the other soft data could improve the integration effect. In particular, the integration accuracy of hard data + Kriging soft data, hard data + GWR soft data and hard data + the soft data except for literature data increased by 85.33%, 85.75% and 91.69%, respectively. The integration of multi-source data through considering the terrain, meteorological factors and etc, can could better reflect the spatial distribution of crop water requirements for winter wheat, and significantly improve the estimation accuracy of crop water requirement for winter wheat. The presented method provided the most important basic data for the precise management and optimization of water and soil resources in the region with sparse monitoring stations. In the paper, we need pay attention to some questions in the soft data processing. The partially missing data of some stations need amend the variance calculation results. In order to avoid the agglomeration phenomenon, the selection of interpolation data need adopt the method of random sampling, and at the same time, the distance between adjacent samplings must be limited, should not be too far or near, and 20 km was advisable. In order to avoid too big error and uncertainty, the literature data must be screened and pretreatment, otherwise, the integration effect may be affected.

data processing; regression; integration; crop water requirement; BME; hard data; soft data; prior information

10.11975/j.issn.1002-6819.2017.09.013

O212.8; S165+.22

A

1002-6819(2017)-09-0099-08

2017-04-09

2017-04-20

水利部公益性行業(yè)科研專項經(jīng)費項目(201501016)；國家自然科學(xué)基金（51609245）；中央級科研院所基本科研業(yè)務(wù)費專項（FIRI2016-09）；河南省基礎(chǔ)與前沿技術(shù)研究（162300410168）

王景雷，男，河南通許人，研究員，主要從事節(jié)水灌溉及作物用水要素的空間變異性研究。新鄉(xiāng) 中國農(nóng)業(yè)科學(xué)院農(nóng)田灌溉研究所，453002。Email：firiwjl@126.com

康紹忠，男，湖南桃源人，中國工程院院士，教授，博士，主要從事農(nóng)業(yè)節(jié)水與水資源可持續(xù)利用研究。北京 中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，100083。Email：kangshaozhong@tom.com

王景雷，康紹忠，孫景生，陳智芳，宋 妮.基于貝葉斯最大熵和多源數(shù)據(jù)的作物需水量空間預(yù)測[J]. 農(nóng)業(yè)工程學(xué)報，2017，33(9)：99－106. doi：10.11975/j.issn.1002-6819.2017.09.013 http://www.tcsae.org

Wang Jinglei, Kang Shaozhong, Sun Jingsheng, Chen Zhifang, Song Ni. Spatial prediction of crop water requirement based on Bayesian maximum entropy and multi-source data[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(9): 99－106. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2017.09.013 http://www.tcsae.org