基于Google Earth Engine的采煤沉陷水體方向變化自動識別

趙艷玲，丁寶亮，何廳廳，肖 武，任 河

(1.中國礦業大學(北京) 地球科學與測繪工程學院，北京 100083；2.浙江大學 公共管理學院，浙江 杭州 310058)

能源是一個國家可持續發展的基礎和重要支柱，對于我國而言，煤炭資源的地位尤為突出，占我國化石能源基礎儲量的94%左右。但隨著煤炭的持續開采，勢必會造成嚴重的環境問題，如地表沉陷、土壤污染等。據統計，平均每開采1萬t煤，就有0.2～0.3 ha土地受到開采沉陷的影響。由于煤礦開采引起的地表擾動一直是一個全球性的問題，在美國、澳大利亞等國家由于地表下沉產生的沉陷水體對當地生態造成了難以修復的破壞。在我國東部高潛水位地區，如安徽，地下水極易由于地表沉降而顯露出來，從而形成地面沉陷積水區，大面積的沉陷水體會損毀大量優質耕地，影響當地的生態系統健康。據安徽省淮南市政府統計，至2017年淮南市采煤塌陷區面積已達3萬ha，近一半塌陷坑被水覆沒，且由于地形平坦，土地復墾和回填的成本較高，大部分塌陷坑長期處于積水狀態。因此，為了制定合理有效的復墾方案，對塌陷區積水的時空變化進行持續監測是有必要的。

遙感技術在監測地表自然資源變化方面已經非常成熟，基于遙感圖像的土地利用變化監測已廣泛用于采礦方向，如基于SAR技術監測礦區地表沉降，其通過跟蹤 SAR圖像像素偏移來監測地表的位移。對沉降邊界進行預測的方法主要有時間函數模型、基于實測數據的預測模型、力學和數值模擬等方法。但是這些方法往往需要采煤相關數據，且處理速度較慢，自動化程度低。對于礦區沉陷積水的變化研究一般是基于3S技術，從水體的“大小”、“形狀”和“變化率”角度出發，在一定的時間和空間尺度進行監測，大量研究采用水體面積年均變化幅度、年均變化速度以及水體向其他土地利用類型轉移的列聯矩陣等方法來描述沉陷積水的時空變化。隨著GIS技術的進一步發展，出現了Digital Shoreline Analysis System(DSAS)方法，其通過繪制基線來對研究區邊界進行監測。STEELE和HEFFERNAN運用海岸線發展系數Shoreline Development Factor(SDF)來描述水體的大小分布、連通性和形狀；BHAGAWAT以城市中心為中點，分別向東、東南、南、西南、西、西北、北、東北共8個方向發射射線，通過這8條射線來監測該城市邊界的空間位置變化，并取得了較高的精度。但在采煤沉陷水體的研究上，對沉陷水體各方向的空間位置變化進行監測的研究還較少。有研究發現，隨著煤炭的持續開采，地面形成的沉陷區通常會隨著開采工作面而呈現規律性的變化，而沉陷積水區也會出現類似的變化。通常，煤炭開采分為橫向單煤層多工作面開采和縱向多煤層工作面開采。在橫向單煤層多工作面開采過程中，隨著工作面推進，地表沿該方向依次沉降，沉陷積水多呈現條帶狀；在縱向多煤層開采過程中，隨著開采程度加深，地表沿垂直工作面方向向四周沉降，沉陷積水多呈現圓狀。在煤炭開采中，2種開采方式往往同時存在。以上的方法與結論為本次研究提供了一條新思路。谷歌大數據云計算平臺GEE能很快地調出大量遙感數據，并且其內部封裝了大量函數可供使用，也逐漸被應用到采礦領域。在進行礦區沉陷水體的研究上，YI等構建了具有空間特征的土壤水分時間序列軌跡，并通過LandTrendr算法，對采煤沉陷區進行時空動態分析；YANG等采用改進歸一化水體指數(mNDWI)和最大類間方差(OTSU)圖像分割算法提取了1988—2018年中國東部高地下水位礦區的塌陷濕地。由此可見，GEE平臺在進行礦區沉陷水體的研究上具有處理速度快、處理數據量大等的特點，大有發展前景。但平臺上還尚未出現類似DSAS等的相關函數，需要使用者通過編寫代碼繪制基線(射線)從而實現監測功能。

