用水總量控制下江蘇省工業廢水排放績效及減排潛力研究

吳鳳平 梁蔓琪

摘要工業生產總值增長對國民經濟發展起主導性作用。伴隨著經濟發展，工業化和城鎮化水平的不斷提高，工業用水總量多、用水效率低下、廢水排放帶來的環境污染嚴重等問題已成為制約工業發展的瓶頸。江蘇省工業生產總值位居全國各省份前列，為了保障江蘇省工業可持續發展，本文在全要素和環境生產技術框架下，利用規模方向距離函數建立工業廢水排放績效測度模型與工業廢水減排潛力測度模型，對江蘇省13個地級市2006—2015年工業廢水排放績效進行測算，將工業灰水足跡作為非期望產出，分別得出用水總量控制下期望產出和非期望產出的改進方向、改進目標值及其變化率，從而測度江蘇省工業廢水排放績效及其減排潛力，并從蘇南、蘇北和蘇中三大區域差異化視角進行綜合分析。然后借助Kernel密度估計方法，分析江蘇省各地級市全要素框架下工業廢水排放績效的空間分布規律和動態演變趨勢。結果表明：江蘇省年均工業廢水排放績效為0.772，工業廢水減排值為263.5億m3，占年均工業灰水排放量的24%。就三大區域而言，年均工業廢水排放績效從高到低依次是蘇南、蘇中和蘇北，其值分別為0.911、0.789和0.624。其中，蘇北地區工業廢水排放績效增長速度最快，在2006—2015年期間由0.529增長至0.79，增長率為49.5%。在工業灰水減排潛力方面，蘇北減排潛力最大，其次是蘇中，蘇南最小。從空間分布和動態演變上看，江蘇省各地區的工業廢水排放績效存在空間分布上的差異，但差異正在逐步減小，并且績效值向高水平聚攏。

關鍵詞工業廢水排放績效；規模方向距離函數；減排潛力；時空變化

中圖分類號X22文獻標識碼A文章編號1002-2104（2018）08-0079-08DOI：10.12062/cpre.20171207

近年來，水資源問題已成為制約世界經濟發展的主要問題之一。我國人均水資源占有量少，水資源時空分布差異顯著。伴隨著經濟的發展，工業化和城鎮化水平的不斷提高，工業用水總量多、用水效率低下、廢水排放帶來的環境污染嚴重已成為制約工業發展的瓶頸。2012年，國務院三號文件《關于實行最嚴格水資源管理制度的意見》，指出當前我國水資源面臨的嚴峻形勢，明確了“用水總量控制、用水效率控制、水功能區限制納污控制”三條紅線的主要目標。本文在用水總量約束下，測算出工業廢水的排放績效與減排潛力，為我國工業的可持續發展提供參考。作為我國的經濟發展大省，江蘇省2015年的地區生產總值突破70 000億元，在全國排名第二，其人均生產總值也位居全國第四。因此，本文選取江蘇省作為實證研究對象具有重要的現實意義。

1文獻綜述

近年來，國內外許多學者圍繞經濟發展的資源消耗和環境影響開展研究。Chung等人[1]在Shepherd距離函數的基礎上提出方向距離函數（DDF），相較于傳統的DEA法，該方法在考慮期望和非期望產出時引入方向矢量，在非期望弱處置原則的基礎上，衡量決策單元在該方向上的效率。

國內外一些學者在二氧化碳排放績效、生態效率等問題上開展了大量的研究。Ramli N A 等[2]在Fre， Grosskorf[3]提出松弛測度（SBM）的基礎上，改進DDF中方向矢量的隨意性問題，考慮期望產出松弛與非期望產出松弛，建立規模方向距離函數（SDDF），以此研究馬來西亞制造業的生態效率。Zaim等[4]和Zofio等[5]采用不同的DEA模型，對經濟合作與發展組織的國家和地區進行二氧化碳排放效率上的分析。Zhou Peng等[6]用非徑向DDF模型來評價生產行業的能源與二氧化碳排放效率。Wang Q W 等[7]基于DDF模型的分析方法，對我國28個省份的工業二氧化碳排放效率進行測度。付麗娜等利用MalmquistDEA模型，對生態效率進行動態對比研究[8]。孫作人等[9]利用非參數距離函數和環境生產技術，對我國工業二氧化碳排放強度進行測算。李健等[10]運用改進后的SDDF模型，對中國省區二氧化碳排放效率及減排潛力進行研究。

