基于SBAS-InSAR技術(shù)的大壩形變探測與監(jiān)測

徐曉新 毛洪孝 蘭麗景

（1.衢州市國土空間規(guī)劃設(shè)計(jì)研究院柯城分院，浙江 衢州 324000；2.浙江振邦地理信息科技有限公司，浙江 衢州 324000）

0 引言

大壩的安全運(yùn)營［1］會影響大壩本身機(jī)構(gòu)與水輪發(fā)電裝置，還會影響下游居民的生命財(cái)產(chǎn)安全［2-3］。大壩投入使用后，受自身與外界環(huán)境的影響，會出現(xiàn)一定程度的形變［4-5］。微小形變不會影響大壩的運(yùn)營狀態(tài)，當(dāng)形變較大情況下，便會影響大壩的安全運(yùn)營［6］。為此，實(shí)時監(jiān)測大壩形變情況意義重大。

例如，楊傳訓(xùn)等［7］根據(jù)黃龍帶水庫大壩兩期高精度三維激光大壩掃描點(diǎn)云數(shù)據(jù)，提出基于正態(tài)分布的點(diǎn)云數(shù)據(jù)配準(zhǔn)算法、改進(jìn)原始漸進(jìn)加密三角網(wǎng)濾波算法，運(yùn)用點(diǎn)到面形變量對比分析模式監(jiān)測大壩形變量，完成大壩形變監(jiān)測。但該方法需要高質(zhì)量的大壩歷史數(shù)據(jù)，當(dāng)大壩歷史數(shù)據(jù)質(zhì)量較低時，則該方法的形變監(jiān)測效果較差。何子鑫［8］基于SBAS-InSAR 技術(shù)的壩體形變獲取方法，根據(jù)獲取的壩體形變信息，分析了大壩不同位置的形變特征及壩體形變與庫區(qū)水位變化的周期性規(guī)律，縮短監(jiān)測時間。該方法雖然能夠有效獲取大壩的形變信息，但易受時空與大氣效應(yīng)的影響，降低大壩形變監(jiān)測精度。

多基線干涉合成孔徑雷達(dá)（small baseline subset-interferometric synthetic aperture radar，SBAS-InSAR）技術(shù)結(jié)合了InSAR 技術(shù)全天候、全天時與高精度等優(yōu)勢，以及SBAS 不受時空與大氣效應(yīng)影響的優(yōu)勢，可提升形變監(jiān)測效果［9-12］，同時該技術(shù)并不需要高質(zhì)量的大壩影像數(shù)據(jù)，便可高效精準(zhǔn)地完成大壩形變監(jiān)測。為此，研究基于SBAS-InSAR 技術(shù)的大壩形變探測與監(jiān)測方法，實(shí)時掌握大壩的運(yùn)營狀態(tài)。

1 大壩形變探測與監(jiān)測

1.1 大壩差分干涉相位圖的濾波處理

經(jīng)過組建的M幅大壩干涉相位圖內(nèi)，包含大量噪聲，不利于后續(xù)大壩形變的探測與監(jiān)測。為此，利用改進(jìn)的Goldstein 濾波算法，濾波處理大壩差分干涉圖，去除干涉圖內(nèi)部噪聲［13］，具體步驟如下。

（1）令大壩差分干涉圖是I(x，y)，其中，大壩差分干涉圖像素位置是(x，y)。通過二維傅里葉變換處理I(x，y)，將空間域信號變更成頻率域［14-15］，得到大壩的干涉小區(qū)域I′(u，v)。

（2）卷積運(yùn)算I(x，y)的矩形平滑窗口與局部功率譜Q(u，v)，并實(shí)施非線性變換處理［16-17］，獲取濾波器Z(u，v)。

（3）傅里葉逆變換，獲取濾波后的大壩差分干涉圖，公式為

式中，s是濾波因子。

大壩差分干涉圖內(nèi)相干值的絕對值|γ|的函數(shù)為

