基于貝葉斯MCMC 方法的城市建筑供熱能耗基準(zhǔn)

那威，邴佳樂(lè)，劉品妍

(1.北京建筑大學(xué) 環(huán)境與能源工程學(xué)院，北京 100044；2.北京建筑大學(xué) 建筑與城市規(guī)劃學(xué)院，北京 100044)

我國(guó)民用建筑終端能耗和碳排放分別占全國(guó)總量的21.2%和21.6%，其中北方城鎮(zhèn)建筑供熱能耗占建筑終端能耗已近1/4，是建筑領(lǐng)域?qū)崿F(xiàn)“雙碳”規(guī)劃目標(biāo)的重點(diǎn)領(lǐng)域.德、英、美、丹麥等國(guó)已經(jīng)應(yīng)用整體能效法對(duì)建筑能耗進(jìn)行管控，而我國(guó)目前仍處于向能耗限額法和整體能效法的過(guò)渡階段[1].建筑能耗限額、整體能效分析是建筑能耗基準(zhǔn)分析的基礎(chǔ)環(huán)節(jié).不同時(shí)空尺度建筑實(shí)際能耗與能耗基準(zhǔn)之間的對(duì)比，可以作為建筑設(shè)計(jì)、建造及運(yùn)行階段用能情況的評(píng)估方法[2].建筑能耗基準(zhǔn)分析已成為現(xiàn)階段建筑全過(guò)程節(jié)能優(yōu)化領(lǐng)域的研究熱點(diǎn).

建筑能耗預(yù)測(cè)分析方法對(duì)能耗基準(zhǔn)指定的支撐顯得尤為重要.大量學(xué)者開展對(duì)建筑能耗預(yù)測(cè)模型、方法和模式的研究及應(yīng)用實(shí)踐.在國(guó)家、省級(jí)區(qū)域尺度上，德國(guó)通過(guò)VDI 3807 標(biāo)準(zhǔn)統(tǒng)計(jì)全國(guó)調(diào)查樣本的分類建筑實(shí)際能耗概率，概率分布中前四分位數(shù)、中位數(shù)作為能耗標(biāo)準(zhǔn)值和指導(dǎo)值[3].英國(guó)、澳大利亞等國(guó)實(shí)施的EEBPP 項(xiàng)目、step to benchmark energy use 項(xiàng)目采用樣本數(shù)據(jù)調(diào)查和統(tǒng)計(jì)分析結(jié)合方法確定建筑基準(zhǔn)能耗[2,4].在單體建筑、街區(qū)等尺度上，部分學(xué)者提出利用物理正演模擬、正向模型仿真來(lái)描述建筑能耗與建筑特征參數(shù)之間的“目標(biāo)-特征”關(guān)系[5]，采用物理仿真模型和計(jì)算流體力學(xué)技術(shù)結(jié)合的能耗模擬引擎和工具，如DOE2[6]、Energy Plus[7]、Design Builder[8]等；另一部分學(xué)者提出數(shù)據(jù)驅(qū)動(dòng)的黑箱或反向模型描述上述數(shù)學(xué)聯(lián)系，建立神經(jīng)網(wǎng)絡(luò)[9]、廣義線性回歸[10-11]、人工智能[12]等黑箱模型，通過(guò)帶有建筑特征和能耗響應(yīng)數(shù)據(jù)的建筑仿真樣本或建筑實(shí)際樣本對(duì)模型進(jìn)行訓(xùn)練及校驗(yàn)，如支持向量機(jī)（support vector machines，SVM）[13]、決策樹和隨機(jī)森林等基于樹的模型[14-15]、人工神經(jīng)網(wǎng)絡(luò)（artificial neutral network，ANN）[16]等.部分學(xué)者結(jié)合物理正演和數(shù)據(jù)驅(qū)動(dòng)算法[17]提出白、黑箱模型結(jié)合的混合模型，如阻容（RC）熱網(wǎng)絡(luò)模型等[18-19].上述建筑能耗預(yù)測(cè)方法用于城市尺度建筑能耗基準(zhǔn)分析，仍存在一些不足之處：一方面我國(guó)研究中仍缺乏自上而下的統(tǒng)計(jì)模型來(lái)計(jì)算能耗基準(zhǔn)所需的數(shù)據(jù)基礎(chǔ),統(tǒng)計(jì)數(shù)據(jù)存在時(shí)空分辨率較低的問(wèn)題，在統(tǒng)計(jì)口徑、范圍、完整性和準(zhǔn)確性方面仍存在諸多爭(zhēng)議[20-22]；另一方面通過(guò)自下而上的模型發(fā)現(xiàn)，若采用物理正演模型，模型構(gòu)建中所需建筑形態(tài)等建?；A(chǔ)信息量較大且通常難以獲得，建模難度高[23].若采用數(shù)據(jù)驅(qū)動(dòng)模型或者混合模型，在訓(xùn)練模型采用仿真樣本時(shí)，數(shù)值仿真建模和求解所需計(jì)算資源和時(shí)間成本通常過(guò)高[5]，采集樣本和建模所需的基礎(chǔ)信息量、建模難度、計(jì)算資源和時(shí)間成本有所下降[24-25].因此，如何在城市尺度建筑樣本數(shù)量和維度有限的情況下，降低校準(zhǔn)模型有限輸入和不確定性對(duì)預(yù)測(cè)性能的影響，仍是目前待解決的問(wèn)題.

