胡鳳鳴, 伍吉倉

(同濟(jì)大學(xué) 測繪與地理信息學(xué)院,上海 200092)

上海磁懸浮全長30 km,是世界第一條商用的高架磁懸浮運營線.磁懸浮工程在建設(shè)和運營期間對沿線的地面沉降控制要求嚴(yán)格,每一跨徑的沉降差異須小于3 mm[1]. 上海浦東地區(qū)第四紀(jì)松散土層厚300余米,極易受地下水開采以及工程活動的影響誘發(fā)地面沉降[2].為保證磁懸浮的安全運營,需要進(jìn)行高頻高密度高精度的沉降監(jiān)測.傳統(tǒng)的精密水準(zhǔn)測量由于空間密度較低,測量時間間隔大,難以滿足磁懸浮沿線沉降監(jiān)測需求.

合成孔徑雷達(dá)干涉測量(interferometric synthetic aperture radar,InSAR)是一種新興的空間對地觀測技術(shù),能夠全天候、大范圍、高精度獲取地球表面地形和形變信息.針對傳統(tǒng)InSAR技術(shù)時空間失相干以及大氣延遲影響[3]提出的多時序InSAR(MT-InSAR)可以有效提高測量精度.研究表明,通過分析處理高質(zhì)量的相干點可以獲得亞毫米級的形變監(jiān)測精度[4].該技術(shù)已經(jīng)被成功應(yīng)用于地殼形變[5]、城市地表[6-7]以及大型結(jié)構(gòu)體的監(jiān)測中[8-9].X波段的TerraSAR-X/TanDEM-X以及COSMO-Skyed系統(tǒng)獲取的高分辨率SAR數(shù)據(jù)使得該技術(shù)可以獲取更加精細(xì)的地表信息[10-11],這對于城市區(qū)域內(nèi)大型結(jié)構(gòu)體健康監(jiān)測具有重要意義.

MT-InSAR技術(shù)僅使用在時間序列上始終保持高相干性的散射體,依據(jù)干涉相位的時空特性分離不同的分量,從而獲取高精度的形變速率.因此選取點的質(zhì)量直接影響最終的解算結(jié)果,選取高相干點的方法可以分為基于點和弧段兩類,前者包括歸一化振幅離差[12-13],記為NAD,以DA表示,信噪比[14],相干性[15],相位分析[16],極大似然法[17],相干矩陣特征值分析[18]等,這些方法通過依據(jù)不同原則分析散射體在時間序列上相位的穩(wěn)定性來選取高相干點.后者首先進(jìn)行點分析得到較多的備選點集合[19-20],之后分析弧段質(zhì)量以及點之間的連通性得到高相干點.選取高質(zhì)量的相干點并保證空間點密度之間保持平衡是MT-InSAR技術(shù)的一大難點.寬松的閾值可以得到高密度的觀測結(jié)果,但是測量精度隨之降低,反之可以保證精度,但是觀測結(jié)果卻過于稀疏.因此需要研究針對不同質(zhì)量的相干點需采用不同的解算方式.目前的解決方案分為兩步,利用較小閾值選取較高質(zhì)量相干點,解算結(jié)果作為一級參考網(wǎng),保證解算精度;之后擴(kuò)大閾值,構(gòu)成二級網(wǎng)求解一般質(zhì)量相干點,以保證空間點密度[21-22].但是現(xiàn)有的空間點分級方法沒有考慮空間點密度的不均勻性,不同級相干點連接方式以及精度評價方法也有待完善,并且計算效率不高,使用原始分辨率數(shù)據(jù)的大范圍處理通常難以實施.

本文首先闡述了基于弧段選取的多級化MT-InSAR技術(shù)處理方法,包括基本觀測量,參考網(wǎng)和附加網(wǎng)的求解和精度評價方法,并對每一級網(wǎng)絡(luò)中相干點的選取以及求解策略作了詳細(xì)討論.利用C波段的ENVISAT、Sentinel-1A以及X波段的TerraSAR-X數(shù)據(jù)對上海磁懸浮及其沿線區(qū)域開展沉降監(jiān)測,分析不同數(shù)據(jù)的適用性,并與精密水準(zhǔn)數(shù)據(jù)驗證,為城市大型結(jié)構(gòu)體健康監(jiān)測提供參考.

