靳 錚，張雪芹

(1：中國科學院地理科學與資源究所,中國科學院陸地表層格局與模擬重點實驗室，北京 100101) (2：中國科學院大學，北京 100049)

渦動相關(eddy covariance, EC)方法在國內外湖泊水熱和成分通量觀測研究中均有廣泛應用[1-6]. 太湖在2013年就建成了由6個通量站組成的中尺度渦動相關通量觀測網[3]，并利用觀測數據分析和模擬了該地區(qū)湖氣間能量、水熱相互作用. 2014年鄱陽湖的渦動相關研究[5]為大型淺水湖湖面能量通量劃分了不同的水文過渡期格局. 洱海[7]和納木錯[8]的渦動相關觀測研究初步揭示了高寒湖泊的水面湍流通量特征. 受制于儀器性能和觀測下墊面的異質性，EC觀測的通量數據時間序列會產生缺測點[9]. 缺測點產生的原因主要包含信號異常、儀器損壞以及湍流平穩(wěn)性統(tǒng)計未能滿足通量計算假設等[10]. 特別是高寒湖泊EC觀測面臨著封凍期維護、強風和低溫等問題. 因此，EC系統(tǒng)一般采取岸基或近岸安裝[11-13]. 這使得EC系統(tǒng)所測來自陸面源區(qū)的通量數據和湖面源區(qū)的混雜一起. 通量源區(qū)劃分并提取目標源區(qū)通量數據后，EC觀測的有效數據比例顯著下降. 湖面源區(qū)通量數據有效數據時間分布存在顯著的不均勻現(xiàn)象，從而嚴重影響EC通量數據的可分析性，因此有必要對EC通量數據進行插補. 湖面水熱通量觀測研究中，學者們把溫度梯度、風速、水汽壓差等作為湍流傳輸過程的基本物理背景信息來估算通量強度(如Bulk Aerodynamic Transfer Model[13]和DYRESM[14]). 這類種估算方式除了考慮動力、熱力機制之外，還需考慮物理量綱的對應，以便將估算模型以等式的形式表達出來. 但是湍流輸送過程的時空尺度決定了這類梯度量綱模型不足以描述湍流傳輸在其宏觀特征背后的微觀過程. 除此之外，通量數據的插補還可以使用查表法[15]、動態(tài)線性回歸方法[16]和非線性回歸方法[17]等基于統(tǒng)計原理的方法. 這類方案是根據觀測經驗，以溫度梯度、水汽壓等關鍵背景物理信息為依據，通過擬合計算得到對通量強度的最優(yōu)統(tǒng)計估計來完成插補的. 相對于上述方案而言，ANN(artificial neural network，人工神經網絡)方法可以更加充分地利用羊卓雍錯觀測實驗中較豐富的輔助觀測系統(tǒng)和較長觀測時間帶來的數據量優(yōu)勢. 由于羊卓雍錯的特殊近地層湍流背景[18]，湍流輸送過程在此處與其背景動力、熱力強迫過程的關系更加復雜，非線性特征更顯著. 故本文在溫度梯度、風速、水汽壓差等核心要素的基礎上，將更多的環(huán)境氣象要素同時引入到湍流傳輸物理背景信息的描述中，并使用ANN模擬當前信息組合與通量強度值的映射關系.