本研究針對(duì)有限樣本和輸入不確定性的處理及計(jì)算效率優(yōu)化，提出基于貝葉斯框架的混合效應(yīng)模型.將熱力站供熱能耗強(qiáng)度監(jiān)測(cè)數(shù)據(jù)作為測(cè)試樣本集，并利用馬爾科夫鏈蒙特卡洛（Markov chain Monte Carlo，MCMC）方法進(jìn)行重要性采樣，對(duì)特征空間實(shí)現(xiàn)全局搜索以推斷模型以及超參數(shù)峰值概率的分布，尋找指定建筑特征的預(yù)測(cè)能耗全局最優(yōu)解鄰域，以最大程度考慮城市尺度全局建筑特征的供熱能耗分布，將模型EUIH 點(diǎn)估計(jì)和區(qū)間估計(jì)結(jié)果作為能耗基準(zhǔn)確定依據(jù).

1 能耗預(yù)測(cè)模型構(gòu)建

對(duì)比典型能耗預(yù)測(cè)模型的所需信息量、建模難度、樣本量、計(jì)算資源和時(shí)間成本，如圖1 所示，圖中圓形面積大小代表相應(yīng)模型所需的樣本量.城市尺度建筑能耗預(yù)測(cè)和校準(zhǔn)的主要挑戰(zhàn)是如何在建筑樣本有限、輸入不確定和建模難度、求解計(jì)算資源和時(shí)間成本條件下構(gòu)建合理模型并求解.

1.1 建筑能耗預(yù)測(cè)模型框架

在建筑使用過(guò)程中，供熱能耗影響因素眾多且龐雜，因素作用效應(yīng)各異，受采集樣本缺乏、反映建筑或熱力站能耗特征的數(shù)據(jù)維度不足的訓(xùn)練數(shù)據(jù)限制，難以充分觀察全部能耗因素及其數(shù)據(jù)水平.將供熱能耗影響因素盡然納入建模過(guò)程，且需考慮建模中難以觀測(cè)的遺漏變量、采集范圍

抽樣偏差、采集數(shù)據(jù)偏誤等數(shù)據(jù)不確定性對(duì)模型準(zhǔn)確性的影響.根據(jù)已有數(shù)據(jù)庫(kù)的數(shù)據(jù)量和數(shù)據(jù)維度情況，參考能耗預(yù)測(cè)[26-30]和循證醫(yī)學(xué)領(lǐng)域[31-35]類似多因素、少樣本、因素異質(zhì)性、不確定性和因素效應(yīng)耦合下建模的探索.本研究構(gòu)建了基于貝葉斯框架建筑能耗預(yù)測(cè)的混合效應(yīng)概率模型，模型框架為一種采用隨機(jī)待估超參數(shù)的分層、變截距、變斜率、廣義線性固定效應(yīng)的數(shù)據(jù)驅(qū)動(dòng)模型,目標(biāo)函數(shù)為

式中：yf(x) 為建筑EUIH 的觀測(cè)數(shù)據(jù),x為建筑特征相關(guān)的輸入變量, θ 為模型調(diào)整超參數(shù)(θm,b,εf)向量,yr(x,θr) 為建筑EUIH 的預(yù)測(cè)值,ym(x,θ) 為模型固定效應(yīng)結(jié)構(gòu)項(xiàng),b(x) 為模型個(gè)體特質(zhì)效應(yīng)項(xiàng),εf為模型預(yù)測(cè)的獨(dú)立隨機(jī)誤差.

