基于廣義回歸神經網絡的全球表層海水1°×1°二氧化碳分壓數據推演

鐘國榮，李學剛，曲寶曉，王彥俊，袁華茂，宋金明

( 1. 中國科學院海洋研究所 海洋生態與環境科學重點實驗室，山東 青島 266071；2. 中國科學院大學，北京 100049；3. 青島海洋科學與技術試點國家實驗室 海洋生態與環境科學功能實驗室，山東 青島 266237；4. 中國科學院 海洋大科學研究中心，山東 青島 266071)

1 引言

當前的研究普遍認為大洋每年可以吸收約2 Pg（以碳計）左右的大氣CO2，這一結果主要是通過海?氣二氧化碳分壓差估算出來的，并且與模式的估計也相一致。但用分壓差估算出的結果仍有很大的不確定性，主要原因是二氧化碳海?氣交換速率的不確定，以及參與計算的表層海水二氧化碳分壓（pCO2）的數據較少且空間分布不均勻[1–3]。盡管pCO2實測數據相對一些其他參數比較容易獲得，可以通過基于非色散紅外法的船舶連續走航觀測獲得[4]，但獲得的實測數據相對于整個大洋來說仍然較少，這使得通過集成多年的數據構建氣候態分布成為過去比較有效的研究方法[5]，但要獲得大范圍區域pCO2的連續時間變化用一般的空間插值法僅依靠這些實測數據是遠遠不夠的。并且實測數據在時空分布上也非常不均勻，特別是早期20世紀70年代前的實測數據幾乎沒有，這極大地限制了基于pCO2演化有關的大洋碳循環研究的時空尺度。雖然美國國家航空航天局(NASA)、歐洲多國合作觀測項目（EPOCA）等一直在著手擴展海洋觀測網絡，但巨大的資金投入換來的也是十分有限的時空覆蓋范圍。在這個背景下，通過大數據技術利用僅有的少量觀測數據和一些輔助參數，構建均勻的大洋pCO2格點數據來研究全球變化成為新的突破方向。有研究者嘗試利用傳統的多元線性回歸來重建二氧化碳分壓變化，但其結果只適用于有限的特定區域[6]，甚至只適用于特定季節[7]。相比之下機器學習算法和人工神經網絡更具有優勢，可以通過實驗建立起大量參數間的實證關聯，來更準確地反映復雜的海水系統中pCO2的變化規律[8]。機器學習算法包括隨機森林算法（Random Forest Algorithm，RFRE）、支持向量機（Support Vector Machine，SVM）等，目前的研究也多局限于單一過程主導的小范圍區域，對復雜區域及全球范圍的預測則顯得比較乏力。人工神經網絡種類很多，現有研究中利用的有前反饋神經網絡（Feed forward Neural Network，FFNN）[9–10]、自組織映射神經網絡（Self-Organizing Map，SOM）[11]等，目前仍存在較大的不確定性，其標準誤差從17.6 μatm到20.2 μatm不等[12–13]。廣義回歸神經網絡（General Regression Neural Network, GRNN）是FFNN中徑向基網絡的一種變形形式，與傳統的前反饋網絡相比，GRNN是非線性擬合能力特化的形式，在各個學科和工程領域應用都更加廣泛。GRNN無需傳統的改變神經元間連接權值的訓練，只需要對一個光滑因子尋優，訓練速度比FFNN快幾十到上百倍，對數據預測的連續性也優于SOM的離散估計。為了獲得誤差更低的高時空分辨率全球表層海水pCO2數據，本文首次嘗試了將GRNN應用于表層海水pCO2格點數據的推演。

2 數據來源