筆者利用遙感大數據平臺GEE，基于1989—2016年潘謝礦區沉陷水體數據，綜合幾何學和統計學知識，對每個獨立的沉陷水體分別繪制射線，通過多組射線來分析潘謝礦區沉陷水體各方向的年際形狀變化，并對其進行擬合，預測其變化情況。此方法可在缺少采煤相關信息的條件下僅依據遙感技術實現，且能達到較高的自動化程度，為礦區內沉陷水體的監測和生態修復方案的制定提供一定的依據。

1 研究區概況及數據源

1.1 研究區概況

潘謝礦區位于安徽省淮南市(圖1)，地理坐標：116.33°E—116.90°E，32.72°N—32.93°N，南北跨度25 km，東西跨度60 km。礦區地形平坦，標高30～40 m，無承壓地下水埋深為1.5 m左右，當礦區地面沉降達1.5 m時可能出現積水。氣候屬溫帶季風氣候，年均氣溫15 ℃，年均降水量970 mm。礦區內共有謝橋礦、張集礦、顧北礦、顧橋礦、丁集礦、潘三礦、朱集礦、潘北礦、潘二礦和潘一礦等10個礦山企業。潘謝礦區的開采歷史多達110余年，煤炭資源豐富，有9～18層可采煤層，平均厚度為20～30 m，屬于近水平厚松散層煤層群。由于長期開采，地面經歷多次沉陷，長期處于不穩定沉陷狀態，地面累計沉陷深度最大可到20～30 m，加之當地較高的地下水位，形成了大面積的沉陷積水區，對當地房屋、耕地造成了嚴重影響。

圖1 研究區位置示意Fig.1 Schematic diagram of the location of the study area

1.2 數據源及數據處理

Google Earth Engine為一個存儲衛星數據以及針對衛星數據進行批量處理的云計算平臺，它存儲了全球40多年的遙感數據，包括Landsat系列、Sentinel,MODIS,NOAA AVHRR等衛星產品。平臺提供Python和JavaScript兩種客戶端庫，用戶可以自由在GEE代碼編輯器進行代碼編寫，之后傳至云端進行大量數據的并行計算，提高了遙感數據的處理效率。相較于傳統的本地計算模式，GEE的高存儲以及云計算能力，使大量遙感數據的快速處理和保存得以實現。

本次研究數據源為GEE平臺提供的1986-01-01—2019-12-30的Landsar SR數據，這些數據已經進行了幾何校正，之后去除50%以上高云像素的數據，最后根據研究區位置共選取了1 460個可用數據。首先使用LandTrendr方法提取變化水體(人工挖掘水體和沉陷水體)，然后采用形態學方法提取出沉陷水體，詳細過程請參考文獻[40-41]，最終得到1989—2017年共29年的沉陷水體數據。選擇1989—2016年數據用于沉陷水體空間位置變化研究，并將2017年沉陷水體數據作為真實值進行精度評價。

由于初始的沉陷水體數據存在細小面積斑塊，且本次研究注重沉陷水體最外層邊界的變化，不考慮區塊內部的其他土地類型，因此選取半徑為2個像素的圓形核對初始數據進行先腐蝕后膨脹的形態學開操作(圖2)，用以消除細小面積斑塊，并且將研究區內由于道路用地等其他土地類型產生的狹長區域一起融入積水斑塊，能減少數據量，處理后的結果如圖3所示。

圖2 形態學開操作前后對比Fig.2 Morphological comparison before and after operation

