樂 穎，夏元平，劉媛媛，胡 昊

(東華理工大學(xué)測(cè)繪工程學(xué)院，江西 南昌 330013)

0 引言

20世紀(jì)70年代左右，合成孔徑雷達(dá)干涉(interferometric synthetic aperture radar，InSAR)技術(shù)出現(xiàn)，Rogers A E等于1969年首次利用 InSAR 技術(shù)研究了金星表面的地形特征[1]。隨后，Gabriel等于1989年首次提出合成孔徑雷達(dá)差分干涉測(cè)量(Differential InSAR，DInSAR)技術(shù)[2]。該技術(shù)是基于InSAR技術(shù)提取地表形變信息的一門極具潛力的新興技術(shù)，由于在一定程度上克服了光照和天氣的影響，可全天時(shí)全天候監(jiān)測(cè)，具有高精度、大區(qū)域、快速等優(yōu)勢(shì)，因此受到了學(xué)者的廣泛關(guān)注[3]。目前，DInSAR技術(shù)精度甚至可達(dá)毫米級(jí)，己經(jīng)成為地震、火山、滑坡、地表形變信息等應(yīng)用領(lǐng)域的主要監(jiān)測(cè)手段[4]。

高分辨率SAR衛(wèi)星的快速發(fā)展，使DInSAR技術(shù)已經(jīng)成為目前獲取地表形變信息的研究熱點(diǎn)[5]。文獻(xiàn)[6]比較了兩種配準(zhǔn)方法的特點(diǎn)，并將DInSAR技術(shù)應(yīng)用于神木礦區(qū)的沉陷監(jiān)測(cè)，獲取最新的礦區(qū)沉陷信息及對(duì)應(yīng)的時(shí)間序列。文獻(xiàn)[7]利用DInSAR技術(shù)監(jiān)測(cè)呂梁山區(qū)某礦區(qū)地表，通過兩個(gè)剖面數(shù)據(jù)分析地表沉陷規(guī)律，表明能夠得到對(duì)應(yīng)的形變信息。文獻(xiàn)[8]將DInSAR技術(shù)運(yùn)用于陜西彬長(zhǎng)礦區(qū)，并與實(shí)際情況相結(jié)合，有效獲取了礦區(qū)的形變信息。文獻(xiàn)[9]基于SAR 影像的振幅信息提出了一種改進(jìn)的PS探測(cè)方法，通過對(duì)上海市局部區(qū)域進(jìn)行研究，準(zhǔn)確地提取了研究區(qū)域的形變信息。文獻(xiàn)[10]基于DInSAR技術(shù)，結(jié)合GPS水平形變結(jié)果，可快速獲得研究區(qū)地表三維形變，為后期的研究提供有效的參考。文獻(xiàn)[11]針對(duì)不均勻地表形變問題，利用TS-DInSAR(time series DInSAR)技術(shù)對(duì)Sentinel-1A數(shù)據(jù)進(jìn)行研究，結(jié)果表明該方法對(duì)大范圍地表形變監(jiān)測(cè)具有可靠的精度。然而，以上文獻(xiàn)為了校正衛(wèi)星軌道和相位偏移，去除處理過程中殘余相位，在運(yùn)用DInSAR技術(shù)處理影像時(shí)憑借人為主觀經(jīng)驗(yàn)選取控制點(diǎn)進(jìn)行軌道精煉，這將影響相位信息轉(zhuǎn)化為形變信息的準(zhǔn)確性。若研究區(qū)域缺少地面控制點(diǎn)信息時(shí)，隨意選取控制點(diǎn)將影響變形信息結(jié)果提取的精度。本文針對(duì)此問題，提出基于相干系數(shù)的DInSAR地表形變信息提取方法。

1 DInSAR技術(shù)

1.1 基本原理

DInSAR形變監(jiān)測(cè)技術(shù)是經(jīng)由InSAR技術(shù)衍生發(fā)展而來，主要包括二軌、三軌與四軌三種。其中，“二軌法”差分干涉技術(shù)運(yùn)用領(lǐng)域最廣，其表達(dá)式如下[12]：

φ=φflat+φtopo+φdefo+φatm+φnoise

(1)

