基于NPP/VIIRS夜間燈光圖像的GDP智能預(yù)測模型研究

王森， 蔡小莉， 鮑云飛， 詹邦成

(1.武漢大學(xué)信息管理學(xué)院， 武漢 430072； 2.河南大學(xué)河南省大數(shù)據(jù)分析與處理重點(diǎn)實(shí)驗(yàn)室， 開封 475004； 3.北京空間機(jī)電研究所， 北京 100094)

近十年來，由于數(shù)據(jù)共享政策和傳感器的快速發(fā)展，夜光遙感成為了遙感領(lǐng)域的熱門分支之一。夜光遙感已經(jīng)被廣泛應(yīng)用于人文地理、區(qū)域經(jīng)濟(jì)、地緣政治等研究領(lǐng)域，甚至金融行業(yè)也開始利用夜光遙感數(shù)據(jù)。由于數(shù)據(jù)豐富和成本低廉，衛(wèi)星遙感是夜光遙感的主流手段。遙感方法是從空間角度對社會經(jīng)濟(jì)發(fā)展?fàn)顩r進(jìn)行測繪和分析的最有效和最廣泛使用的手段之一。遙感影像以其覆蓋面廣、效率高、客觀性強(qiáng)等優(yōu)點(diǎn)，為人類提供了海量的空間信息，在建成區(qū)提取、土地?cái)U(kuò)展探測、土地利用制圖、土地覆蓋變化分析、城市景觀結(jié)構(gòu)檢測、城市空間結(jié)構(gòu)分析方面有豐富的應(yīng)用。在這些研究方向中，利用直接探測到的長時間序列夜間燈光進(jìn)行國內(nèi)生產(chǎn)總值(gross domestic product，GDP)的估算是最具代表性的衡量人類社會活動的發(fā)展方向。

夜間燈光(nighttime light，NTL)遙感影像通過傳感器記錄地球表面的夜間光，可以探測來自城市、城鎮(zhèn)和工業(yè)區(qū)的燈光，在人類社會經(jīng)濟(jì)活動檢測方面顯示出良好的性能。GDP可以有效地反映一個國家的國力和財(cái)力[1]，準(zhǔn)確的像素級GDP數(shù)據(jù)是了解區(qū)域經(jīng)濟(jì)動態(tài)的基礎(chǔ)，是定量分析不同行政區(qū)劃、不同自然屬性和人文屬性地區(qū)經(jīng)濟(jì)差異的基礎(chǔ)。不同地區(qū)的發(fā)展模式呈現(xiàn)出更加突出的地區(qū)差異。當(dāng)今世界上大多數(shù)國家或地區(qū)根據(jù)當(dāng)?shù)匦枨螵?dú)立制定GDP統(tǒng)計(jì)口徑，導(dǎo)致無法直接比較國家/地區(qū)之間的統(tǒng)計(jì)數(shù)據(jù)，這會影響使用效率，很難有效客觀準(zhǔn)確地反映GDP的時空演化規(guī)律和精細(xì)化的區(qū)域發(fā)展研究[2]。當(dāng)使用經(jīng)濟(jì)統(tǒng)計(jì)方法測算GDP時，會存在以下局限性：①統(tǒng)計(jì)尺度不一致、統(tǒng)計(jì)單位內(nèi)部數(shù)據(jù)不統(tǒng)一等問題，難以在精細(xì)尺度上反映區(qū)域經(jīng)濟(jì)發(fā)展的差異；②在一些偏遠(yuǎn)地區(qū)，其經(jīng)濟(jì)資料難以獲取，且數(shù)據(jù)的可靠性不高；③對于少數(shù)地方、少數(shù)企業(yè)、少數(shù)單位在GDP統(tǒng)計(jì)時存在弄虛作假和統(tǒng)計(jì)造假行為、統(tǒng)計(jì)違法違規(guī)的現(xiàn)象，也成“擠水分”等現(xiàn)象。傳統(tǒng)的方法在預(yù)測GDP時，提取模型的建立方法時大多為簡單的線性回歸模型、冪函數(shù)回歸模型、二次函數(shù)模型、多項(xiàng)式法等模型，模型在構(gòu)建時太傳統(tǒng)，因此難以給出準(zhǔn)確地預(yù)測。準(zhǔn)確地計(jì)量GDP，對于了解區(qū)域經(jīng)濟(jì)的動態(tài)變化，滿足跨學(xué)科研究的需要具有重要意義。