模型中固定效應(yīng)結(jié)構(gòu)項(xiàng)利用現(xiàn)有數(shù)據(jù)庫(kù)中的有限數(shù)據(jù)維度信息，反映出建筑共性解釋變量對(duì)能耗的真實(shí)效應(yīng)量，包括城市氣候、建筑分布、建造狀況和樣本的供熱服務(wù)建筑規(guī)模共4 個(gè)維度數(shù)據(jù).1）樣本的建筑建造狀況維度相關(guān)解釋變量.輸入樣本熱力站所服務(wù)居住建筑的建造節(jié)能水平，根據(jù)居住建筑所采用的節(jié)能設(shè)計(jì)標(biāo)準(zhǔn)情況（非節(jié)能或采用30%、50%、65% 居住建筑節(jié)能設(shè)計(jì)標(biāo)準(zhǔn)）采用啞變量分別帶入.2）城市氣候維度相關(guān)解釋變量.本研究根據(jù)樣本供暖期的供熱需求，輸入反映熱力站當(dāng)?shù)毓┡谑彝鈿鉁睾凸崽鞌?shù)的綜合參數(shù)，采用樣本供暖期度日數(shù)為解釋變量之一.供暖期度日數(shù)指標(biāo)是樣本室內(nèi)基準(zhǔn)18 ℃與采暖期室外平均溫度之間的溫差，乘以當(dāng)年采暖期延續(xù)時(shí)間（d）.3）城市建筑分布維度相關(guān)解釋變量.根據(jù)不同建造水平的建筑在城市建筑中的分布狀況，輸入指定建造狀況建筑相關(guān)熱力站樣本占樣本熱力站總數(shù)占比.4）樣本的供熱服務(wù)建筑規(guī)模維度相關(guān)解釋變量.輸入熱力站樣本服務(wù)居住建筑的供熱面積.

模型中個(gè)體特質(zhì)效應(yīng)項(xiàng)反映了難以識(shí)別或觀測(cè)的潛變量對(duì)樣本能耗的效應(yīng)，包括建筑人員活動(dòng)和行為變化、建筑用能設(shè)施設(shè)備配置和使用建筑建造水平等潛變量的相關(guān)效應(yīng).模型中調(diào)整超參數(shù)為模型待估參數(shù)，區(qū)別于頻率學(xué)派極大似然估計(jì)（maximum likelihood estimation，MLE）等算法將模型待估參數(shù)視為固定值.在貝葉斯推斷中，模型固定效應(yīng)結(jié)構(gòu)項(xiàng)及個(gè)體特質(zhì)效應(yīng)項(xiàng)的待估參數(shù)均為隨機(jī)參數(shù).在后續(xù)計(jì)算中，通過(guò)樣本數(shù)據(jù)訓(xùn)練，應(yīng)用貝葉斯推斷和MCMC 采樣推斷出上述參數(shù)的后驗(yàn)概率分布.實(shí)際樣本數(shù)據(jù)提供直接、原始的建筑能耗與不確定性因素之間的關(guān)系，納入建筑實(shí)際運(yùn)行中很多非理想的因素，模型輸出反映輸入不確定性的隨機(jī)效應(yīng)，包括觀測(cè)誤差、采樣偏誤和其他參數(shù)異質(zhì)性因素.該模型通過(guò)尋找具有峰值概率的點(diǎn)估計(jì)表示參數(shù)估計(jì)值的可信度和不確定性，用以控制目標(biāo)預(yù)測(cè)分布的特征.

為了消除原始連續(xù)型變量之間因性質(zhì)、量綱、數(shù)量級(jí)等屬性的差異而帶來(lái)的影響，采用自然對(duì)數(shù)映射方式將模型變量轉(zhuǎn)化為無(wú)量綱的標(biāo)準(zhǔn)化數(shù)值，從而在不改變變量特征的同時(shí)，使得不同變量的效應(yīng)大小具有可比性.上述混合模型中基于貝葉斯框架的模型對(duì)數(shù)映射為

1.2 后驗(yàn)分布構(gòu)建

根據(jù)貝葉斯原理，模型利用觀測(cè)數(shù)據(jù)，并將似然分布和先驗(yàn)分布結(jié)合，推斷、構(gòu)建模型的后驗(yàn)分布函數(shù).將模型調(diào)整超參數(shù)視為隨機(jī)參數(shù)[28,36-37]，觀測(cè)數(shù)據(jù)和建筑特征相關(guān)參數(shù)視為模型輸出和輸入的隨機(jī)變量，基于貝葉斯原理的后驗(yàn)分布函數(shù)為