式(1)中，φflat為平地相位；φtopo為地形起伏引起的地形相位，可以將地形信息恢復(fù)；φdefo為地表發(fā)生形變?cè)斐傻南辔?；φatm為大氣發(fā)生延遲的大氣相位；φnoise為觀測(cè)時(shí)產(chǎn)生的噪聲相位[12]。簡(jiǎn)而言之，就是將式(1)中φflat、φtopo、φatm與φnoise等相位去除，最后只保留φdefo形變相位的技術(shù)[13]。

1.2 數(shù)據(jù)處理流程圖

由于雷達(dá)數(shù)據(jù)與光學(xué)影像不同，屬于主動(dòng)遙感，影像包含的是振幅、相位和極化等多種信息。目前，常用的技術(shù)方法主要有合成孔徑雷達(dá)干涉、極化分析、幅度追蹤、層析建模和立體量測(cè)等五種方法，而其中合成孔徑雷達(dá)干涉方法是最常用也是運(yùn)用范圍最廣的技術(shù)。合成孔徑雷達(dá)差分干涉通過處理覆蓋同一地區(qū)不同時(shí)刻獲取的兩幅SAR影像的相位信息，能夠測(cè)量地表厘米級(jí)的形變。在DInSAR數(shù)據(jù)處理過程中，主要包括：影像配準(zhǔn)、干涉圖生成、干涉圖自適應(yīng)濾波和相干性計(jì)算、最小費(fèi)用流相位解纏、控制點(diǎn)選取、軌道精煉、重去平、相位轉(zhuǎn)形變及地理編碼等，每個(gè)數(shù)據(jù)處理步驟過程中產(chǎn)生的誤差對(duì)干涉結(jié)果的精度和可靠性都具有顯著影響，具體的處理流程見圖1。DInSAR地表形變監(jiān)測(cè)的數(shù)據(jù)準(zhǔn)備主要包括SAR數(shù)據(jù)選取，外部DEM和精密軌道數(shù)據(jù)準(zhǔn)備。其中，外部DEM是為了利用其模擬地形起伏相位分量，精密軌道數(shù)據(jù)則可以有效削弱SAR衛(wèi)星軌道誤差導(dǎo)致的殘余參考相位。

圖1 DInSAR數(shù)據(jù)處理流程圖Fig.1 DInSAR data processing flow chart

2 基于相干性系數(shù)的 DInSAR 地表形變信息提取方法

2.1 干涉相干性測(cè)度

相干圖主要用于體現(xiàn)兩幅影像的相關(guān)程度，可以用于評(píng)價(jià)SAR影像精確配準(zhǔn)后所得到干涉條紋圖的質(zhì)量，主要是通過計(jì)算干涉條紋圖的相干系數(shù)實(shí)現(xiàn)。相干系數(shù)最初由普拉蒂于1993年提出，相干系數(shù)的取值范圍為[0,1]，值越接近1說明相干性越高，其理論模型表達(dá)式如下[14]：

(2)

式(2)中，E[·]表示數(shù)學(xué)期望，u*表示共軛復(fù)數(shù)，u1與u2為主副圖像的信號(hào)。在實(shí)際的InSAR處理過程中，在同一像素中無法得到計(jì)算所需數(shù)量的采樣值。因此，根據(jù)式(2)的理論模型，基于SAR影像復(fù)數(shù)數(shù)據(jù)計(jì)算獲得相干系數(shù)的標(biāo)準(zhǔn)表達(dá)式為[15]：

(3)

式(3)中，u1(n,m)，u2(n,m)分別表示主、副影像數(shù)據(jù)塊內(nèi)某個(gè)坐標(biāo)(n,m)處的復(fù)數(shù)值；|?|2表示對(duì)應(yīng)的二階范數(shù)；M與N表示計(jì)算相干性的數(shù)據(jù)塊尺寸大?。籱與n表示數(shù)據(jù)塊內(nèi)對(duì)應(yīng)的行列號(hào)。

2.2 軌道精煉控制點(diǎn)的選取

利用DInSAR技術(shù)獲取地表形變信息時(shí)，引入穩(wěn)定的地面控制點(diǎn)可以有效估計(jì)殘留相位并去除平地效應(yīng)，從而提高結(jié)果的可靠性。一般而言，選取控制點(diǎn)要遠(yuǎn)離形變區(qū)，認(rèn)為形變?yōu)?，陡峭的地形區(qū)域和有殘余地形相位區(qū)域，地形起伏大的山區(qū)，最好是選擇山谷底部的平地區(qū)域[16]。為了更加快速準(zhǔn)確地選擇軌道精煉的控制點(diǎn)，本文基于相干系數(shù)閾值選取所需的控制點(diǎn)，相干性系數(shù)圖的像元分布情況見圖2。