圖3 潘謝礦區1989—2016年沉陷水體分布Fig.3 Distribution of subsided water bodies in Panxie mining area from 1989 to 2016

2 研究方法

本研究以射線法為核心，提出了采煤沉陷水體方向變化自動識別方法，具體技術流程如圖4所示。

圖4 技術流程Fig.4 Technical flow chart

圖5 射線法示意Fig.5 Schematic diagram of ray method

2.1 射線法基本原理

對于單個沉陷水體區塊而言，射線法實現的基本步驟為：① 選取一個原點，并以正東方向為0°線；② 選擇一定的間隔角度，沿逆時針方向依次發射射線；③ 記錄原點到射線相交于該區塊各年的沉陷水體邊界的距離數據。如圖5所示，以丁集礦某沉陷積水區3 a的變化為例，選取角度為的射線，其中內部的紅色邊界代表沉陷水體某區塊最初產生年(2013年)的邊界，次外層的綠色代表該區塊2014年的邊界，最外層藍色代表該區塊2015年的邊界。圖5中角度為的射線分別交3 a的沉陷水體邊界于,，三點，記錄發射原點到各交點的距離列表，即,,的距離。按照這種規則，將研究區內所有沉陷水體區塊均進行射線法操作，最后得出每一區塊各方向上的距離列表(列表長度=2016-+1，其中為最高產生的年份)。

2.2 射線原點的選取

在采用射線法進行沉陷水體各區塊邊界距離變化提取時，選取的射線原點對結果有很大的影響，因此選取適當的射線原點十分關鍵。由于研究區開采歷史很長，缺少采煤的相關信息，無法以各礦最初的開采位置為原點進行分析，因此，對于每個沉陷水體區塊而言，射線原點的選取約定2個原則：一是在進行射線法擬合時，擬合結果應最接近該區塊的演化過程；二是以區塊的最初產生年為起始年開始研究，能最大程度地接近礦區最初開采位置，原因是在觀察研究區沉陷水體的演化過程中發現，各沉陷水體區塊是由最初產生的多個小面積區塊逐漸演化合并而成的。

以某一區塊為例，首先以2016年最終的區塊為基礎，觀測其是由哪些最初形成的小面積區塊演化合并而成的，并以這些最初的小面積區塊的中心點為射線原點，分別采用射線法構建距離列表，并對每個射線原點采用決定系數來評價，最后選取最大值，即擬合程度最好的點作為該區塊的射線原點。其中決定系數的計算公式為

=

(1)

其中，為預測數據與原始數據均值之差的平方和；為原始數據與其均值之差的平方和。取值為0～1，取值越接近1，表明方程的自變量對因變量的解釋能力越強，這個模型對數據擬合程度越高。和可通過以下公式計算：

(2)

(3)

2.3 射線間隔角度的選取

采用射線法對沉陷水體進行監測時，射線間隔角度與兩條射線間的沉陷水體邊界形狀復雜性成正比，即沉陷水體邊界形狀越復雜，選取的射線間隔角度越小。因此需要選取合適的射線間隔角度。

在采用射線法進行沉陷水體邊界擬合時，由于相鄰2條射線選取的間隔角度較小，可以將射線間曲線型的沉陷水體邊界用圖6中的直線段來進行擬合。其中，距離差越小，說明這2條射線夾角范圍內的沉陷水體邊界變化程度越小，則用一條直線段進行擬合的可靠性越強。研究中發現將1°設定為射線間隔角度時，兩射線間的沉陷水體邊界形狀如圖6(a)所示，即由下面的射線至上面的射線監測過程中，沉陷水體邊界一直保持向外(右)擴張，趨勢不變，此時可用直線擬合沉陷水體邊界；圖6(b)中沉陷水體邊界是先向外(右)擴張至點，再向內(左)擴張的形狀，即趨勢存在變化。當趨勢存在多次變化時，用直線擬合沉陷水體邊界精度較低，就需要在這2條射線間加密一條甚至多條射線。