1 多級化MT-InSAR技術(shù)

1.1 基本觀測量

在MT-InSAR技術(shù)中,基本觀測量為雙差相位觀測值.假設(shè)在時間序列上有m幅已配準(zhǔn)的合成孔徑雷達(dá)(synthetic aperture radar,SAR)數(shù)據(jù),選取其中一幅作為主影像后可以得到m-1幅干涉相位圖.干涉圖中任意兩個像元i與j構(gòu)成弧段,去除平地相位后,第k幅干涉圖中弧段的相位差可以表示為[23]

(1)

式中:Δhi,j和Δvi,j分別為點i與j之間的高程改正之差和形變速率差;ni,j,k是相位模糊度;Bcik和tk為垂直基線和時間基線;Ri為像元點i與衛(wèi)星間的斜距;λ為雷達(dá)波長;θi為入射角;Ci,j為相位常數(shù),與主影像的大氣延遲有關(guān);ei,j,k為相位噪聲.對式(1)相位解纏后使用最小二乘估計即可得到參數(shù)解和對應(yīng)的參數(shù)精度[23].

在基于弧段選取的MT-InSAR技術(shù)中,初始相干點的選取包含大量的噪聲點,最終可靠相干點的篩選依賴于選取高質(zhì)量弧段解.弧段質(zhì)量一般采用時序相干性來衡量[9].由于模型中采用了線性形變速率的假設(shè),因此加入線性速率的假設(shè)檢驗進(jìn)一步篩選弧段.此二者構(gòu)成文獻(xiàn)中的弧段可接受條件[20].

基于本節(jié)得到的弧段解進(jìn)行多級化相干點構(gòu)網(wǎng)求解.由于城市區(qū)域的相干點之間相關(guān)性較低,因此采用歸一化振幅離差方法進(jìn)行相干點選取,該方法在選取相干點時僅考慮相干點的時間序列噪聲,并通過選取不同的閾值來控制選取相干點的質(zhì)量,選點的閾值越低,獲取點的精度越高,但是點的數(shù)目越少.數(shù)據(jù)處理的主要流程見圖1.

圖1 多級化相干點處理流程

Fig.1 Flowchart of processing coherent points with multi-level strategy

1.2 參考網(wǎng)

參考網(wǎng)的求解精度直接決定了最終結(jié)果的精度,構(gòu)成該網(wǎng)絡(luò)的相干點應(yīng)當(dāng)保證有足夠高的精度.因此采用歸一化振幅離差法選取相干點時使用盡可能低的閾值.考慮到大氣效應(yīng)的影響,求解弧段不宜過長[13],參考網(wǎng)中相干點在空間內(nèi)應(yīng)盡可能均勻分布.依據(jù)這兩個準(zhǔn)則,參考網(wǎng)相干點選取的閾值使用0.25,是經(jīng)典永久散射體干涉測量方法中高相干點的選取閾值[12].這一步求解的目標(biāo)是將盡可能多的點連接起來,采用二次構(gòu)網(wǎng)策略.首先使用Delaunay三角網(wǎng)連接所有點,篩選弧段后,對孤立的點集采用漸進(jìn)式星形網(wǎng)構(gòu)成新的弧段.詳細(xì)處理步驟可參見文獻(xiàn)[20].

(2)

式中:B1為構(gòu)成參考網(wǎng)弧段相關(guān)的設(shè)計矩陣;Δv1為求解的所有弧段形變速率之差;(·)+是求解廣義逆運算符;P1為弧段解的權(quán)陣,可采用如下定義:

(3)

參考網(wǎng)的形變速率估計精度計算方式為

(4)

式中:σ1為參考網(wǎng)的單位權(quán)中誤差.

