基于衛(wèi)星熱紅外數(shù)據(jù)的田灣核電溫排水監(jiān)測

王 虎

（1.江蘇核電有限公司，江蘇 連云港 222000）

核電作為清潔能源，越來越受到國家的重視。在《國家核電中長期發(fā)展規(guī)劃（2005—2020年）》的指導(dǎo)下，我國核電事業(yè)得到飛速發(fā)展；截至2020年5月，我國已有47臺機(jī)組投入商業(yè)運(yùn)行。與此同時(shí)，核電站溫排水對周邊海域造成的熱影響也越來越受到人們的關(guān)注[1-2]。溫度作為水體的重要參數(shù)，將給水體中各種生物化學(xué)性質(zhì)帶來影響，從而影響各類生物生長和繁殖活動。因此，核電站周邊海域溫排水監(jiān)測對水環(huán)境保護(hù)至關(guān)重要[3]。核電站溫排水監(jiān)測方法主要包括海上實(shí)測、數(shù)模和物模計(jì)算等[4]。隨著熱紅外遙感技術(shù)的快速發(fā)展，遙感技術(shù)逐漸成為經(jīng)濟(jì)快捷的溫排水監(jiān)測方法[5-6]。

現(xiàn)階段可免費(fèi)獲取的空間分辨率百米級的熱紅外波段航天遙感數(shù)據(jù)源主要包括美國的Landsat5/7/8 數(shù)據(jù)和我國的HJ-1B衛(wèi)星數(shù)據(jù)。經(jīng)過長期研究，國內(nèi)外學(xué)者針對熱紅外數(shù)據(jù)提出了諸多溫度反演算法[7-9]，并將這些反演算法應(yīng)用于溫排水監(jiān)測中[10-15]。本文選用兩組同日過境的Landsat8 TIRS 和HJ-1B IRS 熱紅外數(shù)據(jù)（以下分別簡稱為TIRS和IRS），根據(jù)數(shù)據(jù)特點(diǎn)開展溫度反演，并選取一期數(shù)據(jù)開展海上測量以及與MO?DIS溫度數(shù)據(jù)的交叉驗(yàn)證，研究了兩種熱紅外數(shù)據(jù)溫度反演結(jié)果的可靠性以及同期數(shù)據(jù)溫度分布趨勢的一致性，分析了不同時(shí)相田灣核電站附近海域海水表面溫度的分布狀況，為核電廠溫排水監(jiān)測和管理提供參考。

1 數(shù)據(jù)來源與預(yù)處理

本文采用的數(shù)據(jù)源為2013-11-15 和2017-02-27的Landsat8數(shù)據(jù)、HJ-1B數(shù)據(jù)以及2013-11-15的海面實(shí)測數(shù)據(jù)。為驗(yàn)證數(shù)據(jù)的有效性，還獲取了一景同日過境的MODIS數(shù)據(jù)，衛(wèi)星過境時(shí)田灣核電兩臺機(jī)組在滿功率運(yùn)行，且均處在冬季落急潮態(tài)。Landsat8是美國陸地衛(wèi)星觀測計(jì)劃中的第八顆衛(wèi)星，衛(wèi)星傳感器主要參數(shù)見參考文獻(xiàn)[16]；HJ-1B是我國用于環(huán)境和災(zāi)害監(jiān)測的對地觀測衛(wèi)星，其熱紅外載荷相關(guān)參數(shù)如表1所示。數(shù)據(jù)準(zhǔn)備完成后，對遙感數(shù)據(jù)進(jìn)行輻射定標(biāo)、幾何精校正、濾波、幾何裁剪、水陸分離和去云等預(yù)處理。

2 海面溫度反演

