BP原理及其在林木胸徑模擬中的實現(xiàn)

車少輝 張建國 段愛國 駱昱春

(中國林業(yè)科學(xué)研究院林業(yè)研究所，北京，100091) (江西省林業(yè)科學(xué)院)

人工神經(jīng)網(wǎng)絡(luò)(ANN)是模仿人腦結(jié)構(gòu)及功能的新型信息處理系統(tǒng)，由許多簡單的并行處理單元按某種方式相互連接而成。BP采用誤差反向傳播算法的模型，屬于多層前饋網(wǎng)絡(luò)的一種，也是目前應(yīng)用最廣泛的網(wǎng)絡(luò)類型［1-2］。

近年來，BP在林木生長模型中應(yīng)用越來越多［3-7］，但有關(guān)BP理論和林木生長建模實踐融為一體的介紹很少。盛仲飆［8］介紹BP原理時，主要介紹網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)，而對神經(jīng)元—網(wǎng)絡(luò)結(jié)構(gòu)—算法實現(xiàn)缺少系統(tǒng)連貫的表述。筆者從這3個方面系統(tǒng)剖析BP原理，闡述了其模擬S型曲線的依據(jù)，并借助MATLAB編程語言實現(xiàn)林木胸徑仿真，以期推動BP網(wǎng)絡(luò)的完善和在林業(yè)上的推廣。

1 BP神經(jīng)元的數(shù)學(xué)模型

神經(jīng)網(wǎng)絡(luò)的神經(jīng)元稱節(jié)點或處理單元，具有多輸入單輸出的屬性，是模型的基本組成部分。每個節(jié)點均有相同的結(jié)構(gòu)，其動作在時間和空間上均同步［9-10］。

①神經(jīng)元模型如圖1所示。

②神經(jīng)元模型的數(shù)學(xué)表達(dá)式：

式中：Oj(t)為t時刻神經(jīng)元j的輸出信息;Xi(t)為t時刻神經(jīng)元j接收來自神經(jīng)元i的輸入信息;τij為輸入輸出間的突軸時延;Tj為神經(jīng)元j的閾值;f為神經(jīng)元激活函數(shù)。

τij取單位時間1時：

引入閾值的輸入 X0。令 X0=-1，W0j=Tj，則式中-Tj=X0W0j，合并得

圖1 神經(jīng)元模型

③神經(jīng)元激活函數(shù)。f也稱轉(zhuǎn)移函數(shù)或活化函數(shù)。其基本作用：控制輸入對輸出的激活作用，對輸入、輸出進(jìn)行函數(shù)轉(zhuǎn)換;將可能無限域的輸入變換成指定的有限范圍內(nèi)的輸出。不同的神經(jīng)元數(shù)學(xué)模型主要區(qū)別在于激活函數(shù)不同，使神經(jīng)元具有差異的信息處理特性。

理論上說，任何可微的數(shù)學(xué)函數(shù)都可作為激活函數(shù)，線性的、S型曲線和雙曲線正切是最常用的函數(shù)［11］。

a.線性激活函數(shù)：通常用在輸出層，將輸出結(jié)果映射到任意實數(shù)范圍內(nèi)。

表達(dá)式：

b.Sigmoid函數(shù)：簡稱S型函數(shù)，其特點是函數(shù)本身及其導(dǎo)數(shù)都是連續(xù)的，且是一個嚴(yán)格單增的光滑函數(shù)，并具有漸進(jìn)特性。為實數(shù)域R映射到［0，1］閉集的非減連續(xù)函數(shù)，代表狀態(tài)連續(xù)型神經(jīng)元模型。

Logsigmoid表達(dá)式：

其導(dǎo)數(shù)為：f'(x)=e-x/(1+e-x)2，且 f'(x)=f(x)(1-f(x))。

c.雙曲正切函數(shù)：激活函數(shù)的值域有時在-1～1之間變化，且關(guān)于原點奇對稱。為此，常采用雙極性S型函數(shù)(Tansigmoid)形式：

