基于海面位勢(shì)高度和海表溫度構(gòu)建西太平洋三維溫度場(chǎng)?

范秀梅， 樊 偉， 楊勝龍， 周為峰, 伍玉梅

(中國(guó)水產(chǎn)科學(xué)研究院漁業(yè)資源遙感信息技術(shù)重點(diǎn)開放實(shí)驗(yàn)室,上海200090)

海洋衛(wèi)星提供了連續(xù)、穩(wěn)定和豐富的海面數(shù)據(jù)，如SST和SLA，例如，星載的中分辨率成像光譜輻射計(jì)MODIS (Moderate Resolution Imaging Spectroradiometer)每2天觀測(cè)一次整個(gè)地球表面，因此它可以提供足夠的SST數(shù)據(jù)[1-2]。到2018年1月，世界海洋中大約有3 890個(gè)自由漂流的Argo浮標(biāo)，目前正在對(duì)從海面到2 000 m深度[3]的溫鹽等要素進(jìn)行實(shí)時(shí)觀測(cè)，但是Argo浮標(biāo)在全球大洋中分布不均勻。全球海洋如此廣闊，僅根據(jù)Argo剖面[4]很難獲得高分辨率網(wǎng)格化的三維溫度場(chǎng)。基于豐富的海表數(shù)據(jù)，海洋學(xué)家們已經(jīng)做了大量的研究工作來重建水下海溫?cái)?shù)據(jù)。

Fox等[5]設(shè)計(jì)并介紹了MODAS(Modular Ocean Data Assimilation System)系統(tǒng)，該系統(tǒng)用于重構(gòu)全球海洋的三維溫度場(chǎng)和鹽度場(chǎng)。MODAS使用《1994年世界海洋地圖集》WOA’94(World Ocean Atlas 1994)T/S氣候?qū)W作為背景場(chǎng)，通過最優(yōu)插值方法OI[6](Optimal Interpolation)從海洋觀測(cè)數(shù)據(jù)集( MOODS )中同化吸收歷史T / S剖面，然后可以得到代表海洋平均狀態(tài)的靜態(tài)MOADS三維溫鹽場(chǎng)。利用歷史觀測(cè)記錄的溫度和鹽度廓線求得水下的溫鹽值和海表溫度、海面位勢(shì)高度的回歸方程，借助該方程和近實(shí)時(shí)的衛(wèi)星觀測(cè)的SST、海面高度可以計(jì)算動(dòng)態(tài)MOADS三維溫場(chǎng)。因此，它可以被用作初始的第一猜測(cè)場(chǎng)，通過最優(yōu)插值法同化吸收任何現(xiàn)場(chǎng)數(shù)據(jù)，然后它將提供比靜態(tài)MOADS氣候?qū)W更準(zhǔn)確的“真實(shí)海洋”估計(jì)。這種方法基于長(zhǎng)期歷史觀測(cè)剖面數(shù)據(jù)和氣候態(tài)的溫鹽場(chǎng)。

Nardelli等[7]應(yīng)用m-EOF(Multivariate Empirical Orthogonal Function)方法分析了長(zhǎng)時(shí)間序列的海水溫度、鹽度和海面位勢(shì)高度數(shù)據(jù)，然后基于前三個(gè)EOF模態(tài)獲得了線性回歸方程，它把海表面的性質(zhì)向水下投影。然后基于該回歸方程重建了北太平洋的垂直剖面。這種方法依賴于兩個(gè)假設(shè)：一個(gè)是只需要EOF分析結(jié)果的前幾個(gè)模態(tài)可以解釋矩陣的大部分的值；另一個(gè)是海表溫度和海表位勢(shì)高度、海面鹽度之間存在很強(qiáng)的線性關(guān)系。Nardelli 等[8]和Yosuke等[9]還使用了一種類似于該m-EOF方法的方法來重構(gòu)局部海域的溫度剖面。這種m-EOF方法是基于長(zhǎng)期數(shù)據(jù)的多元EOF分析，因此即使在最壞的情況下，重建溫度也可以與氣候態(tài)的數(shù)據(jù)接近。Guinehut等[10]使用長(zhǎng)期的歷史溫度和鹽度剖面以及月平均數(shù)據(jù)來獲得變量間的線性回歸方程，然后再基于該線性回歸方程和衛(wèi)星觀測(cè)的海表溫度，海面高度異常及月平均數(shù)據(jù)進(jìn)行溫度場(chǎng)的重構(gòu)。