盡管以ANN為基本工具的機器學習方法廣泛應用于地球科學研究中[19]，但其在湖泊通量觀測的數據插補方面仍有待加強. 洱海EC通量觀測實驗[7]使用了一個神經元總數20個以內的ANN模型來模擬和插補通量數據，其ANN模型輸入變量包含風速、水面溫度、氣溫、飽和水汽壓差4項. 但該研究對ANN模型估算通量強度的準確性未作進一步討論. 由于ANN模型的性能取決于數據所包含的有效信息量和噪聲水平[20]，故其模擬性能的檢驗是必要的. 我們從信息利用的角度設計了一種針對EC湍流通量數據的ANN模擬插補及驗證方案. 通量模擬輸入信息中增加變量時，信噪比較低的變量會影響ANN模型的性能，同時大幅增加ANN訓練時間. 為了給輸入ANN模型的大量有效信息和噪聲組合提供充足的概率樣本空間定義域，本文采用了超寬ANN結構[21]. 基于超寬結構的ANN模型神經元總數超過2000個，本研究利用并行計算技術解決了模型訓練的時效問題. 本文對上述超寬結構ANN模型的通量估算性能進行了交叉驗證和誤差分析，提供了客觀且全面的驗證結果和誤差分析結論，并討論了EC湍流通量數據插補問題中ANN模型的泛化，為高寒湖泊湍流通量觀測研究的數據插補提供了一種參考思路.

1 數據與方法

1.1 羊卓雍錯渦動相關觀測

羊卓雍錯湖湖面面積約為643 km2[22]，是喜馬拉雅山脈北麓最大的內陸封閉湖泊. EC設備及輔助觀測系統(tǒng)的設置地點為羊卓雍錯湖靠近浪卡子縣白地水文站的近岸淺灘(29°07′28″N, 90°26′27″E)，海拔高度4420.63 m. 為避免被湖泊封凍期前后的冰凌沖擊摧毀，通量設備于每年的1月上旬進行拆卸，3月下旬重新安裝. 本研究采用通量系統(tǒng)2016年和2017年非封凍期同一觀測時段(4月3日0：00-12月31日23：30)的水熱通量、湍流、氣溫、濕度、水面溫度和水面輻射等同步觀測數據. 此外，2016年3月30日-4月2日以及2018年1月1-4日的數據將用于測試ANN模型的時間擴展性能，不參與ANN模型的訓練及交叉驗證. EC設備信息與數據采集的詳細情況見參考文獻[20].

1.2 神經網絡建模與檢驗方法

1.2.1 同步觀測變量的特征工程 特征工程是對ANN輸入數據進行標準化處理的過程[23]，原始同步觀測數據經過特征工程轉換為可直接輸入ANN的特征向量. 特征向量的處理方式對ANN的訓練時效和模擬性能均有重要影響. 本文使用ANN方法進行通量強度與同步觀測特征間的映射關系擬合，輸入ANN的同步觀測特征變量有12個：氣溫、水面溫度、平均水平風速、平均湍流動能、莫寧-奧布霍夫穩(wěn)定度參數[24]、氣壓、飽和水汽壓、水汽壓和四分量長、短波輻射. 上述加入輸入數據樣本的特征變量除了與渦動通量強度觀測值同一時刻的30 min平均值外，還包含它們前后各30 min的平均值，故ANN的輸入特征向量維度為36. 同步觀測特征變量的標準化處理采用平移縮放法，即根據各個變量的強度概率密度分布剔除異常值后平移并縮放至[0, 1]區(qū)間. 訓練樣本由處理完成后的同步觀測變量特征向量及觀測通量強度序列組成. 由于ANN模型具備較強的非線性映射擬合性能[25]，本研究特征變量選取時并不考慮所選變量和通量強度的統(tǒng)計相關性，而只需要特征變量與湍流輸送過程的熱力、動力過程存在理論關聯(lián)即可.