研究使用的表層海水二氧化碳分壓實測數據來源于表層大洋二氧化碳地圖（Surface Ocean CO2Atlas,SOCATv2019）數據集，該數據集由超過100個成員組成的國際海洋碳研究組織組建，對實測數據進行了質量控制后公開發布。整個數據集包含約2570萬條觀測數據，時間范圍為從1957?2018年。由于受葉綠素濃度數據的時間范圍限制，我們只使用了1998?2018年的數據。數據的總數量分布如圖1所示。實測數據的分布十分不均勻，整體上北半球數據覆蓋率和數據總量都高于南半球，歐洲、日本與美國東部沿岸等少數區域數據總量超過10萬條，而印度洋、南太平洋和一些近岸區域20年間數據總量只有100到1000條左右，不到10條甚至沒有數據的區域也占不小的比例。在時間分布上不均勻的程度更加明顯，如圖2所示的最近20 a獲得的pCO2調查數據，以后10 a數據量多，數據覆蓋范圍更廣，而前10 a數據量少，覆蓋范圍也小。

圖1圖例中n為數據的數量級，代表格點位置中有10n條實測數據，灰色部分代表格點位置中無實測數據。

圖1 1998–2018年間SOCAT二氧化碳分壓數據分布情況Fig. 1 Spatial distribution of pCO2 observations in SOCAT dataset from 1998 to 2018

圖中實測數據覆蓋范圍指有實測數據的網格數占海洋區域總網格數的比例，1°×1°經緯度的分辨率下，海洋區域總網格數約為43000個。

在SOCAT數據集中給出的是二氧化碳逸度（fCO2），在使用時將其轉換成二氧化碳分壓，以便與其他研究或數據集進行對比驗證，二者間的換算關系為[14]

式中，R為氣體常數（8.314 J/（K·mol））；p為大氣壓（單位：Pa）；T為絕對溫度（單位：K）；B和δ為校正系數（單位：m3/mol），計算式為

圖2 SOCAT pCO2實測數據時間分布Fig. 2 Temporal distribution of pCO2 observations in SOCAT dataset

理論上，表層海水二氧化碳分壓主要受海水的熱力學性質、生物活動和物理過程控制。在新構建的方法中，選取了與這些過程緊密相關的溫度、鹽度和葉綠素濃度，加上與時空連續性相關的經緯度、時間等參數作為推演pCO2的輔助參數。這些輔助參數的實測數據時空覆蓋范圍決定了生產出的格點數據的時空覆蓋范圍，因此在相關性高的條件下，應優先選擇實測數據多的參數。本方法中使用的表層海水溫度與葉綠素濃度為衛星遙感數據，具有足夠大的空間范圍和足夠長時間的連續觀測。通過建立與這些參數間的非線性關系來推演表層海水二氧化碳分壓變化：

式中，Lon、Lat為經過三角函數換算的經緯度，經度為0°～360°制，以保證數據在空間上的連續性。Year、Month分別為數據對應的年和月；SST、SSS分別為表層海水的溫度（單位：℃）、鹽度；CHL為葉綠素濃度（單位：mg/m3）。使用的所有參數實測數據來源如表1所示。

表1 數據來源Table 1 Data source

3 格點數據構建方法

3.1 廣義回歸神經網絡原理

廣義回歸網絡是Specht[15]在1991年建立的一種徑向基網絡的變形形式，和徑向基網絡一樣具有良好的非線性問題處理能力，并且訓練更為方便。將訓練樣本作為后驗條件，在Parzen非參數估計的基礎上，廣義回歸網絡計算輸出時遵循最大概率原則[16]。

假設神經網絡的輸入和輸出分別為X和Y，聯合概率密度可表示為f(X,Y)，以X0代表訓練集的觀測值輸入，Y相對X的回歸為

輸入X0時，神經網絡的預測輸出為Y(X0)。給出訓練樣本數據集X0與Y0的情況下，利用Parzen非參數估計對密度函數f(X0,Y)進行估算并化簡可以得到：

式中，n為樣本總數，l為輸入變量X的維數。σ為光滑因子，等同于高斯函數中的標準差。Xi代表第i個計算樣本對應的神經網絡輸入，Yi代表第i個計算樣本對應的神經網絡輸出；X0j代表訓練樣本數據集輸入X0的第j個維度，Xij代表第i個計算樣本對應的輸入Xi的第j個維度。

式（10）即為廣義回歸神經網絡計算出的預測值，其分子為訓練集中所有樣本求出的Yi的加權和，權值等于 e?d(X0,Xi)。