Landsat8 數(shù)據(jù)具有兩個(gè)熱紅外波段，部分學(xué)者針對其熱紅外波段特征開展了分裂窗算法研究[17-18]，但因其B11波段存在不規(guī)律條帶，利用該算法獲取的結(jié)果誤差較大。因此，本文采用輻射傳輸方程算法對其B10波段開展反演。IRS的熱紅外通道(10.4～12.5 μm)與Landsat5 的熱紅外通道（10.5～12.5 μm）接近，本文利用方法簡便、精度較高的單窗算法進(jìn)行溫度反演[19-21]，并根據(jù)其波段設(shè)置的特點(diǎn)對參數(shù)進(jìn)行修訂[22]。

2.1 Landsat8海面溫度反演

首先需模擬出大氣對海面熱輻射的干擾；然后從TIRS觀測到的輻射值總量中減去大氣干擾量，得到海面真實(shí)熱輻射強(qiáng)度；再把這一熱輻射強(qiáng)度轉(zhuǎn)化為相應(yīng)的海面溫度。其計(jì)算公式為：

式中， L10為B10 波段大氣頂層輻射；ε10為海面發(fā)射率；TS為海表溫度，單位為K； L10( )TS為研究對象溫度為TS時(shí)的黑體輻射； L10atm↑為大氣上行輻射；L10atm↓為大氣下行輻射；τ 為B10波段大氣透過率。

由式（1）可知，要獲取海面溫度，需獲取大氣透過率、大氣上行輻射、大氣下行輻射和海面發(fā)射率4個(gè)參數(shù)。

1）大氣頂層輻射主要是將熱紅外波段的DN值通過增益和平移轉(zhuǎn)換為輻射值。其計(jì)算公式為：

式中， DN 為像元灰度值； M10和A10為影像的增益和偏移，這些數(shù)值可從元數(shù)據(jù)文件中獲取。

2）海面發(fā)射率是地物向外輻射電磁波的能力表征，與地物的表面形態(tài)和物理性質(zhì)有關(guān)。因海面與黑體接近（發(fā)射率為1），根據(jù)前人研究，ε10取0.995[19]。

3）大氣上行輻射、大氣下行輻射、大氣透過率均與大氣作用有關(guān)，本文基于NASA網(wǎng)絡(luò)大氣校正平臺，根據(jù)研究區(qū)域海拔，衛(wèi)星過境時(shí)刻的氣壓、溫度、濕度，影像獲取時(shí)間以及中心經(jīng)緯度等參數(shù)進(jìn)行大氣校正[23-24]。

在獲取上述大氣參數(shù)后，可計(jì)算得到海表溫度為Ts時(shí)的實(shí)際輻射值，再求得海面真實(shí)溫度。其計(jì)算公式為：

式中，對于B10，K1=774.89 W/(m2·sr·μm)；K2=1 321.08K。

2.2 HJ-1B海面溫度反演

HJ-1B僅有一個(gè)熱紅外波段。本次反演主要采用覃志豪[19]等推導(dǎo)的單窗算法，并根據(jù)中國資源衛(wèi)星中心韓啟金[22]等的研究對參數(shù)進(jìn)行修訂，大氣透過率計(jì)算部分參考單通道算法。其計(jì)算公式為：

式中，Tsensor為大氣頂部的亮度溫度，單位為K；Ta為大氣等效溫度，單位為K；a、b 為常數(shù)；c=ετ ；d=(1-τ)[1+(1-ε)]； Ts為地表溫度，單位為K； τ 為大氣透過率；ε 為海水發(fā)射率，本文取定值0.995。

1）大氣透過率是影響紅外輻射傳輸?shù)闹匾蛩兀茉谳^短時(shí)間內(nèi)受氣壓、水汽含量影響發(fā)生明顯變化。除水汽含量外，其他氣體混合較均勻，且含量較穩(wěn)定，因此，大氣透過率主要受大氣的影響。本文根據(jù)衛(wèi)星熱紅外波段的實(shí)際設(shè)置情況，建立了水汽含量在0.1~5.0 g/cm2區(qū)間內(nèi)與大氣透過率之間的回歸關(guān)系（R2=0.99），即