Alvera-Azcárate等[11]和王輝贊[12-13]使用EOF分析來恢復(fù)缺失數(shù)據(jù)，如溫度和鹽度，EOF分析首先分解時(shí)間序列數(shù)據(jù)組成的矩陣，然后選擇主成分來計(jì)算缺失數(shù)據(jù)。這種方法要求在該矩陣中缺失的數(shù)據(jù)點(diǎn)附近必須存在觀測(cè)值，即只是少量的數(shù)據(jù)缺失，否則結(jié)果將偏離真實(shí)值。

本文通過以下方式改進(jìn)了上述重構(gòu)三維海洋溫度場(chǎng)的方法。首先，作者沒有應(yīng)用多年觀測(cè)，而是使用了前一年來自CMEMS的3-D月平均溫鹽場(chǎng)，這是根據(jù)Guinehut的方法重建的。使用下一年的衛(wèi)星觀測(cè)的月平均海表數(shù)據(jù) (SLA和SST)來推導(dǎo)下一年的3-D月平均溫度場(chǎng)。使用相鄰或相近年份的數(shù)據(jù)集來獲得海表和水下要素之間的線性回歸關(guān)系可以提高重建的準(zhǔn)確性，因?yàn)閮蓚€(gè)數(shù)據(jù)集越近，它們之間的相關(guān)性就越高，均方根差異越小。其次還利用最有插值方法同化吸收了Argo的實(shí)測(cè)數(shù)據(jù)以提高精度。

2 數(shù)據(jù)

總共需要三種數(shù)據(jù)：三維T/S場(chǎng)用于設(shè)定水下溫度和海表溫度、海表位勢(shì)高度的回歸關(guān)系；SST和SLA等海表的遙感數(shù)據(jù)用于重建三維溫度場(chǎng)；用Argo剖面數(shù)據(jù)同化以提高結(jié)果的準(zhǔn)確性。

2.1 三維溫度和鹽度場(chǎng)

三維溫度和鹽度場(chǎng)是CMEMS提供的MULTIOBS_GLO_PHY_REP_015_002月平均數(shù)據(jù)集,下載地址http://marine.copernicus.eu/services-portfolio/access-to-products/?？臻g分辨率為1/4(°)，垂直范圍從- 5 500～0 m分成33個(gè)等級(jí)；時(shí)間覆蓋范圍是從1993—2016年。這個(gè)數(shù)據(jù)集是使用第1節(jié)引言中提到的Guinehut方法構(gòu)建的。

多年平均海溫和鹽度數(shù)據(jù)來自NOAA(National Oceanic and Atmospheric Administration)提供的WOA2018(World Ocean Atlas 2018)數(shù)據(jù)集。溫鹽數(shù)據(jù)為多年平均(2005—2017年)的各月的三維數(shù)據(jù)，空間分辨率為1/4(°)，下載地址https://www.nodc.noaa.gov/OC5/woa18/woa18data.html。

2.2 遙感數(shù)據(jù)集

SST由NASA的海洋水色網(wǎng)提供，下載地址https://oceandata.sci.gsfc.nasa.gov/，該數(shù)據(jù)通過負(fù)載在Terra衛(wèi)星和Aqua衛(wèi)星上的MODIS觀測(cè)的。Modis-Aqua下午從南向北穿越赤道，而Modis-Terra則是早上從北向南穿越赤道。本文分別從兩顆衛(wèi)星下載月平均9 km的SST，然后將它們合并在一起。

平均動(dòng)力地形MDT(Mean Dynamic Topography)是指參考大地水準(zhǔn)面的平均海面高度。MDT_CNET-CLS13是AVISO發(fā)布的1993—2012年期間海洋MDT的估計(jì)，相應(yīng)的網(wǎng)格分辨率為1/4(°)，如圖1(a)所示。下載地址https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/mdt.html。

