基于卷積神經網絡和狀態時間序列的參數辨識

武 頻，常旭婷，郎佳林，潘凱凱，龔思泉

(1.上海大學 計算機工程與科學學院，上海 200444；2.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，綿陽 621000)

0 引 言

數值模擬技術可以替代部分風洞試驗建立相應流場的數學模型，在降低實驗成本的基礎上更便捷地預測流場結果，相應的對數值模型的可信度也有一定的要求。

數值模型中涉及了一系列參數，有的由觀測儀器間接計算或直接獲得，有的無法測量只能通過經驗推導估計。為了提高數學模型可信度使其能正確再現流場變化，有必要結合流場實測數據對參數進行辨識，目前工程上普遍采用優化算法。

優化算法通過建立目標函數表示數值模型的輸出與實測數據之間的距離，查找目標函數極小時數值模型中的參數組合，實現參數辨識。水文預報領域中陸續引入了遺傳算法[1-3]、下坡單純形法[4]、粒子群算法[5-6]；氣動力領域中廣泛應用了最大似然法[7]。由于最大似然法依賴目標函數梯度以及參數初值的選取，有一定的局限性，張天姣等[8]也將粒子群算法引入到氣動力領域中求解參數辨識問題。優化算法在求解過程中需要對目標函數或目標函數梯度進行大量迭代計算，而目標函數中又包含了流場復雜的數值模型，因此參數辨識過程十分耗時。

隨著深度學習的發展，數值模擬與神經網絡技術的結合逐漸成為一項研究熱點，尤其是在氣動力建模領域之中。Miyanawala等[9]通過卷積神經網絡（Convolutional Neural Networks，CNN）對氣動載荷實現準確預測。陳海等[10]利用CNN和翼型圖像建立氣動模型。王超等[11]基于參數辨識數據利用人工神經網絡建立氣動模型。這些研究都表明深度學習應用在氣動建模中的潛力。

鑒于上述研究背景，本文考慮用深度學習實現參數辨識，克服傳統方法計算量大耗時長這一缺陷。提出利用CNN從數值模擬結果中提取狀態時間序列和模型參數之間的映射關系，建立參數辨識模型，這種參數辨識模型可在實際應用中直接通過一段時間的實測數據對一系列數值模型參數進行快速估計。

1 Lorenz63混沌系統

1.1 數值模型

本文使用Lorenz63對提出的參數辨識方案進行具體的實驗和討論。Lorenz63是Lorenz于1963年建立的一個簡單的大氣對流數學模型[12]，是混沌系統最早的數值表示。混沌系統是指確定性系統中存在看似隨機的不規則運動。混沌的特征在于不確定性，不可重復性和不可預測性。其模式控制方程由非線性系統給出，具體如下：

其中，σ、γ分別表示Prandtl數和Rayleigh數，β表示與對流尺度相聯系的參數。

實驗中采用龍格-庫塔數值積分方法求解方程的解 (x,y,z)。x、y、z在這里就是我們所說的數值模型計算的狀態變量。通常使用一個狀態矢量來表示當前系統的全部狀態變量，即x(x,y,z)。

我們使用三維空間表示Lorenz63系統的狀態變化，每個維度對應每個狀態變量的值，見圖1。圖1是Lorenz63在參數 σ =10.0、 γ =28.0、 β =8/3下的積分結果，積分步長為0.01。參照Lorenz給的原始數據，該系統會在這種參數設置下表現為混沌現象，并且其在其他參數下還會呈現不同的運動狀態。

圖1 σ =10.0， γ =28.0， β =8/3，積分步長為0.01時Lorenz63系統狀態圖Fig.1 A Lorenz63 system with σ=10.0，γ=28.0，β=8/3，and time step 0.01

1.2 參數分析

本文討論的是混沌模式下的Lorenz63系統，根據前人的總結，首先參數 γ要符合條件 2 4.0≤γ≤28.0，且以 σ =10.0、 β =8/3固定不變為前提。如果改變這兩個控制參數，結果不必然。由于本文要實現的是多參數估計，這對于實際工程更有意義，因此只變化一個參數將無法實現任務目標，需要對 σ和 β參數的取值范圍也進行討論。