當(dāng)測區(qū)范圍較大或是處理高分辨率數(shù)據(jù)時,法方程系數(shù)矩陣階數(shù)較大,廣義逆求解緩慢或是由于矩陣過大而難以求解,因此需要減少構(gòu)成參考網(wǎng)的相干點數(shù)目.若使用降低閾值減少相干點的方式會導(dǎo)致相干點在空間分布不均勻,造成最終結(jié)果的精度在空間不均勻分布.此外,附加網(wǎng)在缺少參考網(wǎng)的相干點的區(qū)域無法生成弧段,從而造成相干點的遺漏.因此,本文采用基于點密度的分層抽樣方法來對空間點進(jìn)行稀疏,在保持相干點均勻分布的同時,減少空間點密度較大區(qū)域的相干點數(shù)目.抽樣后得到的相干點組成參考網(wǎng),余下的點使用附加網(wǎng)進(jìn)行求解.

1.3 附加網(wǎng)

采用寬松閾值下的歸一化振幅離差方法選取相干點,使用星形網(wǎng)將這些點與參考網(wǎng)中的點連接構(gòu)成新的弧段.構(gòu)網(wǎng)時仍然采用漸進(jìn)式策略,每次將待求點與若干個參考網(wǎng)中的點連接,循環(huán)處理直到該點與參考網(wǎng)之間有弧段連接,否則剔除該點.聯(lián)合參考網(wǎng)中相干點的參數(shù)對新求解的弧段進(jìn)行積分可以得到附加網(wǎng)中相干點的形變速率:

(5)

式中:B2為構(gòu)成附加網(wǎng)弧段相關(guān)的設(shè)計矩陣;Δv2為求解的所有弧段相對形變速率之差;P2為權(quán)陣,與P1為相同定義方式;k為拉格朗日函數(shù)的系數(shù)矩陣;v2是附加網(wǎng)點待求解的形變速率;G為限制矩陣,定義如下:

Gv2=v1

(6)

附加網(wǎng)的形變速率的精度計算方式為

(7)

式中:σ2為附加網(wǎng)的單位權(quán)中誤差,Ql定義如下:

(8)

其中,B11和B22是B2的子矩陣,定義如下:

B2=[B11|B22]

(9)

求解參考網(wǎng)后,附加網(wǎng)可以使用任意多次以獲取更多相干點的信息.求解附加網(wǎng)時,不同相干點的處理是相互獨立的.當(dāng)區(qū)域較大時,可以將整個區(qū)域分為獨立的若干塊,對每塊區(qū)域內(nèi)的相干點依次處理,從而減少運算的內(nèi)存開銷,適用于大范圍的數(shù)據(jù)處理.

2 數(shù)據(jù)處理結(jié)果

本次實驗采用2007~2018年間的31景C波段ENVISAT、53景X波段TerraSAR-X和51景C波段Sentinel-1A影像,詳細(xì)數(shù)據(jù)參數(shù)見表1.不同衛(wèi)星平臺數(shù)據(jù)均完全覆蓋磁懸浮沿線區(qū)域.由于上海地勢較為平坦,且樓房眾多,處理中不再使用地面高程數(shù)據(jù)進(jìn)行初始高程相位改正,僅改正平地相位.地面點高程引起的干涉相位在求解中通過時間序列方法進(jìn)行分離.

表1 SAR數(shù)據(jù)主要參數(shù)

以高分辨率TerraSAR-X數(shù)據(jù)為例,驗證多級化數(shù)據(jù)方法的有效性.構(gòu)成參考網(wǎng)的相干點采用經(jīng)典的永久散射體[12]判斷方法進(jìn)行選取,DA閾值取為0.25,該閾值可以保證選出的相干點在整個空間均有分布.選取初始相干點數(shù)目為472 929個,點分布參見圖2a,從圖中可以發(fā)現(xiàn)相干點在磁懸浮沿線呈現(xiàn)明顯的線性分布規(guī)律,磁懸浮周邊房屋區(qū)域也存在著高密度的相干點.初始相干點數(shù)目多,難以直接求解,需要降低參考網(wǎng)相干點數(shù)目.

