基于PIE-Engine的太湖水域蓝藻信息反演

1 引言  

太湖是我国第三大淡水湖，地跨江浙两省，水域面积达2338km2，平均水深2m左右。上世纪80年代以来，太湖水体营养化愈加严重，频繁爆发蓝藻水华现象即为一种非常典型的二类水体（Case 2 Waters，光学性质变化受浮游植物及其附属物、外生的粒子和外生有色可溶有机物等其他物质影响的水体）。悬浮物（Suspended Solids，悬浮在水中的固体物质，包括不溶于水的无机物、有机物及泥砂、黏土、微生物等）作为造成水体浑浊的主要影响因素，是衡量水体污染程度的重要指标。利用遥感技术准确获取区域面状水域悬浮物浓度信息，是遥感水质参数监测的一项重要任务。传统的利用遥感手段监测水质的方法大多基于时序遥感影像及地表水质监测数据构建地物信息模型来完成，但受限于本地设备存储、运算能力的限制，当前的监测尺度已逐渐无法满足日益提升的大尺度、长时序遥感应用研究需求，且需要消耗大量人力物力和时间。

当前我们身处大数据时代，云端计算技术蓬勃发展，航天宏图依托行业多年技术积累独立自主研发了安全可控的开放式云计算平台：PIE-Engine（Pixel Information Expert Engine，像素专家引擎），致力于加速我国遥感云计算平台发展进程，实现了遥感数据按需获取、运算以及专题信息聚焦服务，以满足对地观测数据获取能力飞速增长带来的信息高效化处理和服务需求。

2 .反演太湖水域2013-2020年蓝藻密度与浑浊度

基于PIE-Engine获取2013-2020年覆盖太湖水域的Landsat 8 Surface Reflectance遥感影像数据，Landsat 8 SR数据集是OLI/TIRS传感器数据经过大气校正的表面反射率数据，分辨率为30米，重放周期为16天。影像包含5个可见光波段、1个近红外（NIR）波段和2个短波红外（SWIR）波段，经过正射校正后的表面反射率，以及两个经过正射校正后的亮度温度热红外（TIR）波段。

相关研究指出，Landsat 8数据的近红外波段和红光波段分别与蓝藻密度和浑浊度具有较高的相关系数，但相关程度都不高，主要原因有：① 水华与水草区域对水质反演的影响；② 水表面光滑度和微波随时间和空间的干扰；③ 陆地光谱对近岸边水体的干扰；④ 水体中较高浓度悬浮物的高反射干扰。在去除影响区域和异常数据，并利用普通水体与不同波段的组合构建相关性模型后发现，蓝藻密度与B5/B4以及浑浊度与(B4-B3)/(B4+B3)的波段组合具有最高的相关系数。相关性模型如下：

其中，B3、B4、B5分别为绿光、红光以及近红外波段的反射率。

3 太湖水域2013-2020年蓝藻密度 与浑浊度动态变化

3.1蓝藻密度

示例代码：

1. //输入太湖边界矢量  

2. var taihu = pie.FeatureCollection('/PersonalGDB/Vector/TaiHuShuiYu.geojson')  

3.                .first()  

4.                .geometry();  

5. Map.centerObject(taihu, 8);  

6. //加载Landsat 8 SR  

7. var Imgcol = pie.ImageCollection("LC08/01/T1_SR")  

8.             .filterDate("2013-01-01", "2021-01-01")  

9.             .filterBounds(taihu)  

10.             .filter(pie.Filter.lte('cloudCover', 5));  

11. //print(Imgcol.size());  

12. //设置图层显示属性  

13. var visalg = {min: 200, max: 50000,  

14.               palette: ['040274','040281','0502a3','0502b8','0502ce','0502e6',  

15.                         '0602ff','235cb1','307ef3','269db1','30c8e2','32d3ef',  

16.                         '3be285','3ff38f','86e26f','3ae237','b5e22e','d6e21f',  

17.                         'fff705','ffd611','ffb613','ff8b13','ff6e08','ff500d',  

18.                         'ff0000','de0101','c21301','a71001','911003']};  

19. var count = Imgcol.size().getInfo();  

20. var dates = [];  

21. var alg_images = [];  

22. for(var index = 0; index < count; index ++){  

23.     var image = Imgcol.getAt(index).select(["B4", "B5"]).clip(taihu);  

24.     var date =image.get("date").getInfo();  

25.     var b4 = image.select("B4");  

26.     var b5 = image.select("B5");  

27.     //计算蓝藻密度  

28.     var alg = ((((b5.divide(b4)).power(2)).multiply(1352)).subtract  

29.                 ((b5.divide(b4)).multiply(159.08))).add(192.87);  

30.     //print("太湖水域蓝藻密度", date, alg);  

31.     Map.addLayer(alg, visalg, date, true);  

32.     var alg_mean = alg.reduceRegion(pie.Reducer.mean(), taihu, 1);  

33.     //print("太湖水域蓝藻密度平均值", date, alg_mean);  

34.     var alg_max = alg.reduceRegion(pie.Reducer.max(), taihu, 1);  

35.     //print("太湖水域蓝藻密度最大值", date, alg_max);  

36.     var alg_min = alg.reduceRegion(pie.Reducer.min(), taihu, 1);  

37.     //print("太湖水域蓝藻密度最小值", date, alg_min);  

38.     dates.push(date);  

39.     alg_images.push(alg_mean);  

40. };  

41. // 设置图表属性  