SLA (海平面異常)是海面高度減去長(zhǎng)時(shí)間序列的平均海面高度所得的值，平均海面高度由地球重力、海洋環(huán)流等決定。文中使用了CMEMS提供的CLS-DOS-NT-06-034的月平均SLA，這是延遲時(shí)間為4個(gè)月的相對(duì)于1999—2012年的平均海面高度的海面高度異常，海面高度素來自全球海洋很多的衛(wèi)星觀測(cè)。它是分辨率為1/4(°)的網(wǎng)格產(chǎn)品， 2014年6月的海面高度異常如圖1(b)所示。

圖1 平均動(dòng)力地形和月平均海面高度異常Fig.1 Mean dynamic topography and mothly-mean sea level anomalies

2.3 Argo剖面數(shù)據(jù)

Argo是科學(xué)家們于1999年開始布置的全球海洋范圍的溫度/鹽度剖面觀測(cè)浮標(biāo)。下載的Argo剖面數(shù)據(jù)需要剔除不合格的數(shù)據(jù)，例如具有壓力梯度的異常等。計(jì)算的深度是1 000 m，所以選擇剖面深度大于1 000 m，區(qū)域范圍是140°W ～180°W，30°N ～ 50°N，在這個(gè)區(qū)域，2014年Argo觀測(cè)的溫度剖面線數(shù)量是3 512個(gè)，2015年是4 646個(gè)，2016年是4 298個(gè)(見圖2)。

((a)、(b)和(c)分別代表Argo觀測(cè)的2014、2015和2016年的合格溫度剖面線的觀測(cè)地點(diǎn)，不同月份對(duì)應(yīng)不同的顏色。(a), (b) and (c) represent the Argo annual observational points for 2014, 2015 and 2016 separately with different colors corresponding to different months, which have been dealt with by removing the unqualified profiles.)

圖2 Argo 觀測(cè)位置點(diǎn)
Fig.2 Argo observation points

3 重構(gòu)溫度場(chǎng)的流程

重構(gòu)三維溫度場(chǎng)的過程分為2步(見圖3)：(1)通過多元線性回歸方程得到溫度場(chǎng)T1；(2)利用最優(yōu)插值的同化方法，T1吸收同化進(jìn)Argo觀測(cè)的溫度值，從而得到準(zhǔn)確度更高的溫度場(chǎng)T2。下面將更詳細(xì)的說明每一步。

3.1 建立回歸方程

3.1.1 求解多元線性回歸方程的常系數(shù)和平均值 選取西太平洋區(qū)域作研究區(qū)域，該三維水域的水平海域范圍設(shè)定為：30°N～50°N，140°W～180°W，水平分辨率為1°；深度方向?yàn)閺暮１碇?1 000 m，劃分為19層(0， 10， 20， 30， 50， 75， 100， 125， 150， 200， 250， 300， 400， 500， 600， 700， 800， 900， 1 000 m)。在每個(gè)空間網(wǎng)格點(diǎn)上建立海溫和海表位勢(shì)高度、海表溫度的線性方程。在空間位置i，深度k，時(shí)間t處的二元一次線性方程為:

(1)

(2)

(輸入T、S表示第一年的來自CMEMS的再分析溫度和鹽度三維數(shù)據(jù)，輸入SLA和SST是第二年的衛(wèi)星觀測(cè)數(shù)據(jù)。輸入Argo溫度廓線是第二年的衛(wèi)星觀測(cè)數(shù)據(jù)，輸出T1, T2是重構(gòu)過程中每一步輸出第二年的三維溫度場(chǎng)。The inputs T, S are 3-D monthly-mean temperature and salinity in the past year provided by CMEMS, the inputs SST and SLA are satellite data in this year, and the inputs ARGO is the profiles in this year.The outputs T1, T2 are the reconstructed temperature fields in each step.)

圖3 重構(gòu)溫度場(chǎng)的流程圖
Fig.3 The flow chart illustrates the process of temperature reconstruction

(3)

3.1.3 位勢(shì)高度 在3.1.2 節(jié)中，y-1年份的海面位勢(shì)高度可以通y-1年份的三維溫鹽場(chǎng)計(jì)算得到，這節(jié)將顯示如何利用海面高度異常SLA來計(jì)算海面位勢(shì)高度，用到的各個(gè)變量的含義顯示在圖4中。

(海面高度異常SLA(Sea Level Anomaly)，平均動(dòng)力高度 MDT(Mean Dynamic Topography)，絕對(duì)動(dòng)力高度ADT (Absolute Dynamic Height)，位勢(shì)高度GH (Geopotential Height)，海表位勢(shì)高度h(Surface Steric Height).)