從結構上看，廣義神經網絡分為4層：輸入層、隱含層、加和層和輸出層(圖3)。

圖3 廣義回歸神經網絡結構Fig. 3 Structure of general regression neural network

輸入層即負責接收樣本數據的輸入向量X，神經元數量與輸入向量X的維數l相同，以簡單的線性函數作為傳輸函數。其中維數l等于7，即輸入包含7個維度，分別為經度、緯度、年、月、海表溫度、海表鹽度和葉綠素濃度。一些研究中時間僅使用月或者僅使用年作為輔助參數，但同時使用能略微降低整體的誤差，因為增加了輸入向量的維度，而且這對計算時間影響很小。隱含層的神經元數量為馴良樣本數量，通常使用高斯函數作為基礎函數，Φi代表第i個隱含層神經元，其中心向量為對應的輸入向量Xi。加和層的神經元只有兩種，分別為分子單元和分母單元。分子單元將訓練集樣本的輸出期望作為權值，求得隱含層神經元的加權和，即式（12）中的分子部分，分母單元負責的是隱含層神經元的代數和，即式（12）中的分母部分，分子單元和分母單元的輸出在輸出層中相除即得到輸入X對應的預測輸出Y。

3.2 網絡訓練與插值

為快速檢索，原始實測數據根據時間和經緯度存儲在細胞數組中，時間分辨率為12月×21 a，空間分辨率為1°×1°經緯度，對神經網絡進行訓練時需要先將格點化的數據集轉換成向量，再將其輸入到網絡中，過程如圖4所示。

圖4 原始數據向量化過程Fig. 4 Vectorization of original data

訓練使用的樣本數據集為1998?2018年的所有數據中隨機抽取的80%，由于數據總量大，并且格點數據構建的目標最小時間分辨率為1個月，對同1個1°×1°格點里同1個月內的pCO2實測數據進行了算術平均化處理。GRNN程序通過MATLAB自帶的廣義回歸神經網絡函數工具箱實現，網絡的創建、訓練和插值計算均可以通過工具箱函數命令進行。其中訓練過程在創建的同時完成，函數語法為

式中，X0、Y0分別為訓練集的輸入和對應的期望輸出，X0是經度、緯度、年、月、溫度、鹽度、葉綠素濃度實測數據組成的向量，Y0是pCO2實測數據組成的向量；net為網絡名，在同時存在多個網絡時用于辨識；newgrnn為MATLAB自帶的廣義回歸網絡工具箱函數，用于創建并訓練網絡。spread為擴散速度，是人為設定的固定標量，默認為1.0，其值越大擬合出的曲線越平滑，但如果想更精確地接近訓練樣本的期望輸出，應該選擇較小的擴散速度值，經過多次試驗后我們擇優選取的值為1.4。廣義回歸網絡的訓練過程目的是為了求得式（9）中光滑因子σ值的最佳值，這個值很大程度地影響網絡的性能。當σ值非常大時，d(X0,Xi)趨近于0，計算出的輸出Y(X0)近似于所有訓練集樣本輸出的平均值；當σ值趨近于0時，神經網絡會出現過學習的現象，表現為給定的輸入與訓練集中某一數據相同時，計算得到的預測輸出與實測值非常接近，但給出的輸入不在訓練數據集中時，計算出的輸出與實測值偏差較大。避免這各種情況出現的方法是對輸入的各個參數量級進行調整，保證各參數變化范圍的數量級不相差過大。在我們使用的輸入參數中，除年份外的參數均在0～40間變化，因此將所有數據的年份的數量級調整到與其他參數一致：

調整后年份的變化范圍為1～21，這樣就避免了過學習現象的出現。

創建并訓練網絡后，輸入二氧化碳分壓空白區域對應的溫、鹽等參數組成的向量X，即可計算出預測的二氧化碳分壓值Y，函數語法為