綜合來看，現有文獻較少關注廢水排放績效改進對國家經濟建設和生態保護的積極影響，并且其研究對象大多選用省區直轄市或是個別城市，較少涉及到全省范圍內地級市的對比研究。本文基于江蘇省13個地級市2006—2015年的面板數據，在全要素和環境生產技術框架下，利用規模方向距離函數對江蘇省各地級市工業廢水排放績效進行測算，這對提升江蘇省工業廢水排放績效和降低工業廢水排放總量具有一定理論和現實意義，對全國也有借鑒意義。

吳鳳平等：用水總量控制下江蘇省工業廢水排放績效及減排潛力研究中國人口·資源與環境2018年第8期2研究方法

2.1構建生產可能性集

1957年，Farrell[11]提出用生產前沿概念來估計技術效率。1997年，Battese[12]提出隨機前沿生產函數。本文根據Fre等[13]關于環境生產技術的研究，關注工業生產中廢水排放和減排問題，在全要素生產框架下構造生產可能性集如下：設生產系統中有J個決策單元（DMU），第i個x=（x1，x2，…，xk）∈RK+后，產出一種期望產出工業生產總值（y∈RI+）和一種非期望產出工業廢水（b∈RM+）。其生產系統生產可能性集可表示為：

P（x）={（x，y，b）}

∑Jj=1ωjxjk≤xik；k=1，2，…，K；

∑Jj=1ωjyj≥yi；

∑Jj=1ωjbj=bi；

ωj≥0；j=1，2，…，J（1）

其中，下標i表示第i個DMU的投入產出數據，ωj表示相對被評價DMU而重新構造的一個有效DMU組合中第j個DMU的組合比例。生產可能性集P（x）是具有凸性特征的有界閉集，且滿足：①有限的投入只能得到有限的產出，投入要素具有自由可處置性，即當x′≥x時，存在P（x）P（x′）；②在不等式的約束下，投入要素和期望產出具有強可處置性，即若其他條件不變情況，期望產出可多可少；③在等式約束下，非期望產出具有弱可處置性與零結合性特征，即非期望產出總是伴隨著期望產出一起出現。

2.2工業廢水排放績效測度模型

本文采用規模方向距離函數（Scale directional distance functions）來對工業廢水排放績效進行測度，基本思路是考慮在現有水資源投入、勞動力投入和既定經濟產出條件下，假設規模報酬可變，用理論上可達到的最小工業廢水排放量與實際排放量進行比較，進而測度排放績效的高低。通過環境生產技術構建出生產可能性集后，引入方向向量g=（gy，-gb）來構建規模方向距離函數如下：

D0（x，y，b；g）=max{λTβ：（（y，b）+

g·diag（β））∈P（x）}（2）

其中，λ=（λy，λb）T，是期望產出與非期望產出的歸一化權向量，本文設定期望產出與非期望產出同等重要，即λy=λb=0.5。β=（βy，βb）T。其中βy值為擴張因子，其值表征期望產出向生產前沿面擴張的規模。βb值為收縮因子，其值表征非期望產出向生產前沿面收縮的規模，βy，βb∈[0，1]。λTβ值表征了被評價DMU既減少工業廢水排放、又增加經濟產出的能力，值越大表明被評價DMU的工業廢水排放績效越低，也就是說其產生的工業廢水越多，其減排潛力也就越大。本文將1-λTβ值定義為全要素工業廢水排放績效指標，值越大工業廢水排放績效越好。本文參考Fre等[13]的研究，建立線性規劃模型來求解D0（x，y，b；g）。

D0（x，y，b；g）=maxλTβ

s.t.

∑Jj=1ωjxjk≤xik；k=1，2，…，K；

∑Jj=1ωjyj≥yi（1+βy）；

∑Jj=1ωjbj=bi（1-βb）；

λy+λb=1；λy，λb>0；ωj≥0； j=1，2，…，J（3）

2.3工業廢水減排潛力測度模型

根據Ramli等[2]的研究，給出被評價DMU到達生產前沿面的方向向量的計算公式。當λTβ=0，DMU處于生產前沿面上，gy和gb能夠取任意值；當λTβ>0，DMU的期望產出和非期望產出的改進方向向量按公式4計算：

gy=β*yβ*y+β*b；gb=β*bβ*y+β*b（4）