圖4 變量之間的關(guān)系
Fig.4 Diagram illustrates the relation

ADTy-1=MDT+SLAy-1，

(4)

ADTy=MDT+SLAy。

(5)

其中：y和y-1的意義和公式(3)中一樣； SLA是AVISO提供的衛(wèi)星觀測(cè)的海面高度異常。GHy-1可以通過下式計(jì)算：

GHy-1=ADTy-1-hy-1。

(6)

其中,hy-1是y-1年份的海表位勢(shì)高度，是由y-1年份的CMEMS的三維溫鹽場(chǎng)計(jì)算得到。那么y年份的海表位勢(shì)高度：

hy-1=ADTy-GHy-1。

(7)

3.1.4 重構(gòu)三維溫度場(chǎng) 為了測(cè)試多元回歸方程的重構(gòu)三維海溫的效果，選擇利用CMEMS提供的2014年的三維月平均溫鹽場(chǎng)和2015年的衛(wèi)星遙感數(shù)據(jù)來重構(gòu)2015年的三維海溫，并和CMEMS的2015年的三維海溫?cái)?shù)據(jù)做對(duì)比，并檢測(cè)效果。因此公式(3) 中的y=2015 ，y-1 = 2014，公式(3)變成如下格式：

(8)

(其(a)、(b) 為2014和2015年CMEMS的月平均三維溫度場(chǎng)的均方根差異；(c)、(d)為重構(gòu)的2015年的月平均三維溫度場(chǎng)和2015年CMEMS溫度數(shù)據(jù)的均方根差異。(a)，(b) is the RMS difference between the CMEMS monthly-mean temperature of 2015 and 2014. (c),(d) is the RMS difference of the reconstructed temperature field of 2015 and the 2015 CMEMS temperature field.)

圖5 垂直和水平方向的均方根差異
Fig.5 The vertical and horizontal RMS difference

圖5 (a)中, 深度方向的均方根差異從-1 000 m的0.3 ℃增長(zhǎng)到-600 m的1 ℃，在-500 m為 1.5 ℃，最大的深度方向的均方根差異可以達(dá)到2 ℃。而5(c)中，深度方向的均方根差異在-600 m處值為0.5 ℃，且在-200 m以深，都保持在1 ℃以下。圖5(b)中在西南部區(qū)域水平方向的均方根差異基本都大于1.8 ℃，而在5(d)中只有少部分區(qū)域的水平方向的均方根差異超過1.8 ℃。從上面的分析可得重構(gòu)的溫度場(chǎng)的精度還是較高的。

3.2 對(duì)比分析

3.1節(jié)中利用相鄰年份的溫鹽數(shù)據(jù)建立了重構(gòu)海溫的回歸關(guān)系，那么利用多年平均的溫鹽數(shù)據(jù)或者隔一年的溫鹽數(shù)據(jù)建立的回歸關(guān)系來重構(gòu)溫度場(chǎng)的效果如何呢。為了做對(duì)比，利用多年平均的WOA溫鹽數(shù)據(jù)建立了回歸方程，并以此重構(gòu)了2015和2016年的三維溫度場(chǎng)，它們和CMEMS的2015年的溫度場(chǎng)的均方根差異如圖6中的虛線所示。同時(shí)利用2014年的溫鹽場(chǎng)建立的回歸方程重構(gòu)了2016年的三維溫度場(chǎng)，其和CMEMS的2016的溫度場(chǎng)的均方根差異如圖6中點(diǎn)線所示。利用相鄰年份數(shù)據(jù)重構(gòu)的2015和2016年溫度場(chǎng)如圖中實(shí)線所示。從圖6的(a)(b)(c)(d)可見，利用多年平均數(shù)據(jù)建立的回歸方程重構(gòu)的溫度場(chǎng)在-200 m以上比利用相鄰年份重構(gòu)的溫度場(chǎng)精度高，但是在-200 m以下的溫度場(chǎng)，利用相鄰年份重構(gòu)的溫度場(chǎng)精度更高。從圖6的(e)(f)可見，利用隔一年的溫鹽數(shù)據(jù)重構(gòu)的溫度場(chǎng)的精度在-150 m以上比利用相鄰年份重構(gòu)的溫度場(chǎng)精度差，對(duì)-150 m以下的溫度場(chǎng)，二者相似。故得利用相鄰年份的溫鹽數(shù)據(jù)建立的回歸方程重構(gòu)溫度場(chǎng)具有一定的優(yōu)勢(shì)。