式中：P(θ|y) 為模型超參數(shù)后驗(yàn)分布；P(y|θ) 為模型觀測(cè)數(shù)據(jù)的似然分布；P(θ) 為模型超參數(shù)先驗(yàn)分布；P(y) 為觀測(cè)數(shù)據(jù)的邊緣分布，該值是與參數(shù)θ無(wú)關(guān)的標(biāo)準(zhǔn)化常數(shù)，可以忽略P(y) 的具體數(shù)值.

通過(guò)搜索P(θ|y)P(θ) 的峰值概率值，從而確定后驗(yàn)分布的峰值概率值，最大后驗(yàn)（maximum a posteriori, MAP）估計(jì)為

極大似然估計(jì)(maximum likelihood estimation,MLE)或MAP 均使用點(diǎn)估計(jì)來(lái)推導(dǎo)后驗(yàn)分布，但是MLE 的不確定性由漸近特征決定，在計(jì)算過(guò)程中不納入先驗(yàn)信息，完全利用樣本統(tǒng)計(jì)量估計(jì)總體參數(shù).若是MAP 估計(jì)的先驗(yàn)分布設(shè)定為均勻分布，則MAP 的估計(jì)結(jié)果將與MLE 的估計(jì)結(jié)果相同，即

城市尺度建筑能耗模型的預(yù)測(cè)和校準(zhǔn)過(guò)程往往會(huì)有數(shù)據(jù)不足的問(wèn)題[38].MLE 估計(jì)的漸近特性將導(dǎo)致當(dāng)觀測(cè)數(shù)據(jù)不足時(shí)不能很好地利用有限的數(shù)據(jù)，從而容易導(dǎo)致過(guò)擬合和一類誤差問(wèn)題[36-37].因此，提出在貝葉斯框架下與MAP 估計(jì)相結(jié)合的建模方式是提高有限樣本信息利用以及提高預(yù)測(cè)效果的關(guān)鍵改進(jìn).

1.3 先驗(yàn)分布的構(gòu)建

先驗(yàn)分布P(θ) 為模型中超參數(shù)的客觀特征，通常根據(jù)經(jīng)驗(yàn)判斷方式給定，相關(guān)研究認(rèn)為對(duì)先驗(yàn)分布的不同設(shè)定將影響后驗(yàn)分布的計(jì)算過(guò)程，但是否會(huì)造成貝葉斯估計(jì)結(jié)果的顯著差異仍有爭(zhēng)論[39-40].在模型超參數(shù)的先驗(yàn)分布設(shè)定中采用無(wú)信息先驗(yàn).無(wú)信息先驗(yàn)等模糊先驗(yàn)與信息先驗(yàn)相比，未對(duì)后驗(yàn)分布推斷進(jìn)行約束，使得超參數(shù)的推斷處于可以接受的范圍內(nèi)[36-37].通過(guò)模糊先驗(yàn)歸一化，可以適應(yīng)觀測(cè)數(shù)據(jù)存在信息缺失、信息弱或稀疏的情況.模型中超參數(shù) π1、π2、π3、π4和π0的先驗(yàn)分布均設(shè)定為高斯分布G(0,1×10-6) ；超參數(shù)b的先驗(yàn)分布設(shè)定為伽馬分布 Ga(αf,βf) ，其中形狀參數(shù) αb、尺度參數(shù) βb均設(shè)為1 ×10-3；模型誤差 εf為精度因子 σf的函數(shù)，精度因子 σf的先驗(yàn)分布設(shè)定為 Ga(1×10-3,1×10-3) ，具體計(jì)算式為

1.4 MCMC 重要性采樣和推斷

模型超參數(shù)的聯(lián)合后驗(yàn)分布P(θ|y) 整合了先驗(yàn)信息和觀測(cè)數(shù)據(jù)的似然信息，通常具有較復(fù)雜的函數(shù)形式，如式（3）所示，幾乎無(wú)法求出并從其分布中進(jìn)行直接采樣.采取MCMC 方法從后驗(yàn)分布P(θ|y) 中進(jìn)行重要性采樣，當(dāng)生成的馬爾科夫鏈達(dá)到收斂時(shí)，用所獲得平穩(wěn)分布作為所求聯(lián)合后驗(yàn)分布的估計(jì)結(jié)果[28,36]，MCMC 方法流程如圖2 所示.

重要性采樣中采用梅特羅波利斯-黑斯廷斯（Metropolis Hastings，MH）算法生成3 條馬爾可夫鏈[28,36].在MH 算法中采用細(xì)致平衡條件Jt(θ?|θ(t-1))=Jt(θ(t-1)|θ?)推出平穩(wěn)分布，獲取所需重要性樣本的分布即為所求后驗(yàn)分布.該算法的主要流程為

1）獲取 θ(0)，給參數(shù)向量賦初值.