圖6 相鄰射線間沉陷水體邊界變化情況Fig.6 Boundary changes of subsided water bodies between adjacent rays

為了確定合適的射線間隔角度，首先選擇面積較大的5個沉陷水體，從產生年到2016年，以1°為間隔構建射線，統計了相鄰射線間的距離差，共計18 360個數據，并按從小到大的順序排列編號。如圖7所示，當距離差大于103.2 m時發生突變，故認為距離差小于103.2 m的部分是正常的，距離差大于103.2 m的部分視為異常值，即需要進一步加密射線的部分。之后共選取0.25°,0.5°,1°,2°以及4°為備選角度，分別統計其各自總體上相鄰射線距離差之和以及距離差大于103.2 m(異常值)的數量，隨后將0.25°和0.5°組成第1組，0.5°和1°組成第2組，以此類推共4組。以第1組為例，假設間隔角度為0.5°時，異常值的數量為，間隔角度為0.25°時，異常值的數量為，則-代表射線在由間隔角度為0.5°加密到0.25°過程中對異常值的消減程度，最后選擇一個既能最大程度地減小相鄰射線距離差之和又能大量消減異常值的角度。

圖7 相鄰射線距離差統計Fig.7 Statistics of distance difference between adjacent rays

2.4 距離變化數據與年份的相關性檢驗

為了確定沉陷積水區塊是否持續變大，對距離變化與年份之間進行了相關性分析。可直接采用GEE平臺提供的皮爾遜相關性分析(ee.Reducer.pearsonsCorrelation())，對每條射線進行相關性檢驗。其中，在每條射線上，將射線原點到邊界的距離()和相應的年份()作為2個變量。

(4)

其中，為和的皮爾遜相關性;Cov(,)為和的協方差;Var(),Var()分別為和的方差。取值為-1～1，其中正負號分別代表正相關和負相關，的絕對值越接近1，說明和的相關性越強。其中Cov(,)可以通過以下公式計算：

(5)

Var()=Cov(,)

(6)

Var()=Cov(,)

(7)

2.5 預測沉陷水體提取

為了獲得每個沉陷水體各方向的變化規律，利用射線法計算出的各沉陷水體區塊各方向的距離列表數據，采用一元線性最小二乘回歸方法進行沉陷水體的邊界擬合及趨勢分析。GEE提供了線性擬合的ee.Reducer.linearFit()功能。

=+

(8)

其中，為自變量，表示時間(區塊產生年份)；為因變量，表示沉陷水體邊界到原點距離的回歸值；為斜率，表示沉陷水體邊界到原點距離的變化趨勢，>0，越大說明邊界擴張速度越快；為截距，表示水體最初產生年份邊界到原點距離。和可以通過以下公式計算：

(9)

(10)

在每個距離列表中，根據其對應的變化規律來預測下一年即2017年的距離，之后結合該射線的角度，計算出該角度下2017年沉陷水體邊界點，將每個區塊按照其所有射線方向計算出的邊界點依次連接起來，作為預測的2017年沉陷水體邊界結果。

2.6 沉陷水體擴張程度分析

基于得到的2017年預計沉陷水體數據和2016年沉陷水體數據，通過擴張系數來衡量各沉陷水域的擴張速度。

(11)

其中，和分別為每條射線上2016年沉陷水體邊界到射線原點的距離以及預測的2017年沉陷水體邊界到射線原點的距離。其中>0,且的取值越小說明沉陷水體在該方向的擴張速度越慢，沉陷水體越接近穩定。統計研究區內所有沉陷水體區塊各方向射線的擴張系數，并進行分級，分級標準見表1。

表1 擴張系數分級標準Table 1 Expansion coefficient grading standard

3 結果分析

研究以1989—2016年沉陷水體為基礎數據，在對其進行射線法監測時，先將沉陷水體所有區塊從其產生年至2016年間所有沉陷水體邊界數據按年份整合。

