本發(fā)明涉及一種海冰體積的遙感估算方法，屬于遙感地學應用技術領域。

技術背景

兩極海冰作為地球系統(tǒng)的冷源，對全球大氣和海洋環(huán)流有重要和長期的影響，是全球氣候變化的重要指示器。近年來，全球氣候變化的影響越來越明顯，人們越來越需要準確了解兩極海冰正在發(fā)生的變化。海冰厚度和體積的變化日益受到科學界的關注。

隨著空間技術的發(fā)展，利用衛(wèi)星和航空航天對地觀測技術獲取海冰遙感資料，再通過模型算法得出海冰厚度與形態(tài)特征已成為重要的途徑。2010年4月，歐空局發(fā)射的專門為冰凍圈設計的CryoSat-2衛(wèi)星，克服了傳統(tǒng)測高衛(wèi)星ERS-1、ERS-2、ENVISAT軌道傾角的限制，觀測數(shù)據(jù)范圍最高達到北極88°N的區(qū)域。另外，其獨特的多視處理和高分辨率彌補了傳統(tǒng)遙感觀測手段和測高衛(wèi)星的不足。該衛(wèi)星上攜帶的合成孔徑干涉雷達高度計(SIRAL)的垂直測量精度達到1～3cm，相比于先前歐空局雷達高度計(RA)約10km的衛(wèi)星足跡，它采用延遲多普勒雷達高度計(DDA)技術將衛(wèi)星地面足跡減小到沿軌約為0.3km，跨軌約為1.5km。另外，SIRAL對地表點進行多視處理以減少雷達斑點引起的噪聲。其測量海平面高度數(shù)據(jù)精度約為傳統(tǒng)雷達高度計的2倍。

歐洲氣象衛(wèi)星應用組織(EUMETSAT)海洋海冰應用中心(OSI SAF)的核心任務是監(jiān)測南北極冰雪覆蓋的變化。自2005年以來，提供全球海冰類型、密集度、外緣線和反射率等數(shù)據(jù)產品。其海冰類型產品是使用貝葉斯算法和大氣校正模型，基于被動微波散射計(SSMIS)的亮度溫度和主動微波散射計(ASCAT)的后向散射數(shù)據(jù)，反演得到的空間分辨率為10km的極地立體投影網格數(shù)據(jù)。兩者的有效結合，對反演海冰厚度和估算海冰體積有重要意義。

技術實現(xiàn)要素：

本發(fā)明要解決的技術問題是：克服現(xiàn)有技術的上述缺陷，提供一種基于CryoSat-2衛(wèi)星測高數(shù)據(jù)和OSI SAF海冰類型數(shù)據(jù)反演海冰厚度，進而估算海冰體積的方法。本發(fā)明海冰體積估算方法考慮了積雪深度、多年冰和一年冰海冰密度的差別，大大提高了海冰厚度和體積的估算精度。

為了解決上述技術問題，本發(fā)明提出的技術方案是：一種海冰體積的遙感估算方法，包括以下步驟：

第一步、準備數(shù)據(jù)，并對數(shù)據(jù)進行提取和篩選，包括如下兩個方面的內容：

Ⅰ.下載CryoSat-2衛(wèi)星SAR模式L2I級數(shù)據(jù)產品，并提取經緯度坐標、浮冰出水高度、雪深和海冰密集度信息，并根據(jù)指示性參數(shù)剔除不可用的數(shù)據(jù)；

II.下載對應時間和空間范圍的歐洲氣象衛(wèi)星應用組織海洋海冰應用中心的海冰類型數(shù)據(jù)，獲取海冰類型信息以及對應的經緯度信息，所述海冰類型分為一年冰和多年冰；

第二步、分別將上述提取的浮冰出水高度、雪深、海冰密集度以及海冰類型信息與對應的經緯度坐標進行匹配，并將匹配好的數(shù)據(jù)轉換成矢量點數(shù)據(jù)；