夜間燈光數(shù)據(jù)可以彌補(bǔ)上述所說的不足，與傳統(tǒng)的社會經(jīng)濟(jì)普查相比，遙感影像在描述GDP預(yù)測方面具有明顯的優(yōu)勢。首先，夜間燈光數(shù)據(jù)不受區(qū)域間價格因素的困擾。目前，中國僅在地級城市之間計(jì)算各種價格指數(shù)，如果把每個縣的GDP數(shù)據(jù)進(jìn)行之間比較，就無法消除在地區(qū)之間產(chǎn)生的價格因素。其次，夜間燈光數(shù)據(jù)包括的市場經(jīng)濟(jì)的商品和服務(wù)不僅是以GDP衡量的，而且還包括市場不提供的商品和服務(wù)的價值，因此它對于反映一個國家或地區(qū)的經(jīng)濟(jì)發(fā)展是很真實(shí)的。夜間燈光數(shù)據(jù)作為一種全球無差別的夜間觀測量，以其獨(dú)特的方式描述地球經(jīng)濟(jì)活動強(qiáng)度。它可以覆蓋整個世界，并且不由于行政區(qū)域的差異而產(chǎn)生差異，因此，夜間燈光數(shù)據(jù)可以用作獨(dú)特的時空輸入，在預(yù)測GDP方面得到了廣泛的應(yīng)用[3]。Nataliya等[4]利用DMSP/OLS(defense meteorological satellite program/operational linescan system)數(shù)據(jù)研究證實(shí)了利用夜間光衛(wèi)星測量數(shù)據(jù)重建經(jīng)濟(jì)活動地理格局的可行性，且夜間光強(qiáng)度往往隨著人口密度和人均GDP的增加而增加。范強(qiáng)等[5]通過構(gòu)建多項(xiàng)式預(yù)測模型，分析了比較長時間序列燈光集和GDP之間的關(guān)系。Zhou等[6]定量研究了NTL與GDP和人口之間的相關(guān)性，研究結(jié)果表明夜間燈光與GDP和人口具有顯著的正相關(guān)。并通過一系列不平等系數(shù)，得出中國中西部存在明顯的區(qū)域差異。研究使用的夜光遙感數(shù)據(jù)是柵格數(shù)據(jù)，每個圖像中蘊(yùn)藏著大量的數(shù)據(jù)信息。從夜間燈光數(shù)據(jù)中可以看到光亮度的變化和分布，從而提供更多的時間和空間信息。

目前利用Suomi國家極軌伙伴關(guān)系(national polar-orbiting partnership，NPP)衛(wèi)星攜帶的可見紅外成像輻射計(jì)套件(visible infrared imaging radiometer，VIIRS)的數(shù)據(jù)進(jìn)行估算經(jīng)濟(jì)參量的研究相對較多，但很少有研究使用VIIRS數(shù)據(jù)利用機(jī)器學(xué)習(xí)算法進(jìn)行GDP預(yù)測。為此，利用長時序夜間數(shù)據(jù)建立長三角地區(qū)GDP智能預(yù)測模型，以此對該地區(qū)的GDP進(jìn)行及時有效地預(yù)測，并進(jìn)一步分析不同地區(qū)的經(jīng)濟(jì)差異。研究結(jié)果為夜光數(shù)據(jù)的應(yīng)用領(lǐng)域和應(yīng)用方式，以及地區(qū)的可持續(xù)規(guī)劃和發(fā)展提供了決策依據(jù)。

1 基礎(chǔ)理論

隨機(jī)森林(random forest，RF)模型是一種Bagging集成算法，具有靈活度高、不易過擬合、準(zhǔn)確率高的特性，有廣闊的應(yīng)用場景[7]。隨機(jī)森林的有關(guān)理論可參考李倩倩等[8]的研究。隨機(jī)森林可以用來解決回歸問題和分類問題，當(dāng)它被用來做回歸問題時，原理是通過對樣本進(jìn)行隨機(jī)抽取并放回，生成對個決策樹，并把決策樹的結(jié)果進(jìn)行平均作為模型輸出[9]。因此，隨機(jī)森林能夠?qū)Ω呔S特征數(shù)據(jù)進(jìn)行評估及處理的能力[10]。

隨機(jī)森林有兩個重要的參數(shù)：樹的數(shù)目和分割的特征數(shù)。與多元線性回歸、地理擴(kuò)展回歸和人工神經(jīng)網(wǎng)絡(luò)方法相比，隨機(jī)森林算法不需要考慮多重共線性問題，可避免過度擬合問題，運(yùn)算速度高且適合處理高維特征數(shù)據(jù)等優(yōu)點(diǎn)[11-13]，因此在人口模擬，生態(tài)學(xué)，經(jīng)濟(jì)學(xué)，醫(yī)學(xué)等領(lǐng)域均運(yùn)用較廣泛。對于隨機(jī)森林來說其算法的核心步驟如下。

步驟1使用Bootstrap方法選擇可替代的樣本，依次選擇m個集合作為樣本。共生成ntree棵樹的訓(xùn)練集。

步驟2對ntree棵樹進(jìn)行訓(xùn)練并產(chǎn)生相應(yīng)的模型。

步驟3根據(jù)分區(qū)索引為每個拆分選擇最佳特征。

步驟4每個決策樹都是獨(dú)立拆分的，并且在決策樹的拆分過程中不需要修剪，直到該例程的所有訓(xùn)練示例都屬于同一類別為止。生成的多個決策樹形成一個隨機(jī)森林，并且將每個決策樹的分類或預(yù)測結(jié)果進(jìn)行整合。

