利用電參數的工業用戶能效模糊綜合評估

孟金嶺，胡嘉驊，文福拴，林國營，黨三磊

(1.廣東電網有限責任公司電力科學研究院，廣州市 510600;2.浙江大學電氣工程學院，杭州市 310027)

利用電參數的工業用戶能效模糊綜合評估

孟金嶺1，胡嘉驊2，文福拴2，林國營1，黨三磊1

(1.廣東電網有限責任公司電力科學研究院，廣州市 510600;2.浙江大學電氣工程學院，杭州市 310027)

在對工業用戶進行節能改造之前，有必要對其能效水平進行比較準確的評估。現有的一些能效評估方法存在指標選取受限、數據獲取困難、主觀因素較強等問題。在此背景下，針對工業用戶提出基于電參數的能效模糊綜合評估方法。具體地，以計量電表能夠準確測量并記錄的電參數為候選評估指標。 首先，利用主成分分析法(principal component analysis，PCA)識別并刪除次要指標，在此基礎上確定工業用戶綜合能效評估指標體系。之后，利用熵權法計算所選各個指標的權重，并考慮到能效評估結果所呈現的模糊性，采用模糊綜合評估方法對工業用戶的能效進行評估。最后，用廣東省的實際數據說明該方法的可行性與基本特性。

工業用戶；能效評估；電參數；模糊評估；主成分分析法(PCA)；熵權法

0 引 言

近年來，電力系統負荷的峰谷差呈現逐步擴大的趨勢。據統計，2015年我國全社會用電量為55 500億 kW·h，同比增長0.5%。在負荷增長的同時，發電設備累計平均利用小時數卻呈下降趨勢；2015年全國6 000 kW及以上電廠發電設備累計平均利用小時數為3 969 h，同比減少349 h。單純靠增加發電裝機容量來滿足短暫的高峰負荷用電會導致發電設備年利用小時數降低，單位供電成本增加，不利于社會資源的優化配置和合理利用。解決向短暫高峰負荷供電的更經濟的途徑是采用需求側管理[1-4]，通過經濟和技術手段適當調整高峰和低谷負荷時段的用電功率來平滑負荷曲線，改善電力系統運行的安全性和經濟性。能效改造作為需求側管理的重要措施之一，在最近10多年間逐步得到越來越廣泛的應用[5-7]。能效改造通過鼓勵用戶采用先進工藝、先進設備和先進技術，對高耗能低能效的用電設備進行改造，從而提高用電效率，減少電力電量消耗[2]。

能效改造多針對工業用戶，因為其用電負荷和用電量大且高耗能用電設備數量多，具備巨大的能效提升空間。在對工業用戶進行能效改造之前，有必要對其能效水平進行比較準確的評估，以明確對其實施能效改造的必要性。就工業用戶能效評估問題，已有不少研究報道[8-13]。文獻[8]針對配電系統實際運行數據，基于潮流計算和關聯度計算建立了2層節能指標模型，對配電系統節能潛力進行了評估。文獻[9]建立了基于遞階綜合評估框架的能效評估模型，將能效評估過程分解為一級指標評估、二級指標評估和綜合評估3個層次，利用群組判斷法和熵權法對用戶能效進行評估。文獻[10]建立了3層的能效評估指標體系，基于模糊綜合評判原理，對狀態因素集、評語集、權重系數進行分析，在此基礎上建立了能效評估模型。文獻[11]基于層次分析法(analytic hierarchy process，AHP)和逼近理想解排序法(the technique for order of preference by similarity to ideal solution，TOPSIS)，從提升電力系統可靠性與改善負荷輪廓、環境效益、經濟效益、資源利用效率4個方面構建了電力能效項目綜合評估模型。文獻[12]針對高能耗工業用戶，建立了涵蓋經濟、能源、技術、環境的能耗指標評估體系，對能耗水平進行了綜合評估。文獻[13]采用數據包絡分析法，從微觀經濟層面出發，對我國不同地區的工業用戶的能效水平進行了評估。

盡管在能效評估方面已經有不少研究報道，從總體上講，現有的一些能效評估方法存在指標選取受限、數據獲取困難、主觀性程度較高等問題。有鑒于此，本文提出一種基于電參數的工業用戶能效模糊綜合評估方法。這種方法以計量電表能夠準確測量并記錄的電參數為評估指標，采用主成分分析法識別并刪除次要指標，從數據和指標層面保障評估結果的準確性。之后，考慮到能效評估結果一般呈現的模糊性，采用模糊綜合評估方法對用戶能效進行評估。最后，用廣東省的實際數據對所提出的方法做了說明。