其導(dǎo)數(shù)為：f'(x)=1/2［1-f2(x)］。

有研究證明［12］，雙曲正切函數(shù)Tansigmoid計算速度比指數(shù)函數(shù)Logigmoid快，實際中經(jīng)常采用雙曲正切函數(shù)作激活函數(shù)，隱層采用非線性激活函數(shù)才可以實現(xiàn)非線性映射功能。

2 BP的拓?fù)浣Y(jié)構(gòu)

BP拓?fù)浣Y(jié)構(gòu)由大量神經(jīng)元連接構(gòu)成的一個層次型網(wǎng)絡(luò)(如圖2)，包括：①含節(jié)點的輸入層：用來描述問題的自變量;②具有節(jié)點的輸出層：描述因變量;③一個或多個包含節(jié)點的隱層：幫助捕獲數(shù)據(jù)中的非線性特征。

圖2 BP類屋次型網(wǎng)絡(luò)結(jié)構(gòu)

前一層的輸出為下一層的輸入，各神經(jīng)元僅接受前一層的輸入，無反饋;輸入層節(jié)點不具有計算功能，單單接受外來信號，并傳遞給各隱層節(jié)點;隱層是神經(jīng)網(wǎng)絡(luò)的內(nèi)部信息處理層，負(fù)責(zé)信息變換，可設(shè)計多層;最后一個隱層傳遞到輸出層，輸出層節(jié)點具備處理信息能力，可以完成非線性和線性的映射。輸入層節(jié)點不存在激活函數(shù)，隱層通常采用非線性激活函數(shù)(S型函數(shù))，輸出層根據(jù)期望輸出的范圍，選擇采用線性或非線性激活函數(shù)。大多數(shù)輸出層神經(jīng)元采用線性函數(shù)，非線性函數(shù)容易造成預(yù)測結(jié)果失真［13］。

K.Funahashi證明［14］，單隱層的網(wǎng)絡(luò)可以任意精度地逼近任何連續(xù)函數(shù)，只有擬合不連續(xù)函數(shù)(如鋸齒波等)時，才需要兩個隱層來縮小誤差。而誤差精度的提高實際上也可以通過增加隱含層中的神經(jīng)元數(shù)目來獲得，其訓(xùn)練效果也比增加層數(shù)更容易觀察和調(diào)整［15］。輸入層和輸出層節(jié)點數(shù)量取決于采用自變量和因變量的數(shù)量，隱層節(jié)點數(shù)量通常采用“試湊法”和“剪枝法”確定。例如：采用林木密度、林齡、立地質(zhì)量來模擬林木的胸徑和樹高，輸入層應(yīng)選3個節(jié)點，即輸入樣本的3個指標(biāo)值作為一個網(wǎng)絡(luò)輸入向量。輸出層選2個節(jié)點，代表胸徑、樹高。

3 BP標(biāo)準(zhǔn)算法

BP模型須經(jīng)過學(xué)習(xí)和工作兩個階段來實現(xiàn)其功能。學(xué)習(xí)或訓(xùn)練過程中，各神經(jīng)元的連接權(quán)動態(tài)調(diào)整規(guī)則稱為學(xué)習(xí)算法或?qū)W習(xí)規(guī)則或訓(xùn)練規(guī)則。

BP采用一種有導(dǎo)師學(xué)習(xí)模式，學(xué)習(xí)階段分解為2個過程(圖3)：第一，信息正向傳播過程(a)。輸入通過輸入層經(jīng)隱層到達(dá)輸出層，計算每個輸出節(jié)點的實際值，與期望值比較，若誤差未達(dá)到規(guī)定的精度要求，即轉(zhuǎn)入誤差反向傳播。若達(dá)到精度要求，完成此樣本學(xué)習(xí)。第二，誤差反向傳播階段(b)。誤差從輸出層經(jīng)過隱層到輸入層，在此過程中采用梯度下降法調(diào)整權(quán)值和閾值，經(jīng)過數(shù)次迭代使得誤差函數(shù)達(dá)到最小。通過不斷的迭代學(xué)習(xí)，當(dāng)網(wǎng)絡(luò)對于各種給定的輸入均能產(chǎn)生所期望的輸出時，完成學(xué)習(xí)過程。