2）獲取 θ(t)，第t步迭代抽樣中給參數(shù)向量賦值的步驟：從建議概率分布Jt(θ?|θ(t-1)) 中采樣，獲得新候選參數(shù)向量 θ?；從均勻分布中采樣隨機(jī)數(shù)u服從[0,1]上的均勻分布.引入并計(jì)算接受率 αLC為

3）判斷是否接收新候選向量參數(shù) (接受轉(zhuǎn)移的新候選向量參數(shù)時(shí)αLC＞u，拒絕轉(zhuǎn)移的新候選向量參數(shù)時(shí)αLC＜u)得：

該算法通過(guò)完成指定迭代次數(shù)并判斷馬爾科夫鏈達(dá)到收斂后結(jié)束[28,36]，從而使遍歷馬爾可夫鏈的極限分布收斂于平穩(wěn)分布.

2 訓(xùn)練集數(shù)據(jù)來(lái)源

模型訓(xùn)練集包括北京市223 個(gè)居住建筑熱力站隨機(jī)采集樣本，樣本來(lái)自北京某大型熱力公司，該公司采用熱力站供熱的民用建筑面積為950 hm2，共有熱力站總量5 875 個(gè).模型訓(xùn)練集樣本按照樣本EUIH 值區(qū)間、樣本中居住建筑所采用居住建筑節(jié)能設(shè)計(jì)標(biāo)準(zhǔn)的分布情況見表1，其中S為樣本量，Sp為樣本量占比，A為樣本供熱面積，Ap為樣本供熱面積占比.

本研究通過(guò)熱力站能耗監(jiān)測(cè)系統(tǒng)自動(dòng)抄表數(shù)據(jù)庫(kù)中匯總的逐時(shí)能耗數(shù)據(jù)，將年供熱能耗強(qiáng)度的數(shù)據(jù)進(jìn)行整理，包括2017、2018、2019 年共3 年供暖季，同時(shí)采集接入相應(yīng)熱力站的建筑相關(guān)特征數(shù)據(jù).考慮到數(shù)據(jù)采集設(shè)備可能存在運(yùn)行異常等情況，造成所采集的數(shù)據(jù)存在異常值和缺失值，為了保證樣本數(shù)據(jù)質(zhì)量，采用 3σ 準(zhǔn)則[41]對(duì)數(shù)據(jù)異常值進(jìn)行剔除，采用拉格朗日插值法[42]對(duì)數(shù)據(jù)缺失值進(jìn)行插補(bǔ).

3 模型收斂判斷與驗(yàn)證

3.1 模型收斂判斷

通過(guò)馬爾科夫鏈蒙特卡洛（MCMC）方法模擬重要性采樣，獲得的樣本可用前提是馬爾科夫鏈達(dá)到收斂，收斂性可以通過(guò)計(jì)算并檢查潛在尺度縮減因子R? （部分文獻(xiàn)稱為bgr（Brooks-Gelman-Rubin）指標(biāo)）進(jìn)行判斷[43-44].當(dāng)模擬開始時(shí)，M條馬爾科夫鏈在N個(gè)樣本中，進(jìn)行連續(xù)并隨機(jī)采樣，利用每條鏈的軌跡圖來(lái)檢查和設(shè)定退火時(shí)間Nburn-in以去除退火內(nèi)的迭代，跳過(guò)不穩(wěn)定階段.在每條鏈劃分為2 個(gè)階段，M+M鏈中均有(N-Nburn-in)/2的時(shí)間間隔，若是M+M鏈在2 階段中，迭代結(jié)果差異已經(jīng)達(dá)到設(shè)定精度要求，則說(shuō)明該鏈趨于平穩(wěn)，迭代N次已經(jīng)足夠，否則迭代應(yīng)繼續(xù)進(jìn)行.所需迭代次數(shù)根據(jù)因子R? 確定，因子R? 降至1.1 時(shí)停止迭代，當(dāng)達(dá)到1.0 時(shí)的效果最佳，經(jīng)驗(yàn)認(rèn)為鏈中迭代數(shù)應(yīng)至少超過(guò)5×(M+M)[36,45].

本研究模型構(gòu)建和相關(guān)算法采用R 語(yǔ)言和R2winbugs 程序包編寫和實(shí)現(xiàn).計(jì)算中共進(jìn)行30 000 次迭代，其中退火去除初始10 000 次迭代，根據(jù)達(dá)到平穩(wěn)分布后的樣本對(duì)模型中的參數(shù)進(jìn)行估計(jì).以超參數(shù) π4為例，圖3 為3 條馬爾科夫鏈迭代過(guò)程的潛在尺度縮減因子軌跡，t為迭代步數(shù)，左軸為 π4的潛在尺度縮減因子，D為馬爾科夫鏈中采樣80%分布區(qū)間的標(biāo)準(zhǔn)化寬度，可見退火10 000次迭代后該因子迭代中低于1.1 且接近1.0.