從圖2可看出，2015—2020年干涉對(duì)像元個(gè)數(shù)在相干性系數(shù)為0.5時(shí)達(dá)到了峰值，2015、2017和2020年的整體相干系數(shù)數(shù)值比其他年份高。由于相干系數(shù)圖上的像元個(gè)數(shù)數(shù)量眾多，選取軌道精煉時(shí)需要的像元數(shù)量很少。相干性閾值設(shè)置過低會(huì)導(dǎo)致相關(guān)性高像元數(shù)目劇增，大大增加數(shù)據(jù)計(jì)算的時(shí)間，從而增加后續(xù)篩選工作的難度。因此，在保證高相干性又考慮數(shù)據(jù)量大小的基礎(chǔ)上，通過反復(fù)對(duì)比像元個(gè)數(shù)進(jìn)行設(shè)置合適的閾值，最終將相干系數(shù)閾值設(shè)定為0.98，每年干涉相干性系數(shù)的像元分布情況見表1。

圖2 2015—2020年相干性系數(shù)圖像元分布情況Fig. 2 Distribution of image elements of coherence coefficient from 2015 to 2020

表1 每年干涉相干性系數(shù)的像元分布情況

2.3 控制點(diǎn)選取驗(yàn)證

為了更系統(tǒng)地驗(yàn)證本文方法選取控制點(diǎn)的可靠性，首先統(tǒng)計(jì)2015—2020年期間每年干涉相干性系數(shù)圖符合要求的像元數(shù)量，然后根據(jù)相干系數(shù)閾值分別選取25～35個(gè)軌道精煉控制點(diǎn)，最后計(jì)算各點(diǎn)絕對(duì)殘余誤差，其計(jì)算公式為：

(4)

式(4)中，γ是干涉相干性測(cè)度，H是高程模糊度。

以2019年為例，對(duì)比分析DInSAR方法與本文方法用于選取軌道精煉控制點(diǎn)的絕對(duì)殘余誤差，結(jié)果見圖3。本文方法用于選取軌道精煉控制點(diǎn)整體的絕對(duì)殘余誤差較小，只有個(gè)別點(diǎn)的絕對(duì)殘余誤差偏大，誤差符合要求；而DInSAR方法選取控制點(diǎn)的絕對(duì)殘余誤差整體偏大，后續(xù)需要篩選絕對(duì)殘余誤差較大的點(diǎn)。此外，文中為了對(duì)比分析DInSAR方法與本文方法用于選取軌道精煉控制點(diǎn)的準(zhǔn)確性，引入了標(biāo)準(zhǔn)差指標(biāo)，并計(jì)算得到兩種方法絕對(duì)殘余誤差的標(biāo)準(zhǔn)差分別為1.33和0.36，而指標(biāo)值越小說明準(zhǔn)確性越高。因此，在缺少地面控制點(diǎn)信息或研究區(qū)形變先驗(yàn)認(rèn)識(shí)的情況下，若需要選擇可靠的地面控制點(diǎn)，可借助相干性系數(shù)圖進(jìn)行剔除低相干性點(diǎn)。

圖3 2019年選取軌道精煉控制點(diǎn)絕對(duì)殘余誤差Fig.3 Absolute residual error at the control point of orbit refining in 2019

3 實(shí)例驗(yàn)證

3.1 研究區(qū)域概況

贛州市的稀土儲(chǔ)量占全國(guó)最多，而定南縣是贛州市離子型稀土主產(chǎn)縣之一，足以表明定南縣稀土礦礦區(qū)扮演著重要的作用?；诖?，本文以定南縣作為研究區(qū)域，其地理坐標(biāo)為東經(jīng)114°47′49″E～115°22′48″E，北緯24°33′37″N～ 25°03′21″N，總面積為1 321 km2?？h區(qū)內(nèi)分布眾多河流，主要隸屬于定南水和桃江兩水系[17]。同時(shí)，定南縣地處南嶺與武夷成礦帶的中間部位，故而擁有良好的成礦條件，礦產(chǎn)資源非常豐富。在礦山開采過程中，露天礦坑開挖會(huì)直接破壞地表植被，嚴(yán)重影響著當(dāng)?shù)氐目沙掷m(xù)發(fā)展和生態(tài)功能恢復(fù)。因此，研究并提取定南縣的地表形變信息尤為重要[18]。