其中，gy和gb滿足0≤gy，gb≤1且gy+gb=1。同時，根據計算結果，可以得到非有效DMU欲達到效率最優時其期望產出與非期望產出的改進目標值分別為：

yt=∑Jj=1ωjyj；bt=∑Jj=1ωjbj（5）

本文通過計算在用水總量控制下，各DMU趨近生產前沿面時期望產出和非期望產出的改進方向、改進目標值與其變化率，來測算工業廢水的減排潛力。

2.4Kernel密度估計

核密度估計（Kernel density estimation）作為一種估計概率密度函數的非參數方法，在地區分布規律研究中廣泛運用。對于一組隨機變量x1，x2，…，xn，其核密度估計的形式為：

fh（x）=1nh∑ni=1K（xi-xh）（6）

其中，n表示隨機變量的個數，h表示帶寬。本文選用Sliverman[14]證明在高斯核函數的最優帶寬h*，h*=1.364 3δn-0.2s，其中δ=0.776 4，核函數K（x）是一個加權函數，分為高斯核、三角核、四次核等不同類型，根據分組數據的密集程度來選擇運用哪一類型核函數。本文選取高斯核函數來研究江蘇省工業廢水排放績效的分布動態演化趨勢，其函數表達式為：

Gaussian：K（x）=12πe-12x2（7）

因為核密度估計屬于非參數估計，沒有確定的函數表達式，所以可以通過觀察分析核密度估計結果得到的圖形，來了解變量分布的位置、形態等，分析其分布演化。

3實證研究

3.1投入產出變量選取

考慮到用水總量控制下工業廢水排放績效需要反映出環境保護和工業經濟增長兩方面指標，本文選取資本投入（K）、勞動力（L）和工業用水總量（W）為投入變量，以地區工業總產值（G）為期望產出變量，以工業灰水足跡（D）為非期望產出變量。具體投入產出變量選取見表1。

3.2數據來源與處理

由于相關變量統計資料有限，本文將區間定為2005—2015年，以江蘇省13個地級市為研究對象。根據地理位置分布，劃為蘇北、蘇中和蘇南三大地區。蘇北地區5個城市包括徐州、連云港、宿遷、淮安和鹽城；蘇中地區3個城市包括揚州、泰州和南通；蘇南地區5個城市包括南京、常州、蘇州、無錫和鎮江。本文數據來自《中國城市統計年鑒》《江蘇省統計年鑒》《江蘇省水資源公報》及各地市的統計年鑒、環境公報、水資源公報等，根據研究對象的特點對相關變量數據進行相應的選擇和處理。投入產出變量的統計性描述見表2。

（2）勞動力。勞動力投入理論上應包含勞動力的投入數量、素質和勞動時間，但由于江蘇省各地市勞動人員人均受教育水平和勞動時間數據不可得，因此本文選取統計年鑒中各地市各年規模以上工業企業從業人員數作為勞動力投入。

（3）工業用水總量。本文選用各地市水資源公報中歷年工業用水總量數據。

（4）工業總產值。工業總產值作為期望產出，選擇《江蘇省統計年鑒》中公布的各地市工業總產值數據，并按照對應年份的工業品出廠價格指數，折算成以2005年為基期的不變價格。

（5）工業灰水足跡（工業廢水排放總量）。工業生產過程中產生的污染是點源污染，工業廢水中通常包含多種形式的污染物，而化學需氧量（chemical oxygen demand，COD）是工業廢水中含量最大的污染物。工業灰水足跡是指為了稀釋工業生產排放的污染物以達到相關水質標準的水資源需求量，本文選用工業灰水足跡作為工業生產中的非期望產出。

本文參考國際水足跡網絡出版的《水足跡評估手冊》[17]和曾昭等[18]關于灰水足跡的計算方法，給出工業灰水足跡計算方法如下：

WFgrey=LCmax-Cnat（8）

其中，WFgrey為工業灰水足跡（m3/a），L為COD排放負荷（kg/a），Cmax為達到環境水質標準情況下污染物的最高濃度（kg/m3）。本文選擇Ⅲ類水COD標準濃度為計算工業灰水足跡的依據，取Cmax=0.002 kg/m3，Cnat為COD在水體中的原始濃度，常假設為0。為方便進行績效水平的測算，本文將灰水足跡的計量單位更改為億m3/年。

3.3實證結果與分析

3.3.1工業廢水排放績效評估