2）大氣等效溫度。根據(jù)環(huán)境衛(wèi)星熱紅外波段的光譜響應(yīng)特點(diǎn)，利用大氣溫度的等效算法[22]對其進(jìn)行參數(shù)調(diào)整，從而建立中緯度冬季大氣等效溫度與大氣地面溫度的關(guān)系。海表溫度作為估算溫度初始值，以海面亮度溫度為基礎(chǔ)進(jìn)行估算。

2.3 海面溫度反演結(jié)果

本文根據(jù)兩種熱紅外數(shù)據(jù)的溫度反演結(jié)果（圖1、圖2），獲得田灣核電周邊10 km 海域的各個(gè)參數(shù)情況，如表2、3所示。

圖2 2017-02-27的熱紅外溫度場圖

表2 2013-11-15的TIRS與IRS溫度反演結(jié)果對比/℃

圖1 2013-11-15的熱紅外溫度場圖

3 研究結(jié)果與分析

3.1 海面反演溫度可靠性分析

本文在2013-11-15 衛(wèi)星過境期間開展了海面近同步溫度測量，對現(xiàn)場觀測和測量的精確地理位置進(jìn)行定位，以保證測量數(shù)據(jù)與遙感數(shù)據(jù)的位置相對應(yīng)。衛(wèi)星過境前后，在核電溫排水區(qū)域至基準(zhǔn)溫度值海域內(nèi)進(jìn)行反復(fù)測量，如圖3 所示。根據(jù)海面測量獲取的數(shù)據(jù)以及溫度反演獲取的同點(diǎn)反演數(shù)據(jù)，利用最小二乘法進(jìn)行擬合，獲取的擬合結(jié)果如圖4所示。

圖3 2013-11-15的海面實(shí)測點(diǎn)位圖

圖4 2013-11-15的實(shí)測數(shù)據(jù)與反演結(jié)果數(shù)據(jù)擬合結(jié)果

表3 2017-02-27的TIRS與IRS溫度反演結(jié)果對比/℃

2013-11-15 的 Landsat8 反演獲得的 SST 值與實(shí)測數(shù)據(jù)的擬合關(guān)系為y=0.719 1x+4.907 7，相關(guān)系數(shù)R2超過了0.96，標(biāo)準(zhǔn)誤差為0.37；HJ-1B反演獲得的SST值與實(shí)測數(shù)據(jù)的擬合關(guān)系為y=0.511 8x+8.404 6，相關(guān)系數(shù)R2也超過了0.96，標(biāo)準(zhǔn)誤差為0.20。由線性相關(guān)程度和誤差大小可知，溫度反演結(jié)果是準(zhǔn)確可信的。

3.2 與MODIS數(shù)據(jù)的交叉驗(yàn)證

MODIS 是NASA 研制的空間遙感儀器，共設(shè)置了36 個(gè)波段，包含16個(gè)熱紅外波段，每天可獲取全球任意地點(diǎn)影像數(shù)據(jù)。B31波段和B32波段由于對水汽的吸收作用不同，受太陽光反射的影響較微弱，可利用算法消除水汽的影響獲取地表溫度，尤以分裂窗算法最為成熟[23]。NASA 將 MODIS 海表溫度 產(chǎn) 品 免費(fèi)對外公布[24]，經(jīng)過長期驗(yàn)證，其溫度精度介于0.053～0.66℃[23-24]。

以2013-11-15 的Landsa8 熱紅外數(shù)據(jù)反演結(jié)果為例，為進(jìn)一步驗(yàn)證反演結(jié)果，以MODIS數(shù)據(jù)的海表溫度（圖5）為基礎(chǔ)，進(jìn)行交叉驗(yàn)證[25]。衛(wèi)星過境時(shí)間有一定的間隔（表1），主要為落潮末期到漲潮初期，且落潮末期占時(shí)間比例較大（圖6），由于落潮末期海水比較穩(wěn)定，漲潮時(shí)間較短，大量外海海水還未涌入該海域，因此海表溫度變化不大。

圖5 2013-11-15的MODIS熱紅外溫度場圖

圖6 潮汐狀態(tài)變化示意圖

