石海崗 梁春利★ 張建永 張春雷 程旭
(1.核工業(yè)航測遙感中心,河北石家莊 050002;2.東華理工大學(xué),江西南昌 330013)
隨著我國核電事業(yè)的飛速發(fā)展和社會環(huán)境保護(hù)意識的不斷增強(qiáng),核電站溫排水對其受納水域環(huán)境造成的熱影響越來越受到人們的關(guān)注[1-2]。溫排水一方面改變了排水口附近海域的流場,另一方面則會引起局部水域的溫度快速升高。作為重要的水質(zhì)和水域生態(tài)環(huán)境要素,水體溫度幾乎影響水的各種物理、化學(xué)和生物化學(xué)性質(zhì),影響水體的質(zhì)量,從而間接影響到各類水生物的生長和繁殖活動,甚至產(chǎn)生明顯的危害效應(yīng)。因此,核電站溫排水的監(jiān)測評價(jià)對研究水環(huán)境生態(tài)平衡有著重要意義[3]。對于核電站溫排水的監(jiān)測,以往所采用的方法主要包括數(shù)值模擬計(jì)算、物理模擬實(shí)驗(yàn)、現(xiàn)場測量等方法[4]。數(shù)值模擬計(jì)算、物理模擬實(shí)驗(yàn)主要用于核電運(yùn)行前階段。海面實(shí)測是傳統(tǒng)人工測量方法,受時(shí)間、空間的限制較大,同步性較差。隨著科技進(jìn)步,熱紅外遙感技術(shù)因其同步性、周期性和經(jīng)濟(jì)性已成為海溫監(jiān)測有效的調(diào)查手段[5-6]。
田灣核電站廠址所在海域?qū)冱S海北部,漲潮時(shí),潮流從東北方向涌入,自北向南分別進(jìn)入連島鎮(zhèn)南部水域、核電站取水明渠、核電周邊排水區(qū)域,順岸向南東方向形成沿岸漲潮流態(tài);落潮流態(tài)大致與漲潮流相反。排水口附近區(qū)域呈半圓形向東開放,半徑1km 范圍內(nèi),海拔從3米緩慢降低到-2米,屬淤泥分布區(qū),北岸為海水封閉區(qū)域(圖1)。(見中間彩頁,下同。)
近年來,隨著江蘇省沿海開發(fā)國家戰(zhàn)略實(shí)施及其他工程活動,田灣核電周邊海域岸線發(fā)生了顯著的變化:2013-2017 年間,田灣核電站南部的徐圩港建設(shè)初步完工,徐圩港區(qū)工程建設(shè)深入海域;田灣核電排水口導(dǎo)流堤,取水明渠的建設(shè),使核電站周邊海洋環(huán)境發(fā)生持續(xù)變化。國內(nèi)海域岸線變遷關(guān)注的焦點(diǎn)多集中在氣象、漁業(yè)資源、生態(tài)環(huán)境等宏觀方面,海岸線的變化對核電溫排水研究罕有報(bào)道。
Landsat-8 數(shù)據(jù)既具有信息豐富的多光譜波段,又具有信噪比較高的熱紅外波段,本次研究基于Landsat-8 數(shù)據(jù),利用其多光譜數(shù)據(jù)定量獲取了2013、2017 年間2 期田灣核電周邊海域海岸線的數(shù)據(jù),利用其熱紅外數(shù)據(jù)反演獲取了周邊海域的溫度場情況,對田灣核電周邊海域海岸線的時(shí)空變化特征及對溫排水的影響進(jìn)行分析,以查明岸線變化對核電溫排水帶來的影響。因?yàn)槊看芜b感資料對應(yīng)的潮流強(qiáng)度不一樣,只有在相同或相似潮態(tài)的情況下,才能比較不同岸線造成的流速、流向的變化,為保障對比的客觀性,選取的數(shù)據(jù)為同潮態(tài)下兩臺機(jī)組滿功率運(yùn)行條件下的數(shù)據(jù)。
本文數(shù)據(jù)源為:2013年11月15日、2017年2月27 日Lndsat-8 數(shù)據(jù)。為驗(yàn)證反演結(jié)果,獲取了2013年11 月15 日的海面測量數(shù)據(jù)和同日過境的MODIS 數(shù)據(jù)。衛(wèi)星過境時(shí),田灣核電兩臺機(jī)組在滿功率運(yùn)行、且均處在冬季落潮狀態(tài)。
Landsat-8是由美國地質(zhì)調(diào)查局及太空署第八個(gè)陸地衛(wèi)星計(jì)劃,衛(wèi)星上攜帶有OLI 和TIRS 兩個(gè)主要載荷[7],相關(guān)參數(shù)見表1。數(shù)據(jù)準(zhǔn)備完成后,對Landsat-8遙感影像進(jìn)行預(yù)處理,包括輻射定標(biāo)、幾何裁剪、幾何精校正、濾波、水陸分離及去云等處理。