1 工業用戶綜合能效評估體系

1.1 基于電參數的初始能效評估體系

電參數指能夠通過計量表計準確測量并記錄的與用電相關的數據(如電壓、電流、電量、功率等)，或根據這些數據進行簡單數學分析或運算后即可得到的衍生數據(如電流不平衡率、年平均電壓、功率因數、單位產值能耗等)。選用電參數作為評價指標能夠保證評價數據充足且比較準確。

在充分考慮指標與工業用戶能效的關聯性以及指標數據的易獲取性的前提下，通過前期調研和數據特征分析確定初始能效評估體系。假設初選能效評估體系包含m個指標，待評估的對象(工業用戶)數量為n，并假定nm，可以構建維數為nm的初始用戶能效矩陣A(1)：

(1)

1.2 預處理工業用戶原始數據

在利用工業用戶的原始數據進行能效綜合評估之前，需要對其進行包括數據類型統一和去量綱化的預處理。A(1)一般包含效益型、成本型和區間型3種類型的評估指標。效益型指標的數值越大，用戶的能效水平越高；成本型指標的數值越小，用戶的能效水平越高；區間型指標的數值處于某一區間時，用戶的能效水平較高，一般稱這個區間為最佳穩定區間。為便于后續數據處理，首先將效益型指標和區間型指標轉換為成本型指標。

對于成本型指標：

(2)

對于效益型指標：

(3)

對于區間型指標：

(4)

(5)

(6)

對工業用戶原始數據進行預處理后，即可把A(1)轉化為標準初始用戶能效矩陣A(2)：

(7)

1.3 最終能效評估體系

主成分分析(principal component analysis，PCA)[14]也稱主分量分析，利用降維思想把多指標轉化為少數幾個綜合指標(即主成分)，其中每個主成分都能反映原始變量的大部分信息，且所含信息互不重復。

每個主成分可表示為原指標的線性組合，其中每個原指標所對應的系數被稱為該指標對此主成分的構成系數。利用PCA對初始評估指標進行篩選，剔除構成系數較小的指標。首先對A(2)求取相關系數矩陣：

(8)

式中：rij=cov(A(2)(:,i),A(2)(:,j))為第i個指標和第j個指標對應的標準化數據集之間的協方差，易知rij=rji；R為對稱矩陣，且對角元素為各個指標數據集的方差；(:，i)和(:，j)分別表示第i列和第j列。

對于前k個特征根λj(j=1,2,…,k)，求取對應的單位特征向量bj=(b1j,b2j...,bmj)T。這k個特征根對應的特征向量構成主成分分析矩陣B：

(9)

按行對矩陣B進行檢查，如果bij(j=1,2,...,k)與同一列的其他元素相比其絕對值較小，則表示指標l在最能反映能效的前k個主成分的構成系數都很小，則指標l對整個指標體系的影響很小，從而該指標可以被刪除。將A(2)中構成系數小的指標刪除后，可得到用矩陣A表示的工業用戶最終能效評估體系：

(10)

式中：A為維數為np的矩陣；p為經過篩選后剩余的指標數量，p≤m；aij表示第i個用戶的第j個指標經標準化處理后的數值。

2 最優權重向量

在利用多個指標進行綜合評估時，如何合理確定各個指標的權重是一個非常重要的問題，因為其可以明顯影響最終評估結果。在確定指標權重時，需要充分考慮各個指標數據間的客觀差異，盡量避免主觀確定權重所導致的偏差。

在確定了評估對象集和評估指標體系之后，可引入信息論中熵的概念來評估各指標在競爭意義上的相對激烈程度(由于所有指標的權重之和為1，從這個意義上講不同指標之間構成了競爭關系)，進而確定各指標的權重。具體地講，當評估對象在某個指標上的值相差較大時，熵值較小，說明該指標提供的有效信息量較大，該指標的權重應給定得較大，反之亦然。在極端情況下，當所有評估對象在某項指標上的值完全相同時，熵值達到最大，意味著該指標對于決策不能提供任何有用信息，可以從評估指標體系中刪除。