算法流程圖如圖1所示。

C1為學(xué)習(xí)集S1產(chǎn)生的模型;C2為學(xué)習(xí)集S2產(chǎn)生的模型;C3為學(xué)習(xí)集S3產(chǎn)生的模型;X為模型的最佳特征圖1 隨機(jī)森林算法流程圖Fig.1 Algorithm flowchart of random forest

2 研究區(qū)域及數(shù)據(jù)集處理

2.1 研究區(qū)域

中國的長江三角洲(117°29′E～123°01′E，27°64′N～34°33N)地區(qū)是指上海，浙江，江蘇和安徽省(三省一市)。地勢主要包括安徽蘇濱河平原、長江三角洲平原和里下河平原。根據(jù)2019年長三角區(qū)域一體化發(fā)展規(guī)劃，長三角城市群包括214個縣。它的面積為22.5×104km2，總?cè)丝诩s為2.43億，研究區(qū)域如圖2所示。

圖2 研究區(qū)域Fig.2 Study area

2.2 NPP/VIIRS夜間燈光數(shù)據(jù)預(yù)處理

VIIRS原始數(shù)據(jù)采用WGS1984坐標(biāo)系，其夜光影像獲取時間為凌晨1:30，VIIRS空間分辨率為15″[14]，這會造成影像隨著緯度的增大而減小的問題。為了避免影像由于坐標(biāo)系發(fā)生的變形影響及便于方便計(jì)算影像面積，將所有的NPP/VIIRS影像數(shù)據(jù)都轉(zhuǎn)換為Albers等積投影坐標(biāo)系。最后將影像重采樣為0.5 km×0.5 km大小的網(wǎng)格。研究中使用的第一版VIIRS月合成數(shù)據(jù)并沒有處理如極光、火光和其他短暫性光源的影響，因此存在背景噪聲，燈光輻射率值存在負(fù)值、極小值、極高值的特點(diǎn)，這些噪聲和異常值數(shù)據(jù)的存在會影響數(shù)據(jù)模擬精度。針對這種情況，根據(jù)Shi等[15]的方法做了以下處理：把小于零的值設(shè)為0.001(近似為0，對統(tǒng)計(jì)分析沒影響)，把背景值設(shè)為0；將大于235的數(shù)值設(shè)為235，以刪除一些異常值。由于原始獲取的 NPP/VIIRS 夜間燈光數(shù)據(jù)產(chǎn)品時間尺度為月度產(chǎn)品，進(jìn)行空間化時需要使用合成后的全年數(shù)據(jù)，采用平均值法對其進(jìn)行合成獲取得到2012—2018年長三角區(qū)域NPP/VIIRS年度數(shù)據(jù)。研究區(qū)VIIRS影像如圖3所示。

2.3 其他輔助數(shù)據(jù)

除NPP/VIIRS夜間燈光數(shù)據(jù)之外，研究區(qū)還采用建成區(qū)面積、戶籍人口、全社會用電量等其他輔助數(shù)據(jù)。其中經(jīng)濟(jì)統(tǒng)計(jì)數(shù)據(jù)從《中國城市統(tǒng)計(jì)年鑒》中查找并收集了長三角地區(qū)214個縣的2012—2018年GDP數(shù)據(jù)(單位：億元)。表1列出了所使用的數(shù)據(jù)。

DN(digital number)表示像元，DN值越大，亮度越亮圖3 處理后的研究區(qū)夜光數(shù)據(jù)Fig.3 Processed nighttime data of the study area

表1 其他輔助數(shù)據(jù)Table 1 Other auxiliary data

3 貝葉斯優(yōu)化改進(jìn)的隨機(jī)森林模型

3.1 貝葉斯優(yōu)化技術(shù)框架結(jié)構(gòu)

貝葉斯優(yōu)化算法是所有優(yōu)化算法中較高效的一種優(yōu)化算法，許多研究表明貝葉斯優(yōu)化可以替代解決復(fù)雜問題中的其他優(yōu)化算法。它使用先前搜索到的點(diǎn)的信息來確定下一個搜索點(diǎn)。其原理是在不知道目標(biāo)函數(shù)時，充分利用已存在的采樣點(diǎn)來預(yù)測函數(shù)的最大值[16]。

在數(shù)學(xué)上，可以統(tǒng)一替換問題描述為替代目標(biāo)功能的替代最優(yōu)解，可表示為[17]

(1)

式(1)中：x為要進(jìn)行優(yōu)化的參數(shù)；X∈R；x*為要進(jìn)行優(yōu)化的參數(shù)集合；f(x)為要預(yù)測最大值的目標(biāo)函數(shù)。