在交叉驗(yàn)證數(shù)據(jù)上隨機(jī)選取30個(gè)點(diǎn)，利用最小二乘法研究二者之間的關(guān)系，結(jié)果如表4、圖7 所示。由于均處于落潮狀態(tài)，近岸灘涂出露面積較大，根據(jù)比熱容的差異，受太陽輻射影響，灘涂溫度高于海水溫度。MODIS數(shù)據(jù)的空間分辨率為1 km，海陸像元混合比Landsat8 大，點(diǎn)位選取時(shí)優(yōu)先選取遠(yuǎn)離近岸灘涂海域的點(diǎn)位。Landsat8 的最大偏差為1.56℃，最小偏差小于0.1℃。兩組數(shù)據(jù)的擬合關(guān)系為 y=0.844 7x+2.257 6，相關(guān)系數(shù)為0.826 6，標(biāo)準(zhǔn)誤差為0.39。結(jié)果表明，Landsat8的反演結(jié)果與MODIS溫度數(shù)據(jù)之間線性特征明顯，具有較好的線性相關(guān)性和一致性，說明Landsat8 熱紅外波段溫度反演方法獲得的溫度場數(shù)據(jù)是準(zhǔn)確可信的。

圖7 反演結(jié)果與MODIS數(shù)據(jù)的擬合圖

表4 Landsat8反演結(jié)果與MODIS溫度結(jié)果的對比表

3.3 海域溫度特征分析

由圖1可知，2013-11-15核電周邊海域溫度范圍為13.0~23.0℃，溫度分布層次分明；連島北部海域溫度主要集中在13.7~14.5℃；核電南部徐圩港正在建設(shè)過程中，周邊海域溫度主要集中在13.4~14.1℃。核電周邊海域的溫度場明顯受到了溫排水的影響，距離排水口越近溫度越高，在排水口處溫度高達(dá)23℃，外海最低溫度為13.0℃。由于衛(wèi)星過境時(shí)處于落潮潮態(tài)，溫排水沿落潮方向東北擴(kuò)散；而取水口位于排水口東北側(cè)、溫排水?dāng)U散通道上，因此略高于基準(zhǔn)溫度范圍。

與2013-11-15相比，2017-02-27核電周邊海域岸線發(fā)生了顯著的變化，取水口處的明渠由原來的1.9 km增加到4.5 km，排水口處的導(dǎo)流堤向南擴(kuò)建了1.5 km，南側(cè)徐圩港防波堤離岸距離從2.6 km增加到約6.9 km。由于潮汐狀態(tài)相似，2017-02-27的溫度場梯度和空間分布特征與2013-11-15 的溫度場相似（圖2）。由于整體溫度差異，絕對溫度較低：排水口東側(cè)海域溫度場范圍為6.1~14.7℃，主要集中在6.5~13.2℃；徐圩港建成后的內(nèi)部溫度主要集中在5.5~6.6℃；連島北部海域溫度主要集中在6.6~7.2℃。

由于衛(wèi)星過境時(shí)，數(shù)據(jù)同樣處于落潮狀態(tài)，2017-02-27的溫度場明顯受潮態(tài)和岸線影響：從排水口涌出的高溫?zé)崴貙?dǎo)流堤向南流動，越過導(dǎo)流堤后主體沿落潮向東北方向擴(kuò)散，受取水明渠的限制，溫排水被阻擋在明渠南側(cè)，隨著海水混合，溫度逐漸降低，從14.7℃降至6.8℃；沿岸向東南側(cè)擴(kuò)散的溫排水受到南側(cè)徐圩港防波堤的阻擋，沿防波堤向東北展布。

根據(jù)Landsat8和HJ-1B數(shù)據(jù)溫度反演的結(jié)果（圖1、圖2），結(jié)合其他學(xué)者的研究結(jié)果可知，核電周邊海域溫度參數(shù)（表1、表2）基本相同，兩種傳感器反演的海水表面溫度空間分布趨勢一致；而由于空間分辨率的關(guān)系，TIRS數(shù)據(jù)可以反映更精細(xì)的溫度場邊緣[26]。