1.2.2 建模框架與參數調試 本文使用Google?的開源機器學習框架TensorFlow[26]進行ANN模型參數化構建及訓練，并使用Keras[27]開源Python模塊對TensorFlow的功能框架進行調用. 為了保證映射關系擬合計算的時效性，本文采用CUDA?(Nvidia?, Inc.)并行計算方案，硬件計算單元型號為GTX-1080?，峰值性能為每秒9萬億次單精度浮點運算(9 TFLOPS). ANN的初始化模型采用密集連接結構，不加入偏差項神經元，損失函數為訓練樣本的平均絕對誤差，擬合方式為隨機梯度下降算法[28]，并利用隨機參數剔除[29]方案弱化過擬合效應. ANN的激活函數采用PReLU[30](Parametric Rectified Linear Unit)和Tanh[31](Hyperbolic Tangent)，Tanh函數用于激活ANN的第一個隱藏層，PReLU函數則用于激活ANN的輸入層和連接輸出層的隱藏層. 已有研究表明隱藏層寬度于一定的范圍內增加時，ANN可以對非線性過程作出近似與牛頓迭代法的線性求解[21]. 所以，本文參數調試采用隱藏層寬度指數遞增搜索的策略，從初始化的16個隱藏層神經元開始，逐步增加隱藏層神經元至2n個，確定了合適的隱藏層數與神經元數后，再對學習率、剔除率和訓練次數等參數進行調優(yōu).

表1 ANN參數設置的搜索結果

*平均學習率和平均迭代次數為10批次擾動訓練的平均值.

1000 h(計算單元滿載運行時間)左右的參數調試結果表明：淺層超寬密集連接結構的ANN對渦動通量具有較好的模擬效果. 表1為Keras機器學習模塊下ANN設置參數的搜索結果. ANN的兩個隱藏層神經元數量分別達到2048和1024個，這種超寬結構可以為不同質量、不同量綱的輸入物理變量組合提供較大的概率空間，從而在迭代過程中搜索更多的映射關系，實現(xiàn)輸入信息的充分利用. 由于感熱通量、潛熱通量和水汽通量所映射的同步觀測特征是相同的，而通量強度數據的噪聲不同. 故上述3種通量的擬合采用相同的ANN結構，只是利用不同的學習率和迭代次數來適應它們的噪聲水平. 搜索計算所得的ANN模型雖然只有兩個隱藏層，但兩個隱藏層寬度均達到212個神經元后ANN的權重矩陣包含的權重數量已經超過1.6×107個，一次模型擬合所需訓練時間超過12 h. 3109個神經元的ANN模型在對應的訓練迭代次數(表1)上已經接近收斂，每100次迭代的梯度波動率小于0.01%. 考慮到實驗的時效性及模型的收斂情況，我們未繼續(xù)進行更大規(guī)模的ANN擬合計算.

1.2.3 折疊交叉驗證 本文采用折疊交叉驗證方法[32]對機器學習算法的性能進行了評估，并選取了10次折疊的交叉驗證方案. 首先，我們按通量數據的強度概率密度分布將訓練樣本分割為10個折疊子樣本，使各折疊子樣本的通量數據均具有和原始樣本相近的強度概率密度分布和平均值，并且在觀測時段內均勻排列. ANN模型擬合計算時使用其中9個子樣本合并進行訓練，保留1個子樣本用于驗證模擬效果，然后交換用于驗證的子樣本. 在對10個子樣本分別進行效果檢驗后，統(tǒng)計其平均結果和波動情況，以得到模型的誤差期望和穩(wěn)定性估計.

在搜索計算試驗中，雖然10組訓練子樣本數據所驅動的ANN結構及擬合參數是相同的，但是由于ANN擬合計算采用了隨機梯度下降方法，擬合后所得的模型權重參數在每一組交叉驗證模型中存在微小差異，所以有必要分析各交叉樣本模型性能的穩(wěn)定性. 基于湍流輸送過程的高度非線性[33]和ANN模擬輸出可能產生的正負誤差對稱性[34]，本文基于10組經過驗證的ANN模型，額外增加了10批次的擬合計算試驗(共100個ANN擬合擾動模型). 我們對每一批次ANN模型的擬合參數進行了微調擾動，訓練了共計110個ANN模型作為備選集合成員. 而后以平均絕對誤差的方差最低為標準，我們選取了10個成員對通量強度進行集合模擬取平均，以進一步降低最終模擬結果的不確定性.

2 羊卓雍錯通量數據插補及驗證

2.1 湖面水熱通量數據插補

表2 插補前后湖面通量的時間覆蓋率