表1 Landsat-8技術(shù)參數(shù)信息Table 1 Information on technical parameters of Landsat-8
利用ENVI 軟件進(jìn)行預(yù)處理后,為了更準(zhǔn)確進(jìn)行解譯,分別對兩期Landsat-8數(shù)據(jù)進(jìn)行大氣校正,采用Gram-Schmidt Pan Sharpening方法將影像融合為15m,為了突出地物特征波段組合采用6,5,2組合。
經(jīng)過以上處理,將圖像進(jìn)行假彩色合成,制作兩期1∶5萬的影像圖用來解譯分析。
利用Arcgis 制圖軟件,結(jié)合兩期遙感影像,采用遙感動態(tài)監(jiān)測中的常用的目視解譯方法對兩期遙感影像進(jìn)行單獨(dú)解譯,然后通過對各解譯結(jié)果進(jìn)行比較,直接提取變化信息。
解譯結(jié)果見圖3 所示,在2013-2017 年四年間,田灣核電周邊岸線發(fā)生了顯著的變化,進(jìn)水口處的取水明渠由原來的1.9km 增加到4.5km,排水口處導(dǎo)流堤向南擴(kuò)建了1.5km,南側(cè)徐圩港防波堤,則從2.6km,增加到了6.9km,核電排水口處于兩側(cè)防波堤環(huán)抱的人工海灣的灣底(圖2)。
因Landsat-8 TIRS 數(shù)據(jù)與MODIS 熱紅外波段類似,有學(xué)者對其熱紅外波段開展了劈窗算法研究[8-9],但根據(jù)美國地質(zhì)調(diào)查局網(wǎng)站(https://glovis.usgs.gov/#)公布的測試結(jié)果,TIRS 的11 波段由于條帶太突出,反演結(jié)果干擾太大,無法應(yīng)用。本文采用輻射傳輸方程算法,對Landsat-8 TIRS 溫度反演10波段進(jìn)行反演。
衛(wèi)星傳感器TIRS(熱紅外傳感器)接收到的熱紅外輻射值由大氣向上輻射亮度、地面的真實(shí)輻射亮度經(jīng)過大氣層之后到達(dá)衛(wèi)星傳感器的能量組成。公式如下:

式中,Lλ由傳感器接受到的大氣頂層輻射,ελ是地表的比輻射率,TS是地表溫度,Lλ(TS)是溫度為TS時(shí)的黑體輻射,通過普朗克(Planck)定律求得,Lλatm↓是大氣下行輻射,Lλatm↑是大氣上行輻射,τ是地表和傳感器之間的大氣透射率。
由輻射傳輸方程可知,要求算地表溫度,需要知道4個(gè)參數(shù)的值:大氣透射率τ、大氣上行輻射亮度Lλatm↑、大氣下行輻射亮度Lλatm↓,和地表比輻射率ελ。
(1)Lλ的計(jì)算
Lλ的計(jì)算主要是指將傳感器觀測到的圖像灰度值轉(zhuǎn)換成輻射值的過程。公式如下:

其中,Qcal為像元灰度值;ML和AL分別為圖像的增益和偏移。定標(biāo)系數(shù)可以直接從元數(shù)據(jù)中獲取。
(2)比輻射率ελ
物體的比輻射率是物體向外輻射電磁波的能力表征,受很多因素制約,與物體的表面狀態(tài)及物理性質(zhì)有關(guān)。本次反演主要針對海面進(jìn)行,接近于黑體(比輻射率為1),比輻射率取定值0.995。
(3)其他參數(shù)的獲取
大氣下行輻射Lλatm↓,大氣上行輻射Lλatm↑,地表和傳感器之間的大氣透射率τ,與大氣作用有關(guān)。本次研究根據(jù)美國國家環(huán)境預(yù)測中心(NEPC)提供的標(biāo)準(zhǔn)大氣剖面,結(jié)合MODTRAN4.0 模塊建立的大氣校正模型,進(jìn)行大氣校正[10-11],通過輻射傳輸法,消除大氣的影響。根據(jù)田灣核電站提供的衛(wèi)星過境時(shí)刻的氣壓,地表溫度,相對濕度,影像時(shí)間以及中心經(jīng)緯度獲取以上參數(shù)。
在獲取大氣下行輻射Lλatm ↓,大氣上行輻射Lλatm ↑,地表和傳感器之間的大氣透射率τ 參數(shù)后,計(jì)算出海表真實(shí)的輻射亮度值Lλ(TS),根據(jù)普朗克公式的反函數(shù),求得地表真實(shí)溫度Ts:Ts=K2/ln(K1/LT+1)。
對于Band10,K1=774.89W/(m2·sr·μ m),K2=1321.08K。
基于以上的算法,進(jìn)行波段運(yùn)算,獲得兩期數(shù)據(jù)海面的溫度場如圖3、圖4所示。
為了驗(yàn)證溫度反演結(jié)果的可靠性,在2013 年11月15 日衛(wèi)星過境前后一段時(shí)間內(nèi)進(jìn)行了海面溫度測量。測量儀器為JENCO 牌6010 定制版水質(zhì)測量儀,標(biāo)定后儀器測量精度0.1℃。測量時(shí)使用平面定位精度為5~10米的Garmin 60 CSx(GPS)對現(xiàn)場觀測和測量的精確地理位置進(jìn)行定位,以保證測量數(shù)據(jù)和遙感數(shù)據(jù)的位置相對應(yīng)。衛(wèi)星過境前后,在核電溫排水區(qū)域至本底溫度值海域內(nèi)進(jìn)行反復(fù)測量(圖3)。因Landsat-8 熱紅外波段分辨率為100m,測量時(shí)每50~100m間距進(jìn)行一次測溫,共獲取到156個(gè)測溫?cái)?shù)據(jù)。
利用最小二乘法將反演溫度值(SST,sea surface temperature)與實(shí)測值進(jìn)行擬合,驗(yàn)證Landsat-8反演值與海面實(shí)測值之間的關(guān)系。無論是否有相關(guān),都可以用最小二乘法求出最佳的a和b,通過相關(guān)系數(shù)r(通常以其平方值進(jìn)行衡量,0<|r|≤1)衡量線性相關(guān)程度,r 越接近1,線性相關(guān)程度越高,為0時(shí),則不相關(guān);同時(shí),為了衡量實(shí)際值與理論預(yù)測值的偏離程度引進(jìn)了標(biāo)準(zhǔn)誤差σy、殘差(R),和標(biāo)準(zhǔn)殘差R*,這些值偏離越大,相關(guān)性越差,反之,相關(guān)性越好。相關(guān)計(jì)算公式如下:


根據(jù)最小二乘法的計(jì)算方法,將2013年11月15日得到的156 組實(shí)測值分別與Landsat-8 和HJ-1B 數(shù)據(jù)反演的SST值進(jìn)行線性回歸擬合,并對SST值殘差進(jìn)行投點(diǎn),得到了圖4 的擬合結(jié)果。2013 年11 月15日Landsat-8反演獲得的SST值與實(shí)測數(shù)據(jù)擬合關(guān)系式為y=0.7191x+4.9077,擬合后回歸系數(shù)的平方值0.9601,標(biāo)準(zhǔn)誤差為0.37,SST值殘差集中在(-0.8,0.8)的范圍內(nèi),大部分集中在(-0.4,0.4),標(biāo)準(zhǔn)殘差的絕對值也集中在(0.0027,3.3361)區(qū)間內(nèi),大多數(shù)都小于1。從線性關(guān)系的程度和誤差大小上可以反映出溫度反演方法獲得的溫度場數(shù)據(jù)是準(zhǔn)確可信的,海面溫度監(jiān)測結(jié)果是可靠的。
MODIS是美國Terra和Aquar衛(wèi)星的主要傳感器,每天可獲取全球任意地點(diǎn)的影像數(shù)據(jù),含有16 個(gè)熱紅外波段,其中的第31波段和32波段由于對水汽的吸收作用不同,受太陽光反射的影響較微弱,可以用來消除水汽吸收的影響。針對MODIS 數(shù)據(jù)地表溫度的算法,國內(nèi)外很多學(xué)者開展了很多相關(guān)研究,尤以分裂窗算法最為成熟[12],本文不再贅述。經(jīng)過長期的驗(yàn)證,美國NASA 網(wǎng)站針對MODIS 海表溫度二級產(chǎn)品免費(fèi)對外分發(fā)[13],其精度為0.053℃~0.66℃[12-13]。為驗(yàn)證Landsat-8 數(shù)據(jù)反演結(jié)果,以MODIS 數(shù)據(jù)的海表溫度(圖7)為基礎(chǔ),與Landsat-8熱紅外數(shù)據(jù)反演結(jié)果對比分析,進(jìn)行交叉驗(yàn)證[14]。
兩組數(shù)據(jù)過境時(shí)間有一定的時(shí)間間隔(MODIS過境時(shí)間為2013 年11 月15 日13:20,Landsat-8 衛(wèi)星為2013 年11 月15 日10:38),這段時(shí)間主要為落潮末期到漲潮初期,且落潮末期占時(shí)間比例較大(圖8),因?yàn)槁淠╇A段海水比較穩(wěn)定,漲潮時(shí)間較短,大量外海海水還未涌入該海域,海表溫度變化不大。
在兩景溫度場數(shù)據(jù)上隨機(jī)選取30 個(gè)點(diǎn),利用最小二乘法研究兩者之間的關(guān)系。因?yàn)镸ODIS 數(shù)據(jù)空間分辨率為1km,海陸像元混合比Landsat-8 要大,點(diǎn)位選取時(shí)優(yōu)先選取遠(yuǎn)離近岸海域點(diǎn)位。從表中可以看出,偏差最大為1.13℃,最小小于0.1℃。利用最小二乘法進(jìn)行擬合探討兩組數(shù)據(jù)之間的回歸關(guān)系,擬合結(jié)果如圖8 所示。兩組數(shù)據(jù)擬合關(guān)系式為y=0.8447x+2.2576,擬合后回歸系數(shù)的平方為0.8266,標(biāo)準(zhǔn)誤差為0.3918。數(shù)據(jù)擬合結(jié)果表明兩組數(shù)據(jù)之間線性特征非常明顯,具有很好的線性相關(guān)性和一致性,同樣可以反映出Landsat-8 熱紅外波段溫度反演方法獲得的溫度場數(shù)據(jù)是準(zhǔn)確可信的。
2013 年11 月15 日熱紅外溫度場分布圖(圖3、圖4)顯示,核電附近海域溫度分布層次分明,溫度場范圍13.0℃~23.0℃,主要集中在14.5~18.9℃。連島北部海域溫度主要集中在13.7℃~14.5℃。核電南部徐圩港正在建設(shè)過程中,周邊海域溫度主要集中在13.4℃~14.1℃。
2013 年11 月15 日溫度場分布圖(圖3)顯示,核電周邊海域的溫度場明顯受到了溫排水的影響,熱影響強(qiáng)度較高的水體離排水口近,由排水口向外延伸,溫度逐漸降低,到達(dá)環(huán)境本底溫度區(qū)后,變化趨緩。由于兩顆衛(wèi)星過境時(shí)處于落潮潮態(tài),溫排水沿海水落潮方向東北方向擴(kuò)散,取水口處于略高于本底溫度的溫度范圍內(nèi),對核電的冷卻水取水造成了一定的影響,不利于海水置換。
2017 年2 月27 日溫度場分布圖(圖4)顯示,核電周邊海域岸線發(fā)生了顯著的變化,與2013 年11月15 日相比(圖3),取水口處的取水明渠由原來的1.9km 增加到4.5km,排水口處導(dǎo)流堤向南擴(kuò)建了1.5km,南側(cè)徐圩港防波堤,則從2.6km,增加到了6.9km。由于潮汐狀態(tài)相似,2017 年2 月27 日溫度場梯度及空間分布特征與2013 年11 月15 日溫度場相似。由于氣象條件不同,整體溫度較低。核電周邊溫度場范圍6.1℃~14.7℃,主要集中在6.5~13.2℃。徐圩港的半封閉海域溫度主要集中在5.5℃~6.6℃,連島北部海域溫度主要集中在6.6℃~7.2℃。由于衛(wèi)星過境時(shí),獲取的數(shù)據(jù)同樣處于落潮潮態(tài),溫度場明顯受到了潮態(tài)的影響:高溫?zé)崴?2.4℃)漫過排水口導(dǎo)流壩后向東北方向展布,隨著距離增加,海水混合,溫度逐漸降低。因?yàn)榘毒€的阻擋,溫排水被限制在取水明渠南側(cè),此時(shí)溫排水溫度已降至6.8℃;南側(cè)同樣受到了防波堤的影響,溫排水的羽跡被阻斷。高溫海水雖然被阻隔在取水口之外,保護(hù)了取水安全,但岸線的變化,使田灣核電站排水口處于一個(gè)由兩側(cè)海工建筑環(huán)抱下的人工海灣的灣底,改變了海域流場情況,影響了溫排水的展布形態(tài)和規(guī)模。
溫排水影響區(qū)域?yàn)楹穗娬九潘谔幐哂诒镜诇囟?.1℃以上區(qū)域。綜合考慮周邊海域(剔除溫排水影響區(qū)域后海域)平均海面溫度作為本底溫度。在確定了本底溫度之后,將遙感反演的溫度場數(shù)據(jù)整體扣除本底溫度,以獲取核電溫排水形成的溫度場的熱影響數(shù)據(jù)。提取核電站熱影響區(qū),分別劃分出9、10 個(gè)等級,并分別進(jìn)行編碼(見圖9、10),分類統(tǒng)計(jì)各個(gè)等級的面積,采用0.1℃、0.5℃、1.0℃(高于本底水溫溫度)等提取核電站溫排水的熱影響分布信息。根據(jù)各級水溫水體所占的像元數(shù),計(jì)算不同級別水溫分布面積(圖11)。
熱影響編碼圖顯示,核電溫排水明顯影響了核電附近海域的溫度場,由排水口向外延伸,熱影響溫度逐漸降低,到達(dá)本底溫度后,溫度穩(wěn)定。相較于2013 年11 月15 日反演的結(jié)果,2017 年2 月27 日的獲取的結(jié)果,總體上,各級別的面積均有擴(kuò)大,具有溫升級別越高,面積變化越大的趨勢,0.1℃以上的溫升面積增加了18.68%,1℃以上溫升面積增加了約51.21%,并且出現(xiàn)了8℃以上的溫升。因?yàn)閮山M數(shù)據(jù)獲取時(shí)核電的運(yùn)行工況,所處的季節(jié)及潮汐狀態(tài)類似,推斷面積增大的原因?yàn)楣こ探ㄔO(shè)阻礙溫排水?dāng)U散造成的。
周邊海域岸線的變化,改變了海域的流場情況。工程未建成時(shí),漲潮海水從東北方向涌入,沿岸向南東方向流動,落潮時(shí),海水從東北方向流出;工程建成后,徐圩港防波堤阻隔了排水口向南的順岸水流,漲落潮時(shí)形成了向岸、離岸的往復(fù)流,流速減慢,取水明渠延伸,在另一面限制了海水的流動。受此影響,溫排水海域漲落潮流速都有不同程度的減小,流向發(fā)生了偏轉(zhuǎn),不利于溫排水的擴(kuò)散,造成了溫排水影響的面積增大。
本文基于2013年、2017年兩景Landsat-8數(shù)據(jù),完成了田灣核電周邊海域岸線變遷調(diào)查,對其熱紅外波段進(jìn)行了溫度反演,獲取了相似潮汐、氣候條件下核電周邊海域溫度場分布情況,并對2013 年11 月15 日數(shù)據(jù)過境前后進(jìn)行了海溫測量。通過解譯、對比分析得出:
1.在2013-2017年期間,田灣核電周邊岸線發(fā)生了顯著的變化,進(jìn)水口處的取水明渠由原來的1.9km增加到4.5km,排水口處導(dǎo)流堤向南擴(kuò)建了1.5km,南側(cè)徐圩港防波堤,則從2.6km,增加到了6.9km,核電排水口處于兩側(cè)防波堤環(huán)抱的人工海灣的灣底。
2.通過對2013年11月15日Landsat-8熱紅外波段溫度反演結(jié)果與衛(wèi)星過境前后的海上測溫?cái)?shù)據(jù)進(jìn)行擬合、與同日過境的MODIS 數(shù)據(jù)進(jìn)行交叉驗(yàn)證,顯示溫度反演的結(jié)果與海上測溫?cái)?shù)據(jù)和MODIS 數(shù)據(jù)具有很強(qiáng)的相關(guān)性,證明溫度反演結(jié)果是可靠的,但由于測量不是完全同步,數(shù)據(jù)之間存在著誤差。誤差存在的原因除與儀器的精度、人為操作、溫度反演的方法有關(guān)外,還與測量的連續(xù)性和衛(wèi)星數(shù)據(jù)的瞬時(shí)性有很大的關(guān)系。
3.不同時(shí)相遙感數(shù)據(jù)顯示,工程建設(shè)雖然保護(hù)了溫排水的取水安全,但是卻影響了溫排水展布形態(tài)和規(guī)模。核電站需要時(shí)刻關(guān)注,根據(jù)建設(shè)進(jìn)度變化,必要時(shí)進(jìn)行溫排水影響分析,做出相應(yīng)的對策。
4.Landsat-8 數(shù)據(jù)能滿足核電溫排水的監(jiān)測需求,為評估核電站溫排水對其周邊海域溫度環(huán)境的影響提供了迅速便捷的手段。
5.本文僅對冬季落潮落急時(shí)刻的溫度場進(jìn)行了分析,不同季節(jié)、不同潮態(tài)下的核電附近海域的遙感監(jiān)測還有待進(jìn)一步研究。
致謝:文中Landsat-8 數(shù)據(jù)美國USGS網(wǎng)站提供,在此表示誠摯的謝意。