我們可以看Lorenz63在 σ =4.8、 β =8/3、γ=28.0時的結果，見圖2。

圖2 σ =4.8， γ =28.0， β =8/3，積分步長為0.01時Lorenz63系統狀態圖Fig.2 A Lorenz63 system with σ=4.8， γ =28.0，β=8/3，and time step 0.01

在該狀態的Lorenz63系統狀態變量會在剛開始的時間區間中不斷震蕩然后不斷某個點趨近，最后穩定在不動點上。從左圖中可以看到系統狀態只會在最外圈呈現出蝴蝶吸引子，之后不斷收縮到右邊的不動點上。這樣的物理狀態一是和我們需要的混沌模式不一致，二是不符合訓練數據的基本要求。我們需要的數據是一小段連續時間步的狀態變量，用于輸入到CNN中提取運動特征，如果我們輸入上述結果的[60,61]之間的狀態矢量，那么輸入的數據在時間序列上實際是定值，對于神經網絡來說，雖然輸入的原始特征是一個二維矩陣，但其實只有三個特征值，要單從這三個特征中學習到隱藏的運動特征，并預測三個參數，十分困難，再加上狀態變量x和y很接近，則等于要用神經網絡從二個特征值中提取到三個特征。因此，我們需要排除掉會導致這樣結果的參數，設置好參數的取值范圍。

通過手動實驗，我們基本可以確定好三個參數的范圍：

在上述取值范圍下，不同的參數組合還是會導致產生不動點的可能，因此我們對各種參數組合分別進行了Lorenz63數值實驗，實驗結果見表1，其中區間的選取是手動嘗試結果。

表1 Lorenz63各參數取值范圍選定結果Table 1 Selected parameter ranges of each the Lorenz63 system

2 參數辨識模型

本文的參數辨識模型可以依據多個連續時間步的實測系統狀態逆推出流場數值模型的參數。因此，模型的輸入數據是狀態時間序列，輸出數據是參數估計結果。我們用神經網絡提取輸入數據中的隱藏特征，構建參數和狀態時間序列的映射關系。

2.1 狀態時間序列輸入

首先向神經網絡輸入多個時間步的狀態矢量，每個狀態矢量含有多個狀態變量，以Lorenz63系統的3個狀態變量為例，假設輸入5個時間步的數據，后文相同，因此有如下的數據格式，見圖3。

圖3 參數逆推模型輸入數據格式Fig.3 The input format of the parameter estimation model

輸入數據是一個二維矩陣，類似一張單通道圖像。矩陣上的行向量表示一個時刻的狀態矢量，列向量表示特定狀態變量的多個連續時刻的數值。我們將行向量和列向量上的數據分別稱作橫向數據和縱向數據。

從每條縱向數據上看，一條數據表示對應狀態變量在該5個時間步下的數值。我們將數據輸入到CNN中希望網絡能夠根據變量變化的趨勢提取系統的運動特征。然而狀態變量在不同時間區間上會在不同數值范圍上變化，在不同初始場的驅動下也會有不一樣的結果，因此縱向數據需要弱化數值對神經網絡提取特征的影響，使其更關注僅僅在5個時間步下的數值變化趨勢的特征。我們選擇基于原始數據均值以及標準差的標準差標準化方法，對每個輸入數據的每條縱向數據進行單獨標準化，以X列數據為例，將原始值xi標準化到，

其中，δ為一個非常小的常數，是為了防止標準差為0而做的平滑操作，稱作平滑指數。

縱向數據的標準化處理將每列狀態變量進行單獨處理，在數值上僅由單列數據決定。因此，從單條橫向數據上看，每個位置上的狀態變量在數值上有不同的取值范圍，大小不相同。這些狀態變量在神經網絡中也被稱為特征。這些特征如果不處理直接使用原始數據輸入模型中，那么不同特征對模型訓練過程將會產生不同程度的影響，梯度下降過程會更受某些數值大的特征的引導，忽略那些數值小的特征。因此我們需要使不同的特征具有相同的尺度，這樣在訓練神經網絡的時候，不同特征對參數的影響程度就可達成一致。我們通過數據標準化以實現這個目的。同樣采用標準差標準化方法，對橫向數據進行由到變換：