插補前后通量強度數據的時間覆蓋率差異顯著(表2). 湖面通量與陸面通量源區(qū)劃分以及低質量等級數據剔除后，觀測期間通量數據時間覆蓋率均降至40%以下. 由于儀器故障，2016年的潛熱通量和水汽通量數據從10月中下旬開始大量出現(xiàn)異常狀況(如水汽瞬時數濃度測定值為負). 為減緩儀器故障對數據可靠性和整體噪聲造成的影響，本文剔除了2016年9月30日以后的全部潛熱通量和水汽通量數據. 故2016年的潛熱通量與水汽通量時間覆蓋率(26.4%)相對感熱通量而言更低，同時也低于它們2017年同期的水平. 通量數據插補后仍有少量缺測(比例低于2%)，這與輸入ANN模型的同步觀測特征值缺測或異常有關.

2.2 折疊交叉驗證

10組折疊交叉驗證的結果(表3)顯示：ANN模型的通量模擬性能十分穩(wěn)定，各個交叉驗證分組之間的誤差波動較小. 對稱性平均絕對百分比誤差(SMAPE)[35]作為驗證模型性能的關鍵指標，感熱通量10組交叉驗證中的平均值為31.871%,潛熱通量與水汽通量的交叉驗證中的平均值分別為20.8%和20.7%，這表明ANN模型對潛熱通量和水汽通量的模擬效果優(yōu)于感熱通量. 這種差別產生的原因是感熱通量觀測數據的噪聲水平高于潛熱通量和水汽通量. 平均絕對誤差(MAE)在3種通量的各10個交叉驗證組中的平均值為5.4 W/m2、15.7 W/m2和0.35 mmol/(s·m2).MAE與通量觀測平均值的百分比可作為通量模擬的期望誤差. 在感熱通量、潛熱通量和水汽通量的觀測平均值分別為18.8 W/m2、81.5 W/m2和1.84 mmol/(s·m2)的情況下分別為28.7%、19.3%和19.2%. 半小時分辨率下，潛熱通量和水汽通量的模擬期望誤差相對較小. 對于通量強度的觀測期平均值而言，感熱通量、潛熱通量和水汽通量的期望誤差分別為2.0%、1.3%和1.8%. 這由它們在10個驗證組中各組的觀測值平均與模擬值平均計算得出. 通量模擬的整體平均值誤差期望遠小于半小時分辨率單個數值的誤差期望，原因是通量模擬誤差的正負對稱性較好(圖1)，參與平均計算的模擬值越多，整體平均值模擬誤差就于0的附近越穩(wěn)定.

10次交叉折疊驗證表明，MAE差別最大的驗證組之間回歸分布狀況卻十分接近，模型性能在10個交叉驗證組間的波動很小(圖2). 這反映了10組折疊交叉驗證的樣本分割具有良好的均一性. 模擬誤差在不同的通量強度下的分布不均勻，感熱通量對于0～15 W/m2之間時模擬誤差明顯更小，潛熱通量和水汽通量于數值較大時模擬誤差將增大. 這種現(xiàn)象出現(xiàn)的原因是超寬ANN對數據量有較高的敏感性，通量強度數據的概率分布在一定程度上符合正態(tài)分布，高值和低值數據所占比例相對于平均值附近的數據更少. 由對稱絕對百分比誤差(SAPE)的分位數分布(圖3)可見，SAPE的平均值于各個通量的10組交叉驗證中均大于其中位數，表明多數情況下模擬誤差小于模型誤差期望.

表3 ANN模型的10組折疊交叉驗證結果*

*R為Pearson相關系數[36]，SMAPE為對稱性平均誤差百分比，MAE為平均絕對誤差，AVG-Obv為觀測值平均值，AVG-Est 為模擬值平均值，F(xiàn)1～F10為10個折疊組的編號.

圖1 ANN模型模擬誤差的概率密度分布Fig.1 Probability distributions of bias in ANN estimated fluxes