圖3 BP學(xué)習(xí)模式

以單隱層BP(3-2-1)模型為例(圖4)，隱層采用單極性sigmoid函數(shù)，輸出層采用線性函數(shù)，剖析BP算法原理。

圖4 單隱層BP網(wǎng)絡(luò)

輸入向量：X=(x1，x2，x3)T，圖4 中 x0=-1 為隱層神經(jīng)元引入閾值而設(shè)置的;隱層輸出向量Y=(y1，y2)T，圖中y0=-1是為輸出層神經(jīng)元引入閾值而設(shè)置的;輸出層神經(jīng)元O輸出向量為O=(o1);期望輸出向量為d=(d1)。輸入層到隱層之間的權(quán)值矩陣用V表示，V=(V1，V2)，其中列向量V1=(v11，v21，v31)T，V2=(v12，v22，v32)T，vij表示輸入層第 i個神經(jīng)元與隱層神經(jīng)元j的連接權(quán)值。隱層到輸出層之間的權(quán)值矩陣用W表示，W=(w11，w21)T，wjk為第 j個隱層神經(jīng)元與第 k個輸出神經(jīng)元的權(quán)值，由于輸出層僅一個神經(jīng)元，故k=1。v01，v02，v01分別為神經(jīng)元y1，y2，O的閾值。

①信號正向傳播過程中，各層信號之間的數(shù)學(xué)關(guān)系。

接受信號的輸入層：X=(x1，x2，x3)T。

對于隱層：

對于輸出層，有

f(x)=x。將(7)、(8)式及已知量代入(9)中得網(wǎng)絡(luò)模型的輸出：

f(x)=1/(1+exp(-x))。若(10)式中的f采用線性函數(shù)，則輸出向量O永遠(yuǎn)是輸入值加權(quán)的某種線性組合，不可能模擬復(fù)雜的非線性曲線，就無法模擬林木生長的S型曲線。正是由于隱層采用非線性sigmoid類激活函數(shù)，才使得輸出值為不同權(quán)重S型函數(shù)的線性和非線性組合值，適合模擬生長型S型曲線。

②誤差的反向傳播，各層信號之間的數(shù)學(xué)關(guān)系。

當(dāng)網(wǎng)絡(luò)輸出O與期望輸出d存在輸出誤差時，定義誤差能量E。

E=1/2(d-o)2，將(9)、(10)式代入，

由上式可以看出，誤差函數(shù)是各層權(quán)值w，v和輸入x和期望輸出d的函數(shù)，因此在輸入和期望輸出一定的情況下，調(diào)整權(quán)值可改變誤差E。

誤差梯度下降算法又稱最速下降算法，是BP算法中最簡單的一種。誤差函數(shù)E，在某點(wi，vj)的梯度▽E是一個向量，其方向是E(wi，vj)增長最快的方向。顯然，負(fù)梯度方向是誤差減少最快的方向。

對于3層BP網(wǎng)絡(luò)，對輸出層：

no為輸出層神經(jīng)元o的凈輸入。

對隱層：

ny為隱層神經(jīng)元y的凈輸入。

對輸出層和輸入層各定義一個誤差信號，令

綜合式(9)和式(14)，改寫式(12)為：

即輸出層神經(jīng)元k與隱層神經(jīng)元j的連接權(quán)調(diào)整量為輸出層神經(jīng)元k的誤差信號、隱層神經(jīng)元j的輸出、學(xué)習(xí)率η之積的相反數(shù)。