最后將縱向和橫向的結果相加，作為我們雙向數據標準化的最終結果，輸入到神經網絡中去。

2.2 CNN結構

我們的參數辨識模型輸入的是經過上述雙向標準化后的時間序列數據，數據格式見圖3，輸出的是參數估計的結果，格式為一維向量。由于輸入數據格式類似于單通道圖像，每個單元位置可以看作圖像中的像素點，因此選擇常被用作圖像特征提取的CNN[13-14]。典型的CNN在輸入層后由多個交替執行的卷積層、池化層和全連接層構成。故針對本文任務的模型圖如圖4所示。

圖4 CNN結構Fig.4 A CNN structure

1）卷積層

使用多個卷積層從原始特征中提取隱藏特征。由于系統狀態變量之間不像圖像像素，不單單是相近的像素點之間彼此有聯系，頭尾狀態變量之間也由于物理模型的控制方程而存在著相互影響的制約關系，因此進行卷積運算時需要將對應橫向數據上的所有狀態變量包含在內。我們使用“一維卷積”（卷積完后的數據呈一維）進行卷積操作，沿著時間步移動來提取系統在當前時間區內的運動特征，見圖5。

圖5 “一維卷積核”運算示意圖Fig.5 A sketch of the one-dimensional convolution kernel operation

“一維卷積核”設置的維度大小和數據橫向上的大小相同。其中，輸入特征中圈出來的區域為局部感受野。因此“一維卷積”只需要在原始輸入的多通道特征中垂直移動，提取更深層次的特征。

以第一個卷積層為例，輸入特征只有一個通道，因此使用大小為3的卷積核執行卷積運算，通過垂直滑動獲得3×1的特征圖，第一個位置的特征值為：

其中，p為當前局部感受野的位置，b為當前卷積層的偏置，ωij為卷積核 (i,j)位置上的權值。當前卷積層有多個卷積核，每個卷積核的權值不共享，但偏置是共享的，即多個卷積核分別進行卷積運算后都加上當前卷積層的偏置，從而生成多個特征圖，輸入到下一層中。每個特征圖在下一層中表示圖的一個通道，在進行卷積運算時，需要將每個通道上的結果相加，最后加上共同的偏差b。

2）池化層

池化層通過最大池化或平均池化等來減小輸出維的大小，其主要目的是通過壓縮卷積層的輸出來簡化網絡的計算復雜性。每個卷積層之后可以跟隨一個池化層，根據數據特征的維數大小適度添加。

3）全連接層

我們的任務是估計參數，即要輸出一維向量，向量中每個數表示對應參數的值。因此我們需要將卷積操作輸出的二維或三維數據展平，該操作由全連接層進行。直觀來說，就是講原本的三維矩陣展平重新排列變成一個全連接層的各個神經元。

4）輸出層

在多個全連接層的映射之后輸出參數預測結果。由于數值模型的每個參數有不同的物理意義和單位，由真實物理環境決定。因此，為了使得目標輸出數據（也可稱為標簽）符合正態分布，在訓練神經網絡時，我們不僅要對輸入數據進行雙向標準化，標簽也要進行標準化，采用同樣的標準差標準化方式。假設為了得到這個參數辨識模型，生成了a組參數對應的數值模型解，數值模型解用于構造輸入數據，參數值用于構造標簽，以Lorenz63為例對參數 σ 進行如下標準化操作：

γ，β操作相同。

使用標準化后的數據作為標簽訓練網絡，因此輸出層的預測結果 σo需要進行還原，使用公式（6）中的均值和標準差s，即：

構建好網絡結構后采用反向傳播方式對網絡進行訓練。

3 方法驗證

3.1 CNN模型訓練

Lorenz63的CNN參數辨識模型訓練過程由圖6給出。

圖6 參數辨識模型訓練過程Fig.6 The process of training a parameter estimation model