圖2 10次折疊交叉驗證中平均絕對誤差最大和最小的驗證組的回歸分析結果對比(回歸結果均通過了99.9%的置信區(qū)間檢驗，F(xiàn)1、F3、F6、F7和F10為折疊組的編號)Fig.2 Comparisons between the minimum and maximum mean absolute error of regression analysis results in 10-fold cross-validation groups

圖3 10個交叉驗證組的對稱性絕對百分比誤差分位數箱線圖(F1～F10為10個折疊組的編號)Fig.3 Quantile boxes of symmetric absolute percentage error in 10-fold cross-validation groups

本文在計劃修復的觀測期(4月3日0：00-12月31日23：30)之外截取了兩段未參與ANN模型擬合的數據，對比了ANN模型的模擬結果和觀測值(圖4)，結果顯示ANN模型對潛熱通量的模擬于30～120 W/m2的范圍內模擬效果極好，超出該范圍則誤差顯著增加，但變化趨勢的一致性仍然較高，這與上文中ANN模型的折疊交叉驗證結果相符.

圖4 ANN模擬的潛熱通量與觀測值的對比Fig.4 Comparison between ANN estimated and observed latent heat flux

2.3 水熱通量特征分析

2016年與2017年湖面通量數據的缺失狀況較為一致(圖5). 整個觀測期內每天0:00-9:00是湖面通量數據時間覆蓋率最低的時段. 通量源區(qū)定位取決于風速和風向，這表明該時間段內的主風向為陸地指向湖面. 與此相反，9:00-18:00時段湖面通量的時間覆蓋率則遠大于觀測期平均水平. 上述現(xiàn)象表明觀測地點存在顯著的湖陸風循環(huán).

經ANN模型插補，湖面通量于觀測期間的變化特征得以清晰而完整地呈現(xiàn)(圖5). 感熱通量2016年與2017年的狀況總體相似. 4-6月有著穩(wěn)定而顯著的日變化，每天的峰值多于10:00-16:00間出現(xiàn)，其中5月全段和6月上半月的日變化幅度較大. 從7月下旬起日變化減弱，但每天的整體強度顯著增大. 10月下旬-11月底出現(xiàn)了持續(xù)一個月的高峰，且高峰期整體強度上2017年大于2016年. 高峰期過后迅速減小至7月前的水平.

圖5 湖面通量插補前后對比：(a)感熱通量，色標分辨率為2 W/m2；(b)潛熱通量，色標分辨率為5 W/m2；(c)水汽通量，色標分辨率為0.1 mmol/(s·m2)(時間范圍為2016年和2017年的4月3日-12月31日)Fig.5 Comparison between the unpatched and patched lake surface fluxes: (a) sensible heat, with a palette resolution of 2 W/m2; (b) latent heat, with a palette resolution of 5 W/m2; (c) molar water vapor, with a palette resolution of 0.1 mmol/(s·m2) (all included in the period of Apr. 3 and Dec. 31 in 2016 and 2017)

潛熱通量最為顯著的特征是日變化，每天的峰值于12:00-18:00間出現(xiàn). 2016年和2017年的10-12月有著近兩個月的高峰期. 與感熱通量相反，2016年的潛熱通量高峰期整體強度大于2017年. 由于觀測期間氣溫變化幅度小于20℃，此幅度下溫度變化對汽化潛熱的影響低于3%，故水汽通量和潛熱通量的強度變化幾乎一致. 從兩年的觀測其整體狀況來看，湖泊的能量釋放于2016年4月初-10月末是一個漸強的過程，10月末達到頂峰. 而2017年的能力釋放頂峰延后至11月下旬.