3.2 實(shí)驗(yàn)數(shù)據(jù)

由于雷達(dá)數(shù)據(jù)與光學(xué)影像不同，屬于主動(dòng)遙感，影像包含的是振幅、相位和極化等多種信息。目前，合成孔徑雷達(dá)干涉方法是最常用也是運(yùn)用范圍最廣的技術(shù)。本文使用的雷達(dá)影像數(shù)據(jù)為C波段的Sentinel-1A衛(wèi)星的SLC數(shù)據(jù)，選取了時(shí)間跨度為2015—2020年期間每年10月份共12景的影像數(shù)據(jù)，具體雷達(dá)影像數(shù)據(jù)的基本參數(shù)見表2。

表2 雷達(dá)影像數(shù)據(jù)的基本參數(shù)

3.3 結(jié)果分析

通過對(duì)12景覆蓋研究區(qū)的SAR影像數(shù)據(jù)進(jìn)行DInSAR處理，借助相干系數(shù)閾值獲取相干性高的像元，并將其作為軌道精煉的控制點(diǎn)引入后續(xù)的處理中，最終可得到贛州市定南縣鄉(xiāng)鎮(zhèn)在2015—2020年的地表形變信息，其形變結(jié)果圖如圖4所示。監(jiān)測(cè)結(jié)果表明：研究區(qū)在監(jiān)測(cè)時(shí)段內(nèi)地表形變變化相對(duì)穩(wěn)定，整體形變量主要集中在-0.07～0.1 m，部分區(qū)域出現(xiàn)了一定范圍的地表沉降，沉降量集中在-0.11～0.11 m，定南縣在2017、2019和2020年內(nèi)的地表形變較2015、2016和2018年更穩(wěn)定。

定南縣共有7個(gè)鎮(zhèn)，分別為鵝公鎮(zhèn)、黃香鎮(zhèn)、巋美山鎮(zhèn)、老城鎮(zhèn)、歷市鎮(zhèn)、嶺北鎮(zhèn)和天九鎮(zhèn)。根據(jù)雷達(dá)SAR影像數(shù)據(jù)處理得到的結(jié)果，得到定南縣各鎮(zhèn)每年的地表形變變化情況，所得結(jié)果如圖5所示。從圖5的監(jiān)測(cè)結(jié)果可以得出：定南縣各個(gè)鄉(xiāng)鎮(zhèn)在2017、2018、2019和2020年內(nèi)變化趨勢(shì)基本一致，2015年形變量較大的鄉(xiāng)鎮(zhèn)為嶺北鎮(zhèn)，形變量達(dá)到了-25.26 mm；2016年形變量較大的鄉(xiāng)鎮(zhèn)為嶺北鎮(zhèn)、黃香鎮(zhèn)和巋美山鎮(zhèn)，形變量分別為5.26、5.26和26.29 mm；整個(gè)研究區(qū)域在監(jiān)測(cè)時(shí)段內(nèi)最大沉降量和最大上升量分別為-84.82和84.816 mm，分別位于嶺北鎮(zhèn)和鵝公鎮(zhèn)。經(jīng)過查找資料得知，贛州市定南縣內(nèi)具有豐富的離子型吸附性稀土礦礦產(chǎn)資源，且這些礦區(qū)位置大多位于嶺北鎮(zhèn)，其中包含了木子山、甲子背、長(zhǎng)坑尾等稀土礦區(qū)[17]。由于嶺北鎮(zhèn)的平均形變量與其他鄉(xiāng)鎮(zhèn)的平均形變量相差-20.20～-10.10 mm，稀土礦的開采可能會(huì)引起少量的地表形變。

圖4 2015年-2020年研究區(qū)域地表形變結(jié)果圖Fig.4 Results of surface deformation in the study area from 2015 to 2020

圖5 2015—2020年定南縣鄉(xiāng)鎮(zhèn)地表形變信息Fig.5 Surface deformation information of towns in dingnan county from 2015 to 2020