在進行數值實驗時，多組實驗的初始狀態矢量均為(-10,10,20)，積分步長取0.01。每組實驗中的參數是在設置好取值范圍的前提下隨機選擇的。從參數分析可以得到參數組合的取值范圍實際分為兩組，如表2所示。我們在第一組取值范圍下分別對三個參數進行隨機取值，生成1 000組參數組合，并在第二組取值范圍下同樣隨機取值了500個參數組合，每組實驗積分區間為[0，100]。因此一共生成了15 000 000條狀態矢量數據。

表2 訓練數據的參數選取Table 2 Selected parameters for data training

將數值實驗生成的數據處理成模型需要的格式。比如一次實驗中，設定初始特征圖的維度為5×3，采用長度為5的滑動窗口將每次窗口內的5個時間步的狀態變量合成一個二維矩陣，滑動步長為1。格式處理好后進行標準化操作，即準備好了模型訓練需要的樣本數據，進行網絡的訓練。

實驗是在Nvidia Tesla P100 GPU上進行的，使用Keras框架搭建神經網絡。本文的CNN模型包含三個卷積層，每個卷積層分別設置了大小為3的一維卷積核，padding 方式選擇了‘valid’（不在外圈補 0）；由于Lorenz63系統簡單，維數不高，因此當輸入時間步較少時可不添加池化層進行降維操作，如果輸入時間步較大（如 100,1 000···），則適當添加池化層；之后，跟隨著3個全連接層，神經元的激活函數選擇了‘Relu’[15]，

損失函數為平均絕對誤差（Mean Absolute Error，MAE），

其中，Po、Pta分別表示神經網絡輸出結果和目標結果，、的下標i為1、2、3時表示參數 σ 、 γ、 β 。優化器選擇了‘Adam’[16]，學習率為 0.001。

我們將樣本數據隨機打亂，以避免時間序列的影響。取80%的數據作為訓練集，余下的20%為測試集；batch的大小為512，表示每次取512條數據來計算損失函數的平均值，進行梯度下降的訓練；epoch為100，表示最多進行100次重復訓練，其中當損失函數降到最低時停止訓練。

3.2 實驗分析

使用神經網絡對系統參數進行估計的主要目的是減少參數辨識的計算量，縮短計算時間，但前提是要保證參數預測結果的準確度，必須在驗證這種辨識方法的準確可行之后才能進一步討論其對于效率的提升。

因此，實驗結果主要關注準確度和時間兩個方面。同時，為了更好的討論本文方法的優劣，在進行上述實驗的同時，還分別做了一些比較實驗。

1）輸入不同數量時間步的結果

在實驗中我們嘗試改變輸入的時間序列數據的長度，查看了向CNN輸入不同數量特征的結果，見表3。訓練的超參數均與上一節實驗設計的描述相同。

表3 輸入不同長度的時間序列數據的實驗結果比較Table 3 The mean absolute error with different input lengths

表中 M AE直觀地表示了在測試集上模型預測結果與目標結果的整體的平均絕對誤差，我們可以看到這個誤差在0.2之內，且在輸入100、200、500個時間步時僅為0.04。

首先從輸入時間步為10、50、100的結果來看，可以得到一個初步結論：當輸入時間步數量越多時，模型精度越高。這是由于向模型輸入更多時間步的數據，等于使模型能得到更多原始特征，因此我們能認為這種結果是自然現象。

然而當輸入時間步大于100之后，平均誤差并沒有如預期繼續降低，保持在0.04不變。由此我們可以認為神經網絡能從100個時間步的狀態變量中提取到Lorenz63系統的運動特征，再向其輸入更多的狀態變量是沒有必要的。

同時還需要注意一點，由于系統所處的物理場景經常是瞬息萬變的，我們需要參數逆推模型能根據當前的系統狀態變化來估計參數，如果輸入過長時間區間的系統狀態，則會導致參數估計不準確。因此不認為輸入過多時間步的數量是好的選擇。