(其中(b)、(d)、(f)分別為對(duì)應(yīng)于(a)、(c)、(e)的年平均值。(b)、(d)、(f)corresponding to the annual average of (a)、(c)、(e).)

3.3 同化Argo數(shù)據(jù)

Argo觀測(cè)溫度剖面數(shù)據(jù)在該區(qū)域比較豐富，重構(gòu)的溫度場(chǎng)可以同化它以提高精度。本項(xiàng)目選擇最優(yōu)插值(Optimal interpolation，OI)同化方法，該方法簡(jiǎn)單高效，是使用頻率較高的同化方法[14]。從圖7(a)中可見，經(jīng)過同化ARGO觀測(cè)值的2015年重構(gòu)溫度場(chǎng)在次表層的大部分月份的均方根誤差都小于1.5 ℃。在深度200 m以下，垂直方向上的12個(gè)月的重構(gòu)溫度場(chǎng)的均方根誤差都小于1°。從均方根誤差的水平分布圖7(b)來看，大部分區(qū)域都小于1 ℃，只有極少數(shù)區(qū)域超過1.5 ℃。對(duì)比圖5和7，發(fā)現(xiàn)同化后的溫度場(chǎng)確實(shí)提高了精度。

至此共存在三種溫度場(chǎng)：利用回歸方程重構(gòu)得到的2015年的三維溫度場(chǎng)，OI同化吸收Argo觀測(cè)值后的溫度場(chǎng)以及CMEMS的溫度場(chǎng)。把它們分別插值到Argo觀測(cè)點(diǎn)上，并比較它們的值和Argo觀測(cè)值。從圖8可見經(jīng)同化后的重構(gòu)溫度廓線和Argo觀測(cè)廓線較接近，表示重構(gòu)的溫度場(chǎng)的精度是可接受的。

3.4 重構(gòu)的溫度場(chǎng)、CMEMS溫度場(chǎng)和觀測(cè)值比較

本文重構(gòu)月平均三維溫度場(chǎng)依賴于相鄰的年份的溫鹽場(chǎng)，因?yàn)橄噜從攴莸臄?shù)據(jù)之間有很強(qiáng)的相關(guān)性，同時(shí)均方根差異也越小。例如2014和2015年的CMEMS的三維溫度場(chǎng)之間的相關(guān)系數(shù)為97.92%，2015和2016年之間的為 97.97% ，2014和2016年之間的為97.86%；各自之間的均方根差異分別為1.32，1.33和1.37 ℃。圖1中的第一步利用由前一年的溫鹽場(chǎng)得到的回歸方程來重構(gòu)得到的三維溫度場(chǎng)要比利用多年數(shù)據(jù)來重構(gòu)溫度場(chǎng)要更有效。CMEMS的三維溫度場(chǎng)是基于長(zhǎng)時(shí)間序列計(jì)算得到的。把CMEMS的溫度場(chǎng)和本文重構(gòu)的溫度場(chǎng)插值到有Argo觀測(cè)的溫度廓線的地點(diǎn)，并分別計(jì)算它們和Argo觀測(cè)的溫度廓線之間的均方根誤差，結(jié)果顯示在圖9中。從圖9可見，2015年的CMEMS的溫度場(chǎng)與Argo溫度廓線的均方根差異比重構(gòu)的溫度場(chǎng)和Argo溫度廓線的均方根差異要大近0.2 ℃。2016年二者的差距更大，這說明本項(xiàng)目的重構(gòu)的結(jié)果比CMEMS的溫度場(chǎng)略好。

圖7 垂直和水平方向的均方根差異Fig.7 The vertical and horizontal RMS difference

(其中小圖的藍(lán)色點(diǎn)表示argo觀測(cè)的位置，彩色曲線表示插值所得的溫度廓線。The position of the ARGO observation for each figure is shown by the corresponding blue point in the small figure.)

圖8 溫度廓線
Fig.8 Temperature profiles

