一種積雪中黑碳濃度的遙感估算方法
【專利摘要】本發明涉及一種積雪中黑碳濃度的遙感估算方法，利用500米空間分辨率的MODIS每日積雪合成產品（MOD10A1和MYD10A1），結合地面站點實測黑碳濃度數據，進行區域尺度黑碳濃度估算。首先對MODIS每日積雪產品進行影像拼接、地圖投影轉換和數據格式轉換，獲取研究區域積雪反照率數據；然后對影像進行去云合成處理，提取地面站點位置的積雪反照率，并參照相應文獻，建立積雪反照率和黑碳濃度的關系模型，接著對關系模型進行檢驗和反復優化；最后利用關系模型處理積雪反照率影像，利用研究區域矢量邊界數據截取，得到積雪中黑碳濃度的時空分布結果。與傳統方法相比，本發明引入遙感技術手段，可高效獲取大范圍積雪中黑碳濃度分布的時空特征數據。
【專利說明】
一種積雪中黑碳濃度的遙感估算方法
技術領域
[0001]本發明涉及一種積雪中黑碳濃度的遙感估算方法，屬于大氣環境遙感應用技術領域。【背景技術】
[0002]黑碳(Black Carbon)氣溶膠是一種懸浮于空氣中的黑色碳質粒顆粒物，是大氣細顆粒物的重要組成部分，對大氣環境以及人體健康都有著非常重要的影響。積雪中的黑碳顆粒能顯著降低積雪反照率，增強對太陽輻射的吸收，加速積雪融化，造成顯著的氣候效應和水文循環影響。因此進行冰雪中黑碳含量的研究，獲取其時空分布特征，有利于增強對氣候變化及區域水循環的認識。
[0003]目前國內外對積雪中黑碳含量進行了較多的觀測研究，其研究手段主要是在地面建立觀測站點，采集雪樣并進行相應測量，從而獲取該位置積雪中的黑碳含量。這種方法雖然精確度高，但是費工費時，難以做到大范圍積雪中黑碳含量的快速獲取和更新。而大范圍、實時、動態正是遙感技術的特點。考慮到這一特點，發明人考慮基于遙感技術手段，利用 M0DIS積雪數據產品，建立積雪反照率和黑碳濃度的關系模型，從而實現區域尺度黑碳濃度的快速估算。
【發明內容】
[0004]本發明要解決的技術問題是:克服現有技術手段的缺點，建立一種基于遙感光學影像數據和已有地面測站數據的關系模型，利用M0DIS積雪反照率數據，能夠較準確、快速地獲取積雪中黑碳含量及空間分布特征，達到比較好的效果。
[0005]為了解決上述技術問題，本發明提出的技術方案是:一種積雪中黑碳濃度的遙感估算方法，包括以下步驟:第一步、利用地面站點實測積雪中黑碳濃度數據，根據地面站點地理位置及測量時間， 下載對應時間和區域的M0DIS每日積雪合成產品M0D10A1和MYD10A1;第二步、對MOD IS影像數據產品MOD 10A1和MYD10A1分別進行預處理，包括影像拼接、地圖投影轉換和格式轉換，并提取產品中的積雪反照率數據；第三步、對產品MOD 10A1和產品MYD 10A1獲取的積雪反照率數據進行去云合成處理，獲得積雪反照率影像；第四步、依據地面站點地理坐標提取積雪反照率影像上的對應位置積雪反照率； 第五步、參照文獻資料，分析地面站點積雪中的黑碳濃度情況，選取合適站點，初步建立研究區域積雪反照率與黑碳濃度的關系模型；第六步、進行積雪反照率影像向黑碳濃度空間數據的轉換，并對已建立的關系模型及初步得到的黑碳濃度結果，結合樣本數據和相關文獻資料進行檢驗評定，反復優化，直到確立最終的關系模型；第七步、通過研究區域矢量邊界數據，根據最終的關系模型提取出研究區域積雪中黑碳濃度的時空分布結果。
[0006]本方法利用500米空間分辨率的M0DIS每日積雪合成產品M0D10A1和MYD10A1，結合地面站點實測黑碳濃度數據，進行區域尺度黑碳濃度估算。首先對M0DIS每日積雪產品進行影像拼接、地圖投影轉換和數據格式轉換，獲取研究區域積雪反照率數據;然后對影像進行去云合成處理，提取地面站點位置的積雪反照率，并參照相應文獻，建立積雪反照率和黑碳濃度的關系模型，接著對關系模型進行檢驗和反復優化;最后利用關系模型處理積雪反照率影像，利用研究區域矢量邊界數據截取，得到積雪中黑碳濃度的時空分布結果。與傳統方法相比，本發明引入遙感技術手段，可高效獲取大范圍積雪中黑碳含量分布的時空特征數據。
[0007]本發明積雪中黑碳濃度的遙感估算方法，還具有如下改進:1、所述第二步中，影像拼接的方法是，根據研究區范圍，將同一時間相應的幾景影像拼接為一景完整的影像。
[0008]2、所述第三步中，所述去云合成處理的算法如下:a、對于某一坐標的像元，如果產品MOD 10A1對應的反照率為積雪反照率，則把產品 M0D10A1的該坐標像元作為合成后積雪反照率影像的像元；b、對于某一坐標的像元，如果產品M0D10A1對應的反照率為非積雪反照率，而MYD10A1 對應的反照率為積雪反照率，則把產品MYD10A1的該坐標像元作為合成后積雪反照率影像的像元；c、對于某一坐標的像元，如果產品M0D10A1和產品MYD10A1對應的反照率均非積雪反照率，則合成后積雪反照率影像在該坐標的像元為非積雪反照率像元。
[0009]3、所述第四步中，根據地理位置提取積雪反照率，有兩種方法可供選擇:地面站點較少時，在ENVI軟件主窗口中，利用pixel locator提取相應坐標位置的反照率;地面站點較多時，利用ENVI+IDL語言，通過索引到影像上的行列號，批量獲取地面站點對應影像上的反照率。
[0010]4、所述第五步中，為了避免所選樣本中的異常值對最終的擬合結果造成較大的影響，應參考其他已有的研究成果對地面站點進行篩選。
[0011]5、所述第五步中，采用最小二乘法對積雪反照率和黑碳濃度數據進行線性擬合， 從而建立積雪反照率和黑碳濃度的關系模型。
[0012]本發明利用地面站點積雪中黑碳濃度數據，以及M0DIS積雪反照率數據產品，數據獲取比較容易。影像數據的預處理、站點積雪反照率的提取、利用關系模型處理影像等步驟可利用軟件MRT、ENVI及MATLAB處理，簡單易懂、操作方便，大大地減少了工作量，提高了處理的效率，增強了本發明的靈活性和可用性。針對地面站點采樣情況以及該地區黑碳含量的基本狀況，對數據進行了篩選，使建立的關系模型更加符合該地區的自然情況。同時結合 M0DIS積雪反照率產品合成算法，本發明可高效、快速、準確地獲取較大尺度地區積雪黑碳濃度。
[0013]綜上，本發明方法的執行步驟簡單易行，效果較好，解決了黑碳含量研究中樣本采集、測量過程中費時費力的問題，有利于促進對黑碳含量分布變化、黑碳對積雪消融的影響、乃至黑碳的氣候效應等方面的研究。【附圖說明】
[0014]下面結合附圖對本發明作進一步說明。
[0015]圖1是本發明一種積雪中黑碳濃度的遙感估算方法流程圖。
[0016]圖2是使用MRT軟件對M0DIS積雪產品進行預處理示意圖。
[0017]圖3是使用ENVI軟件提取地面站點積雪反照率示意圖。
[0018]圖4是使用MATLAB軟件根據關系模型處理影像示意圖。
[0019]圖5是建立的積雪反照率和黑碳濃度的關系模型。
[0020]圖6是研究區域北疆地區黑碳濃度結果圖。【具體實施方式】[0021 ]下面根據附圖詳細闡述本發明，使本發明的目的和效果變得更加明顯。
[0022]如圖1所示，為本發明一種積雪中黑碳濃度的遙感估算方法流程圖，具體步驟如下:第一步、利用已有文獻(Ye等人2012年在ENVIRONMENTAL RESEARCH LETTERS中的論文 “Black carbon in seasonal snow across northern Xinjiang in northwestern China”)給出的地面站點實測積雪中黑碳濃度數據，根據站點地理位置及測量時間，下載對應時間和區域的MOD IS每日積雪合成產品(MOD 10A1和MYD10A1)。
[0023]本步驟中，已有文獻給出的地面站點共有32個，研究區為新疆北疆地區，測量時間為2012年1月28日至2月27日。研究選取美國國家冰雪數據中心(Nat1nal Snow and Ice Data Center, NSIDC)提供的Terra/MODIS每日積雪合成產品 M0D10A1 和Aqua/MODIS每日積雪合成產品MYD10A1，空間分辨率均為500m JODI S每日積雪合成產品共有四個數據子集: 積雪、積雪反照率、積雪覆蓋比例和質量驗證數據。研究區域覆蓋的軌道號為h23v04/ h24v04,所用影像共80景。
[0024]第二步、對兩種M0DIS影像數據產品分別進行預處理，提取出覆蓋研究區的積雪反照率數據。利用MRT軟件進行處理，如圖2所示，具體包括如下幾個方面的內容:a、根據研究區的范圍，將同一時間相應的兩軌影像拼接為一景完整的影像。
[0025]b、對拼接后的影像數據進行地圖投影轉換，將原始數據產品的正弦曲線投影 (SIN)重新投影為 WGS 1984/Geographic 系統。
[0026]c、對投影轉換后的數據進行格式轉換，將原始數據產品的HDF格式轉化成容易處理的Geotiff格式。[〇〇27]d、M0DIS積雪產品中有四個數據子集，只提取其中的積雪反照率數據即可。[〇〇28]第三步、研究區域地面站點較少，只有32個。讀取M0DIS積雪反照率影像，可以發現不少區域都容易受到云的影響。M0D10A1和MYD10A1的像元編碼方式中，積雪反照率為0? 100，云為150，利用MATLAB軟件對兩種M0DIS反照率產品，運用像元重新處理的方法進行影像的去云合成。具體的算法如下:a、對于某一坐標的像元，如果產品MOD 10A1對應的反照率為積雪反照率，則把產品 M0D10A1的該坐標像元作為合成后積雪反照率影像的像元；b、對于某一坐標的像元，如果產品M0D10A1對應的反照率為非積雪反照率，而MYD10A1對應的反照率為積雪反照率，則把產品MYD10A1的該坐標像元作為合成后積雪反照率影像的像元；c、對于某一坐標的像元，如果產品M0D10A1和產品MYD10A1對應的反照率均非積雪反照率，則合成后積雪反照率影像在該坐標的像元為非積雪反照率像元。[〇〇29]第四步、根據每個地面站點的地理坐標，利用ENVI軟件提取其對應積雪反照率影像上面對應的像元值，如圖3所示。該步驟用到了兩種方法:在ENVI軟件主窗口中，利用 Pixe 1 Locator手動提取影像上某一位置的反照率值;利用ENVI+IDL語言，通過索引到影像上的行列號，批量獲取地面站點對應影像上的反照率值。由于所用到的站點數較少，兩種方法的效率相差不大，但當站點數目很多時，明顯第二種方法效率更高。
[0030]第五步、為了避免所選樣本中的異常值對最終的擬合結果造成較大的影響，地面站點的篩選參照了其它文獻資料關于研究區域的結論，例如Ming等人2009年在 Atmospheric Research中的論文“Black Carbon (BC) in the snow of glaciers in west China and its potential effects on albedos”得到天山地區冰雪中黑碳含量 112 ±27ng g-1;Xu等人2012年在ENVIRONMENTAL RESEARCH LETTERS的論文 “Post-deposit1nal enrichment of black soot in snow-pack and accelerated melting of Tibetan glaciers”得出新疆天山一號冰川冰雪中冬季表層黑碳含量27?31 ng g^1，夏季可達250?500 ng g<。選取積雪表層黑碳濃度<300 ng的地面站點作為合適的樣本點，進行關系模型的構建。[0〇31]對樣本站點數據序列Y(x)，x=l，2,…，100,以線性函數Y(x)=ax+ b來擬合，按最小二乘法可求得a和b，其中，x為積雪反照率，Y為黑碳濃度，a為趨勢項，b為常數項。a值的正負表示黑碳濃度隨積雪反照率變化的方向，a絕對值的大小表示變化的快慢程度。統計樣本數為26,通過顯著性a=0.1的顯著性檢驗。
[0032]第六步、根據初步建立的關系模型，在MATLAB軟件里面編寫程序，完成反照率影像向黑碳濃度空間數據的轉換，如圖4所示。對已建立的關系模型及初步得到的黑碳濃度結果，結合樣本數據、相關文獻資料進行檢驗評定，反復優化，直到確立最終的關系模型。
[0033]對關系模型檢驗評定考慮但不僅限于這些方面:所有樣本點擬合值的均方差要盡可能的小;擬合的線性函數盡可能逼近所有樣本點，相關系數R越接近1越好;整個反照率區間對應的擬合值都要<300 ng g^1;影像處理后黑碳濃度滿足該地區的一些基本情況，如純凈積雪區(積雪深度大、人類活動影響小的地區)黑碳濃度低，城鎮密集的地方(如天山北坡經濟區)黑碳濃度相應要高等。
[0034]由于缺少積雪反照率在20%以下的樣本點數據，所以最終確定的關系模型為Y(x) =-2.4813x + 255.85，R2=0.5243，xG [20,100)。最終的關系模型如圖5所示。
[0035] 第七步、在ArcGIS軟件中，用空間分析(Spatial Analyst)模塊下的Extract1n 工具，利用新疆北疆地區的矢量邊界截取獲得目標區域的黑碳濃度分布。對結果進行重分類，得到最終的結果圖，如圖6所示。
[0036]除上述實施例外，本發明還可以有其他實施方式。凡采用等同替換或等效變換形成的技術方案，均落在本發明要求的保護范圍。
【主權項】
1.一種積雪中黑碳濃度的遙感估算方法，包括以下步驟:第一步、利用地面站點實測積雪中黑碳濃度數據，根據地面站點地理位置及測量時間， 下載對應時間和區域的MODIS每日積雪合成產品M0D10A1和MYD10A1;第二步、對MOD IS影像數據產品MOD 10A1和MYD10A1分別進行預處理，包括影像拼接、地 圖投影轉換和格式轉換，并提取產品中的積雪反照率數據；第三步、對產品M0D10A1和產品MYD10A1獲取的積雪反照率數據進行去云合成處理，獲 得積雪反照率影像；第四步、依據地面站點地理坐標提取積雪反照率影像上的對應位置積雪反照率；第五步、參照文獻資料，分析地面站點積雪中的黑碳濃度情況，選取合適站點，初步建 立研究區域積雪反照率與黑碳濃度的關系模型；第六步、進行積雪反照率影像向黑碳濃度空間數據的轉換，并對已建立的關系模型及 初步得到的黑碳濃度結果，結合樣本數據和相關文獻資料進行檢驗評定，反復優化，直到確 立最終的關系模型；第七步、通過研究區域矢量邊界數據，根據最終的關系模型提取出研究區域積雪中黑 碳濃度的時空分布結果。2.根據權利要求1所述，其特征在于:第二步中，影像拼接的方法是，根據研究區范圍， 將同一時間相應的幾景影像拼接為一景完整的影像。3.根據權利要求1所述，其特征在于:第三步中，所述去云合成處理的算法如下:a、對于某一坐標的像元，如果產品M 0 D10 A1對應的反照率為積雪反照率，則把產品 M0D10A1的該坐標像元作為合成后積雪反照率影像的像元；b、對于某一坐標的像元，如果產品M0D10A1對應的反照率為非積雪反照率，而MYD10A1 對應的反照率為積雪反照率，則把產品MYD10A1的該坐標像元作為合成后積雪反照率影像 的像元；c、對于某一坐標的像元，如果產品M0D10A1和產品MYD10A1對應的反照率均非積雪反照 率，則合成后積雪反照率影像在該坐標的像元為非積雪反照率像元。4.根據權利要求1所述，其特征在于:第四步中，根據地理位置提取積雪反照率，有兩種 方法可供選擇:地面站點較少時，在ENVI軟件主窗口中，利用pixel locator提取相應坐標 位置的反照率;地面站點較多時，利用ENVI+IDL語言，通過索引到影像上的行列號，批量獲 取地面站點對應影像上的反照率。5.根據權利要求1所述，其特征在于:第五步中，為了避免所選樣本中的異常值對最終 的擬合結果造成較大的影響，應參考其他已有的研究成果對地面站點進行篩選。6.根據權利要求1所述，其特征在于:第五步中，采用最小二乘法對積雪反照率和黑碳 濃度數據進行線性擬合，從而建立積雪反照率和黑碳濃度的關系模型。
【文檔編號】G06F17/50GK106021655SQ201610302401
【公開日】2016年10月12日
【申請日】2016年5月9日
【發明人】柯長青, 馬東輝, 蔡宇, 邵珠德, 趙佳曼
【申請人】南京大學