引入熵權法[15]確定指標權重，對各個指標的重要程度進行客觀的量化評估。熵權法可將多個評估對象的同一個評估指標進行綜合考慮來確定權重，計及了多個評估對象間的聯系，可削弱異常值的負面影響，得到的評估結果比較合理。

前文確定的工業用戶最終能效評估矩陣A包含p個指標，n個被評估對象。定義第j個評估指標的熵為

(11)

可以證明0≤Hj≤1。第j個指標的熵權定義為

(12)

3 工業用戶能效模糊綜合評估

所提出的工業用戶能效模糊綜合評估方法主要包括選取指標體系、構造評語集、確定隸屬度函數和進行能效模糊綜合評估4個步驟。所提方法既能得到單項指標的得分，又可綜合所有指標的影響得到每個用戶對不同能效評語等級的模糊隸屬度。通過比較不同用戶的能效評估結果，可為用戶實施節能改造的必要性和迫切性提供參考。

3.1 指標體系的選取

根據現有計量電表的測量水平，同時考慮所選電參數與工業用戶能效的相關性，構造了包含“用電量”、“電壓不合格累計時間”、“電流不平衡率”、“萬元產值能耗”、“主要耗能設備能效等級”、“平均負載率”、“功率因數”、“平均電壓”這8個指標的評價體系。

用電量指工業用戶在考察期內的累計用電量；電壓不合格累計時間指專用變壓器一次側電壓在考察期內超出核定偏差范圍的累計小時數；電流不平衡率反映專用變壓器三相電流的不對稱程度；萬元產值能耗指在考察期內工業用戶生產出來的價值1萬元的產品所消耗的電能；主要耗能設備能效等級采用中國能效標識，按照5個等級劃分，1表示產品能效最好，5表示嚴重能耗產品；平均負載率指在考察期內專用變壓器輸出的平均視在功率與變壓器額定容量的比值；功率因數在一定程度上反映專用變壓器的運行效率，為其有功功率與視在功率的比值；平均電壓指專用線路在考察期內的平均電壓，核定偏差范圍假設為-3%～3%。

3.2 評語集的構造

模糊綜合評估方法能否準確評估工業用戶的能效水平在相當程度上取決于能否構造合理的評語集，進而對能效水平進行分類評估。參考國家發改委、國家質檢總局聯合制定的中國能效標識[16]，將工業用戶的能效狀態由高到低分為5個評語集，依次為：綠、藍、黃、橙、紅。各評語集表示的能效狀態情況大致如下：綠色(G)表示工業用戶能效指標非常好，用電設備運行合理，電能質量高，產品具有很強的競爭力；藍色(B)表示工業用戶能效指標良好，用電合理，產品的市場競爭力強；黃色(Y)表示工業用戶能效一般，用電比較合理，產品具備一定市場競爭力；橙色(O)表示工業用戶能效較低，用電不合理，產品的市場競爭力差；紅色(R)表示工業用戶能效較低，用電很不合理，產品市場競爭力很差。

3.3 隸屬度函數的確定

與1.2節將所有指標轉換為成本型指標的處理方法相對應，本節構造典型的成本型指標的隸屬度函數并采用圖形進行直觀闡述。

針對3.1節構造的評語集，采用常見的半梯形分布函數和線性三角函數構造模糊評估的隸屬度函數：

(13)

(14)

(15)

式中：0

式(13)—(15)所描述的3種隸屬度函數如圖1所示。前已述及，它們都是典型的成本型指標的隸屬度函數，指標值越小，則評估對象屬于高能效水平評語集的隸屬程度就越大，反之亦然。假設歸一化之后的指標值取值x∈[0,1]，對應于這5個隸屬度函數的隸屬度分別為fi(x) (i=1,2,..,5)。分析式(13)—(15)并結合圖1，可知成本型指標的隸屬度函數具有下述特征：

(1)對于任意x∈[0,1]，0≤fi(x)≤1 (i=1,2,...,5)；

圖1 成本型指標的隸屬度函數

(2)對于任意x∈[0,1]，在這5個隸屬度函數中滿足fi(x)>0 (i=1,2,..,5)的隸屬度函數最多只有2個。

3.4 能效模糊綜合評估