算法的思路是首先生成一個初始候選解集合，然后根據(jù)這些點(diǎn)尋找下一個有可能是極值的點(diǎn)，將該點(diǎn)加入集合中，重復(fù)這一步驟，直至迭代終止。最后從這些點(diǎn)中找出極值點(diǎn)作為問題的解。

這里的關(guān)鍵問題是如何根據(jù)已經(jīng)搜索的點(diǎn)確定下一個搜索點(diǎn)。貝葉斯優(yōu)化根據(jù)已經(jīng)搜索的點(diǎn)的函數(shù)值估計(jì)真實(shí)目標(biāo)函數(shù)值的均值和方差(即波動范圍)，如圖4所示。由圖4(a)可以看出，兩條虛線所夾區(qū)域?yàn)樵诿恳稽c(diǎn)處函數(shù)值的變動范圍，在以均值即紅色曲線為中心，與標(biāo)準(zhǔn)差成正比的區(qū)間內(nèi)波動。在搜索點(diǎn)處，紅色曲線經(jīng)過搜索點(diǎn)，且方差最小，在遠(yuǎn)離搜索點(diǎn)處方差更大，這也符合的直觀認(rèn)識，遠(yuǎn)離采樣點(diǎn)處的函數(shù)值估計(jì)的更不可靠。根據(jù)均值和方差可以構(gòu)造出采集函數(shù)，即對每一點(diǎn)是函數(shù)極值點(diǎn)的可能性的估計(jì)，反映了每一個點(diǎn)值得搜索的程度，該函數(shù)的極值點(diǎn)是下一個搜索點(diǎn)，如圖4(b)所示。

算法的核心由兩部分構(gòu)成：對目標(biāo)函數(shù)進(jìn)行建模即計(jì)算每一點(diǎn)處的函數(shù)值的均值和方差，通常用高斯過程回歸實(shí)現(xiàn)；構(gòu)造采集函數(shù)0，用于決定本次迭代時的采樣位置。

紅色曲線為估計(jì)出的目標(biāo)函數(shù)值即在每一點(diǎn)出處的目標(biāo)函數(shù)值的均值；兩條虛線所夾區(qū)域?yàn)樵诿恳稽c(diǎn)處函數(shù)值的變動范圍圖4 貝葉斯優(yōu)化的原理Fig.4 Bayesian optimization principle

3.1.1 高斯過程

高斯過程(Gaussian process，GP)是把多個維度的高斯在無限維上進(jìn)行擴(kuò)展，它是隨機(jī)性的，高斯過程通過均值函數(shù)和協(xié)方差來進(jìn)行定義的。

m(x)=E[f(x)]

(2)

k(x,x′)=E{f(x)-m(x)[f(x′)-m(x′)]}

(3)

式中：m(x)為均值函數(shù)；E為期望，k(x,x′)為協(xié)方差函數(shù)，當(dāng)x=x′時，k(x,x′)=1，x和x′相差越大，k越趨向于0。

則高斯過程可表示為

f(x)～GP[m(x),k(x,x′)]

(4)

式(4)中：GP為高斯過程(Gaussian process)。

在計(jì)算時為了方便，通常將均值初始函數(shù)值設(shè)為0，當(dāng)知道{(xi,fi)|i=1,2,…,n}時，則必然有一個高斯分布存在，其中(xi,fi)為第i個函數(shù)值，且滿足：

f～N[0,K(X,X)]

(5)

式(5)中：N為聯(lián)合分布；高斯過程的協(xié)方差矩陣K(X,X)和協(xié)方差函數(shù)k(x,x′)可分別表示為

(6)

(7)

因此，加入新樣本X*，并利用協(xié)方差矩陣生成一個新的高斯分布，可表示為

f*～N[0,K(X*,X*)]

(8)

由高斯過程的性質(zhì)可得，訓(xùn)練輸出f和測試輸出f*的聯(lián)合分布為

(9)

則f*的聯(lián)合后驗(yàn)分布滿足：

f*|X*,X,f～N[K(X*,X)K(X,X)-1f,

K(X*,X*)-K(X*,X)K(X,X)-1K(X,X*)]

(10)

3.1.2 采樣函數(shù)

通過采樣函數(shù)可以確定下一個要評估的點(diǎn)，減少重復(fù)次數(shù)，大大降低優(yōu)化成本。對采樣點(diǎn)函數(shù)進(jìn)行選擇時，有兩方面需要考慮，一個是利用一個是探索。利用就是利用當(dāng)前的最優(yōu)解，在最優(yōu)解附近尋找下一步的最優(yōu)解，探索是對未評估過的點(diǎn)進(jìn)行探索，避免只在局部探索而造成的不全面的搜索。

對于采樣函數(shù)，常用的有：期望的改進(jìn)函數(shù)EI、改進(jìn)的概率函數(shù)PI以及上限置信區(qū)間UCB。本實(shí)驗(yàn)選擇UCB作為采用函數(shù)，函數(shù)公式為

UCB(x)=μ(x)+εδ(x)

(11)