式中，大壩差分干涉圖的相干系數(shù)是γ；歐拉算法是Li2()；矩形濾波窗口均值是δφ。

1.2 大壩差分干涉相位圖的相位解纏

在區(qū)域增長法內(nèi)，引入最小不連續(xù)閾值，對其進(jìn)行改進(jìn)，可提升大壩差分干涉相位圖的相位解纏效果［18］。因此，利用最小不連續(xù)的區(qū)域增長解纏算法，對1.1 小節(jié)濾波后的大壩差分干涉相位圖進(jìn)行相位解纏。令大壩濾波后大壩差分干涉相位圖內(nèi)隨機(jī)像元（x，y）處的包裹相位是θx，y，展開相位是ψx，y。ψx，y的計(jì)算公式為

式中，cx，y為包裹數(shù)；ρ為修正系數(shù)。

（x，y）處的垂直與水平方向的跳躍數(shù)αx，y、βx，y為

式中，Int(·)代表最接近的整數(shù)。

在式（4）內(nèi)添加式（3）獲取公式為

最小不連續(xù)算法的目標(biāo)是以變更跳躍數(shù)的方式，令η降至最低，獲取最小不連續(xù)閾值。

高相干點(diǎn)相位解纏的具體步驟如下：

（1）在濾波后的大壩差分干涉相位圖內(nèi)，選取質(zhì)量較優(yōu)的種子，以其為起點(diǎn)［19-20］。

（2）求解已展開的像素點(diǎn)附近的最小不連續(xù)閾值，通過式（3）對最小閾值對應(yīng)的像素點(diǎn)實(shí)施相位展開。

（3）通過多個方向的信息，展開各像素點(diǎn)，避免因單個方向信息展開像素點(diǎn)出現(xiàn)偶然誤差情況。

（4）通過最小不連續(xù)閾值，確定像素點(diǎn)相位展開過程，確保每次展開的像素點(diǎn)均是最小不連續(xù)閾值點(diǎn)。

（5）在區(qū)域生長重合情況下，需以路徑連接的方式，融合重合的生長區(qū)域，確保像素點(diǎn)解纏相位的連續(xù)性，完成大壩差分干涉相位圖的相位解纏。

1.3 大壩差分干涉圖的高相關(guān)點(diǎn)選取

對于1.2 小節(jié)相位解纏的大壩相位解纏圖，通過相干系數(shù)閾值法，可有效剔除時間去相干導(dǎo)致的隨機(jī)散射相位問題。在大壩相位解纏圖內(nèi)，選擇高相干點(diǎn)是指自身的散射特性波動性較小，同時不會根據(jù)時間的變化而變化。

通過相干性衡量大壩相位解纏圖內(nèi)的主輔影像間相關(guān)程度的物理量。大壩相位解纏圖的相干系數(shù)計(jì)算公式為

其中，大壩相位解纏圖的空間矩形窗大小是(2m+1)(2n+1)；代表取共軛（·）*；大壩相位解纏圖內(nèi)主圖像數(shù)據(jù)是g1；輔圖像數(shù)據(jù)是g2。

將全部符合式（8）的像元定義為被選擇的高相干點(diǎn)，公式為

式中，εa、εb是相干系數(shù)閾值；γ的最小值與均值為γmin、γmean。

1.4 大壩形變探測與監(jiān)測的高程測量估計(jì)數(shù)學(xué)模型構(gòu)建

依據(jù)1.3 小節(jié)選取的高相干點(diǎn)，構(gòu)建大壩形變探測與監(jiān)測的高程測量估計(jì)數(shù)學(xué)模型。令大壩相位解纏圖內(nèi)高相干點(diǎn)o的干涉相位是ξ(o)，計(jì)算公式為

大壩形變探測與監(jiān)測的高程測量估計(jì)數(shù)學(xué)模型為