根據算法插補后的水汽通量數據及ANN的模擬誤差期望，羊卓雍錯湖面蒸發(fā)量的EC觀測結果于4-12月為740±9 mm(2016年)和703±7 mm(2017年). 觀測期湖面的感熱通量平均值為14.9±0.2 W/m2(2016年)和14.3±0.2 W/m2(2017年)，潛熱通量平均值為78.9±1.5 W/m2(2016年)和74.4±1.5 W/m2(2017年)，水面凈輻射的平均值為126.7 W/m2(2016年)和139.5 W/m2(2017年). 由此可知湖面在2016年和2017年觀測期間的能量平衡狀況均為凈流入，凈流入能量為7.737×108J/m2(2016年)和1.198×109J/m2(2017年). 在不考慮湖水熱力層結動態(tài)機制的情況下，凈流入的熱量可使湖面下方20 m深的湖水平均每月上升約1～1.5℃，在水體的垂直動態(tài)熱平衡下，不同深度的水溫變化可能各不相同.

3 結論與討論

針對高寒湖泊渦動通量觀測有效數據比例偏低的問題，本研究利用了近年來機器學習研究中GPU并行計算和算法優(yōu)化方面的發(fā)展成果，通過基于信息利用原則的超寬ANN模型，有效優(yōu)化了EC通量數據的連續(xù)性，并利用10次折疊交叉驗證方法檢驗了ANN模擬的效果. 主要結論如下：

在ANN模型插補通量數據中有效地利用了同步觀測特征信息與通量強度的熱力、動力關聯(lián)性. 從數據插補效果看，該方案揭示了湍流通量回歸中跨量綱、跨維度映射擬合的可行性和有效性. 羊卓雍錯湖面感熱通量、潛熱通量和水汽通量的有效觀測平均值分別為18.8 W/m2、81.5 W/m2和1.84 mmol/(s·m2). 根據TensorFlow機器學習框架下超寬結構ANN模型交叉驗證結果，模擬的AME分別僅有5.4 W/m2、15.7 W/m2和0.35 mmol/(s·m2). 10個交叉驗證分組之間的誤差波動幅度分別不超過1 W/m2、2 W/m2和0.05 mmol/(s·m2). 交叉驗證組中全體觀測值平均與模擬值平均的期望誤差分別為2.0%、1.3%和1.8%. 這表明超寬ANN模型的通量模擬性能十分穩(wěn)定，在各個交叉驗證組中波動較小，且模擬誤差具有較好的正負對稱性.

高寒湖泊環(huán)境帶來了湍流通量觀測的缺測問題，通量數據插值是觀測研究的必要部分. 湍流通量過程具有高復雜度、高非線性的特征，觀測數據量大、種類多. ANN模型對輸入數據的原始形態(tài)沒有任何限制，無論所測同步觀測數據的時間同步性、空間尺度、量綱、精度、采樣率、噪聲水平存在怎樣的差別，都可以通過特征工程的標準化處理輸入ANN. 超寬結構和GPU并行計算同時保證了ANN模型模擬通量時輸入數據的充足和模型訓練的時效，是ANN模型在通量數據插補問題上泛化的必要條件. 本研究對ANN模擬通量的驗證分析表明，超寬ANN結構飽和輸入同步觀測要素的方式利用了更多有效信息. 本文將并行計算技術作為輔助工具引入通量觀測研究之中，這顯然是湍流通量觀測研究所需進一步發(fā)展的方向. 基于大數據理論，通量模擬在數據量更大的區(qū)間效果更好，表明隨著數據量增加通量模擬的效果還具備提升潛力. 本文驗證ANN模型性能時發(fā)現(xiàn)，同步觀測特征與通量強度通過超寬結構ANN模型完成了較好的映射關系擬合，這實現(xiàn)了同步觀測特征變量自相關信息的有效利用. 因此，在今后的EC湍流通量觀測實驗中，通過增加同步觀測要素，可以利用此插補方法改善高寒湖泊環(huán)境所帶來的通量數據質量問題. 本文提出基于同步觀測信息充分利用的數據插補思路為高寒湖泊等特殊環(huán)境下的EC通量觀測實驗提供了提升數據質量的參考.