式(11)中：μ(x)和δ(x)分別為利用高斯過程并結(jié)合后驗(yàn)分布函數(shù)得到的均值和協(xié)方差；利用式(11)可以對參數(shù)ε進(jìn)行修剪和調(diào)節(jié)。

3.2 實(shí)驗(yàn)驗(yàn)證

將實(shí)驗(yàn)數(shù)據(jù)集中的2012—2017年長時間序列數(shù)據(jù)為訓(xùn)練集，2018年數(shù)據(jù)為測試集，將訓(xùn)練集進(jìn)行訓(xùn)練。并且通過交叉驗(yàn)證確定參數(shù)以獲得最佳的模型效果。

為了提高預(yù)測精度，利用貝葉斯優(yōu)化對參數(shù)進(jìn)行優(yōu)化，貝葉斯優(yōu)化作為調(diào)參利器的一種，它是一種黑盒優(yōu)化算法，用于求解表達(dá)式未知的函數(shù)的極值問題。算法根據(jù)一組采樣點(diǎn)處的函數(shù)值預(yù)測出任意點(diǎn)處函數(shù)值的概率分布，這通過高斯過程回歸而實(shí)現(xiàn)。根據(jù)高斯過程回歸的結(jié)果構(gòu)造采集函數(shù)，用于衡量每一個點(diǎn)值得探索的程度，求解采集函數(shù)的極值從而確定下一個采樣點(diǎn)。最后返回這組采樣點(diǎn)的極值作為函數(shù)的極值,優(yōu)化流程如圖5所示。

超參數(shù)是進(jìn)行模型訓(xùn)練時用于控制算法行為的參數(shù)。由于這些參數(shù)并不能從常規(guī)訓(xùn)練中得到。因此在進(jìn)行訓(xùn)練模型前，需要為它們分配值。

對于隨機(jī)森林，當(dāng)“子模型數(shù)”的值增加時，它對減少模型的總體方差具有顯著作用，但對子模型的偏差和方差沒有影響。但是，由于改變是為了減少模型方差公式的第二項(xiàng)，因此提高精度存在上限。在不同條件下，“分裂準(zhǔn)則”對模型的準(zhǔn)確性有不同的影響，并且在實(shí)際使用中需要靈活調(diào)整此參數(shù)。如果希望模型具有較低的偏差和較高的方差，則需要制作更多的葉節(jié)點(diǎn)或樹的深度變的更深；同時，通過調(diào)整“min_samples_leaf”和“min_weight fraction_leaf”，這意味著子模型結(jié)構(gòu)較為復(fù)雜，一般情況下，使用bootstrap進(jìn)行子采集可以減少樣本之間的相關(guān)性。在分割過程中，適當(dāng)減少考慮的最大特征數(shù)量“max_features”，或者將其他的隨機(jī)性考慮到子模型中也能使子模型之間相關(guān)程度降低。但是，如果盲目地對此參數(shù)進(jìn)行減少是不行的，因?yàn)樵诓鸱诌^程中的可選功能比較少，因此會導(dǎo)致模型的偏差逐漸增大。圖6中，可以看到上述參數(shù)對隨機(jī)森林模型的整體性能的影響。

本實(shí)驗(yàn)主要對max_depth、min_samples_leaf、n_estimators、min_samples_split這4個主要參數(shù)進(jìn)行貝葉斯優(yōu)化以提高預(yù)測精度。其中，將max_depth的值設(shè)置為None。貝葉斯優(yōu)化過程主要包括先驗(yàn)函數(shù)和采集函數(shù)，PF是利用高斯過程來實(shí)現(xiàn)的，采集函數(shù)包括EI、PI、UCB[18-19]

構(gòu)建貝葉斯優(yōu)化改進(jìn)的隨機(jī)森林模型使用sklearn庫以及scikit-optimization的BayesSearchCV，模型的超參數(shù)范圍設(shè)max_depth[5,15],n_estimators[10,250],min_samples_split[2,25],min_sam-ples_leaf[1,10]。使用sklearn 中的貝葉斯搜索函數(shù)對上述函數(shù)設(shè)置的區(qū)間進(jìn)行搜索，同樣以RMSE作為精度評價指標(biāo)，進(jìn)行最優(yōu)參數(shù)組合。偽代碼如下。

輸入:f,x,s,m

輸出:一組最近的超參數(shù)

1.begin

2.D←initsamples(f,x)

3.fori←|D|toTdo

4.p(y|x,D)←fitmodel(m,D)

6.yi←f(xi)

7.D←D∪(xi,yi)

8.end for