在由規模方向距離函數構建的工業廢水排放績效模型的基礎上，本文變量值代入，得到江蘇省13個地級市2006—2015年工業廢水排放績效表，如表3所示。就各地年平均值來看，無錫和鎮江的績效值都在0.9以上，為最佳績效。然而，宿遷的年均績效值最低，僅為0.393，且與江蘇省其它地市的績效存在較大的差距。就蘇南、蘇中和蘇北三個地區而言，蘇南績效最高，為0.911；其次為蘇中地區，績效為0.789；蘇北最差，績效值僅有0.624，且低于江蘇省平均績效水平0.772。導致江蘇各地區工業廢水排放績效差異的主要原因是各地區間的發展差異，影響因素大致分為：經濟發展水平、工業技術水平、工業產業結構以及江蘇各地市地方環保政策的差異。

為了能夠更加直觀地觀察江蘇省三個地區工業廢水排放績效的變化，本文將蘇南、蘇中、蘇北2006—2015年的年均績效值做成趨勢走向圖，并且加入整個江蘇省的年度平均值作為參考（見圖1）。就江蘇省的趨勢曲線來看，全省工業廢水績效呈上升趨勢，且在2014年后加速上升，這與2014年江蘇省響應國家號召發布的節能減排戰略目標，以及最嚴格水資源管理制度實施過程中“三條紅線”控制體系的全面建立有關。由圖1可知，蘇南、蘇中和蘇北地區的工業廢水排放績效時間趨勢走向存在差異。蘇南地區工業廢水排放績效與其他地區相比而言，一直處于0.9左右，具有較少的提升空間。在經濟發展速度上，蘇中和蘇北地區的經濟發展潛力更大，增長速度也比蘇南地區更快。因此，蘇南地區在工業化規模較大、高用水產業集聚的背景下，工業廢水減排潛力在相對平穩中略有降低。

在三個趨勢曲線中，蘇中地區的變化頻率和幅度最大。在2012年之前，蘇中地區的工業廢水排放績效值均高于全省平均水平，尤其在2007—2010年期間穩步提升。而在2012年后，蘇中地區的績效值跌落在全省平均水平之下，且成為三大地區中績效最低的地區。這一變化的原因在于蘇中地區的工業生產結構單一，高新技術產業占比較小。以揚州市為例，揚州市的主要工業產業為制造業，并且傳統制造業如紡織業、普通機械制造業、塑料制品業等占比較大。泰州市的工業生產結構中，高污染、高耗水的工業占比較大，并且隨著電力行業和鋼鐵行業的發展，阻礙了當地工業廢水減排的進程。

以蘇北地區為例，該地區在研究初期的工業廢水排放績效處于較低的水平，且具有很大的提升空間。因為在2006—2015年期間，整個江蘇省處于工業化、城鎮化加速發展的階段，尤其是蘇北地區吸納了來自蘇南地區和其他發達地區的先進生產技術、資金和人才。在這樣的社會大環境下，蘇北地區的經濟發展水平迅猛增長，尤其在2011年后，伴隨著江蘇省“十二五”規劃的推進，工業廢水排放績效的增長速度大幅提升。在2006—2015年期間，工業廢水排放績效由0.529增長至0.79，增長率為49.5%。高新技術產業和第三產業的引入改變了原先的生產結構，并且重污染、高能耗的工業企業在技術調整后減少了工業廢水的排放。

3.3.2工業廢水排放績效時空分布分析

為了更直觀地展現江蘇省工業廢水排放績效的空間分布，方便在研究過程中對績效的動態演變進行分析，本文將2006年、2009年、2012年和2015年江蘇省工業廢水排放績效地區分布進行對比，詳情見圖2。從圖中可知，江蘇省各地區的工業廢水排放績效分布是不均衡的，總體而言，績效水平大體在逐步提升。2006年，江蘇省各地市的工業廢水排放績效呈現出地區分布不均衡的特點，且由南至北績效值不斷下降，績效水平高的地區主要集中在蘇南和部分蘇中地區。而在此之后，伴隨著蘇北地區和部分蘇中地區績效水平的提升，江蘇省內南北差距在不斷的縮小。2006年，績效值在0.65以下的地級市數量為5個，2009年、2012年都減少為3個。而2015年該類別地級市數量減少為2個，并且該年不存在績效值為0.5以下的地級市。在變化幅度上看，宿遷和連云港三個地市擁有最明顯的變化，年均績效從2006年至2015年提升0.4以上。相反，南京、蘇州和南通三個地市的年均績效值卻在十年內有著較為明顯的降低。其中南通市變化幅度最大，2015年較2006年績效值降低0.315。

