1?引言??

太湖是我國(guó)第三大淡水湖，地跨江浙兩省，水域面積達(dá)2338km2，平均水深2m左右。上世紀(jì)80年代以來(lái)，太湖水體營(yíng)養(yǎng)化愈加嚴(yán)重，頻繁爆發(fā)藍(lán)藻水華現(xiàn)象即為一種非常典型的二類水體（Case 2 Waters，光學(xué)性質(zhì)變化受浮游植物及其附屬物、外生的粒子和外生有色可溶有機(jī)物等其他物質(zhì)影響的水體）。懸浮物（Suspended Solids，懸浮在水中的固體物質(zhì)，包括不溶于水的無(wú)機(jī)物、有機(jī)物及泥砂、黏土、微生物等）作為造成水體渾濁的主要影響因素，是衡量水體污染程度的重要指標(biāo)。利用遙感技術(shù)準(zhǔn)確獲取區(qū)域面狀水域懸浮物濃度信息，是遙感水質(zhì)參數(shù)監(jiān)測(cè)的一項(xiàng)重要任務(wù)。傳統(tǒng)的利用遙感手段監(jiān)測(cè)水質(zhì)的方法大多基于時(shí)序遙感影像及地表水質(zhì)監(jiān)測(cè)數(shù)據(jù)構(gòu)建地物信息模型來(lái)完成，但受限于本地設(shè)備存儲(chǔ)、運(yùn)算能力的限制，當(dāng)前的監(jiān)測(cè)尺度已逐漸無(wú)法滿足日益提升的大尺度、長(zhǎng)時(shí)序遙感應(yīng)用研究需求，且需要消耗大量人力物力和時(shí)間。

當(dāng)前我們身處大數(shù)據(jù)時(shí)代，云端計(jì)算技術(shù)蓬勃發(fā)展，航天宏圖依托行業(yè)多年技術(shù)積累獨(dú)立自主研發(fā)了安全可控的開放式云計(jì)算平臺(tái)：PIE-Engine（Pixel Information Expert Engine，像素專家引擎），致力于加速我國(guó)遙感云計(jì)算平臺(tái)發(fā)展進(jìn)程，實(shí)現(xiàn)了遙感數(shù)據(jù)按需獲取、運(yùn)算以及專題信息聚焦服務(wù)，以滿足對(duì)地觀測(cè)數(shù)據(jù)獲取能力飛速增長(zhǎng)帶來(lái)的信息高效化處理和服務(wù)需求。

2 .反演太湖水域2013-2020年藍(lán)藻密度與渾濁度

基于PIE-Engine獲取2013-2020年覆蓋太湖水域的Landsat 8 Surface Reflectance遙感影像數(shù)據(jù)，Landsat 8 SR數(shù)據(jù)集是OLI/TIRS傳感器數(shù)據(jù)經(jīng)過大氣校正的表面反射率數(shù)據(jù)，分辨率為30米，重放周期為16天。影像包含5個(gè)可見光波段、1個(gè)近紅外（NIR）波段和2個(gè)短波紅外（SWIR）波段，經(jīng)過正射校正后的表面反射率，以及兩個(gè)經(jīng)過正射校正后的亮度溫度熱紅外（TIR）波段。

相關(guān)研究指出，Landsat 8數(shù)據(jù)的近紅外波段和紅光波段分別與藍(lán)藻密度和渾濁度具有較高的相關(guān)系數(shù)，但相關(guān)程度都不高，主要原因有：① 水華與水草區(qū)域?qū)λ|(zhì)反演的影響；② 水表面光滑度和微波隨時(shí)間和空間的干擾；③ 陸地光譜對(duì)近岸邊水體的干擾；④ 水體中較高濃度懸浮物的高反射干擾。在去除影響區(qū)域和異常數(shù)據(jù)，并利用普通水體與不同波段的組合構(gòu)建相關(guān)性模型后發(fā)現(xiàn)，藍(lán)藻密度與B5/B4以及渾濁度與(B4-B3)/(B4+B3)的波段組合具有最高的相關(guān)系數(shù)。相關(guān)性模型如下：

其中，B3、B4、B5分別為綠光、紅光以及近紅外波段的反射率。

3 太湖水域2013-2020年藍(lán)藻密度?與渾濁度動(dòng)態(tài)變化

3.1藍(lán)藻密度

示例代碼：

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.?//設(shè)置圖層顯示屬性??

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.?????//計(jì)算藍(lán)藻密度??

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

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

30.?????//print("太湖水域藍(lán)藻密度",?date,?alg);??

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

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

33.?????//print("太湖水域藍(lán)藻密度平均值",?date,?alg_mean);??

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

35.?????//print("太湖水域藍(lán)藻密度最大值",?date,?alg_max);??

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

37.?????//print("太湖水域藍(lán)藻密度最小值",?date,?alg_min);??

38.?????dates.push(date);??

39.?????alg_images.push(alg_mean);??

40.?};??

41.?//?設(shè)置圖表屬性??

42.?var?line_a?=?{title:?'太湖水域藍(lán)藻密度動(dòng)態(tài)變化',??

43.???????????????legend:?['藍(lán)藻密度'],??

44.???????????????xAxisName:?"日期",??

45.???????????????yAxisName:?"藍(lán)藻密度（萬(wàn)/L）",??

46.???????????????chartType:?"line",??

47.???????????????yMin:?0,??

48.???????????????yMax:?2500,??

49.???????????????smooth:?true};??

50.?//?顯示折線圖??

51.?ChartImage(alg_images,?dates,?line_a);??

52.?//動(dòng)畫顯示??

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

54.?//?圖例??

55.?var?data?=?{title:?"藍(lán)藻密度",??

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（萬(wàn)個(gè)/L）",?"2000（萬(wàn)個(gè)/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.//設(shè)置水域預(yù)覽顏色組合??

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.? ? // 計(jì)算太湖水域渾濁度??

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.//設(shè)置圖表屬性??

43.var line_t = {title: '太湖水域渾濁度動(dòng)態(tài)變化',??

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.//動(dòng)畫顯示??

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);??