圖3 潛在尺度縮減因子軌跡Fig.3 Trace of potential scale reduction factor

圖4為3 條馬爾可夫鏈的迭代軌跡，可見3 條馬爾科夫鏈的軌跡歷史明顯重疊.上述結(jié)果顯示3 條馬爾可夫鏈的重要性采樣已收斂，模型中其他超參數(shù)均具有相似收斂特征.

圖4 3 條馬爾可夫鏈迭代軌跡圖Fig.4 Trace plot of 3 Markov chains in iterations

3.2 模型驗(yàn)證

采用預(yù)測(cè)殘差自相關(guān)性檢驗(yàn)驗(yàn)證模型，以評(píng)估所構(gòu)建模型描述輸入和輸出關(guān)系的可靠性.圖5為模型參數(shù)估計(jì)中3 條馬爾可夫鏈殘差自相關(guān)函數(shù)， ACF 為自相關(guān)函數(shù)，L為迭代滯后階數(shù).可見迭代過(guò)程中3 條馬爾可夫鏈的自相關(guān)函數(shù)均呈現(xiàn)明顯向零值震蕩衰減趨勢(shì)，說(shuō)明殘差已與歷史輸入獨(dú)立.

圖5 模型參數(shù)估計(jì)中3 條馬爾可夫鏈殘差自相關(guān)函數(shù)Fig.5 ACF of residuals of model from 3 Markov chains

進(jìn)一步利用新測(cè)試集數(shù)據(jù)對(duì)模型進(jìn)行泛化能力測(cè)試.通過(guò)獲取模型訓(xùn)練集數(shù)據(jù)相同數(shù)據(jù)渠道，隨機(jī)采集2017—2019 年北京市48 座居住建筑熱力站能耗監(jiān)測(cè)數(shù)據(jù)為測(cè)試集樣本.圖6 為測(cè)試集熱力站能耗EUIH 的觀測(cè)值和預(yù)測(cè)值對(duì)比，圖中下軸為EUIH 的觀測(cè)值，左軸為EUIH 的預(yù)測(cè)值；短虛線和長(zhǎng)虛線分別為模型5%和10%的絕對(duì)百分比誤差.

圖6 測(cè)試集熱力站能耗EUIH 的觀測(cè)值和預(yù)測(cè)值對(duì)比Fig.6 Comparison between metered EUIH with calibrated EUIH for heating substation samples in validation dataset of model

參考關(guān)于廣泛用于評(píng)估能耗校準(zhǔn)模型和預(yù)測(cè)模型的ASHRAE 指南14-2004 標(biāo)準(zhǔn)[46]，采用絕對(duì)百分比誤差（percentage error, PE）、歸一化平均偏差（normalized mean bias error, NMBE）、均方根誤差（root mean sqnare error, RMSE）和均方根誤差變異系數(shù)（cofficient of variation root mean squaren error, CVRMSE）等預(yù)測(cè)精度統(tǒng)計(jì)評(píng)估指標(biāo)進(jìn)行誤差分析，即

式中：RMSE 為指標(biāo)衡量觀測(cè)值與預(yù)測(cè)值之間的偏差；CVRMSE 為衡量觀測(cè)值與預(yù)測(cè)值之間的無(wú)量綱偏差；NMBE 為衡量觀測(cè)值與預(yù)測(cè)值之間的歸一化偏差，上述3 個(gè)指標(biāo)數(shù)值越小說(shuō)明預(yù)測(cè)精度更高;y?j為樣本EUIH 預(yù)測(cè)值, 單位為GJ/m2；yj為樣本EUIH 觀測(cè)值, 單位為GJ/m2；Nt為測(cè)試集樣本總數(shù)；yˉj為測(cè)試集樣本EUIH 觀測(cè)值均值, 單位為GJ/m2.

預(yù)測(cè)精度指標(biāo)公式結(jié)果顯示NMBE 為0.4%，RMSE 為0.034，CVRMSE 為13.1%，模型的CVRMSE和NMBE 指標(biāo)已經(jīng)明顯低于ASHRAE 指南14-2014標(biāo)準(zhǔn)中建議CVRMSE 限值（15%）和NMBE 限值（5%），說(shuō)明本研究模型進(jìn)行居住建筑熱力站EUIH預(yù)測(cè)具有較強(qiáng)魯棒性和泛化能力.