式中，系數(shù)矩陣是D(M×N)，D內(nèi)各行數(shù)值為每個大壩解纏相位圖，各列數(shù)值為原始大壩影像；N幅大壩影像內(nèi)待求形變相位建立的矩陣是ξ；M幅大壩相位解纏圖內(nèi)；高相干點(diǎn)干涉相位組建的矩陣是Δξ。

因?yàn)镸≥N，所以D的秩是N，利用最小二乘法求解式（11）得到

式中，D的廣義逆矩陣是D*；轉(zhuǎn)置符號是T；（·）-1代表求逆運(yùn)算。

利用奇異值分解法，求解式（10），公式為

式中，DTD內(nèi)M×M的正交矩陣是U；DTD內(nèi)M×N的正交矩陣是R；DTD內(nèi)對角線元素的對角矩陣是C。令D的秩是H，則DTD的前H個特征值不是0，其余特征值均是0。因此式（12）可變更為

式中，C內(nèi)第l個元素是?l；U與R內(nèi)第l個元素是ul、rl。

式中，G(j′，B)是0 的矩陣是G。通過奇異值分解法計(jì)算G，可獲取υ的最小范數(shù)解，在求解大壩形變探測與監(jiān)測的高程測量估計(jì)數(shù)學(xué)模型時，需考慮DEM誤差，Δ?的求解結(jié)果為

2 實(shí)驗(yàn)分析

以某地區(qū)的水電站大壩為實(shí)驗(yàn)對象，該水電站大壩屬于一級建筑物，裝機(jī)容量是430萬千瓦，發(fā)電量接近186.5萬千瓦。該水電站大壩的主要功能是發(fā)電，還具備防洪與灌溉等功能。該大壩的拱壩壩頂高程約1 255 m，最大壩高約295 m，拱冠梁寬約15 m，最低建基基面高程約955 m。該大壩的水庫容量在1.55×1010m3左右，調(diào)節(jié)庫容量在1.02×1010m3左右。利用Sentinel-1A 衛(wèi)星搭載C 波段合成孔徑雷達(dá)，采集該大壩的影像數(shù)據(jù)，該大壩的部分影像數(shù)據(jù)如圖1所示。

圖1 大壩的部分影像數(shù)據(jù)

根據(jù)圖1 可知，C 波段合成孔徑雷達(dá)采集的大壩影像數(shù)據(jù)清晰度較高。

設(shè)置時間基線閾值是100 d，垂直有效基線閾值是100 m，依據(jù)設(shè)置的時間基線閾值，在采集的大壩影像數(shù)據(jù)內(nèi)，選擇符合要求的大壩影像，符合要求大壩影像數(shù)據(jù)的時空基線分布情況如圖2所示。

圖2 時空基線分布情況

根據(jù)圖2 可知，本文方法可有效按照設(shè)置的時間基線閾值，選擇符合要求的大壩影像數(shù)據(jù)，選擇的大壩影像數(shù)據(jù)內(nèi)，時間基線控制在100 d，符合設(shè)置的時間基線閾值要求，垂直有效基線均在100 m 以內(nèi)，符合垂直有效基線閾值要求。按照選擇的大壩影像數(shù)據(jù)時間基線點(diǎn)，對其進(jìn)行兩兩相連，組成大壩差分干涉相位圖對，大壩差分干涉相位圖對組成情況如圖3所示。

圖3 大壩差分干涉相位圖對組成情況

由圖3 可知，依據(jù)符合要求的大壩影像數(shù)據(jù)時間基線，共組成27 幅大壩差分干涉相位圖對。實(shí)驗(yàn)證明：本文方法可有效生成大壩差分干涉相位圖對。

在27幅大壩差分干涉相位圖對，隨機(jī)選擇一幅大壩差分干涉相位圖，利用本文方法對其進(jìn)行濾波處理，本文方法的濾波效果如圖4所示。

圖4 濾波前后的大壩差分干涉相位圖