3.4 熱影響統(tǒng)計(jì)分析

對于田灣核電附近海域，基準(zhǔn)溫度采用核電附近不受溫排水影響區(qū)域的連島北部海域的平均溫度[27]。在確定基準(zhǔn)溫度后，將海域溫度場數(shù)據(jù)整體減去基準(zhǔn)溫度；再根據(jù)數(shù)據(jù)熱影響實(shí)際情況，劃分熱影響等級，并進(jìn)行編碼（圖8、9）；最后分類統(tǒng)計(jì)各等級像元面積（圖10、11），并提取核電溫排水熱影響分布信息[26]。

由圖8、9 可知，核電溫排水由排水口向外擴(kuò)散，溫度迅速降低，到達(dá)外海基準(zhǔn)溫度后，溫度基本保持不變，與同日過境數(shù)據(jù)海水表面溫度空間分布基本一致，IRS數(shù)據(jù)與TIRS數(shù)據(jù)反演得到的熱影響擴(kuò)散形狀分布趨勢一致。由于氣象條件、過境時(shí)間、潮態(tài)方面不完全相同，兩種傳感器獲取的結(jié)果存在一定差別。由圖10、11可知，同日過境數(shù)據(jù)各級別的熱影響區(qū)面積差別很小；相較于2013-11-15 的反演結(jié)果，2017-02-27獲取的結(jié)果各級別面積均有擴(kuò)大，但取水口處于基準(zhǔn)溫度范圍內(nèi)。岸線的變化，客觀上使核電排水口位于一個(gè)由兩側(cè)海工建筑環(huán)抱的人工海灣的灣底，影響了溫排水的擴(kuò)散。由于兩組數(shù)據(jù)獲取時(shí)的核電運(yùn)行工況、所處季節(jié)和潮汐狀態(tài)類似，因此面積增大的原因?yàn)楣こ探ㄔO(shè)阻礙了溫排水的擴(kuò)散。

圖8 2013-11-15的熱影響編碼圖

圖9 2017-02-27的熱影響編碼圖

圖10 2013-11-15的IRS與TIRS熱影響面積對比圖

兩個(gè)時(shí)間TIRS監(jiān)測到的最高熱影響級別均比IRS高1℃，說明TIRS數(shù)據(jù)可更準(zhǔn)確地描述溫排水對周邊海域熱影響的分布狀況；同時(shí)，熱影響面積統(tǒng)計(jì)結(jié)果顯示，7℃以上、6℃以上的熱影響面積分別小于總面積的1%和2%，兩種數(shù)據(jù)反演的溫排水對環(huán)境的熱影響主要集中在前9 個(gè)級別。因此，IRS 和TIRS 對核電溫排水環(huán)境影響監(jiān)測的結(jié)果基本相同。

圖11 2017-02-27的IRS與TIRS熱影響面積對比圖

4 結(jié) 語

1）本文對TIRS 數(shù)據(jù)和IRS 數(shù)據(jù)獲取的溫度反演結(jié)果進(jìn)行實(shí)測數(shù)據(jù)擬合和交叉驗(yàn)證。結(jié)果表明，溫度反演的結(jié)果可靠，兩種數(shù)據(jù)在數(shù)值和空間分布上具有很好的一致性。

2）不同時(shí)期遙感數(shù)據(jù)調(diào)查結(jié)果顯示，田灣核電周邊工程建設(shè)改變了岸線分布，溫排水分布形態(tài)和規(guī)模受到工程影響。

3）TIRS 數(shù)據(jù)和IRS 數(shù)據(jù)獲取的溫排水溫度場在數(shù)值和空間上基本一致，TIRS可監(jiān)測到更高的熱影響級別。由于溫升主要集中在6℃以下的熱影響級別，因此HJ-1B 可達(dá)到與Landsat8 基本等同的調(diào)查效果。兩種熱紅外數(shù)據(jù)相互補(bǔ)充，為評估核電站溫排水對周邊海域溫度環(huán)境的影響提供了迅速便捷的手段。