圖1 田灣核電站周邊海域遙感影像圖(時(shí)相:2017-02-27,波段組合:6 5 2)Fig.1 Remote sensing image of the sea area around Tianwan Nuclear Power Station(Time:2017-02-27,band combination:652)

圖2 田灣核電周邊岸線解譯結(jié)果圖Fig.2 Interpretation results of the coastline around Tianwan Nuclear Power Station

圖3 2013年11月15日熱紅外溫度場圖Fig.3 Distribution of Thermal Infrared Temperature on Nov.15,2013

圖4 2017年2月27日熱紅外溫度場圖Fig.4 Temperature Distribution of Thermal Infrared on Feb.27,2017

圖5 2013年11月15日海面實(shí)測點(diǎn)位圖Fig.5 Distribution of Sea Surface Temperature on November 15,2013

圖6 2013年11月15日海上實(shí)測值與反演SST值線性擬合圖和殘差投點(diǎn)圖Fig.6 Linear Fitting and Residual Point Map of Measured values and Inverted SST Values at Sea on Nov.15,2013

圖7 2013年11月15日MODIS熱紅外溫度場圖Fig.7 Temperature Distribution Map of MODIS Thermal Infrared on November 15,2013

圖8 2013年11月15日潮汐狀態(tài)變化示意圖(為Landsat-8過境時(shí)潮汐狀態(tài),為Aquar過境時(shí)潮汐狀態(tài))Fig.8 Diagram of tidal state changes on November 15,2013(Tidal state of Landsat-8,Tidal state of Aquar)

圖8 2013年11月15日Landsat-8與MODIS反演結(jié)果數(shù)據(jù)擬合圖Fig.8 Fitting of Landsat-8 and MODIS Inversion Results on November 15,2013

圖9 2013年11月15日熱影響編碼圖Fig.9 Heat Impact Coding Chart,15 November 2013

圖10 2017年2月27日熱影響編碼圖Fig.10 Heat Impact Coding Chart,27 February 2017

圖11 不同時(shí)相熱影響面積對比圖(單位:km2)Fig.11 Contrast chart of heat affected area at different time phases