42. var line_a = {title: '太湖水域蓝藻密度动态变化',  

43.               legend: ['蓝藻密度'],  

44.               xAxisName: "日期",  

45.               yAxisName: "蓝藻密度（万/L）",  

46.               chartType: "line",  

47.               yMin: 0,  

48.               yMax: 2500,  

49.               smooth: true};  

50. // 显示折线图  

51. ChartImage(alg_images, dates, line_a);  

52. //动画显示  

53. Map.playLayersAnimation(dates, 0.5, 100);  

54. // 图例  

55. var data = {title: "蓝藻密度",  

56.             colors: ['#040274','#040281','#0502a3','#0502b8','#0502ce','#0502e6',

57.                      '#0602ff','#235cb1','#307ef3','#269db1','#30c8e2','#32d3ef',

58.                      '#3be285','#3ff38f','#86e26f','#3ae237','#b5e22e','#d6e21f',

59.                      '#fff705','#ffd611','#ffb613','#ff8b13','#ff6e08','#ff500d',

60.                      '#ff0000','#de0101','#c21301','#a71001','#911003'],  

61.             labels: ["200（万个/L）", "2000（万个/L）"],  

62.             step: 30};  

63. var style = {  

64.     top: "85%",  

65.     left: "40%",  

66.     width: "350px",  

67.     height: "70px"  

68. };  

69. var legend = ui.Legend(data, style);  

70. Map.addUI(legend);  

3.2 浑浊度

示例代码：

1.//输入太湖边界矢量  

2.var taihu = pie.FeatureCollection('/PersonalGDB/Vector/TaiHuShuiYu.geojson')  

3.               .first()  

4.               .geometry();  

5.Map.centerObject(taihu, 8);  

6.//加载Landsat 8 SR  

7.var Imgcol = pie.ImageCollection("LC08/01/T1_SR")  

8.            .filterDate("2013-01-01", "2021-01-01")  

9.            .filterBounds(taihu)  

10.            .filter(pie.Filter.lte('cloudCover', 5));  

11.//print(Imgcol.size());  

12.//设置水域预览颜色组合  

13.var vistbd = {min: 100, max: 800,  

14.              palette: ['040274','040281','0502a3','0502b8','0502ce','0502e6',  

15.                        '0602ff','235cb1','307ef3','269db1','30c8e2','32d3ef',  

16.                        '3be285','3ff38f','86e26f','3ae237','b5e22e','d6e21f',  

17.                        'fff705','ffd611','ffb613','ff8b13','ff6e08','ff500d',  

18.                        'ff0000','de0101','c21301','a71001','911003']};  

19.var count = Imgcol.size().getInfo();  

20.var dates = [];  

21.var tbd_images = [];  

22.for(var index = 0; index < count; index ++){  

23.    var image = Imgcol.getAt(index).select(["B4", "B3"]).clip(taihu);  

24.    var date = image.get("date").getInfo();  

25.    var b3 = image.select("B3");  

26.    var b4 = image.select("B4");  

27.    // 计算太湖水域浑浊度  

28.    var tbd = (((b4.subtract(b3)).divide(b4.add(b3)).power(2)).multiply(3117.4))

29.           .add(((b4.subtract(b3)).divide(b4.add(b3)).power(2)).multiply(1083.6))

30.           .add(106.17);  

31.    //print("太湖水域浑浊度", date, tbd);  

32.    Map.addLayer(tbd, vistbd, date, true);  

33.    var tbd_mean = tbd.reduceRegion(pie.Reducer.mean(), taihu, 1);  

34.    //print("太湖水域浑浊度平均值", date, tbd_mean);  

35.    var tbd_max = tbd.reduceRegion(pie.Reducer.max(), taihu, 1);  

36.    //print("太湖水域浑浊度最大值", date, tbd_max);  

37.    var tbd_min = tbd.reduceRegion(pie.Reducer.min(), taihu, 1);  

38.    //print("太湖水域浑浊度最小值", date, tbd_min);  

39.    dates.push(date);  

40.    tbd_images.push(tbd_mean);  

41.};  

42.//设置图表属性  

43.var line_t = {title: '太湖水域浑浊度动态变化',  

44.              legend: ['浑浊度'],  

45.              xAxisName: "日期",  

46.              yAxisName: "浑浊度",  

47.              chartType: "line",  

48.              yMin: 0,  

49.              yMax: 400,  

50.              smooth: true};  

51.//显示折线图  

52.ChartImage(tbd_images, dates, line_t);  

53.//动画显示  

54.Map.playLayersAnimation(dates, 0.5, 100);  

55.//图例  

56.var data = {title: "浑浊度",  

57.            colors: ['#040274','#040281','#0502a3','#0502b8','#0502ce','#0502e6',

58.                     '#0602ff','#235cb1','#307ef3','#269db1','#30c8e2','#32d3ef',

59.                     '#3be285','#3ff38f','#86e26f','#3ae237','#b5e22e','#d6e21f',

60.                     '#fff705','#ffd611','#ffb613','#ff8b13','#ff6e08','#ff500d',

61.                     '#ff0000', '#de0101', '#c21301', '#a71001', '#911003'],  

62.            labels: ["清澈", "轻度浑浊", "中度浑浊", "重度浑浊"],  

63.            step: 30};  

64.var style = {top: "85%",  

65.             left: "40%",  

66.             width: "350px",  

67.             height: "70px"};  

68.var legend = ui.Legend(data, style);  

69.Map.addUI(legend);