根據(jù)圖4（a）可知，濾波器前大壩差分干涉相位圖內(nèi)，包含大量噪聲點(diǎn)，無法為后續(xù)大壩形變探測與監(jiān)測，提供精準(zhǔn)的數(shù)據(jù)支持；根據(jù)圖4（b）可知，濾波后的大壩差分干涉相位圖內(nèi)，噪聲點(diǎn)明顯減少，可清晰呈現(xiàn)大壩的細(xì)節(jié)信息。實(shí)驗(yàn)證明：本文方法可有效濾波處理大壩差分干涉相位圖，去除差分干涉相位圖的內(nèi)部噪聲，提升其清晰度。

利用本文方法對濾波后的大壩差分干涉相位圖進(jìn)行相位解纏，相位解纏結(jié)果如圖5所示。

圖5 大壩差分干涉相位圖的相位解纏結(jié)果

根據(jù)圖5（a）可知，改進(jìn)前的大壩差分干涉相位圖的相位解纏結(jié)果中，包含大量大氣誤差；根據(jù)圖5（b）可知，本文方法改進(jìn)后的大壩差分干涉相位圖的相位解纏結(jié)果中，僅包含少量的大氣誤差，可明顯提升相位解纏效果。實(shí)驗(yàn)證明：本文方法可有效相位解纏大壩差分干涉相位圖，且相位解纏效果較優(yōu)，僅存在少量的大氣誤差。

利用本文方法在大壩相位解纏圖內(nèi)，選擇高相關(guān)點(diǎn)，用于建立大壩形變探測與監(jiān)測的高程測量估計(jì)數(shù)學(xué)模型，高相干點(diǎn)選擇結(jié)果如圖6所示。

圖6 高相干點(diǎn)選擇結(jié)果

由圖6 可知，本文方法可有效在大壩相位解纏圖內(nèi)，選擇高相干點(diǎn)，數(shù)量為56 個。本文方法選擇高相關(guān)點(diǎn)，主要集中在河流附近，高相關(guān)點(diǎn)的分布密度較高；其他地方的高相干點(diǎn)較少，分布密度較小。

在該大壩的左邊與右邊各選擇一個監(jiān)測點(diǎn)，利用本文方法在選擇的高相關(guān)點(diǎn)上，建立大壩形變探測與監(jiān)測的高程測量估計(jì)數(shù)學(xué)模型，分析這兩個監(jiān)測點(diǎn)的大壩形變情況，結(jié)果如圖7所示。

圖7 不同水位時的大壩形變探測與監(jiān)測結(jié)果

由圖7（a）與圖7（b）可知，不同水庫水位時，隨著時間的延長，大壩左右兩邊的形變量均呈上升趨勢。根據(jù)圖7（a）可知，當(dāng)水庫水位較高時，大壩左邊的形變量較小，大壩右邊的形變量略大于左邊。根據(jù)圖7（b）可知當(dāng)水庫水位較低時，大壩整體形變趨勢與水庫水位較高時的變化趨勢非常接近，但整體形變量明顯低于水位較高時的形變量；原因是大壩長時間受到上游較高水位的壓力，導(dǎo)致大壩出現(xiàn)向下游的位移，即大壩形變量較大。實(shí)驗(yàn)證明：本文方法可有效探測與監(jiān)測大壩形變。

3 結(jié)束語

傳統(tǒng)大壩形變監(jiān)測方式是操作人員通過測距儀與經(jīng)緯儀，觀測并記錄各監(jiān)測點(diǎn)的位置，完成形變監(jiān)測，該方式需要操作人員具有豐富的工作經(jīng)驗(yàn)，且監(jiān)測成本較高，無法大面積監(jiān)測大壩形變。為此，利用SBAS-InSAR 技術(shù)的大區(qū)域、高精度等監(jiān)測優(yōu)勢，研究基于SBAS-InSAR 技術(shù)的大壩形變探測與監(jiān)測方法，精準(zhǔn)探測與監(jiān)測大壩形變量，為維護(hù)大壩安全運(yùn)營提供參考。