根據(jù)1.2節(jié)中的處理策略,首先使用降低閾值法,將初始選點的閾值更改為0.1,可以得到17 680個相干點,點分布可參見圖2b.從該圖可以發(fā)現(xiàn),相干點在空間的分布很不均勻,磁懸浮沿線幾乎不存在相干點.在多級化數(shù)據(jù)處理中,估算精度隨著分級的進(jìn)行而逐步遞減,若磁懸浮沿線的相干點只存在于附加網(wǎng)中,那么這些點的精度相對于整個測區(qū)是偏低的.因此降低閾值選取參考網(wǎng)的相干點是不合理的.

使用基于空間點密度的分層抽樣方法來合理選擇參考網(wǎng)的規(guī)模.基于圖2a的初始相干點集,先計算相干點的空間密度.取100 m為搜索半徑,對每個點統(tǒng)計搜索窗口內(nèi)的相干點數(shù)目,從而得到各點的空間相干點密度.為簡化運算,將空間點密度量化為等間距的5級密度,點密度分布見圖3.由圖3可以發(fā)現(xiàn),磁懸浮沿線相干點密度遠(yuǎn)低于周邊房屋區(qū)域.獲取空間點密度后進(jìn)行分層抽樣,依據(jù)空間點密度將空間點分為5層,每層抽樣時采用等間距隨機選點,空間密度大的區(qū)域采用較大的采樣間隔,反之則采用較小的采樣間隔.在本次試驗中,5層抽樣的采樣間隔比為1∶2∶3∶4∶5.設(shè)置目標(biāo)抽樣點數(shù)目為15 000,可以得到5層的采樣間隔依次是14,28,42,56,70.進(jìn)而從472 929個初始相干點中抽樣得到17 187個相干點,這與降低閾值法得到的相干點數(shù)目相當(dāng),分層抽樣得到的相干點分布見圖2c.對比圖2b可以發(fā)現(xiàn),相干點在空間的分布更加均勻,并且較好地保留了磁懸浮沿線的相干點.

a 初始選取相干點

b 降低選點閾值法

c 基于點密度的分層抽樣法

圖2 參考網(wǎng)相干點分布

Fig.2 Distribution of coherent points in the reference network

圖3 空間點五級密度分布

求得參考網(wǎng)中相干點的參數(shù)后,擴(kuò)大DA閾值至0.4,選取中等質(zhì)量相干點,并利用附加網(wǎng)求解分層抽樣余下的相干點以及中等質(zhì)量的相干點,為保證高質(zhì)量相干點的精度,將附加網(wǎng)分為兩次進(jìn)行,分層抽樣余下的相干點構(gòu)成一級附加網(wǎng),中等質(zhì)量相干點構(gòu)成二級附加網(wǎng),各級網(wǎng)獲取的相干點數(shù)目見表2.

表2 各級網(wǎng)絡(luò)獲取的相干點數(shù)目

Tab.2 The number of accepted coherent points in each network

參考網(wǎng)一級附加網(wǎng)二級附加網(wǎng)初選點數(shù)17 155470 1392 210 867接受點數(shù)15 813435 2341 711 388通過率/%92.292.677.4

基于弧段解的時間序列分析方法與選點策略無關(guān),在數(shù)據(jù)處理的流程中通過網(wǎng)絡(luò)連通性分析來去除噪聲點,以保證最終結(jié)果的可靠性.初始選點質(zhì)量越高,經(jīng)網(wǎng)絡(luò)連通性分析后可接受的相干點越多.由表2可以發(fā)現(xiàn),參考網(wǎng)與一級附加網(wǎng)中的相干點是通過同一DA閾值選取的,因此有著相似的通過率,二級附加網(wǎng)的相干點由于選點閾值較為寬松,噪聲點多,因此可接受相干點的通過率低于前者.

Sentinel-1A與ENVISAT數(shù)據(jù)的處理使用類似的多級化處理方式.不同數(shù)據(jù)集得到的磁懸浮沿線及其周邊區(qū)域的視線向形變速率如圖4所示.為了盡可能突出沉降結(jié)果的細(xì)節(jié)信息,所有結(jié)果均保持SAR圖像的原始分辨率,未做多視處理.