結合以上，基于CNN-時間序列的參數辨識方法，在針對具體任務時需要綜合考量以確定輸入數據格式。針對Lorenz63實驗，輸入100時間步是較好的選擇。

2）實驗預測結果

用圖7抽樣展示實驗結果，更直觀的體現模型預測的準確度。圖中[0,200]區間的參數為 σ =10.0，γ=28.0，β = 8/3，[200,400]的參數為 σ =14.0，γ = 24.0，β=1.5。在這兩組參數下另外進行Lorenz63實驗得到了兩段狀態變量時間序列，均未參與模型的訓練，選擇的是[0,200]區間的數據輸入到模型進行參數估計，估計結果為圖7中虛線。我們可以看到這個結果近似真實值。其中用到的模型為輸入時間步為100的模型（稱為CNN-100）。

圖7 CNN-100的參數估計結果圖Fig.7 Estimated parameters by the CNN-100

由于實際問題中觀測中存在噪聲，因此向這兩段狀態變量數據添加標準差為對應量值的1%的噪聲擾動，再次查看結果，見圖8。結果表明，模型依舊能近似估計真實參數。

圖8 添加噪聲的CNN-100參數估計結果圖Fig.8 Estimated parameters by the CNN-100 with noise.

3）同類實驗比較

有研究工作者也在進行通過神經網絡估計Lorenz 63模型參數的實驗（https://github.com/Yajing-Zhao/Lorenz_MLP）。

本課題進行了同類實驗的對比分析。

該實驗使用到的神經網絡為多層感知機（Multi-layer Perceptron，MLP），輸入數據為展平成一維向量的狀態時間序列。本文對MLP的實驗方案和模型進行了改進：采用我們CNN實驗相同的數據，輸入50/100個時間步的狀態變量，輸入層的神經元數量為150/300，添加了三個隱藏層，神經元數量分別為128/256、64/128、16/64，激活函數選擇‘Relu’。結果見表4。

表4 同類實驗效果比較Table 4 The comparison between the MLP and the CNN

4）與PSO參數辨識方法的時間比較

本文選用粒子群算法[17]也進行了實驗，展示了基于神經網絡的參數辨識在時間上帶來的優勢。

我們設置了粒子群法的慣性權值因子為0.8，學習因子c1、c2為2，群體大小為40，在合適范圍內隨機初始化，最后輸入上一時間步的系統狀態變量（數值模型輸出）以及當前的實測值（添加了1%噪聲擾動的狀態變量），對目標函數進行全局最優化，最大迭代次數為200，并同時計算了10個結果取平均。

選用CNN-100與粒子群算法的結果進行比較，兩個程序均運行在Inter(R) Core(TM) i5-6200U CPU上，結果見表5。

表5 CNN-100與粒子群算法實驗結果的比較Table 5 The comparison between the CNN-100 and the Particle Swarm Optimization

從準確度上看，兩種方法都能估計出系統參數的近似解，傳統的粒子群算法并沒有更勝一籌。

從時間上看，我們的模型在通過逆推參數實現參數辨識的過程中只耗費了0.09 s（包括了對輸入數據進行雙向標準化并對輸出數據進行還原的時間），比粒子群算法2.23 s少了96%。從該點上看，使用神經網絡的參數辨識方法能夠大大的縮短參數辨識過程的時間，從而使參數辨識工作能夠更快的完成。

另外，神經網絡模型在進行參數辨識之前需要耗費一定時間進行數據生成和網絡學習，不過這些工作并不發生在參數辨識階段，不會耽誤原本工程的開展；且一旦構建好模型，便可以反復使用，無需在每次參數辨識時都進行神經網絡訓練。因此從整體上看，本文方法耗費的時間依舊較少。

4 結 論

本文基于卷積神經網絡對數值模擬的參數進行估計，結合雙向數據標準化的方式提取系統變量在時間上的運動特征。通過Lorenz63混沌系統的驗證，本文的參數辨識模型在具有高精度的同時，比粒子群參數辨識法的計算時間降低了96%。

總的來說，本文的研究可以節省參數辨識過程耗費的大量計算時間，實現快速估計參數的目標。