4 模型結(jié)果分析

以北京市223 個(gè)居住建筑熱力站隨機(jī)采集樣本的2017、2018、2019 年共3 年供暖季逐時(shí)EUIH數(shù)據(jù)為基礎(chǔ)，結(jié)合本研究構(gòu)建的模型和能耗預(yù)測(cè)方法，使用重要性采樣達(dá)到平穩(wěn)分布后的迭代樣本，推斷模型中參數(shù)和EUIH 的聯(lián)合后驗(yàn)分布.圖7為模型中參數(shù)后驗(yàn)分布計(jì)算結(jié)果，其中左軸為模型參數(shù)的后驗(yàn)分布概率密度ρ，下軸 π0～π4為模型結(jié)構(gòu)項(xiàng)解釋變量固定效應(yīng)彈性系數(shù)計(jì)算結(jié)果， π0～π4分別為-6.852、0.495、0.170、-0.512、0.484；b1～b3為個(gè)體特質(zhì)效應(yīng)項(xiàng)潛變量固定效應(yīng)彈性系數(shù)計(jì)算結(jié)果，b1～b3分別為0.047、0.039、0.037，分別對(duì)應(yīng)居住建筑熱力站樣本2017、2018、2019 年潛變量固定效應(yīng)情況.π0～ π4、b1～b3的計(jì)算結(jié)果與式（2）聯(lián)合，進(jìn)一步推斷出樣本所在北京市居住建筑EUIH的聯(lián)合后驗(yàn)分布替代模型.北京市不同居住建筑群采用該聯(lián)合后驗(yàn)分布替代模型快速預(yù)測(cè)EUIH值，通過(guò)圖7 的參數(shù)計(jì)算結(jié)果與式（2）同時(shí)獲得EUIH 的峰值概率期望值和不同可信區(qū)間閾值.

圖7 混合效應(yīng)模型參數(shù)后驗(yàn)分布Fig.7 Posterior distributions of parameters in mixed effects model

圖8 為北京市居住建筑EUIH 的聯(lián)合后驗(yàn)分布計(jì)算結(jié)果，左軸為EUIH 基準(zhǔn)值EUIHB（GJ/m2），右軸為蒙特卡洛標(biāo)準(zhǔn)誤差E（GJ/m2），下軸為年份.聯(lián)合后驗(yàn)分布具有峰值概率的EUIH 期望值，可以設(shè)定作為北京市居住建筑EUIH 的基準(zhǔn)值，2017、2018、2019 年EUIH 基準(zhǔn)值分別為0.194、0.192、0.193 GJ/m2.若根據(jù)城市管理不同實(shí)際需要，將不同可信區(qū)間上下限EUIH 數(shù)值設(shè)定作為能耗準(zhǔn)入值或限定值、先進(jìn)值或引導(dǎo)值[47]，如95%可信區(qū)間上下限、上下四分位數(shù)等.以2018 年為例，95%可信區(qū)間上下限、上下四分位數(shù)由高至低分別為0.163、0.180、0.201、0.233 GJ/m2，即2018 年EUIH 引導(dǎo)值為0.163 GJ/m2（反映優(yōu)秀、更低的能耗強(qiáng)度），約束值為0.180 GJ/m2（反映良好、較低的能耗強(qiáng)度），基線值為0.192 GJ/m2（反映一般、平均的能耗強(qiáng)度）.

圖8 北京市居住建筑EUIH 的聯(lián)合后驗(yàn)分布Fig.8 Joint posterior distribution for EUIH of residential buildings in Beijing

同一年度、同一城市的居住建筑區(qū)域樣本的EUIH 值與上述基準(zhǔn)值進(jìn)行對(duì)比，差值反映出該建筑區(qū)域的能耗水平與所在城市正常能耗水平的差距.該基準(zhǔn)值與同一年度、其他國(guó)家或城市的居住建筑EUIH 相關(guān)基準(zhǔn)進(jìn)行對(duì)比，差值反映出北京市居住建筑供熱節(jié)能潛力.圖9 為北京市與其冬季環(huán)境溫度相近的國(guó)際能源署成員國(guó)（包括美國(guó)、韓國(guó)、日本、法國(guó)、德國(guó)、加拿大、丹麥、芬蘭、瑞典、烏克蘭、比利時(shí)、匈牙利等國(guó)家）城市居住建筑EUIH 基準(zhǔn)值對(duì)比，從2010—2018 年國(guó)際能源署成員國(guó)城市居住建筑EUIH 基準(zhǔn)值呈明顯下降趨勢(shì)，說(shuō)明供熱節(jié)能水平和節(jié)能效果提升；與北京市2018 年居住建筑EUIH 基準(zhǔn)值相比，北京市基準(zhǔn)值明顯低于絕大多數(shù)國(guó)際能源署成員國(guó)，說(shuō)明北京市居住建筑供熱節(jié)能潛力較小.