其中，f為所謂的黑盒子，即輸入一組超參數(shù)，得到一個輸出值;T為目標(biāo)函數(shù)；x為超參數(shù)搜索空間；s為采集函數(shù)；m為對數(shù)據(jù)集D進(jìn)行擬合得到的模型；D←initsamples(f,x)這一步驟就是初始化獲取數(shù)據(jù)集D=(xi,yi)，(xi+1,yi+1)，…，(xn,yn)，其中yi=f(xi)，這些都是已知的；每次選出參數(shù)x后都需要計(jì)算f(x)，因此一般需要固定選參次數(shù)(或者是函數(shù)評估次數(shù))，p(y|x,D)←fitmodel(M,D)，預(yù)先假設(shè)M服從高斯分布，且已知了數(shù)據(jù)集D，因此可通過模型函數(shù)表示；xi為一組超參數(shù)；yi←f(xi)表示由xi去得到超參數(shù)yi；D←D∪(xi,yi)表示更新數(shù)據(jù)集。

貝葉斯參數(shù)優(yōu)化過程如表2所示。

a表示對模型性能的影響；p表示參數(shù)的遞增；o為原點(diǎn)坐標(biāo)圖6 隨機(jī)森林參數(shù)對整體模型性能的影響Fig.6 The influence of random forest parameters on overall model performance

可以看出，僅迭代了30次，貝葉斯優(yōu)化就得到了最優(yōu)參數(shù)組合。

從圖7可以看到，在前面15次迭代時，R2變化較大，但是迭代15次后，R2就很穩(wěn)定，這說明貝葉斯優(yōu)化已經(jīng)找到相對較好的參數(shù)組合。最后得到最優(yōu)的參數(shù)組合如表3所示。

優(yōu)化結(jié)果顯示當(dāng)4個超參數(shù)的值分別取以上值時可以取得最小均方根誤差(root mean squared error, RMSE)，將參數(shù) max_depth、min_samples_leaf、min_samples_split、n_estimators取14.43、1.289、2.828、241.3時預(yù)測精度最高，編譯函數(shù)采用fit(·)。

表2 貝葉斯參數(shù)優(yōu)化Table 2 Bayesian parameter optimization

4 網(wǎng)格搜索改進(jìn)的隨機(jī)森林模型

4.1 網(wǎng)格搜索技術(shù)框架結(jié)構(gòu)

網(wǎng)格搜索即通過對所有的參數(shù)進(jìn)行循環(huán)遍歷，最終選擇使得機(jī)器學(xué)習(xí)模型表現(xiàn)最好的參數(shù)或參數(shù)組合，也稱為“窮舉法”，無確定的函數(shù)公式與之相照應(yīng)。它的關(guān)鍵是調(diào)優(yōu)搜索。通過構(gòu)建參數(shù)候選集合，對所有候選參數(shù)進(jìn)行選擇時，通過循環(huán)嘗試各種可能性，具有最佳性能的參數(shù)是最終結(jié)果，它類似于手動調(diào)整。對于連續(xù)型的超參數(shù)，對其可行域進(jìn)行網(wǎng)格劃分，選取一些典型值進(jìn)行計(jì)算。假設(shè)需要確定的超參數(shù)有2個，第1個的取值為[0,1]之間的實(shí)數(shù)，第2個的取值為[1,2]之間的實(shí)數(shù)。則可以按照如下的方案得到若干離散的取值，以這些值運(yùn)行算法：將第1個參數(shù)均勻的取3個典型值，將第2個參數(shù)均勻的取3個典型值。對于所有的取值組合運(yùn)行算法，將性能最優(yōu)的取值作為超參數(shù)的最終取值，這種方法如圖8所示。

網(wǎng)格搜索隨著參數(shù)數(shù)量的增加呈指數(shù)級增長，因此對于超參數(shù)較多的情況，該方法面臨性能上的問題。著名的支持向量機(jī)(support vector machine，SVM)開源庫libsvm使用了網(wǎng)格搜索算法確定SVM的超參數(shù)。

圖7 迭代結(jié)果曲線Fig.7 Curve of iteration results

表3 最優(yōu)參數(shù)組合Table 3 Optimal parameter combination

圖8 網(wǎng)格搜索的原理Fig.8 Principle of grid search

網(wǎng)格搜索的基本原理是，首先設(shè)置參數(shù)的值范圍，然后通過特定步長將參數(shù)的可行區(qū)間劃分為網(wǎng)格，然后通過搜索每個網(wǎng)格將參數(shù)的值控制在一定范圍內(nèi)。并根據(jù)特定規(guī)則搜索所有交點(diǎn)，并以此方式找到每個最佳目標(biāo)函數(shù)的交點(diǎn)。此時，最佳參數(shù)是交點(diǎn)處的參數(shù)[20-21]。

當(dāng)徑向基函數(shù)(rodial bosis function, RBF)的核函數(shù)選擇后，懲罰因子C和核函數(shù)參數(shù)σ也要確定，以提高故障識別率。在網(wǎng)格搜索方法中，懲罰因子C的間隔為C∈[C1,C2](其中，C1為懲罰因子最小值，C2為懲罰因子最大值)，其改變步長為CS，核函數(shù)參數(shù)的間隔為σ∈[σ1,σ2](其中，σ1為核函數(shù)參數(shù)最小值，σ2為核函數(shù)參數(shù)最大值)，其改變步長為σS。網(wǎng)格搜索方法對網(wǎng)格區(qū)間中的每對參數(shù)(C′,σ′)執(zhí)行模型訓(xùn)練和故障識別率計(jì)算，最后把識別率最高的參數(shù)作為模型的參數(shù)。網(wǎng)格搜索實(shí)際上是一個n層循環(huán)(n為參數(shù)的數(shù)量)，其計(jì)算步驟如下。