不同衛(wèi)星數(shù)據(jù)在磁懸浮軌道及其周邊區(qū)域得到的相干點數(shù)目依次為:ENVISAT 129 182個,Sentinel-1A 356 663個,TerraSAR-X 2 165 435個.由圖4b和4c可知,Sentinel-1A與TerraSAR-X得到的結(jié)果中,磁懸浮沿線均存在大量呈線性排列的相干點,而在ENVISAT的結(jié)果中(圖4a),磁懸浮沿線相干點極少.Sentinel-1A與ENVISAT均為中等分辨率C波段數(shù)據(jù),且Sentinel-1A數(shù)據(jù)的方位向分辨率低于ENVISAT數(shù)據(jù).由此可以得知磁懸浮沿線上相干點的數(shù)目與距離向分辨率有關(guān),而與方位向分辨率關(guān)系較弱.ENVISAT數(shù)據(jù)由于距離向分辨率較低,在磁懸浮沿線難以形成有效的點目標(biāo),不能獲取磁懸浮軌道上足夠的形變信息.

對比表1中不同傳感器的測繪幅寬,盡管TerraSAR-X數(shù)據(jù)分辨率較高,但是單景數(shù)據(jù)覆蓋范圍僅為Sentinel-1A數(shù)據(jù)的1/30.ENVISAT數(shù)據(jù)與Sentinel-1A數(shù)據(jù)有著相似的空間分辨率,但是Sentinel-1A數(shù)據(jù)的覆蓋范圍是ENVISAT數(shù)據(jù)的5倍,且相干點的密度更高.因此Sentinel-1A數(shù)據(jù)的可用性更好.

此外,選取一塊形變顯著的區(qū)域(圖4中白色方框區(qū)域),詳細(xì)研究磁懸浮周邊區(qū)域沉降對軌道的影響.不同衛(wèi)星數(shù)據(jù)獲取的該區(qū)域沉降速率見圖5.圖5中表示沉降速率的色條與圖4一致,白線表示磁懸浮軌道位置.從圖5中可以發(fā)現(xiàn),在2007~2011年間,磁懸浮周邊區(qū)域沉降較大,最大年沉降速率大于20 mm·年-1(圖5a);2011~2015年間,該區(qū)域趨于穩(wěn)定,但是仍然存在沉降速率大于15 mm·年-1的點(圖5b);2016~2018年間,該區(qū)域大多數(shù)點保持穩(wěn)定,沉降速率小于10 mm·年-1(圖5c).不同數(shù)據(jù)得到的結(jié)果中沿軌道的線相干點密度有很大差異,

a ENVISAT

b TerraSAR-X

c Sentinel-1A

Fig.4 Subsidence velocity along the line of sight within the maglev zone

但是這些點在不同時段均保持穩(wěn)定,年平均沉降速率小于3 mm·年-1.因此周邊區(qū)域的沉降并未影響磁懸浮軌道的安全.

a ENVISAT

b TerraSAR-X

c Sentinel-1A

圖5 沉降較大區(qū)域放大圖

Fig.5 Magnified view of the selected area with large subsidence

3 結(jié)果驗證

利用實測水準(zhǔn)數(shù)據(jù)驗證本文多級MT-InSAR結(jié)果的可靠性.精密水準(zhǔn)包含2009~2014年間5期數(shù)據(jù),與TerraSAR-X數(shù)據(jù)獲取的時間段一致.該區(qū)域共26個水準(zhǔn)點(BM1~BM26),位于磁懸浮軌道以及周邊區(qū)域,具體位置見圖6.

圖6 磁懸浮區(qū)域水準(zhǔn)點分布

Fig.6 Distribution of levelling points within the maglev zone

水準(zhǔn)測量基于地面單個點,而SAR圖像中的單個像元反映了地面某一塊區(qū)域的綜合信息,二者不能直接對比.為了獲取MT-InSAR結(jié)果中與水準(zhǔn)點對應(yīng)的沉降信息,以每個水準(zhǔn)點為中心,選取30 m×30 m的搜索窗口,以搜索窗口內(nèi)所有相干點的加權(quán)平均值作為與水準(zhǔn)點對應(yīng)的結(jié)果.MT-InSAR結(jié)果表明磁懸浮軌道與周邊區(qū)域呈現(xiàn)不同的沉降信息,在使用搜索窗口時必須將軌道與周邊區(qū)域相干點分開處理.因此水準(zhǔn)數(shù)據(jù)的驗證分為兩組,磁懸浮軌道上點以及軌道周邊區(qū)域的點.