為了探尋稀土礦的開采與地表形變之間的關(guān)系，并獲取稀土礦開采區(qū)域的地表形變細(xì)節(jié)信息和地表形變發(fā)展趨勢(shì)信息，本文結(jié)合稀土礦權(quán)在已有礦權(quán)范圍內(nèi)尋找了一些相關(guān)性高且能代表形變規(guī)律的特征點(diǎn)，特征點(diǎn)選取分布見圖6，最后對(duì)特征點(diǎn)進(jìn)行時(shí)間序列分析。

分別計(jì)算選取的30個(gè)特征點(diǎn)在2015—2020年間的地表形變量，圖7繪制了地表形變細(xì)節(jié)信息和發(fā)展趨勢(shì)。

圖6 30個(gè)特征點(diǎn)選取分布圖Fig.6 Distribution of 30 feature points

圖7 30個(gè)特征點(diǎn)的地表形變情況Fig.7 Surface Deformation of 30 feature points

從圖7中不難發(fā)現(xiàn)，選取的30個(gè)特征點(diǎn)在同年份的變化趨勢(shì)基本一致，在2015年的地表形變發(fā)展整體為下降趨勢(shì)，沉降量范圍在-25.5～-15.5 mm，結(jié)果表明這些稀土礦區(qū)均存在開采現(xiàn)象；2016年和2017年的地表形變?cè)?附近波動(dòng)，地表變化幅度分別為-5.26～15.77 mm和-12.82～4.27 mm，由于這段時(shí)間國(guó)家開始加大對(duì)非法開采行為的打擊力度，使得一些稀土私礦被查封；2018年整體呈下降趨勢(shì)，部分點(diǎn)為向上抬升，地表形變?cè)?32.984～14.136 mm之間，贛州地區(qū)對(duì)稀土礦山環(huán)境進(jìn)行了綜合治理以解決歷史遺留問題，將稀土礦區(qū)平整為建設(shè)用地投入到當(dāng)?shù)亟?jīng)濟(jì)發(fā)展中去；2019年和2020年的地表變化幅度分別為-26.03～18.07 mm和-3.81～3.81 mm，說明贛州定南地區(qū)開展的廢棄稀土礦山治理取得效果顯著，對(duì)當(dāng)?shù)氐纳鷳B(tài)環(huán)境有了明顯的修復(fù)。監(jiān)測(cè)結(jié)果表明：由于離子型吸附性稀土礦的礦區(qū)范圍不大，而且分布稀疏，從而導(dǎo)致稀土礦的開采引起的地表形變較緩慢。由于稀土礦的開采嚴(yán)重破壞了生態(tài)環(huán)境，江西贛州自2013年來不斷推進(jìn)稀土產(chǎn)業(yè)整頓，從加大稀土礦非法開采行為打擊力度到對(duì)生態(tài)環(huán)境進(jìn)行綜合治理，很大程度上促進(jìn)了礦山生態(tài)環(huán)境保護(hù)與稀土資源開發(fā)利用的協(xié)調(diào)發(fā)展。

4 結(jié)論

本文提出了基于相干性系數(shù)的 DInSAR地表形變信息提取方法。該方法通過設(shè)置一定的相干系數(shù)閾值，篩選出適合的高相干性像元作為控制點(diǎn)進(jìn)行軌道精煉，最后提取出地表形變信息。以贛州市定南縣作為研究區(qū)域，借助 2015—2020 年間的Sentinel-1A衛(wèi)星數(shù)據(jù)，利用相干性系數(shù)形變信息提取方法，監(jiān)測(cè)出研究區(qū)范圍內(nèi)-25.5～18.07 mm的地表形變信息。實(shí)驗(yàn)結(jié)果表明，該方法相較于常規(guī)的DInSAR方法，用于提取控制點(diǎn)的標(biāo)準(zhǔn)差提高了72.9%，定量表明了基于相干性系數(shù)的DInSAR方法用于選取控制點(diǎn)的準(zhǔn)確性更高，能夠?qū)崿F(xiàn)對(duì)研究區(qū)域地表形變信息的動(dòng)態(tài)監(jiān)測(cè)。由于本文采用的數(shù)據(jù)量不足，研究過程中沒有考慮大氣因素的影響，導(dǎo)致得到的結(jié)果仍表現(xiàn)出一定的偏差。將雷達(dá)影像數(shù)據(jù)與光學(xué)遙感影像相結(jié)合，并與外業(yè)水準(zhǔn)測(cè)量結(jié)果進(jìn)行對(duì)比，將是下一步要進(jìn)行的工作。