大连港位于辽东半岛南端的大连境内，始建于1899年，距今已有百余年的历史。作为中国“一带一路”重点建设的15个港口城市之一，长期以来承担着东北原油和成品油转运的主要任务。在原油装卸的过程中，由于原油泄漏、大型船舶运输等因素，导致码头区域内的水域遭到严重的石油污染，特别是2010年7月16日(简称“7.16”溢油事故)，中石油大连大孤山新港码头的一个储油罐输油管线发生起火爆炸事故，据估算有1500 t原油进入海洋，随着海流和潮汐的运动，不可避免地为大连港海域带来一定的影响[1]。另外据不完全统计，2010—2014年期间，大连港发生油泄漏事件6次以上，加之历史上的各种油泄漏事件的发生，对大连港水环境已经造成了危害。掌握大连港海域油含量空间分布特征，可为海洋环境监管与污染防治、海洋生态保护与修复提供重要的参考依据。
遥感具有大面积、快速、动态、低成本获取区域信息的优势，目前，利用遥感图像来提取水环境石油类污染信息的相关研究已经取得了长足的进展，主要集中在水体石油类物质固有光学特性[2-6]、表观光学特征[7-9]、荧光特性[10-13]、石油含量提取[14-19]、海洋溢油监测[20-21]等方面。黄妙芬等[22]利用现场实测数据和Landsat 8/OLI遥感数据，采用归一化处理方法，发现含油水样与无油水样在绿色、红色和近红外波段的遥感反射比存在明显的差异，由此提出了水体石油污染归一化遥感反射比指数NDPRI(Normalized difference petroleum remote sensing reflectance index)，该指数可用来判断水体是否存在油污染，并进一步建立了利用NDPRI来反演水中石油含量的遥感模式，为利用遥感技术进行海域油含量空间分布特征分析奠定了基础。
本研究中，利用2018年8月25—27日在辽宁大连港海域连续3 d的现场测量数据，以及Landsat 8卫星数据，确定油含量随时间变化的参数化模型，修正利用NDPRI指数反演石油污染浓度的遥感模型，使其适用于大连港海域，并对大连港海域油含量的时空特征进行了分析，旨在为海洋环境油污染的监管与防治提供参考依据。
1 方法
1.1 试验区及站位设置
图1为大连港海域空间分布图。2008年8月，《大连港总体规划》中提出了大连港以“一岛三湾”(大孤山半岛、大窑湾、鲇鱼湾、大连湾)综合运输港区的发展布局。目前，石化产业主要分布在大孤山半岛，有西太平洋石油化工公司(简称西太公司)、逸盛大化石化有限公司、福佳大化石油化工公司、大连港油品码头、昆仑能源大连LNG码头等。石油化工生产及油轮运输不可避免带来油污染，这使得大连港成为含油水体光学特性和油含量变化特征研究的天然试验基地。
试验时间为2018年8月25—27日，共3 d，试验设3个试验点。8月25日的观测点位于“7.16”石油管道爆炸事故点(图1的A点)，当时事故发生的处理办法之一是用油毡吸附原油成块沉底，2018年8月20号正好是台风“温比亚”经过大连，狂浪把海底的油物质搅动起来，搅动后水体会受到油的影响，加上天气热水温高，沉底的油块挥发，使得水体含油，该站点水深21.5 m；8月26日的观测站点位于大连港的航道上(图1的B点)，船舶过往对水体油含量影响较大，该站点水深24.8 m；8月27日的观测点位于大山岛和小山岛之间水域的东侧(图1的C点)，与上两个站点相比，船舶经过该站点的频率小，受到轮船扰动小，随着涨潮和落潮，流速和流向有明显改变，该站点水深33.0 m。
表1为3个站点各观测时刻对应的时间和涨落潮过程。对于油含量分析水样，3 d观测时刻总数为30个，每个观测时刻观测深度分别为0、3、5、10、15 m，一共获取150个数据，其中A站点9:00由于其他配套仪器出现问题缺测，B站点15:00由于流速过大缺测。
1.2 方法
1.2.1 现场数据的采集 采用便携式紫外荧光测油仪(美国特纳TD-500D)测量油含量，所使用的监测标准与中国的《海洋监测规范》(GB17378.3—1998)相符。使用仪器测量之前，需先将正己烷和石油类标准物进行混合制备标准油样，以对设备的固体标样进行校准，然后在现场使用该固体标样进行标定，油标样来自国家海洋环境监测中心。
表1 观测时间及潮汐变化记录
Tab.1 Record sheet of observation time and tidal changes
观测时间time潮汐过程tidalprocess2018-08-252018-08-262018-08-27siteAsiteBsiteC7:00涨潮涨潮涨潮8:00涨潮(缺测)涨潮涨潮9:00涨潮涨潮涨潮10:00落潮平潮涨潮11:00落潮落潮平潮12:00落潮落潮落潮13:00落潮落潮落潮14:00落潮落潮落潮15:00落潮落潮(缺测)落潮16:00停潮落潮落潮17:00涨潮涨潮
采用水面之上法测量水体光谱，测量仪器为可见光-近红外地物光谱仪(美国ASD FieldSpec 3，350～2500 nm)，参考板采用30%反射率的标准板。
采用数字直读式海流仪(SLC9-2DV)测量流速、流向。测量船的吃水深度为3 m，为了保证数据的准确性，海流计布放深度为6 m。
本研究中研究区域没有潮汐观测站，试验期间也没有相应的潮汐测量装置，故采用由海事服务网http://ocean.cnss.com.cn提供的潮汐数据。试验站点的位置与辽宁省大连市大窑湾(南大圈)港口距离最短(图1的黄色圆圈点)，因而采用大窑湾(南大圈)的潮汐数据。
1.2.2 卫星数据 本研究中主要用到的Landsat 8卫星数据，数据来源于该卫星的官方数据发布网站(https://earthexplorer.usgs.gov/)。Landsat 8搭载的陆地成像仪(Operational Land Imager，OLI)在可见光近红外及短波红外设置了7个波段，波长范围分别为B1，433～453 nm；B2，450～515 nm；B3，525～600 nm；B4，630～680 nm；B5，845～885 nm；B6 SWIR1，1570～1650 nm；B7 SWIR2，2110～2290 nm；另外设置了针对大气校正的波段9，波长范围为B9，1360～1380 nm；空间分辨率均为30 m，轨道周期/重访周期为16 d。
2 结果与分析
2.1 潮汐特征分析
从图2可见：根据日变化特征来看，大连港主要为不正规半日混合潮港，一天内有两个涨落潮过程；一天中两次涨落潮的涨潮时时长和最高潮高及最低潮高不等，涨潮时要长于落潮时，两次高潮的潮高位置不同，但两个低潮位置相差不是很明显。
图1 大连港空间分布图及试验点位置
Fig.1 Spatial distribution map and the sampling location at Dalian port
图2 大连港潮汐表(基于平均海平面)
Fig.2 Tidal table in Dalian port based on average sea level
分析图2还可以看出，A站点、B站点和C站点的最大高潮与次大高潮潮高分别相差0.67、0.63、0.58 m，最低潮与次低潮潮高分别相差为0.05、0.17、0.22 m，可见，连续3 d，最大高潮与次大高潮潮高之差在逐渐减小，而最低潮与次低潮潮高之差在逐渐增大。
2.2 海流特征分析
图3为观测期的流速和流向，图中0°表示正北方向，线条的长度表示流速的大小，箭头所指方向为流向。
分析图3可以看出，所有站点处于不同涨落潮时刻，海流方向和流速不同。A站点受潮汐影响明显，涨潮为西向流，落潮为东向流，流速为27～65 cm/s。B站点涨潮为东南向流，落潮为西南向流，总体看南向流为主，该站点位于航路上，海流的测量受到过往船舶的影响比较大，流速为22～75 cm/s。C站点位于大三山岛和二三山岛之间的偏东方向，两岛之间形成的水道对海流产生一定程度的影响，总体看涨潮时为南向流，落潮时为西向流，流速为21～85 cm/s。
图3 流速与流向(观测深度6 m)
Fig.3 Current velocity and direction(observation at water depth of 6 m)
2.3 水体油含量随深度的变化特征
图4为8月25—27日3 d油含量观测值。分析图4可以看出：
(1)A站点呈U型变化特征，涨潮时刻油含量呈现下降的趋势，因为涨潮时海水从东边过来，外海水会冲淡该站点本身具有的油含量；落潮后海水从西边过来，带来西面的高油含量水，使得油含量有明显的增加。
(2)B站点在涨潮时刻不同深度的变化比较平稳，11:00落潮后，0、3、5、10 m深度的油含量急剧上升，15 m深度变化不大，到16:00停潮时油含量又开始下降。
(3)C站点从7：00—11:00平稳变化，油含量相对比较高，这是由于涨潮时海水为从西北大连湾海域过来，带来一定的油污水，使得油含量相对较高，11:00—13:00出现一个陡降后，13:00—16:00油含量又是平稳变化，但处于一个低值范围，这主要是因落潮时为西向流，海水从外海过来，加上这个区域远离油污水的影响，故油含量比较低。
(4)从浓度的变化值来看，A站点油含量为0.7～4.1 mg/L，B站点油含量为1.0～9.3 mg/L，C站点油含量为0.2～7.8 mg/L，3个观测站点中，B站点油含量全天处于高值状态。
对所有油含量随时间变化曲线进行拟合，结果如图5所示，发现可以用六次方的多项式来拟合曲线,即
Cp=az6+bz5+cz4+dz3+ez2+fz+g。
(1)
其中：Cp为水中石油物质浓度(mg/L)；z为深度(m)；a(L/mg·m-6)、b(L/mg·m-5)、c(L/mg·m-4)、d(L/mg·m-3)、e(L/mg·m-2)、f(L/mg·m-1)、g(L/mg)为系数。
图4 不同深度的油含量随时间的变化
Fig.4 Variation in petroleum concentration at different depths with time
图5 5 m深度油含量随时间的变化
Fig.5 Curve fitting of petroleum concentration changing with time at 5 m depth
表2为六次方多项式拟合的R2汇总。分析表2可以看出：A站点拟合的R2为0.70～0.92，B站点拟合的R2为0.89～0.98，C站点拟合的R2为0.73～0.94，均有较好的拟合效果；5 m和10 m水深拟合效果都较好，均在0.83～0.98之间。这表明，油含量随时间的变化可以用6次方多项式函数来描述。
2.4 油含量遥感反演
从网站https://earthexplorer.usgs.gov/下载2018年8月25日在大连港过境的卫星数据Landsat-8/OLI，所下载的遥感图像属于一级产品，其像元数值用灰度值(Digital Number，DN)表示，不具有物理意义，因而需要对其进行辐射定标。将像元DN值转为入瞳光谱辐亮度，转换公式为[23]。
表2 六次方多项式拟合的R2汇总表
Tab.2 Summary of R2 fitted by the sixth power polynomial
深度depth/mA站点siteAB站点siteBC站点siteC00.820.930.7330.700.910.7550.860.980.94100.830.890.85150.920.940.86
Lλ=ML×Qcal+AL
(2)
其中：Lλ为入瞳光谱辐亮度[W/(m2·sr·μm)]；ML为各波段的增益量[W/(m2·sr·μm)]；AL为各波段的偏移量[W/(m2·sr·μm)]；Qcal为Landsat 8/OLI遥感数据产品中的灰度值(DN值，无量纲)；ML和AL可从相应的遥感数据产品的元数据文件中获取。
Lλ可由下式表示[24]:
Lλ=Lr(λ)+La(λ)+t0(λ,θ0)Lw(λ)。
(3)
其中：Lλ含义同上；λ为波长(μm)；Lr(λ)为空气分子瑞利散射的辐亮度[W/(m2·sr·μm)]；La(λ)为大气气溶胶米氏散射的辐亮度[W/(m2·sr·μm)]；Lw(λ)为离水辐亮度[W/(m2·sr·μm)]；t0(λ,θ0)为太阳方向的大气漫射透射率(无量纲)；θ0为太阳天顶角(°)。
ENVI软件提供了FLAASH大气校正模块，该模块适用于对高光谱遥感数据和多光谱数据(包括陆地资源卫星等)进行大气校正。基于FLAASH大气校正模块，结合Landsat 8/OLI波段9数据进行大气校正，可得到Lr(λ)、La(λ)和t0(λ,θ0)。将公式(2)计算得到的Lλ和这些大气参数代入公式(3)，就可得到经过大气校正后的离水辐亮度Lw(λ)。
在海洋遥感中，遥感反射率定义为[25]
Rrs=Lw(λ)/Ed(0+)。
(4)
其中：Rrs为遥感反射比(sr-1)；Ed(0+)为水表面上总的入射辐照度,可由大气辅助数据或模型反演而得。根据公式(4)可计算出Landsat 8/OLI传感器过境时的遥感反射比。
水体边界提取采用二次改进的归一化水体指数(Second modified normalized difference water index，SMNDWI)[26]，其计算公式如下：
SMNDWI=(B7-B3)/(B7+B3)。
(5)
其中：B7和 B3分别对应的是短波红外波段(SWIR 2)和绿光波段，波段范围分别为2110～2290 nm和525～600 nm，SMNDWI值大于-0.1时作为水体处理。
在得到遥感反射比Rrs后需要进一步进行归一化处理，归一化处理的表达式为
(6)
其中：Rrs′为归一化后的遥感反射率(无量纲)；Rrs,i为第i波段对应的遥感反射比(sr-1)；本研究中采用的卫星数据是Landsat 8/OLI，所以n取值5。将在每个波长处的遥感反射比除以所研究波长范围的所有遥感反射比的总和，这样将原始光谱范围变为(0，1)范围，同时通过比值，把有量纲表达式的遥感反射比(单位sr-1)变为无量纲表达式。
水体石油污染归一化遥感反射比指数NDPRI(Normalized difference petroleum remote senising reflectance index)，其计算公式[22]为
NDPRI=Rrs, G′-Rrs,NIR′/Rrs, G′+Rrs,NIR′。
(7)
其中：NDPRI为归一化遥感反射比指数(无量纲)；Rrs, G′为绿光波段G(green band)的归一化遥感反射比(无量纲)，即Landsat 8/OLI对应的波段3(B3)的归一化遥感反射比；Rrs,NIR′为近红外波段NIR(near infrared band)的归一化遥感反射比(无量纲)，即Landsat 8/OLI对应的波段5(B5)的归一化遥感反射比。试验表明，在本研究区域，对应像元的NDPRI值小于0.49时，该像元有油污染。
文献[22]中给出了NDPRI反演水体油含量的遥感模式，即
Cp=-7.4629ln(x)-4.6031。
(8)
其中：Cp为水体油含量(mg/L)；x为NDPRI值(无量纲)。直接利用公式(8)反演大连港水体油含量时发现反演的值偏高，主要原因是公式(8)的参数是基于盘锦辽河油田污水处理厂和配比试验油含量数据而确定，油含量值属于高值范围。因而在将该模式用于本试验区时需要对参数进行修正，利用现场光谱数据模拟Landsat 8/OLI波段，对参数进行了修正(篇幅所限，不在此赘述过程)，得到适合大连港的油含量反演模型为
Cp=-6.019ln(x)-9.89。
(9)
2.5 大连港油含量的空间分布特征
将公式(9)应用于2018年8月25日过境的LANDSAT 8卫星数据，反演出大连港海水表层的油含量，结果如图6所示。对比图4和图6可以看出：在8月25日10:34(Landsat-8卫星过境)的时候，A站点的遥感反演值在2.7 mg/L左右，与当天11：00的实测值(2.2 mg/L)非常接近；此时B站点的遥感反演值在3.9 mg/L左右，与8月26日的实测值4.2 mg/L比较接近；C站点遥感反演值在3.5 mg/L左右，与8月27日的实测值3.8 mg/L比较接近。这表明，公式(9)适用于大连港海域，有较好的油含量反演精度。
从空间分布上来看，大连港海域油含量分布具有明显的差异，浓度值的变化范围为0.4～5.0 mg/L；紧贴沿岸的海域油含量相对较低，鲇鱼湾东边及三山岛的东边、北边呈现高值区，而大窑湾、三山岛西侧及A站点的南边为低值区。在大窑湾湾内没有石化工业，油含量相对比较低，在1.5 mg/L以下；鲇鱼湾东边及三山岛的东边、北边都属于油含量高值区，可达3.5 mg/L以上；在三山岛西侧及A站点的南边出现低值区域。A站点在卫星过境时，处于涨潮阶段，为西向流，海水从外海过来，油含量相对较低；B站点为偏西南向流，该站点位于航道上，观测期内油浓度一直都比较高；C站点处于大连湾湾口，北边受油码头和石化生产污水影响较大，而且大量的倾废物质会流经这个站点，在观测阶段经常能看到成片的油花。
图6 大连港海域油含量空间分布图
Fig.6 Spatial distribution of petroleum concentration in the sea area of Dalian port
3 讨论
掌握水体油含量随潮汐的时间变化及水平空间分布特征，对于利用水色遥感反演水体油含量模型精度的提高、研究含油水体辐射传射特性具有重要的意义。大连港周边由于石油化工生产企业的存在，加上历年所发生的溢油事故，不可避免地给港口水域环境带来一定油污染。
通过分析大连港海域油含量随潮汐变化数据，初步掌握了油含量的时间分布特征并确定了参数化模型。大连港油含量随时间变化比较复杂，与潮汐和海流有着密切的关系，但可以用六次方多项式来描述，这为利用潮汐数据推断油含量的时间变化特征提供了一种有效的方式。通过修正水体石油污染归一化遥感反射比指数NDPRI反演水中石油含量遥感模型的参数，确定了适用于大连港的油含量遥感反演模型，并计算出大连港海域的油含量。通过对油含量的空间分布图及实测数值进行比对分析，表明修正后的模式所反演的结果与实际情况比较吻合。
本研究中提出的油含量随潮汐的时间变化模型及利用遥感技术提取油含量空间分布信息的模型，有助于实时掌握大连港油含量变化特征，研究成果可为“陆海统筹、以海定陆”的环境治理思想提供科学、及时、有效的数据基础与技术支撑，还可为水体石油类污染监测提供技术手段。
本研究中数据取自大连港海域，所得结论和所建立的遥感模型是否具有普适性，还需要更多的观测数据来支撑。发展普适模型而无需根据不同区域的测量数据来修正相关系数，这是今后需要进一步研究的工作。
致谢：感谢所有参加实验的人员！
[1] 兰国新,马龙,陈冠文.环境小卫星在大连新港溢油事故监测中的应用[J].中国航海,2011,34(4):89-93.
[2] 黄妙芬,宋庆君,唐军武,等.石油类污染水体后向散射特性分析—以辽宁省盘锦市双台子河和绕阳河为例[J].海洋学报,2009,31(3):12-20.
[3] 黄妙芬,唐军武,宋庆君.石油类污染水体吸收特性分析[J].遥感学报,2010,14(1):131-147.
[4] 宋庆君,黄妙芬,唐军武,等.水中石油类含量对后向散射系数光谱的影响[J].光谱学与光谱分析,2010,30(9):2438-2442.
[5] Otremba Z,Król T.Modeling of the crude oil suspension impact on inherent optical parameters of coastal seawater[J].Polish Journal of Environmental Studies,2002,11(4):407-411.
[6] 黄妙芬,邢旭峰,宋庆君,等.油砂混合水体后向散射系数光谱贡献分离算法Ⅰ:实验理论[J].光谱学与光谱分析,2017,37(1):205-211.
[7] Król T,Stelmaszewski A,Freda W.Variability in the optical properties of a crude oil - seawater emulsion[J].Oceanologia,2006,48(S1):203-211.
[8] Otremba Z,Piskozub J.Modeling the remotely sensed optical contrast caused by oil suspended in the sea water column[J].Optics Express,2003,11(1):2-6.
[9] Huang Miaofen,Song Qingjun,Xing Xufeng,et al.An experimental study on the feasibility of monitoring petroleum-polluted waters with remote sensing technology[C]//Proceedings of the Sixth International Symposium on Digital Earth:Data Processing and Applications.Beijing:SPIE,2010:78411O.
[10] 慈兴华,向巧玲,陈方鸿,等.定量荧光分析技术在原油性质判别方面的应用探讨—以胜利油区为例[J].石油实验地质,2004,26(1):100-102.
[11] 黄妙芬,宋庆君,邢旭峰,等.石油污染水体荧光图谱特征分析[J].光谱学与光谱分析,2014,34(9):2466-2471.
[12] 吴静,谢超波,曹知平,等.炼油废水的荧光指纹特征[J].光谱学与光谱分析,2012,32(2):415-419.
[13] 黄妙芬,王忠林,刘远,等.辽东湾海域有色可溶性有机物荧光特征分析[J].大连海洋大学学报,2018,33(2):137-144.
[14] Huang Miaofen,Xing Xufeng,Song Qingjun,et al.A new algorithm of retrieving a petroleum substances absorption coefficient in sea water based on a remote sensing image[J].Acta Oceanologica Sinica,2016,35(11):97-104.
[15] Huang Miaofen,Song Qingjun,Xing Xufeng,et al.A biooptical model of retrieving petroleum concentration in seawater[J].Acta Oceanologica Sinica,2014,33(5):81-85.
[16] 黄妙芬,宋庆君,简伟军,等.水体石油类信息遥感提取模式研究[J].海洋学报,2012,34(5):74-80.
[17] Otremba Z.Influence of oil dispersed in seawater on the bi-directional reflectance distribution function(BRDF)[J].Optica Applicata,2005,35(1):99-109.
[18] Sun Shaojie,Hu Chuanmin.Sun glint requirement for the remote detection of surface oil films[J].Geophysical Research Letters,2016,43(1):309-316.
[19] Hu Chuanmin,Li Xiaofeng,Pichel W G,et al.Detection of natural oil slicks in the NW Gulf of Mexico using MODIS imagery[J].Geophysical Research Letters,2009,36(1):L01604.
[20] 陆应诚,胡传民,孙绍杰,等.海洋溢油与烃渗漏的光学遥感研究进展[J].遥感学报,2016,20(5):1259-1269.
[21] Zhong Zhixia,You Fengqi.Oil spill response planning with consideration of physicochemical evolution of the oil slick:a multiobjective optimization approach[J].Computers & Chemical Engineering,2011,35(8):1614-1630.
[22] 黄妙芬,宋庆君,陈利搏,等.基于归一化遥感反射比反演水中石油含量模式研究[J].海洋技术学报,2015,34(1):1-9.
[23] Using the USGS landsat 8 product[EB/OL].(2013-07-13).https://www.scoop.it/topic/remote-sensing.
[24] 黄妙芬,牛生丽,孙中平,等.环境一号卫星CCD相机水体信息采集特性分析[J].遥感信息,2010(4):68-75.
[25] 刘良明.卫星海洋遥感导论[M].武汉:武汉大学出版社,2005:69.
[26] 李明,郑小慎.二次改进归一化差异水体指数(SMNDWI)提取海岸线的研究(英文)[J].海洋通报,2016,18(2):15-27.