根據前文確定的工業用戶最終能效評估體系A以及最優權重向量w，對所有待評估的工業用戶依次進行能效模糊綜合評估，步驟如下所述。

(1)建立每個工業用戶的模糊綜合評估矩陣Spq。p為最終采用的能效評估體系所包含的指標數量；q為所構造的評語集中包含的評語數量(q=5)；Spq的第i行Si=(si1，si2,…，si5)表示針對第i個指標的單指標評估向量；sij為第i個指標對評語j的隸屬度。

(2)將最優權重向量與用戶的模糊綜合評估矩陣相乘，求取用戶對不同能效評語集的綜合隸屬度，得到用戶能效評估結果。

(3)在對所有用戶進行綜合評估后，根據最大隸屬原則對不同用戶的能效水平進行比較分析。

在Matlab環境下對上述方法進行了實現。

4 算例

針對廣東省2015年10 kV工業專線用戶的相關統計數據，采用所提出的方法對其中10個用戶進行綜合能效評估。首先對應3.1節選取的8個電參數指標，建立工業用戶初始能效評估體系。在所給出的評估指標體系下，這10個用戶在考察期內(2015年整年)的相關數據如表1所示。按照1.2節的方法對初始數據進行預處理，得到標準化后的用戶能效指標數據，如表2所示。

表1中列出的前5個指標為成本型指標，平均負載率和功率因數為效益型指標，年平均電壓為區間型指標，假定其最佳穩定區間為[374.2，383.8]。

表1 10個待評估工業用戶的指標數據

Table 1 Index data of 10 industrial customers to be evaluated

表2 標準化后的用戶能效指標數據

對標準化后的用戶能效指標數據進行主成分分析，這里選用第一主成分和第二主成分進行指標篩選。求得的2個最大特征根所對應的單位特征向量(即第一主成分和第二主成分的構成系數)分別為b1=(-0.022 7，0.489 0，0.076 5，0.145 4，0.519 2，-0.156 0，-0.191 9，0.634 6)T和b2=(-0.225 2，0.139 6，0.482 2，0.420 4，-0.036 5，0.578 9，0.428 0，0.031 4)T。第一主成分的主要代表指標為電壓不合格累計時間、主要耗能設備能效等級和年平均電壓；第二主成分的主要代表指標為電流不平衡率、萬元產值能耗、平均負載率和功率因素。年用電量在第一主成分和第二主成分的構成系數均較小，作為次要指標予以剔除。這樣，電壓不合格累計時間等余下的7項指標構成了最終的能效評估指標體系。

利用熵權法求得的最終能效評估體系中各指標的最優權重向量為w=(0.274 4，0.070 7，0.099 8，0.135 7，0.092 6，0.199 4，0.127 4)。

假定歸一化之后的指標值對5個評語集的分布是均勻的，即式(13)—(15)中的參數K1～K5依次分別為0.1，0.3，0.5，0.7和0.9，據此求得的用戶1的模糊綜合評估矩陣為

將最優權重向量與模糊綜合評估矩陣相乘，得到用戶1的模糊評估結果為E1=(0.448 3，0.289 4，0.042 2，0.098 2，0.121 8)。

采用同樣方法對用戶2—9進行能效模糊評估，最終得到所有用戶的能效模糊綜合評估結果，如圖2所示。

基于上述評估結果，可得到：

(1)根據最大隸屬度原則，用戶1—6和用戶9的能效評估等級為“綠”，用戶7、8和10的分別為“紅”、“橙”和“藍”；

(2)用戶2、3、5、6對評語“綠”的模糊隸屬度都在0.69以上，對評語“橙”和評語“紅”則均小于0.08(其中用戶2、5、6對這2個評語的模糊隸屬度均為0)，這表明這4個用戶的能效水平較高，不需要進行能效改造；

圖2 10個工業用戶能效模糊綜合評估結果

(3)用戶1、4、9的能效評估等級雖然為“綠”，但對其他4個評語的隸屬度也相對較高，尤其是用戶9，其對評語“紅”的隸屬度大于0.35，這表明這3個用戶的總體能效水平較高，但在個別指標上得分較低，相應地就需要對這些指標進行完善；

(4)用戶8對評語“黃”和評語“橙”的模糊隸屬度較高，表明其能效處于中下水平，有必要對其用電設備和用電習慣進行更深入研究，提高用電效率；

