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

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

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

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

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

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

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況

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

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

1.2 數(shù)據(jù)源及數(shù)據(jù)處理

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

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

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

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

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

2 研究方法

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

圖4 技術(shù)流程Fig.4 Technical flow chart

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

2.1 射線法基本原理

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

2.2 射線原點的選取

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

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

=

(1)

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

(2)

(3)

2.3 射線間隔角度的選取

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

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

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

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

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

2.4 距離變化數(shù)據(jù)與年份的相關(guān)性檢驗

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

(4)

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

(5)

Var()=Cov(,)

(6)

Var()=Cov(,)

(7)

2.5 預測沉陷水體提取

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

=+

(8)

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

(9)

(10)

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

2.6 沉陷水體擴張程度分析

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

(11)

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

表1 擴張系數(shù)分級標準Table 1 Expansion coefficient grading standard

3 結(jié)果分析

研究以1989—2016年沉陷水體為基礎(chǔ)數(shù)據(jù)，在對其進行射線法監(jiān)測時，先將沉陷水體所有區(qū)塊從其產(chǎn)生年至2016年間所有沉陷水體邊界數(shù)據(jù)按年份整合。

3.1 射線原點選取結(jié)果

至2016年，潘謝礦區(qū)各個小面積沉陷水體已經(jīng)演化成面積可觀的共41塊沉陷水體，各區(qū)塊選取的射線原點結(jié)果及相應(yīng)的決定系數(shù)見表2。

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

由表2可看出，41塊沉陷水體每個區(qū)塊的決定系數(shù)均大于0.6。統(tǒng)計每個區(qū)塊的面積并將其作為權(quán)重，將決定系數(shù)進行加權(quán)平均操作，計算公式為

(12)

計算出加權(quán)平均決定系數(shù)為84.56%，回歸方程擬合程度較好。

3.2 射線間隔角度選取結(jié)果

由上所述，選取0.25°，0.5°，1°，2°，4°為備選角度，分別統(tǒng)計其各自總體上相鄰射線距離差之和以及距離差大于103.2 m(異常值)的數(shù)量，見表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 皮爾遜相關(guān)性檢驗結(jié)果

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

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

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

=

(13)

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

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

圖8 沉陷積水區(qū)邊界變化與相應(yīng)年份的相關(guān)性統(tǒng)計Fig.8 Correlation statistics between the boundary change of the subsidence water area and the corresponding year

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

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

3.5 沉陷水體擴張性分析

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

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

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

4 結(jié) 論

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

(2)通過對構(gòu)建出的射線統(tǒng)計出的年際距離變化數(shù)據(jù)進行皮爾遜相關(guān)性分析，發(fā)現(xiàn)射線原點距該沉陷水體邊界距離與其相對應(yīng)年份具有較強的相關(guān)性，其中86.33%以上的射線相關(guān)性系數(shù)達到0.8以上，67.45%的射線相關(guān)系數(shù)達0.9以上。

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

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

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

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