3.3.3工業廢水排放績效分布演變趨勢

圖3是描述2006—2015年江蘇省工業廢水排放績效的核密度曲線，采用2006、2009、2012和2015四個年度數據。分布圖中的橫軸表示工業廢水排放的績效水平，縱軸表示核密度[19]。總體而言，2006—2015年間，績效呈單峰分布的總體態勢保持不變。就密度的曲線位置和平移來看，2015年相對于2006年，峰值大幅上升，變化范圍大幅縮小，這說明十年間績效的地區差異大大縮小。在2006—2012年間，分布曲線總體呈現右移的趨勢，說明在這段時間內，江蘇省工業廢水排放績效總體提升。但是，在2015年分布曲線較2012年相比左移，但峰值急劇增加，變化范圍小幅縮小，這說明2015年江蘇省各地市的工業廢水排放績效向著高水平聚攏。

3.3.4工業廢水減排潛力測度

改進方向的基礎上，得到投入產出模型中非期望產出改進的目標值。在表4顯示的是江蘇各地市和三個地區的年均工業灰水足跡、年均工業灰水減排值、減排潛力和減排貢獻率，以及各地市在減排潛力和貢獻率上的排名。就城市而言，蘇州工業生產所產生的年均工業灰水足跡最多，為256.75億m3，并且同時占據最高的減排貢獻率

整個江蘇省工業廢水中污染物排放所需的年均灰水足跡為1 191.58億m3，其中蘇南的年均灰水足跡最多，為623.4億m3。這與蘇南地區經濟發展水平高，工業生產規模大，工業用水量和排污量大等因素有關。其次是蘇北地區，平均工業灰水足跡為308.13億m3，這與蘇北地區重工業和化工工業等高污染工業為主的產業發展結構有關，也與蘇北地區2006年以來經濟的迅猛發展速度息息相關。蘇中地區的年均灰水足跡最少，為260.05億m3。就減排潛力而言，蘇北地區0.32的減排潛力大于蘇中地區的0.27和蘇南地區的0.13。并且在減排貢獻率上，蘇北地區貢獻最大，為38.7%；蘇南地區貢獻35.44%，位列其次；蘇中地區貢獻最少，為25.8%。

4結論

在生產前沿函數和規模方向距離函數的基礎上，本文對江蘇省13個地級市的2006—2015年期間年度工業廢水排放績效進行研究，對研究期間績效的時空分布和分布的動態演變趨勢進行分析，得到各地市的工業廢水減排潛力和減排值，并進行排序。結論如下：

（1）研究期間，江蘇省各地市的工業廢水排放績效水平都有所提升。就蘇南、蘇北和蘇中三個地區而言，績效水平與經濟發展水平呈正相關，蘇南地區績效值0.911，為最高；蘇中地區和蘇北地區分居二三位，分別為0.789和0.624；蘇北地區提升速度最快，由0.529增長至0.79，增長率為49.5%。蘇中地區2010—2014年期間績效下降最明顯。

（2）從空間分布和動態演變上看，江蘇省各地區的工業廢水排放績效存在空間分布上的差異，但2006—2015年其趨勢顯示績效差異正在逐步減小，并且績效值向高水平聚攏。

（3）在減排潛力方面，蘇北地區擁有最大的減排潛力0.32，并且提供了整個江蘇省38.77%的減排貢獻值。整個江蘇省的年均減排潛力值為263.5億m3，占年均工業灰水實際排放值的24%。說明江蘇省的工業廢水減排仍存在較大的提升空間。

（編輯：劉照勝）

參考文獻（References）

[1]CHUNG Y H， FRE R， GROSSKOPF S. Productivity and undesirable outputs： a directional distance function approach[J]. Microeconomics， 1995， 51（3）：229-240.

[2]RAMLI N A， MUNISAMY S， ARABI B. Scale directional distance function and its application to the measurement of ecoefficiency in the manufacturing sector[J]. Annals of operations research， 2013， 211（1）：381-398.

[3]FRE R， GROSSKOPF S. Directional distance functions and slacksbased measures of efficiency[J]. European journal of operational research， 2010， 206（1）：320-322.