步驟1初始化精度為0和參數(shù)C,σ的初始值C1、σ1。

步驟2固定步長C增長CS,C=C+CS。

步驟3固定步長σ增長σS，σ=σ+σS。

步驟4引入?yún)?shù)(C,σ)建立支持后量機(jī)(support vector machine,SVM)模型，使用測試樣本計(jì)算預(yù)測準(zhǔn)確率，與上一步準(zhǔn)確率的進(jìn)行比較，如果準(zhǔn)確率得到提高，則將該參數(shù)C、σ替換上一步C、σ，反之亦然。

步驟5判斷σ是否增加到σ2,如果是，則跳至步驟6；如果不是，則執(zhí)行步驟3。

步驟6判斷C是否增長到C2,如果是則執(zhí)行第7步,如果不是，則執(zhí)行第2步。

步驟7當(dāng)前C和σ作為為最優(yōu)參數(shù),計(jì)算結(jié)束。

通過上述循環(huán)找到RBF函數(shù)的最佳參數(shù)之和就足夠了。

4.2 實(shí)驗(yàn)驗(yàn)證

為了對預(yù)測算法的精度進(jìn)行提高，需要對參數(shù)進(jìn)行優(yōu)化，提出一種基于袋外數(shù)據(jù)估計(jì)的回歸誤差，采用改進(jìn)的網(wǎng)格搜索算法對隨機(jī)森林(random forest,RF)模型進(jìn)行參數(shù)調(diào)優(yōu)，流程圖如圖9所示。設(shè)隨機(jī)森林樹的數(shù)量為n_estimators, 范圍設(shè)置為[50,100,150,200],樹的高度為max_depth, 范圍設(shè)置為[3,4,5,6],拆分內(nèi)部節(jié)點(diǎn)所需的最小樣本數(shù)為min_samples_split,范圍設(shè)置為[2,3,4,5]。葉子節(jié)點(diǎn)所需的最小樣本數(shù)為min_samples_leaf,范圍設(shè)置為[1,2,3,4]。采用網(wǎng)格搜索算法根據(jù)步長對參數(shù)進(jìn)行依次調(diào)整，對參數(shù)范圍依次遍歷找到使精度最高的參數(shù)，并利用5折交叉驗(yàn)證對其進(jìn)行評估。5折交叉驗(yàn)證是通過將數(shù)據(jù)集平均分為5份，隨機(jī)選擇4份作訓(xùn)練集，1份作測試集，一共進(jìn)行5次訓(xùn)練和測試，采用精度等評分方法計(jì)算平均值，然后找到評分最大的參數(shù)組合，得到預(yù)測精度最高的一組最優(yōu)參數(shù)。

考慮到隨機(jī)森林算法中決策樹數(shù)k′和候選分裂屬性數(shù)mtry為離散值，提出的基于網(wǎng)格搜索改進(jìn)的隨機(jī)森林模型進(jìn)行參數(shù)優(yōu)化時目標(biāo)函數(shù)值選用袋外數(shù)據(jù)估計(jì)誤差。隨機(jī)森林模型在構(gòu)建過程時存在隨機(jī)性，回歸誤差也會在一定范圍內(nèi)存在波動，因此為減小不確定性對參數(shù)選擇的影響，在計(jì)算回歸誤差時選取多個結(jié)果的平均值作為隨機(jī)森林模型回歸誤差。當(dāng)搜索完畢時輸出最優(yōu)參數(shù)，其搜索過程用流程圖(圖10)來表示。

構(gòu)建網(wǎng)格優(yōu)化改進(jìn)的隨機(jī)森林模型采用sklearn庫中的sklearn.ensemble模塊，從模塊中調(diào)用RandomForestRegressor函數(shù)對參數(shù)進(jìn)行網(wǎng)格優(yōu)化，模型的超參數(shù)及范圍設(shè)置如表4所示。在本次實(shí)驗(yàn)中，它總共構(gòu)建了256個不同的模型。

圖9 優(yōu)化的隨機(jī)森林算法整體流程Fig.9 Overall flow of the optimized random forest algorithm

圖10 基于改進(jìn)的網(wǎng)格搜索算法的隨機(jī)森林參數(shù)尋優(yōu)流程圖Fig.10 Flow chart of random forest parameter optimization based on improved grid search algorithm