綜合式(7)、式(8)和式(15)改寫式(13)為：

即隱層神經(jīng)元j與輸入神經(jīng)元i的連接權(quán)調(diào)整量為隱層神經(jīng)元j的誤差信號、輸入神經(jīng)元i的輸入、學(xué)習(xí)率之積的相反數(shù)。

如果已知輸出層和輸入層的誤差信號就可以求出權(quán)值調(diào)整量的值，展開式(14)、式(15)，輸出層

隱層

利用E=1/2(d-o)式計算式(18)、式(19)中網(wǎng)絡(luò)誤差對輸出層的偏導(dǎo)利用式(11)中網(wǎng)絡(luò)誤差對隱層輸出的偏導(dǎo)：

將以上兩個偏導(dǎo)結(jié)果代入式(18)和式(19)，應(yīng)用sigmoid函數(shù)f'(x)=f(x)(1-f(x))，得：

其中，輸出層神經(jīng)元激活函數(shù)為f(x)=x;

至此兩個誤差信號的推導(dǎo)已經(jīng)完成，將式(20)，(21)代回式(16)和式(17)推出此3層BP權(quán)值閾值調(diào)整公式為

從權(quán)值調(diào)整過程來看，BP算法有自身的缺陷：容易陷入誤差局部最小值。而且標(biāo)準(zhǔn)算法收斂結(jié)果與初始點的選擇有關(guān)，最終不一定能夠收斂到全局誤差最小點。因而各種改進(jìn)BP算法應(yīng)運而生，例如：建立在確定梯度下降的最優(yōu)化方向上的增加動量項法，或者基于最優(yōu)化理論的訓(xùn)練算法(共軛梯度算法、擬牛頓法、Levenberg-Marquardt算法)。

在單個處理單元層次，無論采用哪種學(xué)習(xí)規(guī)則進(jìn)行調(diào)整，其算法都十分簡單，但當(dāng)大量處理單元集體進(jìn)行權(quán)值調(diào)整時，網(wǎng)絡(luò)就出現(xiàn)智能的特性，其中有意義的信息就分布地存儲在調(diào)節(jié)后的權(quán)值矩陣中。

經(jīng)過學(xué)習(xí)過程的網(wǎng)絡(luò)，才是一個成熟的網(wǎng)絡(luò)模型。給網(wǎng)絡(luò)一個輸入向量，網(wǎng)絡(luò)根據(jù)儲存的信息輸出模擬結(jié)果。BP網(wǎng)絡(luò)的工作過程很簡單，不需要迭代運算。

4 BP網(wǎng)絡(luò)林木平均胸徑仿真

試驗數(shù)據(jù)取自大崗山年珠林場杉木密度實驗林。密度試驗林于1982年使用1年生的實生苗造林，每個區(qū)組分5個水平：2 m×3 m(A 密度，1 667 株/hm2)、2 m×1.5 m(B 密度，3 333株/hm2)、2 m×1 m(C 密度，5 000 株/hm2)、1 m×1.5 m(D 密度，6 667株/hm2)、1 m×1 m(E 密度，10 000株/hm2)，3 個區(qū)組，共15個小區(qū)。采用隨機區(qū)組排列，每個小區(qū)外設(shè)置保護行，四角設(shè)置固定水泥樁為界。

調(diào)查方法：對每株林木進(jìn)行掛牌，各小區(qū)沿等高線蛇形順序編號，并標(biāo)注胸高位置。每次調(diào)查均于林分停止生長后或下一年開始生長前進(jìn)行。樹齡在10 a前逐年記錄，之后隔年調(diào)查1次。2～5 a時，測定每木樹高、冠幅及樹高達(dá)到1.3 m以上的胸徑。6年生開始，在每個小區(qū)的上、中、下各選2株優(yōu)勢木，用以求算立地指數(shù)。截止2008年，連續(xù)觀測14次。根據(jù)時間間隔一致的原則，剔除第7、9年觀測數(shù)據(jù)，選用12次測量數(shù)據(jù)，共計180個樣本對建立神經(jīng)網(wǎng)絡(luò)模型。