式中，net和創建時的net為網絡名，可替換為其他名稱，但兩個過程的名稱必須保證一致。在這個計算過程中，輸入的向量X包括了所有待插值的格點，不需要每個格點單獨計算。最后再將輸出的二氧化碳分壓預測值向量Y還原成180°×360°大小的矩陣，存儲到細胞數組中，插值方法結束。

4 構建數據的準確性分析

由于在插值方法訓練時僅使用80%的實測數據，剩余的20%實測數據就可用于對方法進行準確性評估。這20%實測數據在訓練完成后輸入到神經網絡中，比較實測值Y0和神經網絡計算出的預測值Y的差異來評估所構建數據的準確度。通常用標準誤差（RMSE）和平均相對誤差(MRE)來描述方法的精確度。

其中標準誤差的計算公式為

式中，Yi代表第i個樣本的pCO2神經網絡預測值，Y0i代表第i個樣本的pCO2實測值，n為進行誤差評估的測試樣本總數，參與神經網絡訓練的數據不用于誤差評估。標準誤差對一次預測中的特大或特小誤差十分敏感，反映出預測值相對于實際值的偏離程度，是用于評價擬合效果的指標中最常用的。

平均相對誤差的計算公式為對每個點誤差占實測值的比重求平均：

對比廣義回歸網絡計算出的預測值Y和實測值Y0可以發現，預測結果與實測數據基本一致。以Y為x軸，Y0為y軸作圖，絕大部分數據點聚集在y=x直線處，部分偏離較遠但仍均勻地靠近直線并分布在直線兩側（圖5），回歸線也十分逼近y=x直線。全球大洋的預測值相較于實測值的平均相對誤差為2.97%，標準誤差為16.93 μatm，相關系數為0.8847。實測數據多的區域，如亞熱帶太平洋、赤道太平洋和亞熱帶大西洋，插值預測的表現最好，標準誤差為10.45～13.87 μatm，平均相對誤差為1.93%～2.66%。南太平洋數據量較少，誤差卻也很低，可能是因為數值變化范圍較其他區域小，實測pCO2值均在250～450 μatm之間，而不管哪個區域在這一區間內的預測值與實測值都很接近。表2給出了與其他方法的標準誤差對比[11–13,17–21]，在整體上，GRNN略優于FFNN與SOM，具體到特定區域范圍時，一些機器學習算法的表現可能更加精確，例如Chen等[17]使用隨機森林算法重建了墨西哥灣的pCO2變化，標準誤差僅9.10 μatm，然而僅適用于主導因素較為單一的小范圍區域。也有研究將SOM和FFNN結合在一起，利用SOM將大西洋劃分成多個區域，對每個區域訓練一個FFNN來進行插值預測[20]，但精確度并沒有顯著提升，并且由于同時使用了兩種神經網絡，應用起來更加繁瑣。近岸區域由于受到陸地徑流和人類活動等因素影響，規律十分復雜，廣義回歸網絡做出的預測表現與大洋區域相比較差，但相近于Laruelle等[21]使用SOM-FFNN法對近岸區域的預測表現。如果包含近岸區域，整體的標準誤差將上升到21.60 μatm，但其他的研究也只關注大洋區域，或者只關注近岸區域，并沒有統一研究。

圖5 廣義回歸神經網絡預測值與實測值對比Fig. 5 Comparation of GRNN predict pCO2 and in situ measurements

表2 GRNN與其他方法誤差對比Table 2 Comparation of errors between GRNN and other approaches

圖6 GRNN與實測值及其他神經網絡方法同時間點pCO2數據結果對比Fig. 6 Comparation of pCO2 distribution between in situ measurements, GRNN and other approaches