(5)用戶10對評語“綠”和評語“藍”的模糊隸屬度較大，表明其能效處于中上水平，只需要進行適當程度的能效改造；

(6)用戶7對評語“紅”的模糊隸屬度在0.5以上，且對其他4個評語的隸屬度較低，表明其能效水平較低，用電很不合理，需要進行大幅度的能效改造。

5 結 語

針對工業用戶能效評估問題，提出了一種基于電參數的模糊綜合評估方法，具有候選指標充分、所需數據易于獲取且比較準確等優點。首先建立了以電參數為基礎的評估指標體系并利用主成分分析法剔除次要指標，之后求取各指標最優權重向量，進而采用模糊綜合評估方法確定各個用戶的能效水平。最后，采用廣東省工業專線用戶2015年的實際數據對所提出的方法做了說明。通過對工業用戶的能效進行綜合評估，可以讓用戶清楚了解其總體能效水平和各項指標的評估結果，從而對主要設備和用電習慣進行適當調整，達到提高能效水平進而改善運營效益的目的。

[1]BOSHELL F, VELOZA O P. Review of developed demand side management programs including different concepts and their results[C]// Proceedings of IEEE Transmission and Distribution Conference and Exposition. Bogota, Colombia, Latin America：IEEE，2008：1-7.

[2]王蓓蓓, 李揚, 高賜威. 智能電網框架下的需求側管理展望與思考[J]. 電力系統自動化, 2009, 33(20):17-22. WANG Beibei, LI Yang, GAO Ciwei. Demand side management outlook under smart grid infrastructure[J]. Automation of Electric Power Systems, 2009, 33(20): 17-22.

[3]劉道新, 郭萬祝, 王世成, 等. 面向能源互聯網的新型需求側管理模式研究[J]. 電力建設, 2016, 37(6):10-16. LIU Daoxin, GUO Wanzhu, WANG Shicheng, et al. New demand-side management pattern for energy internet[J]. Electric Power Construction, 2016, 37(6):10-16.

[4]于長廣, 劉穎. 提高技術水平實現電力需求側有效管理[J]. 電力建設, 2007, 28(5):80-83. YU Changguang, LIU ying. Improving technology to achieve effective electric power demand side management[J]. Electric Power Construction, 2007, 28(5):80-83.

[5]WORRELL E, BERNSTEIN L, ROY J, et al. Industrial energy efficiency and climate change mitigation[J]. Energy Efficiency, 2009, 2(2): 109-123.

[6]王丹, 范孟華, 賈宏杰. 考慮用戶舒適約束的家居溫控負荷需求響應和能效電廠建模[J]. 中國電機工程學報, 2014, 34(13): 2071-2077. WANG Dan, FAN Menghua, JIA Hongjie. User comfort constrain demand response for residential thermostatically-controlled loads and efficient power plant modeling[J]. Proceedings of the CSEE, 2014, 34(13): 2071-2077.

[7]于建成, 王旭東, 張東,等. 面向多元能源互聯的天津生態城智能電網創新示范區建設[J]. 電力建設, 2015, 36(11):58-63. YU Jiancheng, WANG Xudong, ZHANG Dong, et al. Innovation demonstration area construction of Tianjin eco-city smart grid facing multi energy interconnection[J]. Electric Power Construction, 2015, 36(11): 58-63.

[8]王煜晗, 范瑞祥, 董泰青, 等. 基于兩層指標模型的配電網節能潛力分析[J]. 電力建設, 2014, 35(12): 94-99. WANG Yuhan, FANG Ruixiang, DONG Taiqing, et al. Energy-saving potential analysis of distribution network based on two layered index models[J]. Electric Power Construction, 2014, 35(12): 94-99.

[9]羅耀明, 毛李帆, 姚建剛, 等. 電力用戶綜合能效評估模型[J]. 電力系統及其自動化學報, 2011, 23(5): 104-109. LUO Yaoming, MAO Lifan, YAO Jiangang, et al. Evaluationmodel of integrated energy efficiency for power users[J]. Proceedings of the CSU-EPSA, 2011, 23(5): 104-109.

[10]吳劍飛, 姚建剛, 陳華林, 等. 電力客戶能效狀態模糊綜合評估[J]. 電力系統保護與控制, 2010, 38(13): 94-98. WU Jianfei, YAO Jiangang, CHEN Hualin, et al. Fussy comprehensive evaluation of power customer energy efficiency assessment[J]. Power System Protection and Control, 2010, 38(13): 94-98.