3.1 射線原點選取結果

至2016年，潘謝礦區各個小面積沉陷水體已經演化成面積可觀的共41塊沉陷水體，各區塊選取的射線原點結果及相應的決定系數見表2。

表2 各沉陷水體區塊選取的射線原點情況及決定系數Table 2 Ray origin situation and determination coefficient selected in each mining area

由表2可看出，41塊沉陷水體每個區塊的決定系數均大于0.6。統計每個區塊的面積并將其作為權重，將決定系數進行加權平均操作，計算公式為

(12)

計算出加權平均決定系數為84.56%，回歸方程擬合程度較好。

3.2 射線間隔角度選取結果

由上所述，選取0.25°，0.5°，1°，2°，4°為備選角度，分別統計其各自總體上相鄰射線距離差之和以及距離差大于103.2 m(異常值)的數量，見表3。

表3 射線間隔角度選取Table 3 Ray interval angle selection

可見，從0.5°劃分到0.25°時，相鄰射線距離差之和減少了4 087.8 m，消除了7個異常值；從1°劃分到0.5°時，相鄰射線距離差之和減少了3 299.8 m，消除了16個異常值；從2°劃分到1°時，相鄰射線距離差之和減少了7 522.5 m，消除了76個異常值；從4°劃分到2°時，相鄰射線距離差之和減少了3 019.3 m，消除了20個異常值。因此，當射線間隔劃分到1°時，相鄰射線的距離差之和減少得最多，且能最大程度地消減異常值，故最終選取1°為射線間隔角度。

3.3 皮爾遜相關性檢驗結果

經過對沉陷水體邊界至原點的距離與其相應的年份進行皮爾遜相關性分析，結果如圖8所示，其中86.33%以上的射線相關性系數達到0.8以上，67.45%的相關系數達0.9以上，說明每條射線上，原點距該沉陷水體邊界距離與相對應年份具有較強的相關性，之后進行的回歸分析可靠性強。

3.4 預測沉陷水體提取及其對比分析

基于前面的線性擬合提取出2017年預測沉陷水體，并將其與通過遙感影像提取出的2017年沉陷水體進行疊加(圖9)，統計出2者相交的面積和各自的面積，并計算精度。精度計算公式為

=

(13)

式中，為預測沉陷水體的精度；為真實2017年沉陷水體中預測正確部分的面積；為真實2017年沉陷水體的面積。

精度評價的結果見表4。其中面積較大的礦區，如謝橋礦、張集礦、顧橋礦、潘三礦和潘一礦，其預測精度均在80%以上，潘謝礦區總體預測與遙感影像提取出的相交面積為8 836.68 hm，遙感影像提取面積為10 466.56 hm，總體預測精度為84.43%。

圖8 沉陷積水區邊界變化與相應年份的相關性統計Fig.8 Correlation statistics between the boundary change of the subsidence water area and the corresponding year

圖9 潘謝礦區2017年沉陷水體邊界Fig.9 Boundary of subsided water body in Panxie mining area in 2017

表4 潘謝礦區沉陷水體預測與遙感提取面積比較Table 4 Comparison of subsidence water body prediction and real area in panxie mining area

3.5 沉陷水體擴張性分析

沉陷水體各方向的擴張性分析結果如圖10所示。圖10(a)中A區域整體呈現藍色，沉陷水體擴張速度慢，沉陷水體發育較完全，穩定性強；B區域大面積呈現紅色和黃色，擴張速度中等偏慢，處于較慢的發育階段，整體較穩定；C區域擴張速度較慢，穩定性強；D區域大部分呈現棕色，處于急劇擴張狀態，活躍程度高；E區域大面積呈現綠色，擴張速度較快，活躍性較強。圖10(b)中A區域大面積呈現紫色和棕色，擴張速度快，活躍性強；B,C，E三個區域擴張速度均較慢，整體較穩定；D區域大部分為紫色，擴張速度快，活躍性強。圖10(c)中A,C，D,E，I,J六個區域擴張速度較慢，其中C,D，E,I均處于較穩定的狀態，A,J擴張速度緩慢，有較弱的穩定性；B區域從東北方向沿逆時針至西南方向沉陷水體擴張速度較慢，其余的部分擴張速度較快；F區域左邊擴張速度快，右邊擴張速度中等偏慢；G區域大面積呈綠色，擴張速度較快，活躍程度高；H區域處于急劇擴張狀態，極不穩定，活躍性強；K區域從東北方向沿逆時針至西方向區間內，擴張速度中等偏快，剩余區間擴張速度慢。