在利用網(wǎng)格搜索算法對隨機(jī)森林模型進(jìn)行改進(jìn)時采用的函數(shù)是sklearn中的grid search函數(shù)[sklearn.Model_selection.GridSearchCV()]，通過對上述設(shè)置的區(qū)間范圍進(jìn)行搜索，以尋求模型精度最大時的最優(yōu)參數(shù)組合。各參數(shù)與結(jié)果表現(xiàn)關(guān)系如圖11所示。圖11(a)為樹的最大深度與擬合優(yōu)度R2的關(guān)系，圖11(b)為葉子節(jié)點(diǎn)所需的最小樣本數(shù)與擬合優(yōu)度R2的關(guān)系，圖11(c)為拆分內(nèi)部節(jié)點(diǎn)所需的最小樣本數(shù)與擬合優(yōu)度R2的關(guān)系，圖11(d)為隨機(jī)森林樹的數(shù)目與擬合優(yōu)度R2的關(guān)系。

表4 設(shè)置需要調(diào)整的參數(shù)及范圍Table 4 Set the parameters and ranges that need to be adjusted

test_score表示模型精度即擬合優(yōu)度R2圖11 各參數(shù)與結(jié)果關(guān)系圖Fig.11 Relationship between parameters and nesults

調(diào)參結(jié)果顯示當(dāng)n_estimators，min_samples_split，max_depth，min_samples_leaf′的取值分別為150、2、6、1時可以得到最小均方根誤差(root mean squared error,RMSE)和最大的擬合度R2值，將參數(shù)調(diào)整后的最優(yōu)參數(shù)組合輸入模型，使用predict函數(shù)來預(yù)測。此時得到的最優(yōu)R2值為0.961 8，RMSE值為0.039 5，平均絕對誤差(mean absolute error，MAE)值為0.021 6。

5 結(jié)果與分析

模型預(yù)測性能對比分析如表5所示。經(jīng)過貝葉斯優(yōu)化和網(wǎng)格優(yōu)化后的R2分別從0.951 0提高至0.975 3和0.961 8，同時參考其他幾個評價指標(biāo)，都顯示預(yù)測效果有明顯提高。其中基于貝葉斯優(yōu)化改進(jìn)隨機(jī)森林模型預(yù)測精度最高，且有明顯提高，均方根誤差值也最小。預(yù)測結(jié)果對應(yīng)的柱狀圖如圖12所示。

表5 模型預(yù)測性能對比分析Table 5 Comparative analysis of model prediction performance

圖12 改進(jìn)前后模型預(yù)測對比柱狀圖Fig.12 Histogram of comparison of model predictions before and after improvement

可以看出，網(wǎng)格優(yōu)化改進(jìn)隨機(jī)森林模型后運(yùn)行時間T由原來的4.95 s提高至3.23 s。貝葉斯優(yōu)化改進(jìn)隨機(jī)森林模型后運(yùn)行時間由原來的4.95 s提高至3.12 s，比網(wǎng)格優(yōu)化提升更高且參數(shù)搜索范圍大大增加。

為了更形象生動的比較上述算法的預(yù)測精度，將改進(jìn)前后的3種模型的預(yù)測結(jié)果與真實(shí)值進(jìn)行比較，如圖13所示，可以看出，基于貝葉斯優(yōu)化和網(wǎng)格優(yōu)化調(diào)整超參數(shù)改進(jìn)后的模型的預(yù)測結(jié)果基本和零誤差線y=x重合，說明貝葉斯優(yōu)化改進(jìn)的隨機(jī)森林模型和網(wǎng)格優(yōu)化改進(jìn)的隨機(jī)森林模型預(yù)測的結(jié)果和真實(shí)值基本一致；其預(yù)測精度在所有算法中達(dá)到最大。

圖13 不同模型測試集真實(shí)值與預(yù)測值相關(guān)性Fig.13 Correlation between the true value and the predicted value of the test set of different models

6 結(jié)論

采用機(jī)器學(xué)習(xí)算法并結(jié)合NPP/VIIRS夜間燈光數(shù)據(jù)進(jìn)行長三角地區(qū)GDP的智能預(yù)測研究，通過對數(shù)據(jù)集的一系列處理，對隨機(jī)森林算法進(jìn)行了貝葉斯優(yōu)化和網(wǎng)格搜索自動調(diào)參，對算法中的超參數(shù)調(diào)優(yōu)，該優(yōu)化方法將預(yù)設(shè)的模型參數(shù)自由組合，對所有參數(shù)組合進(jìn)行遍歷搜索，通過對3種指標(biāo)的分析，比較不同模型的特征并對模型進(jìn)行評估。結(jié)果表明，預(yù)測效果最好的是貝葉斯優(yōu)化改進(jìn)的隨機(jī)森林模型，R2達(dá)到了0.97，優(yōu)化后的隨機(jī)森林模型有明顯提高。后續(xù)的研究可以圍繞著擴(kuò)展模型的訓(xùn)練樣本集、采用更長年限的長時間序列數(shù)據(jù)集、改進(jìn)目前使用的模型等方面進(jìn)一步開展，促使模型更加健壯和可靠。