圖9 國(guó)際能源署成員國(guó)與北京的居住建筑EUIH 基準(zhǔn)值對(duì)比Fig.9 Comparison of benchmark of EUIH for residential building between city of IEA number countries and Beijing

5 結(jié) 論

本研究提出一種基于自下而上數(shù)據(jù)驅(qū)動(dòng)模型的建筑能耗預(yù)測(cè)和校準(zhǔn)方法，并在建筑樣本數(shù)據(jù)有限、輸入不確定性強(qiáng)、建模復(fù)雜度和計(jì)算成本限制下，推斷城市尺度建筑能耗高概率分布點(diǎn)、區(qū)間及聯(lián)合后驗(yàn)分布，以確定城市尺度建筑供熱能耗強(qiáng)度基準(zhǔn).

1）系統(tǒng)綜述對(duì)比目前建筑能耗預(yù)測(cè)模型構(gòu)建、求解和應(yīng)用中所需建模信息量、建模難度、樣本量、計(jì)算資源和時(shí)間成本，指出城市尺度數(shù)據(jù)驅(qū)動(dòng)建筑能耗預(yù)測(cè)領(lǐng)域現(xiàn)有的研究差距和適合于解決建筑樣本有限、輸入不確定等問(wèn)題的模型形式.

2）開發(fā)貝葉斯框架的混合效應(yīng)模型，采用城市氣候、形態(tài)、建筑功能等特征指標(biāo)為模型輸入，采集2017—2019 年北京市223 座居住建筑熱力站能耗監(jiān)測(cè)數(shù)據(jù)為測(cè)試集樣本訓(xùn)練模型.通過(guò)MCMC方法進(jìn)行重要性采樣，并推斷出模型超參數(shù)峰值概率分布，預(yù)測(cè)城市尺度建筑供熱能耗強(qiáng)度（EUIH）的精度指標(biāo)NMBE 為0.400%，RMSE 為0.034，CVRMSE 為13.100%，驗(yàn)證所提模型具有較強(qiáng)魯棒性和泛化能力.

3）基于貝葉斯框架的混合效應(yīng)模型和MCMC方法平穩(wěn)分布后迭代樣本，推斷城市尺度建筑供熱能耗強(qiáng)度（EUIH）的聯(lián)合后驗(yàn)分布，EUIH 點(diǎn)估計(jì)和區(qū)間估計(jì)結(jié)果為確定城市建筑供熱能耗基準(zhǔn)提供依據(jù).

4）北京市2017、2018、2019 年EUIH 基準(zhǔn)值分別為0.194、0.192、0.193 GJ/m2.與其冬季環(huán)境溫度相近的國(guó)際能源署成員國(guó)EUIH 基準(zhǔn)值對(duì)比，北京市基準(zhǔn)值明顯低于絕大多數(shù)國(guó)際能源署成員國(guó)，說(shuō)明了北京市居住建筑供熱節(jié)能潛力較小.

本研究中的貝葉斯框架混合效應(yīng)模型描述了輸入不確定性對(duì)模型輸出的影響，模型輸入、輸出、模型中的超參數(shù)等均以貝葉斯學(xué)派概率方式而不是固定值的方式進(jìn)行描述，允許在一定值域范圍內(nèi)變化，從而數(shù)值上包含了模型輸入在觀測(cè)誤差影響、設(shè)備使用狀態(tài)、人員活動(dòng)行為變化、建筑建造水平和功能差異等不確定因素對(duì)模型輸出影響.在所提模型中，有以下幾個(gè)方面存在局限性：共軛先驗(yàn)函數(shù)的形式和初值設(shè)定、似然分布函數(shù)的形式和初值設(shè)定、貝葉斯框架優(yōu)化、模型輸入優(yōu)化、實(shí)際數(shù)據(jù)顆粒度差異優(yōu)化、重要性采樣的方式效度等.在建模過(guò)程中，若是將建筑人員活動(dòng)等不確定因素的數(shù)據(jù)作為相關(guān)解釋變量，則需要使用全貝葉斯、經(jīng)驗(yàn)貝葉斯模型和推理方法來(lái)改善模型結(jié)構(gòu)和結(jié)果.