由林木生長理論可知胸徑生長是一個非線性過程，適合于BP模擬。K.Funahashi證明單隱層BP模型能夠以任意精度逼近任意函數(shù)，筆者選取單隱層BP模型。輸入向量為林分3個屬性值：林木密度D(株/hm2);立地指數(shù)S(m);林齡A(a)。輸出：林木平均胸徑Dg(cm)。首先，將輸入輸出向量分別歸一化到［-1，1］區(qū)間，以消除量綱差異和加速網(wǎng)絡(luò)收斂過程。隱層和輸出層預(yù)設(shè)激活函數(shù)分別為tansigmoid函數(shù)和purelin函數(shù);算法采用最速梯度下降法;隱層節(jié)點數(shù)用“試湊法”確定。輸入向量格式為一列一個樣本，包含3個屬性值，同時此樣本有對應(yīng)的輸出向量，兩者構(gòu)成輸入輸出樣本對。輸入的列數(shù)與輸出的列數(shù)一致，即樣本數(shù)量一定;而兩個向量的維數(shù)(行數(shù))可以不同。

利用MATLAB神經(jīng)網(wǎng)絡(luò)工具箱內(nèi)置函數(shù)newff創(chuàng)建BP模型，格式為：newff(P，T，［S1S2…S(N-1)］，{F1F2…F(N+1)})，P，T分別為輸入和期望響應(yīng)矩陣;Si為各隱層神經(jīng)元的數(shù)目，不包括輸入層和輸出層;Fi為網(wǎng)路各層的傳遞函數(shù)類型，包括隱層和輸出層的傳遞函數(shù);數(shù)據(jù)集隨機分配，且訓(xùn)練集∶驗證集 ∶測試集=0.7 ∶0.15 ∶0.15，未給出的參數(shù)均取 MATLAB默認(rèn)值。

P=［D;S;A］;網(wǎng)絡(luò)模型輸入三維向量;

T=［Dg］;期望輸出一維向量;

＞＞W(wǎng)=newff(p，t，5，{‘tansig’，‘purelin’}，‘traingd’);創(chuàng)建BP網(wǎng)絡(luò)模型W;

＞＞W(wǎng).trainparam.lr=0.05;設(shè)置學(xué)習(xí)率;

＞＞W(wǎng).trainparam.epochs=1000;最大迭代次數(shù);

＞＞W(wǎng).trainparam.goal=0.01;最小誤差;

＞＞［W，tr］=train(W，p，t);訓(xùn)練網(wǎng)絡(luò)，存儲訓(xùn)練參數(shù)。

＞＞Y=sim(W，p);仿真輸出;

最后，利用postmnmx()函數(shù)對Y反歸一化還原為胸徑值，完成模擬過程。

為了直觀理解模擬情況，筆者選擇MATLAB輸出的兩幅圖說明。圖5為不同數(shù)據(jù)集實測值(x軸)與模擬值(y軸)的關(guān)系，其中訓(xùn)練集、測試集、驗證集和全部建模數(shù)據(jù)集的相關(guān)系數(shù) R 分別為：0.977、0.962、0.985、0.977。圖 6 展示了各種數(shù)據(jù)均方誤差隨迭代次數(shù)的變化，圖中從細(xì)到粗依次為測試數(shù)據(jù)、訓(xùn)練數(shù)據(jù)、驗證數(shù)據(jù)的模擬誤差變化曲線。其中迭代935次時，模型收斂到Ems=0.01。模型參數(shù)存儲在W中，利用MATLAB相關(guān)函數(shù)可以顯示具體值，在此不再累述。

圖5 模擬值與實測值的相關(guān)關(guān)系

5 結(jié)論與討論