使用緩沖區(qū)法分離磁懸浮軌道上點,選取10 m緩沖區(qū)半徑,在磁懸浮沿線建立緩沖區(qū),位于緩沖區(qū)內(nèi)部的點視為磁懸浮軌道相干點.圖7表示磁懸浮軌道上水準(zhǔn)點與MT-InSAR相干點的形變速率對比結(jié)果,二者結(jié)果差異的平均值是-0.020 4 mm·年-1,標(biāo)準(zhǔn)差是0.628 3 mm·年-1.圖8表示磁懸浮周邊區(qū)域內(nèi)水準(zhǔn)點與MT-InSAR相干點的形變速率對比結(jié)果,二者差異的均值為-0.414 4 mm·年-1,標(biāo)準(zhǔn)差是1.170 8 mm·年-1.MT-InSAR的結(jié)果中,磁懸浮軌道上點與水準(zhǔn)測量結(jié)果符合度較高.而磁懸浮周邊區(qū)域中相干點難以與水準(zhǔn)點直接對應(yīng),因此符合度較低.

圖7 磁懸浮軌道上水準(zhǔn)點與MT-InSAR形變速率對比

Fig.7 Comparison between deformation velocities derived from levelling and from MT-InSAR along the maglev track

圖8 磁懸浮周邊區(qū)域水準(zhǔn)點與MT-InSAR形變速率對比

Fig.8 Comparison of deformation velocity between levelling and MT-InSAR nearby the maglev track

選取4個水準(zhǔn)點BM6、7、11、12,對比水準(zhǔn)與MT-InSAR獲取的時間序列累積形變量,如圖9所示.由圖可以發(fā)現(xiàn),MT-InSAR得到的形變時間序列與水準(zhǔn)結(jié)果有著良好的一致性.并且水準(zhǔn)結(jié)果表明BM11與BM6在2013年10月至2014年期間出現(xiàn)了明顯的抬升現(xiàn)象,這與MT-InSAR獲取的結(jié)果有著相似的形變趨勢.2013年開始出現(xiàn)的抬升現(xiàn)象與上海市實行的限制地下水開采以及地下水回灌有關(guān)[24].

a 水準(zhǔn)點11

b 水準(zhǔn)點6

c 水準(zhǔn)點12

d 水準(zhǔn)點7

圖9 水準(zhǔn)點與MT-InSAR時間序列形變量對比

Fig.9 Comparison between time series displacement derived from levelling and from MT-InSAR

4 結(jié)語

本文介紹了一種多級化的MT-InSAR數(shù)據(jù)處理方法,提出了一種基于相干點空間密度的抽樣方法,綜合考慮選取相干點的密度以及空間分布,得到優(yōu)化的分級MT-InSAR處理方法.并利用該方法獲取了2007~2018年間的上海磁懸浮及其周邊區(qū)域沉降信息,不同數(shù)據(jù)結(jié)果均顯示磁懸浮軌道保持穩(wěn)定,而周邊區(qū)域存在不同程度的沉降,隨著時間的累積,沉降速率是逐漸變小的.采用精密水準(zhǔn)數(shù)據(jù)對沉降監(jiān)測結(jié)果進(jìn)行了精度驗證,兩者在時間以及空間上吻合度較高.MT-InSAR技術(shù)得到的形變速率可以滿足高精度的地面形變監(jiān)測要求.

此外,本文還分析了不同SAR影像數(shù)據(jù)分辨率對于磁懸浮沿線相干點獲取的影響,發(fā)現(xiàn)相干點的密度與SAR數(shù)據(jù)的距離向分辨率有著密切關(guān)系,而方位向分辨率影響較小.對比上一代C波段中等分辨率ENVISAT數(shù)據(jù),開源的Sentinel-1A數(shù)據(jù)有著更大的覆蓋范圍,并且憑借較高的距離向分辨率,可以獲取更多的相干點,在城市地表以及大型結(jié)構(gòu)體形變檢測中有著突出的優(yōu)勢.