[11]李慧玲, 蘆新波, 劉大川, 等. 基于AHP-TOPSIS的電力能效項目綜合評價[J]. 現代電力, 2014, 31(4): 88-94. LI Huiling, LU Xinbo, LIU Dachuan, et al. A comprehensive evaluation method for power energy efficiency project based on AHP-TOPSIS[J]. Modern Electric Power, 2014, 31(4): 88-94.

[12]徐勇, 王堅. 基于AHP和PCA的企業能耗水平評價[J]. 價值工程, 2008, 27(8): 9-12. XU Yong, WANG Jian. To evaluate on the level of enterprise energy consumption based on AHP and PCA[J]. Value Engineering, 2008, 27(8): 9-12.

[13]SHI G M, BI J, WANG J N. Chinese regional industrial energy efficiency evaluation based on a DEA model of fixing non-energy inputs[J]. Energy Policy, 2010, 38(10): 6172-6179.

[14]BRO R, SMILDE A K. Principal component analysis[J]. Analytical Methods, 2014, 6(9): 2812-2831.

[15]聶宏展, 呂盼, 喬怡, 等. 基于熵權法的輸電網規劃方案模糊綜合評價[J]. 電網技術, 2009, 33(11): 60-64. NIE Hongzhan, LYU Pan, QIAO Yi, et al. Comprehensive fuzzy evaluation for transmission network planning scheme based on entropy weight method[J]. Power System Technology, 2009, 33(11): 60-64.

[16]國家發展改革委, 國家質檢總局. 能源效率標識管理辦法[EB/OL]. (2016-2-29) [2016-11-01]. http://www.sdpc.gov.cn/gzdt/201603/ W020160309388699982349.pdf.

(編輯 張媛媛)

Fuzzy Comprehensive Evaluation for Energy Efficiency of Industrial Consumers Employing Electric Parameters

MENG Jinling1, HU Jiahua2, WEN Fushuan2, LIN Guoying1, DANG Sanlei1

(1.Electric Power Research Institute of Guangdong Power Grid Co., Ltd., Guangzhou 510600, China; 2.School of Electrical Engineering, Zhejiang University, Hangzhou 310027, China)

It is necessary to accurately evaluate the energy efficiency of an industrial consumer before implementing an energy-saving scheme. Several existing evaluation methods for energy efficiency have some problems such as limited selection of indices, limited access to data and strong subjectivity. Given this background, this paper presents a fuzzy comprehensive evaluation method for the energy efficiency of industrial consumers based on electric parameters. Specifically, the electric parameters are candidate evaluation index, which can be accurately measured and recorded by smart meters. Firstly, we adopt principal component analysis (PCA) to identify and remove secondary indicators, on this basis determine the comprehensive energy efficiency evaluation index system for industrial consumers. Then, we calculate the weights of all selected indices with using entropy weight method and evaluate the energy efficiency of industrial consumer with using fuzzy evaluation method, considering the exhibited fuzziness of evaluation results. Finally, the actual data from Guangdong province are served for demonstrating the feasibility and features of the proposed method.

industrial consumer; energy efficiency evaluation; electric parameters; fuzzy evaluation; principal component analysis (PCA), entropy weight method

國家自然科學基金項目(51477151)；南方電網公司科研項目(K-GD2014-192)

TM 72

A

1000-7229(2016)11-0016-07

10.3969/j.issn.1000-7229.2016.11.003

2016-09-05

孟金嶺(1986)，男，碩士，工程師，研究方向為電力系統諧波治理、電力需求側綜合節能技術；

胡嘉驊(1992)，男，博士研究生，主要從事電能交易、需求側與能效管理、可再生能源等方面的研究工作；

文福拴(1965)，男，教授，博士生導師，本文通信作者，主要從事電力系統故障診斷與系統恢復、電力經濟與電力市場、智能電網與電動汽車等方面的研究工作；

林國營(1983)，男，碩士，高級工程師，研究方向為智能用電、數字化電能計量技術；

黨三磊(1982)，男，碩士，高級工程師，研究方向為電力市場與智能用電技術。

Project supported by National Natural Science Foundation of China (51477151)