第三步、獲取浮冰出水高度柵格數(shù)據(jù)、雪深柵格數(shù)據(jù)、海冰密集度柵格數(shù)據(jù)和海冰類型柵格數(shù)據(jù)，并將上述柵格數(shù)據(jù)轉換為25km×25km的空間分辨率；

第四步、以浮冰出水高度柵格數(shù)據(jù)為掩膜對海冰類型柵格數(shù)據(jù)進行處理，獲得CryoSat-2數(shù)據(jù)覆蓋范圍內浮冰的海冰類型；

第五步、通過流體靜力平衡公式計算每個像元的海冰厚度值H_i，并生成海冰厚度的柵格數(shù)據(jù)，公式表達如下，

$H_i=\frac{{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{F}_i+\frac{{\mathrm{\ρ}}_s}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{H}_s$

其中，ρ_w、ρ_i和ρ_s分別是海水密度、海冰密度和海冰上覆蓋積雪的密度，F(xiàn)_i為浮冰出水高度，H_s為海冰上的積雪深度，一年冰和多年冰具有各自的海冰密度值；

第六步、每個像元所代表的實際面積乘以該像元對應的海冰密集度及海冰厚度值，得到該像元對應的海冰體積，所有像元的海冰體積累加，得到整個海域的海冰體積，公式表達如下：

$V=\underset{j=1}{\overset{n}{\Σ}}(S\×C_j\×H_j)$

其中V表示整個海域的海冰體積，C_j、H_j分別代表第j個像元的海冰密集度和海冰厚度，S為像元的實際面積，即像元的空間分辨率，n為像元總個數(shù)。

本方法數(shù)據(jù)準備階段利用MATLAB軟件對數(shù)據(jù)進行了提取、篩選、匹配和矢量轉換，以便利用ArcGIS進行后續(xù)的處理，使得處理過程更加高效準確。

本發(fā)明海冰體積的遙感估算方法，還具有如下特征：

1、所述第一步中，CryoSat-2衛(wèi)星SAR模式L2I級數(shù)據(jù)產品為DBL格式，讀取DBL格式的原始文件獲得數(shù)據(jù)集，從數(shù)據(jù)集中提取經緯度、浮冰出水高度、雪深和海冰密集度信息；歐洲氣象衛(wèi)星應用組織海洋海冰應用中心海冰類型數(shù)據(jù)格式為NetCDF4格式。提取數(shù)據(jù)后，剔除異常值，減小了數(shù)據(jù)量，提高了數(shù)據(jù)處理效率。

2、所述第二步中，使用MATLAB軟件中shapewrite函數(shù)將匹配好的數(shù)據(jù)轉換成帶經緯度坐標的矢量點數(shù)據(jù)，存儲為shapefile文件。將提取的數(shù)據(jù)轉換成帶經緯度坐標的矢量點，有利于ArcGIS的后續(xù)處理。

3、所述第三步具體包括如下步驟：

a.對第二步中獲得的矢量點數(shù)據(jù)在WGS_1984地理坐標系下進行極地立體方位投影；

b.將極地立體方位投影后的矢量點數(shù)據(jù)進行柵格轉換，獲得浮冰出水高度柵格數(shù)據(jù)、積雪深度柵格數(shù)據(jù)、海冰密集度柵格數(shù)據(jù)和海冰類型柵格數(shù)據(jù)；

c.通過空間重采樣將上述柵格數(shù)據(jù)的空間分辨率采樣至相同分辨率。

獲取同一坐標系下統(tǒng)一空間分辨率的浮冰出水高度、雪深、海冰類型和海冰密集度柵格數(shù)據(jù)，有利于海冰類型和浮冰信息相匹配，從而確定浮冰類型。

4、所述第四步中，海水密度ρ_w取1023.8kg/m³，一年海冰和多年海冰的密度分別取916.7kg/m³和882.0kg/m³，積雪密度取400.0kg/m³。

海冰密度誤差在海冰厚度反演中占據(jù)主要地位，一年冰密度大于多年冰。對一年冰和多年冰賦予不同的密度參數(shù)，提高了厚度反演的精度。海水和海冰的密度參考了其他文獻，例如Wadhams等人2010年發(fā)表于The Cryosphere的“The relation between sea ice thickness and freeboard in the Arctic”通過689個觀測站點的數(shù)據(jù)得到一年冰的平均密度是916.7±35.7kg/m³，多年冰的平均密度是882±23kg/m³，當測量值較多時，將一年冰和多年冰的密度設定為916.7kg/m³和882.0kg/m³是可行的。Rachel等人2015年發(fā)表于Nature的“Increased Arctic sea ice volume after anomalously low melting in 2013”將海水密度設定為1023.8kg/m³，一年冰和多年冰的密度分別為916.7kg/m³和882.0kg/m³。

5、所述第五步中，通過ArcGIS軟件中的Raster Calculator工具，根據(jù)流體靜力學平衡模型反演得到海冰厚度。

本發(fā)明利用了CryoSat-2衛(wèi)星數(shù)據(jù)和OSI SAF海冰類型數(shù)據(jù)，數(shù)據(jù)獲取簡單。CryoSat-2SIRAL測高數(shù)據(jù)能提供精確的浮冰出水高度信息，另外其輔助測量數(shù)據(jù)中還包含了積雪深度、海冰密集度的信息。OSI SAF海冰類型數(shù)據(jù)是通過貝葉斯算法和大氣校正模型反演得到的空間分辨率為10km的極地立體投影網格數(shù)據(jù)。兩者的有效結合提高了海冰厚度的反演精度，為海冰體積的估算提供了數(shù)據(jù)基礎。

本發(fā)明數(shù)據(jù)提取和矢量轉換過程均通過MATLAB編程自動實現(xiàn)，大大減少了工作量，同時減少了人工誤差。后續(xù)柵格計算通過ArcGIS中的Raster Calculator工具完成。

綜上，本發(fā)明方法的執(zhí)行步驟簡單易行，效果較好。目前海冰厚度大面積觀測資料較少，傳統(tǒng)的經驗方法估算海冰厚度精度不高。本發(fā)明使用高精度大范圍的SIRAL浮冰出水高度數(shù)據(jù)和OSI SAF海冰類型產品相結合，根據(jù)流體靜力學平衡原理反演海冰厚度，并估算海冰體積。對精確量化南北極海冰物質平衡變化有深遠的科學意義。

附圖說明

下面結合附圖對本發(fā)明作進一步的說明。

圖1是本發(fā)明海冰體積遙感估算方法流程圖。

圖2(a)為提取的浮冰出水高度、雪深、海冰密集度信息。

圖2(b)為海冰類型信息。

圖3(a)、(b)、(c)、(d)分別為25km空間分辨率的浮冰出水高度柵格數(shù)據(jù)、雪深柵格數(shù)據(jù)、海冰密集度柵格數(shù)據(jù)、海冰類型柵格數(shù)據(jù)的圖像。

圖4為海冰厚度結果顯示圖。

圖5為一年冰和多年冰海冰體積的估算結果。

具體實施方式

下面根據(jù)附圖詳細闡述本發(fā)明，使本發(fā)明的技術路線和操作步驟更加清晰。本發(fā)明實例采用的數(shù)據(jù)是CryoSat-2衛(wèi)星L2I級SAR模式基線C數(shù)據(jù)和OSI SAF10km分辨率網格數(shù)據(jù)。CryoSat-2數(shù)據(jù)獲取時間為2016年1月1日——2016年1月31日，OSI SAF數(shù)據(jù)獲取時間為2016年1月15日。

圖1為海冰體積遙感估算方法流程圖，具體步驟如下:

第一步、準備數(shù)據(jù)階段，具體操作步驟如下：

Ⅰ.依次讀取CryoSat-2數(shù)據(jù)3597個DBL格式原始文件，獲得數(shù)據(jù)集。從數(shù)據(jù)集中提取每個經緯度坐標點的浮冰出水高度、雪深和海冰密集度信息，并提出不可用數(shù)據(jù)；

II.下載對應時間的歐洲氣象衛(wèi)星應用組織(EUMETSAT)海洋海冰應用中心(OSI SAF)海冰類型數(shù)據(jù)。讀取OSI SAF數(shù)據(jù)NetCDF4格式文件，獲取海冰類型以及對應的經緯度信息。

第二步、分別將上述提取的浮冰、雪深、密集度以及海冰類型信息與對應的經緯度匹配。使用shapewrite函數(shù)將上述匹配好的數(shù)據(jù)轉換成帶經緯度坐標的矢量點數(shù)據(jù)。分別為矢量文件a和b。其中矢量文件a包含了浮冰出水高度、積雪深度、海冰密集度信息。矢量文件b包含了海冰類型信息。矢量數(shù)據(jù)的屬性值如圖2(a)、圖2(b)所示。其中圖2(b)中屬性值Ice_type為2代表一年冰，為3代表多年冰。

第三步、使用ArcGIS軟件獲取浮冰出水高度、雪深、海冰類型以及海冰密集度柵格數(shù)據(jù)，并將上述數(shù)據(jù)統(tǒng)一為25km空間分辨率的柵格數(shù)據(jù)集，結果如圖3所示。用海冰類型鑒別浮冰。

a.定義投影。將矢量文件a和b加載到ArcGIS中，定義WGS_1984地理坐標系，進行極地立體方位投影。

b.獲得柵格數(shù)據(jù)集。利用ArcGIS工具Conversion Tools->To Raster->Point to Raster，計算方式選擇mean，柵格大小為默認值。獲得浮冰出水高度、積雪深度、海冰密集度和海冰類型柵格數(shù)據(jù)。

c.統(tǒng)一空間分辨率。利用ArcGIS工具Data Management Tools->Raster->Raster Processing->Resample，選取最鄰近分配法，將上述柵格數(shù)據(jù)統(tǒng)一為25km×25km空間分辨率。其中浮冰出水高度、雪深和海冰密集度柵格數(shù)據(jù)分別如圖3(a)、3(b)、3(c)。

第四步、掩膜提取。以浮冰出水高度數(shù)據(jù)掩膜提取海冰類型，從而獲得CryoSat-2數(shù)據(jù)覆蓋范圍浮冰的海冰類型，如圖3(d)所示。本例中，利用ArcGIS工具Spatial Analyst Tools->Extraction->Extract by Mask。

第五步、通過流體靜力平衡公式計算每個像元的海冰厚度值H_i，并生成柵格圖層，如圖4所示。公式表達如下，

$H_i=\frac{{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{F}_i+\frac{{\mathrm{\ρ}}_s}{{\mathrm{\ρ}}_w-{\mathrm{\ρ}}_i}{H}_s$

式中，ρ_w、ρ_i和ρ_s分別是海水密度、海冰密度和海冰上覆蓋積雪的密度，F(xiàn)_i為浮冰出水高度，H_s為海冰上的積雪深度。本例中，海水密度ρ_w取1023.8kg/m³，一年海冰和多年海冰的密度分別取916.7kg/m³和882.0kg/m³，積雪密度取400.0kg/m³。

第六步、估算海冰體積。每個像元所代表的實際面積乘以該像元對應的海冰密集度及海冰厚度值，得到該像元對應的海冰體積，所有像元對應的海冰體積相加，得到整個海域的海冰體積，估算結果如圖5所示，其中一年冰體積約為16.33×10³km³，多年冰體積約為7.22×10³km³。具體公式表達如下：

$V=\underset{j=1}{\overset{n}{\Σ}}(S\×C_j\×H_j)$

其中V表示整個海域的海冰體積，C_j、H_j分別代表第j個像元的海冰密集度和海冰厚度，S為像元的實際面積，即像元的空間分辨率，n為像元總個數(shù)。

除上述實施例外，本發(fā)明還可以有其他實施方式。凡采用等同替換或等效變換形成的技術方案，均落在本發(fā)明要求的保護范圍。