BP網(wǎng)絡(luò)是一種有效的建模工具，廣泛適應(yīng)于各個領(lǐng)域。林木生長建模上應(yīng)用BP模型的文章也越來愈多，筆者利用MATLAB神經(jīng)網(wǎng)絡(luò)工具箱實現(xiàn)了胸徑仿真功能(R=0.977)，結(jié)果說明，采用BP靜態(tài)網(wǎng)絡(luò)逼近動態(tài)數(shù)據(jù)是一種可行的方法，為BP模型在林木生長模型中的應(yīng)用提供支持。

圖6 均方誤差隨迭代次數(shù)的變化

BP標(biāo)準(zhǔn)算法(梯度最速下降算法)進(jìn)行訓(xùn)練時，模型收斂速度慢(迭代935次)，甚至有多次不收斂的情況出現(xiàn);筆者側(cè)重于研究BP原理及其在MATLAB中實現(xiàn)的過程，未對采用改進(jìn)算法的模型進(jìn)行比較。下一步的工作討論各種算法對模型的收斂情況及其精度、泛化能力的影響。不同版本的MATLAB構(gòu)建BP模型格式有些細(xì)微差別，但都能滿足建模需求。

［1］Basheer I，Hajmeer M.Artificial neural networks：fundamentals，computing，design，and application［J］.Journal of Microbiological Methods，2000，43(1)：3-31.

［2］Rumelhart D，Hinton G，Williams R.Learning internal representations by error propagation，Parallel distributed processing：explorations in the microstructure of cognition，vol 1：foundations［M］.Cambridge：MIT Press，1986.

［3］陳晨，郭芳，黃家榮，等.杉木人工林直徑分布BP模型的研究［J］.河南農(nóng)業(yè)大學(xué)學(xué)報，2005，39(4)：390-393.

［4］吳建華，盧炎生.基于神經(jīng)網(wǎng)絡(luò)的樹木生長預(yù)測［J］.長沙電力學(xué)院學(xué)報：自然科學(xué)版，2002，17(2)：28-30.

［5］黃家榮，任誼群，高光芹.森林生長的人工神經(jīng)網(wǎng)絡(luò)建模［M］.北京：中國農(nóng)業(yè)出版社，2006.

［6］金星姬，賈煒瑋，李鳳日.基于BP人工神經(jīng)網(wǎng)絡(luò)的興安落葉松天然林全林分生長模型的研究［J］.植物研究，2008，28(3)：370-374，384.

［7］洪偉，吳承禎.基于人工神經(jīng)網(wǎng)絡(luò)的森林資源管理模式研究［J］.自然資源學(xué)報，1998，13(1)：69-72.

［8］盛仲飆.BP神經(jīng)網(wǎng)絡(luò)原理及MATLAB仿真［J］.渭南師范學(xué)院學(xué)報：綜合版，2008，23(5)：65-67.

［9］韓力群.人工神經(jīng)網(wǎng)絡(luò)教程［M］.北京：北京郵電大學(xué)出版社，2006.

［10］田雨波.混合神經(jīng)網(wǎng)絡(luò)技術(shù)［M］.北京：科學(xué)出版社，2009.

［11］Diamantopoulou M J，Milios E.Modelling total volume of dominant pine trees in reforestations via multivariate analysis and artificial neural network models［J］.Biosystems Engineering，2010，105(3)：306-315.

［12］Swingler K.Applying neural networks：a practical guide［M］.London：Academic Press，1996.

［13］Khashei M，Bijari M.An artificial neural network(p，ád，áq)model for timeseries forecasting［J］.Expert Systems with Applications，2010，37(1)：479-489.

［14］Funahashi K.On the approximate realization of continuous mappings by neural networks［J］.Neural Networks，1989，2(3)：183-192.

［15］韓震，姜照華.基于神經(jīng)網(wǎng)絡(luò)的大連經(jīng)濟增長模擬與預(yù)測［J］.大連理工大學(xué)學(xué)報，2001，41(6)：752-756.