圖9 2015和2016年重構(gòu)的溫度、CMEMS溫度數(shù)據(jù)與Argo觀測(cè)值的均方根差異Fig.9 The RMS difference between Argo observation and reconstructed temperature,temperature of CMEMS during 2015 and 2016

4 結(jié)語(yǔ)

本文建立了一種重構(gòu)三維溫度場(chǎng)的方法，重構(gòu)了流程。該方法基于相鄰年份的3維溫度場(chǎng)來重構(gòu)，因?yàn)橄噜從攴莸娜S溫度場(chǎng)之間有很強(qiáng)的相關(guān)性。例如，2014和2015年的CMEMS三維溫度場(chǎng)的相關(guān)系數(shù)為97.92%，均方根差異為1.32 ℃；2015和2016年為97.97%，均方根差異為1.33 ℃；2014和2016年的為97.86%，均方根差異為1.37 ℃。按照?qǐng)D1的流程，利用2014年的CMEMS三維溫度場(chǎng)重構(gòu)2015年的溫度場(chǎng)，第一步的結(jié)果即通過線性回歸方程重構(gòu)后所得結(jié)果和CMEMS溫度場(chǎng)的相關(guān)系數(shù)為99.09%，均方根差異為0.88 ℃。第二步加入Argo同化后的重構(gòu)溫度場(chǎng)與與CMEMS溫度場(chǎng)的相關(guān)系數(shù)為99.29%，均方根差異為0.79 ℃。利用2015年的CMEMS三維溫度場(chǎng)重構(gòu)2016的溫度場(chǎng)的這二步的相關(guān)系數(shù)分別為：98.94%、99.12%，均方根差異分別為0.99和0.89 ℃。所以利用前一年的三維溫度場(chǎng)來逐月重構(gòu)當(dāng)年的溫度場(chǎng)是可行的。

重構(gòu)的三維溫度場(chǎng)在200 m以下具有較好的精度，垂直均方根誤差小于1 ℃，在次表層垂直均方根誤差超過1°C。從7 —9月，重構(gòu)溫度場(chǎng)的均方根誤差比其他月份大。水平的均方根誤差分布不均勻，在大部分區(qū)域都小于1 ℃，在少數(shù)區(qū)域超過1.5°。所以如果在那些均方根誤差較大的地方的次表層，存在從可能少數(shù)字月至9月的觀測(cè)數(shù)據(jù)，那么重構(gòu)的溫度場(chǎng)就可以同化吸收這些觀測(cè)值從而可以進(jìn)一步提高精度。

致謝：Argo觀測(cè)的溫度剖面數(shù)據(jù)由國(guó)際Argo項(xiàng)目組提供，下載地址：ftp://ftp.ifremer.fr/ifremer/argo/dac/。三維月平均溫/鹽場(chǎng)由CMEMS 負(fù)責(zé)制作和發(fā)布，下載地址：http://marine.copernicus.eu。

附錄

a．位勢(shì)高度

在溫度T，鹽度S和壓力p處，單位重量的海水體積為：

(1)

δ(S,T,p)=α(S,T,p)-α(35,0,p)。

(2)

叫做特征體積異常。動(dòng)力高度φ(相當(dāng)于大氣中的重力勢(shì))表示從給定的位置垂直移動(dòng)單位體積的水體至海面所要做的功，單位為J/kg或者m2/s2。相對(duì)于海面高度的φ表達(dá)式如下：

(3)

從靜力學(xué)公式

p=-ρgz

(4)

可以導(dǎo)出如下關(guān)系式：

(5)

其中：g表示重力加速度；z向上增長(zhǎng)。動(dòng)力高度φ和動(dòng)力高度異常ψ的關(guān)系能夠被推導(dǎo)出，如下

(6)

該式表明動(dòng)力高度異常等于重力勢(shì)φ減去一個(gè)壓力函數(shù)。動(dòng)力高度異常可通過下式計(jì)算：

(7)

重力加速度可以設(shè)為在45°N的中的值9.806 65 m/s2，因此位勢(shì)高度可以定義為：

(9)

它是具有長(zhǎng)度單位米的地轉(zhuǎn)流函數(shù)。位勢(shì)高度可以和衛(wèi)星觀測(cè)的海面高度異常集合，來計(jì)算某一深度至大地水準(zhǔn)面之間的重力勢(shì)高度。