[4]ZAIM O， FATMA T. Environmental efficiency in carbon dioxide emission in the OECD a nonparametric approach[J]. Journal of environmental management， 2000， 58（2）： 95-143.

[5]ZOFIO J L， PRIETO A M. Environmental efficiency and regulatory standards： the case of CO2 emission from OECD industries[J]. Resource and energy economics， 2001， 23（1）： 63-83.

[6]ZHOU P， WANG H. Energy and CO2 emission performance in electricity generation： a nonradial directional distance function approach[J]. European journal of operational research， 2012，221（3）： 625-635.

[7]WANG Q W， ZHOU P， SHEN N， et al. Measuring carbon dioxide emission performance in Chinese provinces： a parametric approach[J]. Renewable and sustainable energy review， 2013， 21：324-330.

[8]付麗娜，陳曉紅，冷智花.基于超效率DEA模型的城市群生態效率研究——以長株潭“3+5”城市群為例[J].中國人口·資源與環境，2013， 23（4）： 169-175. [FU Lina， CHEN Xiaohong， LENG Zhihua. Urban agglomerations ecoefficiency analysis based on superefficiency DEA Model： case study of ChangZhuTan ‘3+5 urban agglomeration[J]. China population， resources and environment， 2013， 23（4）： 169-175.]

[9]孫作人，周德群，周鵬.工業碳排放驅動因素研究：一種生產分解分析新方法[J].數量經濟技術經濟研究，2012（5） ： 63-74. [SUN Zuoren， ZHOU Dequn， ZHOU Peng. Driving forces of industrial CO2 emissions： a new productionbased decomposition analysis method[J]. The journal of quantitative & technical economics， 2012（5）： 63-74.]

[10]李健，鄧傳霞.基于SDDF的中國省區二氧化碳排放效率及減排潛力測度[J].軟科學，2015， 29（3）：70-73. [LI Jian， DENG Chuanxia. Study on carbon dioxide emission performance and reduction potential of China provinces based on SDDF [J]. Soft science， 2015， 29（3）： 70-73.]

[11]FARRELL M J. The measurement of productive efficiency[J]. Journal of the royal statistical society， 1957： 253-290.

[12]BATTESE G E， CORRA G S. Estimation of production frontier model： with application to the pastoral zone of Eastern Australia[J]. Australian journal of agricultural and economics， 2012， 21（3）： 169-179.

[13]FRE R， GROSSKOPF S， LOVELL C A K， et al. Multilateral productivity comparison when some outputs are undesirable： a nonparametric approach[J]. The review of economics and statistics， 1989， 71（1）：90-98.

[14]王惠，王樹喬.中國工業CO2排放績效的動態演化與空間外溢效應[J].中國人口·資源與環境， 2015， 25（9）：29-36. [WANG Hui， WANG Shuqiao. A study on dynamic evolution of industrial carbon emissions performance and its factors in China [J]. China population， resources and environment，2015， 25（9）： 29-36.]

[15]張軍，章元.對中國資本存量K的再估計[J].經濟研究， 2013（7）： 35-43. [ZHANG Jun， ZHANG Yuan. Reestimation of Chinas stocks of capital[J]. Economic research journal， 2013（7）： 35-43.]

[16]周五七，聶鳴.中國碳排放效率的區域差異研究[J].數量經濟技術經濟研究， 2012（9）： 58-70. [ZHOU Wuqi， NIE Ming. Regional differences in the efficiency of industrial carbon emission in China[J]. The journal of quantitative & technical economics， 2012 （9）： 58-70.]

[17]HOEKSTRA A Y， CHAPAGAIN A K， ALDAYA M M， et al. The water footprint assessment manual： setting the global standard [M]. London， UK： Earthscan， 2011： 27-35.

[18]曾昭，劉俊國.北京市灰水足跡評價[J].自然資源學報， 2013， 28（7）： 1169-1178. [ZENG Zhao， LIU Junguo. Historical trend of grey water footprint of Beijing， China.[J]. Journal of natural resources， 2013， 28（7）： 1169-1178.]

[19]施祖麟，黃治華.基于核密度估計法的中國省區經濟增長動態分析[J].經濟經緯， 2009（4）：60-63. [SHI Zulin， HUANG Zhihua. A dynamic analysis of the economic growth on the provincial level in China based on kernel density estimation method.[J]. Economy survey， 2009 （4）： 60-63.]