傳遞函數模型研究張柏發摘要：本文對時間序列理論模型進行了研究。從ARMA模型的結構與建模方法出發，重點介紹了傳遞函數模型的定義、原理、識別、擬合、診斷。傳遞函數模型是時間序列分析中一種典型動態響應模型，其傳遞函數部分能體現一個輸出序列與輸入序列的關系，噪聲部分又能反映了輸出序列自身的時序關系。并進一步介紹了傳遞函數模型在土壤溶質運移模擬中的應用和傳遞函數模型的構建方法之一——概率神經網絡模型。本章主要介紹了理論模型和理論知識，包括時間序列平穩性檢驗，ARMA模型識別、擬合、預測，再重點介紹了傳遞函數模型的原理、傳遞函數的性質、互相關函數、互相關函數與傳遞函數的關系，以及模型的建立過程等內容。傳遞函數模型是時間序列分析中一種基于一個時間序列當前觀測值不僅與本身過去觀察值有關，還受到其它時間序列的影響，而提出的動態響應模型，在經濟、工業生產、生物、醫學等領域均有廣泛的應用。傳遞函數模型是時間序列分析中一種典型動態響應模型，其傳遞函數部分能體現一個輸出序列與輸入序列的關系，噪聲部分又能反映了輸出序列自身的時序關系。向陽等(2006)通過分析降水和地下水埋深之間的關系，利用傳遞函數理論和回歸分析方法建立降水估計地下水埋深統計模型，提高了擬合精度。LiChunyan和ChenJun(2009)分析了交通事故的宏觀影響因素建立了傳遞函數模型，XuYajing和LiXianglu(2009)應用傳遞函數模型進行經濟預測，Kinley、Pratop(2010)也應用傳遞函數模型在生物方面做研究。傳遞函數模型可以看成是帶有輸入序列的ARIMA模型,它是ARIMA模型與普通回歸模型的結合。如果多個動態數據所成的時間序列都是平穩的,或經過適當差分后平穩,即可建立傳遞函數模型.傳遞函數模型的建模過程非常復雜,必須借助于計算機才能完成。當前最為流行的統計分析軟件—SAS在進行多元時間序列分析中時,具有其他統計軟件無可比擬的優勢。1時間序列理論模型介紹1.1時間序列基本概念時間序列分析可以從數量上揭示某種系統或現象的發展變化規律，或者從動態的角度刻畫某現象與其它現象的依存關系以及發展的變化規律，運用時間序列模型可以預測系統的未來情況。時間序列模型已經成為現代各個領域的重要研究工具，廣泛應用于各個領域的分析與預測，包括經濟、金融、工業、商業等方面。1.1.1時間序列的平穩性⑴平穩序列如果時間序列{yt,t=1,2,…}滿足下面①均值為常數，(t=1,2,…)；②方差(t=1,2,…)；③協方差不依賴于t，(t=1,2,…)；則稱時間序列{yt,t=1,2,…}為寬平穩時間序列，也叫廣義平穩時間序列，簡稱平穩序列。對于平穩的時間序列at，如果對于所有的t和k≠0成立，則稱時間序列at為白噪聲序列。⑵非平穩序列所謂一個時間序列是非平穩的，是指時間序列的統計規律會隨著時間的平移而發生變化，也就是說生成變量時間序列數據的隨機過程的統計特征隨時間的變化而變化。只要弱平穩的三個條件不全滿足，就稱該時間序列是非平穩的。通常我們所說的平穩性就是弱平穩性。因為平穩的時間序列有穩定的發展趨勢(均值)、波動性(方差)和橫向聯系(協方差)，因此可以用時間序列的樣本均值和方差推斷各時點隨機變量的分布特征。對于平穩序列，任何震蕩的影響都是暫時的。隨著時間的推移，這些影響將會逐漸消失，時間序列將回復到長期的平均水平。因此，用平穩的時間序列數據的回歸才是有效的。因此，在分析之前，判斷我們所采用的時間序列的平穩性是非常必要的。1.2時間序列平穩性檢驗時間序列的平穩性是時間序列計量分析有效性的基礎，因此，時間序列的平穩性檢驗具有重要的實際意義。下面主要介紹單位根檢驗方法。在式子(1)中，若=0，則有(2)其中，ut為白噪聲序列，上式稱為一階自回歸過程，記為AR(1)。(3)當時，有①；②；③當時，yt滿足平穩性的三個條件。所以，接受表明yt是非平穩序列，而拒絕原假設則就表明序列yt是平穩的。式稱為p階自回歸過程，記為AR(p)。1.3ARMA模型1.3對于平穩序列{yt},{at}為零均值白噪聲，滿足：(4)其中,B為向后推移算子，稱模型為自回歸移動平均模型(Auto-Regressive-Moving-AverageModel)，簡記為ARMA(p,q)。當q=0時，稱為p階自回歸模型，記為AR(p)；當p=0時，稱為q階滑動平均過程，記為MA(q)。1.模型定階的方法主要有殘差方差圖定階法、F檢驗定階法和準則函數定階法。常用的方法是準則函數定階法，下面介紹幾種準則函數定階法。準則定階法，即確定一個準則函數，建模時按照該準則函數的取值大小確定模型的優劣，使準則函數達到極小的是最佳模型階數。⑴AIC準則AIC準則(A-InformationCriterion最小信息準則)最先是日本學者Akaike提出，也稱為赤池信息準則。它適用于AR、MA、ARMA三類模型的定階。設{yt：1≤t≤n}為一隨機時間序列，對其擬合ARMA(p,q)模型，采用極大似然方法對模型估計參數，L為模型的極大似然值，AIC準則函數定義如下：(5)其中,r=p+q為模型中獨立參數的個數；是殘差方差的極大似然估計。在實際中，常用以下定義的AIC準則函數(用樣本大小n標準化)(6)對于事先給定的最高階數M(n)，如果(7)就認為p0和q0是最佳模型階數。⑵BIC準則BIC信息準則與AIC類似，是由Schwarz提出來的，也稱為Schwarz信息準則(記為SIC)，對于ARMA模型定階，BIC準則函數定義如下;(8)對于事先給定的最高階數M(n)，若某一階數(9)就認為p0和q0是最佳模型階數。⑶CAT準則parzen在1977年建議采用如下模型選擇準則，稱為CAT：(10)其中，是當用AR(j)模型對序列擬合時的無偏估計，n是觀測個數。當CAT(p)達到極小時，稱p為最優階數。1.模型識別結束后，需要進行參數的估計，得出模型才是我們要的目的。參數估計方法主要有矩估計法、最小二乘估計法和極大似然估計法等。這里主要介紹矩估計法。矩估計法的基本思想是：ARMA模型的自相關函數(矩函數)可以表示為未知的模型參數的函數，反過來，模型參數原則上也可由自相關函數(矩函數)來表示。理論自相關函數由計算出的樣本自相關函數代替，就得到了參數的估計值。如AR(p)過程：(11)均值用估計。為了估計，通過(12)得到如下的Yule-Walder方程組(13)接著，用代替，通過計算前面的線性方程組得到矩估計，…,,即(14)該估計量我們通常稱為Yule-Walder估計。在計算出，…,后，利用下面結果(15)就可以得到的矩估計為(16)1.模型的適應性是指一個時間序列模型解釋系統動態性(即數據序列的相關性)的程度。一個適合的時間序列模型應該是完全或者基本上解釋了系統的動態性。因此，模型中殘差序列{at}應該是白噪聲序列。故{at}序列的獨立性檢驗就是模型的適應性檢驗。相關函數法是先通過計算殘序差列{}的自相關函數，然后進行判斷殘差序列的獨立性。設表示{}序列的自相關函數，即(17)如果殘差序列{at}是白噪聲序列，那么是互不相關的，且~N(0,1/n)。1.4ARIMA模型在實際中，許多時間序列過程是非平穩的，這些非平穩時間序列可能有隨時間變化的均值，和隨時間變化的二階矩(如方差)。此時就需要對時間序列做適當的變換使其成為平穩序列。對于線性趨勢的時間序列，可以對時間序列進行差分，即可消去線性趨勢，差分后的時間序列就是平穩的。然后在對差分后得到的平穩時間序列進行建模建立ARMA模型。利用差分來建立的時間序列模型，全稱為差分自回歸移動平均模型(AutoregressiveIntegratedMovingAverageModel,簡記ARIMA)。其中ARIMA(p,d,q)稱為差分自回歸移動平均模型，AR是自回歸,p為自回歸項；MA為移動平均，q為移動平均項數，d為時間序列成為平穩時所做的差分次數。即ARIMA模型為(18)對于隨著時間推移，方差不斷增加的時間序列，可以通過其它適當的變換消除方差的不平穩性對時間序列的影響。常用的變量替換有：對數變換、平方根變換等等，具體方法視實際情況而定。ARIMA模型的基本思想是:將預測對象隨時間推移而形成的數據序列視為一個隨機序列,用數學模型來描述,模型被識別后可以從時間序列的過去值及現在值來預測未來值。ARIMA模型預測的基本程序如下:1.根據時間序列的散點圖、自相關函數和偏自相關函數圖以ADF單位根檢驗其方差、趨勢及其季節性變化規律,對序列的平穩性進行識別。2.對非平穩序列進行平穩化處理。如果數據序列是非平穩的,并存在一定的增長或下降趨勢,則需要對數據進行差分處理,如果數據存在異方差,則需對數據進行技術處理,直到處理后的數據的自相關函數值和偏相關函數值無顯著地異于零。3.根據時間序列模型的識別規則,建立相應的模型。若平穩序列的偏相關函數是截尾的,而自相關函數是拖尾的,可斷定序列適合AR模型;若平穩序列的偏相關函數是拖尾的,而自相關函數是截尾的,則可斷定序列適合MA模型;若平穩序列的偏相關函數和自相關函數均是拖尾的,則序列適合ARMA模型。4.進行參數估計,檢驗是否具有統計意義。5.進行假設檢驗,診斷殘差序列是否為白噪聲。6.利用己通過檢驗的模型進行預測分析。可用SPSS統計分析軟件來進行ARIMA模型預測,可以省去繁瑣的計算環節,簡便、快捷,精度高。ARIMA模型是單變量模型中最為完善的預測模型,但是它要求時間序列自身的變化呈現出某種規律性,而且這種規律還會持續到未來,另外,使用ARIMA模型進行預測時,未來各期的預測往往是單調上升或下降的,除非序列具有明顯的季節性變化規律。普通多元回歸分析在許多情況下并不適用,原因是自變量對因變量的作用和影響方式并非同步變化,一些自變量可能會領先因變量的變化,一些會對若干期的因變量持續產生影響。因此上述兩類模型的應用都有一定的局限性。如果能將上述兩種模型的優點結合在一起,會得到更佳的預測效果,模型的應用范圍會更加廣泛。這種模型就是傳遞函數模型,也稱多變量ARIMA模型或ARIMAX模型。1.5傳遞函數模型一般形式傳遞函數模型可分為無噪聲的傳遞函數模型和附加噪聲的傳遞函數模型，實際問題中，輸出序列往往會受到輸入序列之外的各種干擾，因此附加噪聲的傳遞函數模型應用更廣。假設xt，yt都是可以通過適當變換使之成為平穩的序列。在單輸入單輸出系統中，輸出序列yt與輸入序列xt通過一個線性濾波相關聯：(19)其中，作為濾波的傳遞函數，xt將輸入序列過濾傳遞到輸出序列yt，B為向后推移算子，nt是與輸入序列xt獨立的系統噪聲序列。Box，Jenkins和Reinsel稱方程(19)為傳遞函數模型。若xt和nt服從ARMA模型時，式子(19)就是我們所熟知的ARMAX模型。系數vj稱為系統的脈沖響應權重或脈沖響應函數，是j的函數。如果這些脈沖響應權重序列是絕對可和的，即，則傳遞函數模型被稱為是穩定的。實際應用中，為了降低難度，我們用有理函數形式給出傳遞函數：(20)因此，實際應用的單變量傳遞函數模型一般表達式：(21)其中at為白噪聲序列,b是一個滯后參數，表示輸入變量對輸出變量產生影響前需要的實際時間滯后期數。式(21)稱為傳遞函數。傳遞函數模型一般要求觀測值至少應有30個。建模前要求序列均應為平穩時間序列,平穩是指該序列的概率結構不隨時間而變,即均數不隨時間變化,方差不隨時間變化,自相關系數只與時間間隔有關而與所處的時間點無關。最簡單的平穩序列就是白噪聲。通常,我們可以通過自相關圖檢驗法或單位根檢驗法等方法來檢驗序列的平穩性。如果序列不平穩,可利用對數變換、平方根轉換、濾波、差分等手段來使之平穩。傳遞函數模型的建模一般分為3步:第1步,利用預白噪化進行識別（確定參數b,s,r）,使用的分析工具是輸入和輸出之間的互相關函數。具體做法就是用一個ARIMA模型擬合輸入序列,該模型把殘差降為白噪聲,然后再用此模型過濾輸入序列得到白躁聲殘差序列。隨后用同樣的模型過濾輸出序列,并計算過濾后輸出序列和過濾后輸入序列的互相關系數。第2步,運用極大似然法、最小二乘法或非線性最小二乘法進行模型參數估計（確定）。第3步,利用Q統計量對殘差進行自相關檢驗,利用S統計量對殘差和白噪化輸入進行互相關檢驗。如果殘差無自相關,同時殘差和白噪化輸入也無相關,模型即為所求。1.ARMA模型定階通過自相關函數和偏自相關函數初步判定，那么傳遞函數模型的統計性質也為傳遞函數的三個參數r，s和b的判定提供了依據。由傳遞函數式(20)，即可以得到r，s和b的階數以及它們與脈沖響應權重vj的關系因此，可得(22)r個脈沖響應權重vb+s，vb+s-1，…，vb+s-r+1為差分方程(23)提供了初值。因此，傳遞函數(20)的脈沖響應權重應具有下述性質：(1)b個零權重，v0=v1=…=vb-1=0；(2)s-r+1個權重vb，vb+1，…vb+s-r不固定形式；若s<r，不會出現這樣的值；(3)r個初始脈沖響應權重vb+s-r+1，vb+s-r+2，…vb+s；(4)對于j>b+s，vj服從等式(23)給出的模型。這里b是由j<b時vj=0,vb≠0確定的。r的值是由脈沖響應權重的形式來確定的。對于一個給定的b值，若r=0，則s的值可以通過j>b+s時vj=0計算得出；若r≠0，則s的值可以通過觀測脈沖響應權重模式在何時開始衰減得出。1.穩定性意味著當一個有界的輸入就會產生一個有界的輸出。傳遞函數模型的穩定性要求與ARMA模型平穩性要求類似。有所不同的是，除了要求傳遞函數部分的穩定性，同時也要求噪聲部分的平穩性。如果脈沖響應權重是絕對收斂的，即,且特征方程的根在單位圓內，此時系統稱為穩定系統。對于噪聲部分平穩性的要求與對ARMA模型平穩性要求一樣，要求特征方程的根在單位圓內。對于非穩定系統，可以通過適當的差分或其它變換將系統轉化為穩定系統。1.⑴互相關函數(CCF)①互協方差函數和互相關函數互相關函數是測度兩個隨機變量之間相關的強度和方向的常用函數。給定兩個隨機過程xt和yt，t=0，±1，±2，…當xt和yt都是一唯平穩過程，并且xt和yt的互協方差函數cov(xt,ys)僅與時間差(s-t)有關時，就稱序列xt和yt是聯合平穩序列。這種情況下，計算xt和yt的互協方差函數：(24)其中，,通過標準化，可得互相關函數（CCF）：(25)其中，和是xt和yt的標準差。當且僅當k>0，=0，或者k＜0，=0時，序列xt和yt存在一個因果傳遞函數模型。若對某些k>0，≠0時，而對所有的k＜0，=0時，則稱序列xt導致yt。若對某些k＜0，≠0時，而對所有的k>0，=0，則稱序列序列yt導致xt。若對某些k>0，≠0，且對某些k＜0，≠0，則稱序列xt和yt存在互相反饋關系。若≠0，則稱xt和yt存在同期關系。②樣本互相關函數在實際應用中，因為總體的互相函數是未知的，所以通常用樣本互相關函數作為總體的互相函數估計。對于給定平穩時間序列xt和yt(1≤t≤n)，互相關函數可以通過樣本互相關函數(26)估計出來，其中(27)其中,、、Sx和Sy分別是兩個序列的均值和標準差。③互相關函數與傳遞函數關系互相關函數是識別傳遞函數的工具，而脈沖響應權重與互相關函數有關，故只要確定脈沖響應權重，傳遞函數隨之確定。對于時間t+k,轉換函數模型可以寫為：(28)不失一般性，假設在方程(28)兩邊同時乘以xt,兩邊求期望得：(29)因為所以因此(30)從(30)可以看出，互相關函數受到輸入序列xt自相關函數和脈沖響應權重vj的影響。當輸入序列xt是白噪聲序列時，即對k≠0，=0，方程(30)可以化簡成：(31)此時，脈沖響應權重vk與互相關函數成比例關系。⑵模型初步識別①預白化輸入序列xt如果輸入序列xt是白噪聲，則可以得到如式(31)脈沖響應權重與互相關函數的關系式，因此先對輸入序列xt做預白化處理。設傳遞函數為(32)設輸入序列xt是一個平穩序列，其適應的模型為即其中，at是期望為0,方差為的白噪聲序列。將看成一個濾波器，對yt進行濾波變換，得(33)將(32)代入式(33)，化簡得令,則,計算和之間的樣本的CCF，,估計vk：②傳遞函數部分識別由式(22)的第一式可知前面b個脈沖響應權重為零，即由傳遞函數的常見形式，知可由脈沖響應權重的呈現形式得到階數r的值，又由式(31)，顯然可以直接由互相關函數圖判斷得到階數r的值。對于階數s，若衰減從vb開始，則s=0；若衰減從vb+1開始，則s=1；若衰減從vb+2開始，則s=2。③噪聲部分識別由模型，有(34)傳遞函數部分估計后，計算噪聲序列估計（殘差序列）噪聲部分模型可以通過檢驗其樣本的ACF和PACF進行識別：最終即可得到傳遞函數模型：(35)1.⑴模型估計傳遞函數模型識別之后，就要進行參數估計，和ARMA模型一樣，主要的估計方法有矩估計法、最小二乘估計法和極大似然估計法等。以最小二乘估計法為例，記待估計參數為：由傳遞函數模型(35)，可以改寫模型為(36)由樣本序列（xt，yt）(t=1,2,…,N)求得樣本殘差序列(t=1,2,…,N),是未知參數的函數，即,當(37)達到最小值時，求出的就是參數的最小二乘估計，且是白噪聲at的方差估計量。⑵模型診斷檢驗在模型識別和參數估計之后，還需對模型的進行適當性檢驗。檢驗內容分兩部分，第一，傳遞函數模型是否擬合得當；第二，殘差序列與輸入序列xt派生出的預白化序列是否互相關。①殘差序列自相關檢驗傳遞函數模型包括傳遞函數部分和噪聲部分。如果傳遞函數模型擬合不當，則殘差序列表現出自相關，如果是適當的模型，則殘差序列自相關和偏自相關函數應無任何表現形式，此時殘差序列可以看成白噪聲序列的一個樣本,當樣本容量足夠大時，的自相關函數相互獨立地服從N(0,1/m)，這里，m=n-u-p為可計算的觀測值個數，其中，n是有效樣本個數,u=max(r,s+b)檢驗統計量其中，p和q是噪聲部分模型的參數個數；K一般取足夠大。②殘差序列互相關檢驗在傳遞函數模型(36)中，假設at是噪聲序列，而且和輸入序列xt是互相獨立，因此它就和輸入序列xt派生出的預白化序列獨立。如果傳遞函數適當，殘差和輸入序列xt派生出的預白化序列之間樣本互相關函數獨立漸近服從N(0,1/m)。檢驗統計量S(38)因為互相關函數有同期相關，所以式(38)從k=0開始。r+s+1是傳遞函數部分的參數個數。1.6多變量傳遞函數模型輸出序列可能被多個輸入序列影響，所以推廣到多變量傳遞函數模型或者其中，是輸入序列xjt的j次傳遞函數，且at與每一個輸入序列xjt，j=1，2，…，k獨立。1.7變量帶干預的時間序列模型時間序列常受突發事件或特殊事件影響，諸如經濟政策的改變、新環境法規的采納、罷工以及廣告促銷活動。我們稱這類外部事件為干預事件。干預分析模型是傳遞函數模型的一個特殊情況，使用傳遞函數模型可以解釋干預事件對時間序列的影響。這里的輸入序列是簡單的指示性變量形式，一個僅僅取0或1的輸入變量。1.單干預輸入情形，干預分析模型有如下形式其中，表示干預事件對時間序列yt的影響，It是輸入干預變量，nt是噪聲序列，假設nt服從ARIMA模型，。噪聲模式通常是時間序列yt在受到干預事件影響之前，用單變量模型識別方法加以識別。⑴干預變量常用的干預變量有兩種類型：一類表示干預在時刻T發生，但其影響仍然保持，即干預是一個階躍函數：另一類表示是干預只在一個時點T發生，即干預是脈沖函數：這里脈沖函數可以通過對階躍函數差分而得到，即因此，一個干預模型可以等價地用階躍函數或脈沖函數表示。⑵干預響應形式對于階躍或脈沖干預有各種可能的響應，下面是常見的響應形式。①突然發生，影響持續時間長久的輸入干預影響的形式為這里輸入干預變量是，輸出變量b期才做出反映，且影響強度為，之后繼續保留下去，再不回到以前的狀況。②緩慢發生，影響持續時間長久的輸入干預影響形式為。③突然發生，影響持續時間短暫的輸入干預影響形式為，這里輸入干預變量是,輸出變量b期才做出反映，且影響強度為,之后干預影響消失，回到以前的狀況。④緩慢發生，影響持續時間短暫的輸入干預影響形式為。⑤各種響應都可以用階梯和脈沖函數的不同組合來生成，如一般地，一個響應可以表示為有理函數。b是干預影響的時間延遲，多項式的權重是用來表示干預影響的初始期值，多項式的根假設在單位圓上或者在圓外。單位圓上的根刻畫額線性增長的響應現象，而單位圓外的根則表明具有漸變的影響。1.干預分析模型雖然是傳遞函數模型中的特殊情況，但是在模型識別上有差異，這里的干預傳遞算子能以預白化技術為基礎，這里應該通過考慮那些有可能造成變化或者影響的機制以及預期變化的隱含形式，來假定干預模型的形式。實際中，我們經常通過直接觀察數據，以及數據序列圖，以初步判定干預事件對時間序列的影響形式，初步得出干預算子形式。模型初步識別后，再利用傳遞函數模型參數估計方法估計參數值，通過參數顯著性檢驗、殘差序列的自相關檢驗，檢驗模型是否是適應的。1.對于多干預輸入情形，有一般模型這里，Ijt，j=1，2，…，k是干預變量，這些干預變量或者是階梯函數或者是脈沖函數。一般地，可以是特定的示性變量。第j個響應的式子是根據對給定干預時的預期響應而設定的。對于一個非平穩過程，模型通常不包含常數項。1.8變量帶干預的傳遞函數模型當時間序列受到外界突發事件影響時，需要考慮帶干預變量，做干預分析模型。那么當一個時間序列不僅受到其它時間序列影響，而且還受到外界突發事件影響時，例如股市上的成交量與收盤價有影響，同時收盤價是個波動較大的序列，一旦有突然事件發生，收盤價受到較大的波動，此時不能再簡單只是對該時間序列建立干預分析模型，也不能簡單用輸入序列建立輸出序列的傳遞函數模型。針對這類情況的出現，提出了變量帶干預的傳遞函數模型。1.⑴輸出變量帶干預的傳遞函數模型當一個時間序列受到其它序列影響，同時受到外界干預時，并且這里的突發事件只對該時間序列有影響。那么我們考慮對輸出變量建立帶有輸入變量和干預變量的傳遞函數模型。簡稱輸出變量帶干預的傳遞函數模型，模型一般形式如下：yt是輸出變量，xt是輸入變量，b(x)是xt的延遲參數，It(y)是yt的干預變量，b(I)是It(y)的延遲參數，這些干預變量或者是階梯函數或者是脈沖函數,at(y)是yt白噪聲序列。推廣到更多變量多干預情況：yt是輸出變量，xjt(j=1,2,…,k)是輸入變量，bj(x)(j=1,2,…,k)是xjt的延遲參數，Ijt(y)(j=1,2,…,m)是yt的干預變量，bj(I)(j=1,2,…,m)是Ijt(y)的延遲參數，這些干預變量或者是階梯函數或者是脈沖函數,at(y)是白噪聲序列。⑵輸入變量帶干預的傳遞函數模型當輸入變量受到外界影響時，且影響是十分顯著的。由于輸出變量與輸入變量之間存在單向關系，所以在事件發生后，輸出變量也受到影響，在這種情況下可以考慮事件發生時，輸入序列建立干預分析模型，而不再是ARMA模型，再對輸出序列建立傳遞函數模型。簡稱輸入變量帶干預的傳遞函數模型。模型一般形式如下：yt是輸出變量，xt是輸入變量，b(x)是xt的延遲參數，It(x)是xt的干預變量，b(I)是It(x)的延遲參數，這些干預變量或者是階梯函數或者是脈沖函數,at(y)和at(x)是yt和xt白噪聲序列。推廣到多變量多干預情況：yt是輸出變量，xjt(j=1,2,…,k)是輸入變量，bj(x)(j=1,2,…,k)是xjt的延遲參數，Ijt(x)(j=1,2,…,n)是xjt的干預變量，bj(I)(j=1,2,…,n)是Ijt(x)的延遲參數，這些干預變量或者是階梯函數或者是脈沖函數,at(y)和at(x)是yt和xt白噪聲序列。⑶輸入輸出變量帶干預傳遞函數模型當輸出變量和輸入變量同時受到不同的干預事件影響時，當這些的干預事件不相影響，就需要分別對輸出變量和輸入變量都建立帶有干預變量的模型，再建立傳遞函數模型，簡稱輸入輸出變量帶干預傳遞函數模型。模型一般形式如下：yt是輸出變量，xt是輸入變量，b(x)是xt的延遲參數，It(y)是yt的干預變量，是It(y)的延遲參數，It(x)是xt的干預變量，是It(x)的延遲參數，這些干預變量或者是階梯函數或者是脈沖函數,at(y)和at(x)是yt和xt白噪聲序列。推廣到多變量多干預情況：yt是輸出變量，xjt(j=1,2,…,k)是輸入變量，bj(x)(j=1,2,…,k)是xjt的延遲參數，Ijt(y)(j=1,2,…,m)和Ijt(x)(j=1,2,…,n)分別是yt和xt的干預變量，和分別是Ijt(y)和Ijt(x)的延遲參數，這些干預變量或者是階梯函數或者是脈沖函數,at(y)和at(x)是yt和xt白噪聲序列。2傳遞函數模型在土壤溶質運移中的模擬Jury(1982)等提出了模擬田間穩定流條件下保守溶質Br-運移的TFM（TromsferFunctionModel）。該模型通過研究溶質從土壤表面運移到土壤某一深度歷經時間的概率密度函數(probabilitydensityfunction,簡稱pdf),來推求溶質平均濃度與時間和土壤深度的關系,并以此來模擬溶質在土壤中的運移.在以后的研究中,Jury等(1986)、White等(1986)、Jury等(1990)又把該模型推廣運用于模擬非穩定流場中田間土壤示蹤溶質Br-的運移。TFM把土壤視為時不變系統,把溶質輸入和輸出過程與溶質分子在土壤中的運移與轉化過程的概率聯系起來,其輸出可通過卷積方程得到。它主要考慮溶質在土體“兩端”的輸入和輸出過程,利用系統辨識的方法,求出溶質分子在土體中歷經時間的概率密度函數從而對土壤溶質運移過程作出模擬和預測。一般而言,溶質進入土體的方式可概化為階躍輸入和脈沖輸入,無論對于哪種方式,傳遞函數方程的一般形式均為:(39)式中Qout(t)為溶質質量在觀測時間t通過下邊界,或吸附進入固相,或揮發進入氣相,或轉化為另一不同形式而消失等方式離開土壤的損失速率;Qin(t′)為較早時間t＇溶質質量進入土壤的速率;g[(t-t＇)′]為溶質壽命關于溶質進入土壤時間的條件概率密度函數或傳遞函數。式(39)表明:在時間t溶質質量從土壤中的累積損失速率等于累積的加權溶質質量輸入速率,該權重因子是溶質壽命的條件概率密度,溶質的壽命等于時間t和較早時間t＇之差。傳遞函數方程將溶質的輸入過程Qin(t′)與輸出過程Qout(t)聯系起來。某一特定溶質運移現象的g[(t-t′)t′]的求解問題,即為給定有關溶質輸入和輸出速率的實驗數據將式(39)作為時間的方程來求解。脈沖輸入(如降雨或灌溉歷時與對溶質運移過程的監測時間相比很短時,可視為此種輸入),其對土壤的物理擾動最輕微,因此,系統對溶質的反應完全可由隱含于溶質壽命pdf的內在的土壤過程所決定,溶質質量輸入速率在決定溶質質量損失速率方面不起主要作用。故可假設溶質從土壤表層的輸入Qin(t′)與從出流表面的輸出Qout(t)是兩個獨立的隨機過程,即假定溶質壽命的概率密度函數與輸入時間t′無關,則式(39)可簡化為:(40)又脈沖輸入速率可近似表為:(41)這里,k為歸一化常數,為Dirac函數,由其性質,可將方程(40)進一步簡化為:(42)若溶質僅僅通過一個出流面流失,而忽略通過其它方式溶質的損失,式(42)中g(t)便成為溶質歷經時間而非方程(39)所定義的溶質壽命的pdf.這樣,式(42)就可寫成:(43)式中Qex(t)為溶質的出流速率.另一方面,顯然有(44)比較(43)、(44)兩式便知:式(44)即為通過實驗數據由計算間接獲得pdf的公式。經田間和室內土柱實驗的驗證,表明運用歷經時間pdf呈如下對數正態分布的TFM能夠較好地模擬溶質在土壤中的運移:(45)其中μ和σ2分別為lnt的均值和方差。若已知g(t;λ),μ和σ2∈參數空間λ,則用極大似然法估算μ和σ2的似然函數L形如:(46)式中N為實驗取樣觀察次數,CK為不同觀測時段以式(44)表示的規一化濃度,ΔtK為觀測時段。式(46)取對數并對其中參數求導得:(47)進一步地可導出:(48)(49)在研究排水條件下溶質歷經時間的pdf時,時間中值tm和時間均值是描述pdf特性十分重要的時間變量。tm是指施于土壤表面的一個溶質分子有50%概率出現在所感興趣的的排水深度時所用的時間,其計算公式為:tm=exp(μ)(50)而定義為表施溶質的一半已經到達所感興趣的排水深度時所用的時間,計算公式如下:=exp(μ+σ2/2)(51)White(1987)定義:參與溶質運移的土壤體積含水量為:θst=q0tl-1,其中q0為平均降雨強度;t為中值或均值歷經時間;l為從土壤表面到地下排水工程的距離。上式適用于雨強較大、排水工程埋深較淺的濕潤氣候區。對于脈沖灌水條件下,土壤非穩定流場中的溶質運移而言,入滲水和原有土壤水是一個不穩定的混合過程,它們對土壤溶質運移的貢獻不同。一般用出流速率d(t,l)代替q0則更能反映參與輸運土壤溶質出流的那部分水體。于是有(52)和(53)以氮素為例，White(1987)研究了降雨條件下氮素在土壤中的淋失,給出在時段[t1,t2]內淋失量的表達式:(54)而對于給定的觀測時段[0,t],以脈沖注入的溶質通過出流口的累積淋失量(忽略上式中表示淋失量是減少意義的負號)為:(55)令聯立和兩式解得:(57)則(58)3傳遞函數模型在神經網絡中的應用目前，傳遞函數模型的構建方法主要有線性回歸、非線性回歸、神經網絡等，其中，基于前饋神經網絡理論建立的傳遞函數由于不需要先驗假設，可通過迭代運算獲得目標變量的最優函數，具有較大的優越性。下面簡要介紹了概率神經網絡（PNN）的結構模型、功能、基本學習算法。3.1概率神經網絡模型概率神經網絡(ProbabilisticNeuralNetworks，PNN)是由D.F.Specht在1990年提出的。主要思想是用貝葉斯決策規則，即錯誤分類的期望風險最小，在多維輸入空間內分離決策空間。它是一種基于統計原理的人工神經網絡，它是以Parzen窗口函數為激活函數的一種前饋網絡模型。PNN吸收了徑向基神經網絡與經典的概率密度估計原理的優點，與傳統的前饋神經網絡相比，在模式分類方面尤其具有較為顯著的優勢。3.1.1由貝葉斯決策理論：(59)其中。一般情況下，類的概率密度函數是未知的，用高斯核的Parzen估計如下：(60)其中，是屬于第類的第k個訓練樣本，是樣本向量的維數，是平滑參數，是第類的訓練樣本總數。去掉共有元素，判別函數可簡化為：(61)3.1.2PNN的結構以及各層的輸入輸出關系量如圖1所示，共由四層組成，當進行并行處理時，能有效地進行上式的計算。圖1概率神經網絡結構（以3個類別為例）如圖1-1所示，PNN網絡由四部分組成：輸入層、樣本層、求和層和競爭層。PNN的工作過程：首先將輸入向量輸入到輸入層，在輸入層中，網絡計算輸入向量與訓練樣本向量之間的差值，差值絕對值的大小代表這兩個向量之間的距離，所得的向量由輸入層輸出，該向量反映了向量間的接近程度；接著，輸入層的輸出向量送入到樣本層中，樣本層結點的數目等于訓練樣本數目的總和，，其中M是類的總數。樣本層的主要工作是：先判斷哪些類別與輸入向量有關，再將相關度高的類別集中起來，樣本層的輸出值就代表相識度；然后，將樣本層的輸出值送入到求和層，求和層的結點個數是M，每個結點對應一個類，通過求和層的競爭傳遞函數進行判決；最后，判決的結果由競爭層輸出，輸出結果中只有一個1，其余結果都是0，概率值最大的那一類輸出結果為1。3.2基本學習算法以下幾步構成了PNN神經網絡的算法：第一步：首先必須對輸入矩陣進行歸一化處理，這樣可以減小誤差，避免較小的值被較大的值“吃掉”。設原始輸入矩陣為：(62)從樣本的矩陣如式（62）中可以看出，該矩陣的學習樣本有m個，每一個樣本的特征屬性有n個。在求歸一化因子之前，必須先計算矩陣：然后計算：=(63)式中，則歸一化后的學習矩陣為C。在式（63）中，符號“”表示矩陣在做乘法運算時，相應元素之間的乘積。第二步：將歸一化好的m個樣本送入到網絡樣本層中。因為是有監督的學習算法，所以很容易就知道每個樣本屬于那種類型。假設樣本有m個，那么一共可以分為c類，并且各類樣本的數目相同，設為k，于是m=k*c。第三步：模式距離的計算，該距離是指樣本矩陣與學習矩陣中相應元素之間的距離。假設將由P個n維向量組成的矩陣稱為待識別樣本矩陣，則經歸一化后，需要待識別的輸入樣本矩陣為：(64)計算歐式距離：就是需要識別的樣本向量，樣本層中各個網絡節點的中心向量，這兩個向量相應量之間的距離：(65)第四步：樣本層徑向基函數的神經元被激活。學習樣本與待識別樣本被歸一化后，通常取標準差=0.1的高斯型函數。激活后得到初始概率矩陣：(66)第五步：假設樣本有m個，那么一共可以分為c類，并且各類樣本的數目相同，設為k，則可以在網絡的求和層求得各個樣本屬于各類的初始概率和：(67)在上式中，代表的意思是：將要被識別的樣本中，第i個樣本屬于第j類的初始概率和。第六步：計算概率，即第i個樣本屬于第j類的概率。(68)參考文獻：[1]倪雅茜,張文華,郭生練.流量過程線分割方法的分析探討[J].水文,2005,25(3):10-14,19.[2]詹道江,葉守澤.工程水文學[M].3版.北京:中國水利電力出版社,2000.[3]劉清華,史紅亮,安書全.人工神經網絡在地表徑流預報中的應用[J].三峽大學學報:自然科學版,2007,29(4):292-294.[4]顧慰祖,謝民.同位素示蹤劃分藤橋流域流量過程線的試驗研究[J].水文,1997(1):29-32,23.[5]MazvimaviD,MeijerinkAMJ,SteinA.PredictionofBaseFlowsfromBasinCharacteristics:aCaseStudyfromZimbabwe[J].HydrologicalSciencesJournal,2004,49(4):703-715.[6]黃國如.流量過程線的自動分割方法探討[J].灌溉排水學報,2007,26(1):73-78.[7]SlotoRA,CrouseMY.HYSEP:HYSEP:AComputerProgramforStreamflowHydrographSeparationandAnalysis:U.S.GeologicalSurveyWater-ResourcesInvestigationsReport96-4040[M].Lemoyne,Pennsylvania:U.S.[8]LyneVD,HollickM.StochasticTimeVariableRainfallRunoffModeling:ProceedingsofHydrologyandWaterResourcesSymposium[M].Perth,Australia:NationalCommitteeonHydrologyandWaterResources[9]InstituteofHydrology.Lowflowstudies:ResearchReport[M].Wallingford,UK.:[10]EckhardtK.HowtoConstructRecursiveDigitalFiltersforBaseflowSeparation[J].HydrologicalProcesses,2005,19:507-515.[11]NathanRJandMcMahonTA.EvaluationofAutomatedTechniquesforBaseFlowandRecessionAnalyses[J].WaterResourcesResearch,1990,26(7):1465-1473.[12]黃真理,李玉粱,李錦秀,等.三峽水庫水環境容量計算[J].水利學報,2004,3:7-14.(HUANGZhenli,LIYuliang,LIJinxiu,etal.ThecalculationofThreeGorgesreservoirwatercapacity[J].JournalofHydraulicEngineering,2004,3:7-14.(inChinese))[13]陳強,李德靖,趙興明.淺析青海省地下水資源評價中的基流分割[J].水資源與水工程學報,21(1):164-167.（CHENQiang,LIDejing,ZHAOXingming.SimplyanalysisontheevaluationofgroundwaterresourcesinbaseflowseparationinQinghaiProvince[J].JournalofWaterResourcesandWaterEngineering,21(1):164-167.(inChinese)）[14]成小松,劉宗義.淺析貴州省地下水資源評價中的基流分割[J].貴州水力發電,2006,20(5):11-14.(CHENGXiaosong,