盡管存在一定的誤差，GRNN法的結果與pCO2實測值的分布在高值和低值區域的位置上基本一致（圖6a，圖6b）。與同樣使用SOCAT數據集的其他神經網絡方法的數據產品進行對比(圖6b至圖6d)，圖6c為SOM法，標準誤差為23.3 μatm，圖6d為SOM?FFNN聯用法，標準誤差為22.8 μatm，幾種方法整體的季節趨勢表現出高度的一致。在1月份，南半球中緯度海域pCO2高，南太平洋東部高于西部；北半球中緯度海域整體pCO2低，而北太平洋和北大西洋近極地區域與中緯度區域相反，pCO2高。7月中緯度海域整體分布規律與1月大致相反，南大洋pCO2高，北極區域低，這些特征都與其他方法給出的結果相一致。盡管使用的實測pCO2數據集一致，不同方法使用的輔助參數種類也不同，例如圖6中另外兩種方法中，圖6c使用的參數中沒有經度，使用了混合層深度，圖6d沒有經緯度和時間，這可能也是特定區域的pCO2分布規律上各有差異的原因，特別是印度洋等實測數據少的區域。此外，不同研究使用的輔助參數空間覆蓋范圍不同，導致構建出的pCO2空間范圍存在差異。不同類型的神經網絡本身的特性也存在差異，由于SOM給出的是離散的估計，數據的空間連續性最差，存在明顯的斑塊狀分布。SOM-FFNN雖然是連續的估計，但是數據過渡也不太平滑，銳利的邊界仍存在。而表層海水并不是相互獨立的，由于物理混合過程的影響，高分辨率的情況下相鄰網格間月平均pCO2不應該相差太大。比起其他方法，GRNN法推演出的數據平滑程度更高，不需再進行人為的二次處理來達到平滑過渡效果，有潛力應用于更高分辨率的數據構建上，如0.25°×0.25°甚至更高。

除了神經網絡法外，與Takahashi等[24]通過將數十年的實測數據校正到同一年，構建氣候態分布的方法對比，pCO2的整體分布規律也存在較高的一致性。如圖7b是Takahashi等[24]的研究中給出的校正到2005年的1月pCO2全球分布，圖7a是同時間GRNN法給出的結果，北太平洋的高值帶和低值帶、南大洋的低值區域等非常相似。盡管使用的源數據和方法本身上存在一些差異，使不同研究的結果在具體的區域分布各有不同，但整體的分布規律高度相似，結合標準誤差和平均相對誤差來看，有足夠的理由相信廣義回歸神經網絡在二氧化碳分壓的插值預測上是可靠的。

圖7 GRNN法與Takahashi等[24]氣候態pCO2數據對比（2005年1月）Fig. 7 Comparation of pCO2 distributions between GRNN output and Takahashi[24] climatological mean (January, 2005)

5 結論

基于廣義回歸神經網絡，建立了表層海水二氧化碳分壓與經緯度、時間、溫度、鹽度和葉綠素濃度間的非線性關系，并據此重建了近20年來表層海水二氧化碳分壓的全球分布，標準誤差為16.93 μatm，平均相對誤差為2.97%。與其他方法的對比證實了本插值方法的可靠性，并且廣義回歸神經網絡法的適用性更強，對近岸區域也能做出預測，其表現與只關注近岸區域的其他研究相近，網絡訓練的速度也遠高于其他神經網絡。使用廣義回歸網絡進行插值預測時，由于不需要設定擴散速度外的參數，插值結果的表現主要受實測數據本身影響。在參數選擇方面，輸入參數對神經網絡的預測表現影響很大，但增加相關性低的參數并不能提高精確度，反而會降低輸出數據的平滑度，導致分布呈現斑塊狀，并且會增加計算時間。雖然增加相關性高的參數可以顯著地提高精確度，但受該參數可獲得性的極大限制。例如本研究在pCO2參數的構建中選用的葉綠素濃度，該參數與pCO2有較高的相關性，但由于僅能獲取該參數在1998年以后的大量衛星遙感數據，也僅能用該參數重建1998年以后的pCO2數據，而由于無法獲得足夠的1998年之前的數據，就無法用本研究建立的插值方法重建1998年之前的pCO2數據。現有研究也大部分是通過使用衛星遙感數據來解決參數可獲得性的問題，但滿足條件的衛星遙感數據也只有最近幾十年的，這也是大部分現有研究都只能重建近幾十年pCO2數據的原因。而更早期的觀測數據過少，很難支撐大范圍的預測插值，因此如何重建早期的pCO2數據成為待解決的下一個難題。