謝橋礦、潘三礦、潘北礦、潘二礦、潘一礦建成投產時間分別為1997年、1992年、2004年、1989年、1983年，開采歷史較長，且近幾年開采量較小，沉陷水體擴張速度慢。張集礦、顧北礦、丁集礦、朱集礦建成投產時間分別為2001年、2007年、2007年、2007年，且其中張集礦開采時間較晚，這些礦區近幾年仍處在煤炭大量開采階段，沉陷水體擴張速度較快。顧橋礦于2007年建成投產，且處在煤炭大量開采階段，其西北方向擴張速度較快、東南方向擴張速度慢，但圖3中顯示，該礦區沉陷水體是由多個平行“條帶”狀沉陷水體匯聚而成，射線法對連續擴張型沉陷水體的監測效果較好，而對前面這種沉陷水體的演化方式監測效果一般。

圖10 潘謝礦區沉陷水體擴張速度Fig.10 Expansion speed of subsidence water body in Panxie mining area

4 結 論

(1)提出了射線法應用于采煤沉陷水體空間位置變化監測時原點和射線角度間隔確定方法。其中原點的選取是以2016年獨立的沉陷水體為基礎，選取出各獨立沉陷水體最初形成的小面積區塊，再從這些小面積區塊中基于決定系數篩選出最能代表(擬合)這一獨立沉陷水體演化過程的區塊，將這一區塊中心作為射線原點；射線間隔角度選擇既能最大程度地減小相鄰射線距離差之和又能大量消減異常值的角度。據此，本研究對潘謝礦區41個沉陷水體區塊選取了各自的射線原點，并以1°為射線間隔，構建了共14 760條射線。

(2)通過對構建出的射線統計出的年際距離變化數據進行皮爾遜相關性分析，發現射線原點距該沉陷水體邊界距離與其相對應年份具有較強的相關性，其中86.33%以上的射線相關性系數達到0.8以上，67.45%的射線相關系數達0.9以上。

(3)構建一元線性最小二乘法回歸方程擬合了2017年沉陷水體邊界，總體決定系數為84.56%，擬合程度良好。并將預測出的2017年沉陷水體數據與遙感影像提取的2017年沉陷水體數據進行對比，預測精度為84.43%。

(4)利用擴張系數評價了沉陷水體各方向的擴張速度。其中，謝橋礦、潘三礦、潘北礦、潘二礦、潘一礦沉陷水體擴張速度慢，張集礦、顧北礦、丁集礦、朱集礦沉陷水體擴張速度較快，顧橋礦西北方向擴張速度較快、東南方向擴張速度慢，與礦山企業的開采情況基本對應。

在缺少礦區采煤相關信息的條件下，本研究將射線法的思想引入采煤沉陷積水的研究中，提供了一個利用遙感影像數據，對采煤沉陷水體的空間位置進行擬合和預測的新思路。實現了在GEE平臺上射線的自動繪制及應用，從沉陷水體的提取到射線法監測，整個研究基本實現自動化，只需使用者選取合適的遙感影像數據集、劃定研究區即可。

由于影響采煤沉陷水體變化的因素不僅僅是地下采煤，比如自然降水、區域水資源調配等，對預測精度有一定影響，后續可考慮加入改進。