本發(fā)明涉及遺產(chǎn)地監(jiān)測領(lǐng)域，特別是涉及基于分布式散射體時序干涉sar技術(shù)的遺產(chǎn)地形變監(jiān)測方法。
背景技術(shù)：
遺產(chǎn)地，如文物古跡、自然奇觀，是人類歷史上具有重要意義和巨大價值的無以替代的財富。5000余年的歷史文明以及遼闊的地域使得中國成為遺產(chǎn)大國，截止2015年7月列入《世界遺產(chǎn)名錄》的有48項，數(shù)量位居世界第二，僅次于意大利。不幸的是，遺產(chǎn)地在自然災害和人類活動的侵蝕下正在變得日益不穩(wěn)定。自然災害包括滑坡、地震、洪水、惡劣的天氣以及氣候突變。人類活動包括戰(zhàn)爭、資源的過度開采、城鎮(zhèn)化建設(shè)和失控的旅游業(yè)。這些外界因素使得很多遺產(chǎn)地岌岌可危，監(jiān)測并實施有效的保護勢在必行。
遺產(chǎn)地不穩(wěn)定的程度可以用形變來衡量，大致可以分為兩類，一是長期的緩慢形變，尺度在毫米級或厘米級，主要由自然的侵蝕、地殼運動、人類活動引起的地表沉陷等引起；二是突發(fā)的劇烈形變，尺度在厘米級以上，主要由地震、滑坡等自然災害以及人類破壞等引起。傳統(tǒng)的形變測量方法如水準、gps或物探技術(shù)均存在著以下不足：1)效率低、成本高、需要耗費大量的人力物力，受野外作業(yè)環(huán)境的嚴重限制；2)在空間維度上無法獲得大面積形變場，時間維度上難以獲取長時間序列形變，時空變化規(guī)律難以展現(xiàn)。
現(xiàn)有的時間序列差分干涉sar技術(shù)(multi-temporaldifferentialinterferometricsyntheticapertureradar,mtinsar)聯(lián)合多基線重軌sar影像測量地表微小形變，克服了傳統(tǒng)方法的局限性，使得快速、大面積、時序監(jiān)測成為了可能，而且理論上精度可以達到毫米級，因此更加適合于遺產(chǎn)地的監(jiān)測。常用的mtinsar技術(shù)包括：永久散射體技術(shù)(persistentscatterers,psinsar),小基線集技術(shù)(smallbaselinessubset,sbas)，stamps以及squeesar技術(shù)。
但是現(xiàn)有的技術(shù)在遺產(chǎn)地監(jiān)測中仍存在一些亟待解決的問題：1)對sar影像數(shù)量的要求較高(通常要求為25景以上，如psinsar和squeesar)；2)郊區(qū)干涉對的時間失相干現(xiàn)象嚴重，永久散射體(persistentscatterers,ps)稀少，可用于參數(shù)反演的點數(shù)量不足，影響相位解纏和形變估計的精度。sbas方法雖然利用短時間/空間基線的干涉對來減少失相干的影響，以挖掘短時間內(nèi)能夠保持相干的像元信息來增加監(jiān)測點的密度，但是多視操作降低了空間分辨率，相干性估計也容易出現(xiàn)偏差。
mtinsar技術(shù)目前在郊區(qū)的研究主要集中在分布式散射體(distributedscatterers，ds)上，它指的是占據(jù)一定面積均勻地物，在sar影像上對應的像元呈現(xiàn)出相似的后向散射特征，裸地和低矮的植被都屬于ds點。引入ds點能夠顯著提升遺產(chǎn)地可監(jiān)測目標的數(shù)量，改進形變反演的結(jié)果。但是分布式散射體相干性差、信噪比低，需要進行復雜特殊的濾波處理。
由此可見，現(xiàn)有的mtinsar技術(shù)在遺產(chǎn)地形變監(jiān)測中仍然存在著缺陷，如何能創(chuàng)設(shè)一種更加適用的基于分布式散射體時序干涉sar技術(shù)的遺產(chǎn)地形變監(jiān)測方法，成為業(yè)界追求的目標。
技術(shù)實現(xiàn)要素：
本發(fā)明要解決的技術(shù)問題是提供一種基于分布式散射體時序干涉sar技術(shù)的遺產(chǎn)地形變監(jiān)測方法，使其實現(xiàn)大范圍、長時間、高精度的監(jiān)測遺產(chǎn)地不穩(wěn)定現(xiàn)象，從而克服現(xiàn)有的遺產(chǎn)地形變監(jiān)測方法的不足。
為解決上述技術(shù)問題，本發(fā)明提供1、基于分布式散射體時序干涉sar技術(shù)的遺產(chǎn)地形變監(jiān)測方法，包括如下步驟：
(一)確定遺產(chǎn)地的監(jiān)測區(qū)域及資料收集
(11)確定遺產(chǎn)地的監(jiān)測范圍，包括遺產(chǎn)地和所述遺產(chǎn)地周邊的緩沖帶；
(12)調(diào)查所述遺產(chǎn)地的背景知識，并收集所述遺產(chǎn)地的時間序列sar單視復數(shù)影像；
(二)sar數(shù)據(jù)處理與參數(shù)反演
(21)對所述步驟(12)收集的時間序列sar單視復數(shù)影像進行互配準、輻射校正和影像裁剪預處理；
(22)提取永久散射體目標和分布式散射體目標；
(23)去除分布式散射體目標的地形與平地相位，采用最大似然估計和lbfgs技術(shù)進行優(yōu)化相位估計；
(24)融合初始相位的永久散射體目標和經(jīng)過優(yōu)化相位的分布式散射體目標，設(shè)置短時間/空間基線閾值生成小基線集干涉對，對所述小基線集干涉對進行相位解纏并去除地形、大氣延遲造成的誤差；
(25)求解平均年形變速率及形變歷史變化值。
作為本發(fā)明的一種改進，所述方法還包括結(jié)果分析與解譯步驟，具體為：利用野外實地調(diào)查測量結(jié)果對步驟(25)所得結(jié)果進行精度評定；并結(jié)合野外實地調(diào)查測量結(jié)果與步驟(25)所得結(jié)果共同分析所述遺產(chǎn)地的形變特征及引發(fā)不穩(wěn)定的因素。
進一步改進，所述步驟(12)中所述遺產(chǎn)地的背景知識包括所處的環(huán)境、遭受的威脅、陸地景觀、地表特征、空間尺度、遺產(chǎn)單體的幾何形狀和材料屬性。
進一步改進，所述步驟(22)中提取永久散射體目標的具體方法為：根據(jù)單幅影像的光譜特征與多幅影像的振幅離差閾值法提取。
進一步改進，所述步驟(22)中提取分布式散射體目標的具體方法包括：
a、對所述時間序列sar單視復數(shù)影像中像元的時序幅度值進行排序，轉(zhuǎn)換得到像元幅度值累積分布函數(shù)，其像元幅度值累積分布函數(shù)的無偏估計表示為:
其中，fn(x)為像元幅度值累積分布函數(shù)的無偏估計，x代表像元幅度值，n為配準后的sar影像數(shù)量，x1為時序幅度值排序中的第1個元素的幅度值，xk為時序幅度值排序中的第k個元素的幅度值；xk+1為時序幅度值排序中的第k+1個元素的幅度值；k取值范圍為0…n-1。
b、再用ks檢驗或ad檢驗方法對相鄰像元的無偏估計值進行判定，滿足一定閾值的即為同質(zhì)像元；
c、對于每個像元可獲得其估計窗口內(nèi)連通的所有同質(zhì)像元，同質(zhì)數(shù)量大于設(shè)定閾值的像元即為分布式散射體目標。
進一步改進，所述ad檢驗方法定義為：
其中，p,q代表兩個像元，
進一步改進，所述步驟(24)中對所述小基線集干涉對的相位解纏方法采用時空3d相位解纏方法。
進一步改進，針對第一代sar數(shù)據(jù)進行軌道、地形相關(guān)誤差去除的模型如下：
其中，x和y為距離向和方位向的坐標，h為高程；ε為隨機相位誤差；ai為待估計參數(shù)。
進一步改進，所述步驟(25)中平均年形變速率是針對差分干涉對的相位模型采用最小二乘法求得；所述形變歷史變化值是利用svd算法求得，具體算法為：
假設(shè)第i幅干涉圖，主輔影像獲取時間分別為ta和tb，且tb＞ta，像元x處的差分干涉相位表示為：
式中φ(ta,x)和φ(ta,x)為在ta和tb時刻影像上的形變相位值，d(ta,x)和d(tb,x)是相對于參考時間t0的雷達視線向形變量，參考時間形變量d(t0,x)≡0，λ為雷達波長，用φ表示n個時刻sar圖像上的形變相位矩陣，δφ表示m個差分干涉圖上的相位矩陣，其矩陣形式表示為：
δφ＝aφ
其中，系數(shù)矩陣a[m×n]的每行對應一幅干涉圖，每列對應一個時間上的sar影像，用svd算法求解得每個時間sar影像上形變對應的相位，經(jīng)轉(zhuǎn)換可得每個時間上的形變值。
采用這樣的設(shè)計后，本發(fā)明至少具有以下優(yōu)點：
本發(fā)明基于分布式散射體的時序干涉sar技術(shù)用于遺產(chǎn)地形變監(jiān)測，使其能利用較少的sar影像數(shù)據(jù)實現(xiàn)更大區(qū)域、更多地物的形變測量，從而推進遺產(chǎn)地的穩(wěn)定監(jiān)測與保護。
本發(fā)明屬于一種適用于大范圍、長時間、高精度監(jiān)測遺產(chǎn)地不穩(wěn)定現(xiàn)象的方法，特別適合于建筑物稀少的遺郊區(qū)。
附圖說明
上述僅是本發(fā)明技術(shù)方案的概述，為了能夠更清楚了解本發(fā)明的技術(shù)手段，以下結(jié)合附圖與具體實施方式對本發(fā)明作進一步的詳細說明。
圖1是本發(fā)明基于分布式散射體時序干涉sar技術(shù)的遺產(chǎn)地形變監(jiān)測方法示意圖。
具體實施方式
基于遺產(chǎn)地的監(jiān)測需求和地表類型，本發(fā)明融合了sbas和squeesar算法，識別、提取ds點和ps點并進行聯(lián)合分析，采用短時間/空間基線生成干涉對集合，最后去除地形、大氣等誤差生成形變結(jié)果，具體方法如下。
參照附圖1所示，本發(fā)明基于分布式散射體時序干涉sar技術(shù)的遺產(chǎn)地形變監(jiān)測方法，包括如下步驟：
一、確定遺產(chǎn)地的監(jiān)測區(qū)域及資料收集
1.1、確定監(jiān)測范圍，應包括遺產(chǎn)地以及遺產(chǎn)地周邊一定范圍的緩沖帶，該緩沖帶有助于了解遺產(chǎn)地所處的環(huán)境，解釋遺產(chǎn)地發(fā)生不穩(wěn)定的原因和分析其發(fā)展規(guī)律。
1.2、調(diào)查確定遺產(chǎn)地的背景先驗知識，包括其所處的環(huán)境、遭受的威脅(如滑坡等自然災害或地面沉降等)，陸地景觀、地表特征、空間尺度、遺產(chǎn)單體的幾何形狀和材料屬性等數(shù)據(jù)。這些數(shù)據(jù)決定了要采取的監(jiān)測方法和數(shù)據(jù)處理流程。
1.3、收集與該遺產(chǎn)地相關(guān)的時間序列sar單視復數(shù)影像；
二、sar數(shù)據(jù)處理與參數(shù)反演
2.1、對收集的時間序列sar單視復數(shù)影像進行數(shù)據(jù)預處理，包括輻射校正、互配準、影像裁剪和主影像選擇等。
2.2、提取永久散射體目標和分布式散射體目標；
根據(jù)sar影像的后向散射特征和研究區(qū)域的地物類型提取永久散射體目標(ps點)和分布式散射體目標(ds點)。
由于ps點在sar影像上表現(xiàn)為十分穩(wěn)定的強散射，能夠保持高相干性，對應的地物為遺產(chǎn)區(qū)的房屋、橋梁、道路以及裸露的巖石等。ps點具體的提取方法為：根據(jù)單幅影像的光譜特征與多幅影像的振幅離差閾值法提取。
ds點對應著分辨單元內(nèi)所有散射體的后向散射系數(shù)大致相同的地物，多為郊區(qū)普遍存在的裸地和稀疏植被等，它們僅在部分干涉圖保持一定的相干特性。
ds點的提取則需要對sar影像上像元的后向散射特征進行統(tǒng)計分析，具體方法為：
1)對像元的時序幅度值進行排序，轉(zhuǎn)換得到其累積分布函數(shù)；
2)用kolmogorov-smirnov(ks)檢驗或anderson-darling(ad)檢驗方法對相鄰像元的無偏估計值進行判定，滿足一定閾值的即為同質(zhì)像元；
3)對于每個像元可獲得其估計窗口內(nèi)連通的所有同質(zhì)像元，同質(zhì)數(shù)量大于設(shè)定閾值的像元即為ds點。
其中，像元幅度值累積分布函數(shù)的無偏估計可表示為:
其中，fn(x)為像元幅度值累積分布函數(shù)的無偏估計，x代表像元幅度值，n為配準后的sar影像數(shù)量，x1為時序幅度值排序中的第1個元素的幅度值，xk為時序幅度值排序中的第k個元素的幅度值；xk+1為時序幅度值排序中的第k+1個元素的幅度值；k取值范圍為0…n-1。
無參數(shù)雙尾ad檢驗方法可定義為：
p,q代表兩個像元，
2.3、去除分布式像元的地形與平地相位，然后進行時空三維的優(yōu)化濾波；
因為ds點相干性差、信噪比低，需要進行時空濾波。首先像元p基于同質(zhì)像元集為ω的相干系數(shù)矩陣γ(p)可改進為：
相干系數(shù)矩陣第k行、j列的元素γkj表示為第k景和第j景sar影像組成干涉對的相干系數(shù)。主對角線為自身影像干涉且值為1。
基于相干系數(shù)矩陣和最大似然估計準則，假設(shè)第一幅sar影像相位
其中，λ＝exp(iφ)為n維向量，
2.4、融合永久散射體目標和分布式目標；
對ps點初始相位和ds點的優(yōu)化相位進行差分干涉對重新組合。
2.5、設(shè)置短時間/空間基線閾值生成小基線集干涉對；
好的相干性是差分干涉技術(shù)測量形變的前提，而引起干涉對失相干最嚴重的因素是時間基線和空間基線。以單一主影像形成的干涉對集合因為時間/空間基線大、而且對影像數(shù)據(jù)量要求多，難以滿足遺產(chǎn)地監(jiān)測需求。以短時間/空間基線為標準形成的多基線干涉對能夠克服失相干嚴重的問題，且能生成大量的干涉對降低對數(shù)據(jù)量的需求。
假設(shè)有n景影像，則小基線干涉對數(shù)量m滿足：
差分干涉對的相位模型可表示為：
其中，a為常數(shù)項，ρξ、ρη是大氣效應和基線引起的干涉相位分別沿方位向ξ和距離向η的線性變化系數(shù)，δh為高程誤差，δt為輔影像到主影像的時間基線，v表示線性形變速率，ε為殘差，包含了大氣殘余、非線性形變等。β⊥為干涉對的垂直基線，λ為雷達波長，r為衛(wèi)星傳感器到目標的斜距，θ為主圖像的入射角。其中
2.6、相位解纏并去除地形、大氣延遲等誤差；
纏繞干涉圖的相位值范圍是(-π,π)，解決2π整數(shù)倍問題的過程稱為相位解纏。這個過程是至關(guān)重要的，點密度過稀、地形起伏和形變較大等都能引起解纏錯誤?；谏鲜霾罘指缮鎸Φ南辔荒Ｐ?，對于大量的小基線集干涉對可采用時空3d相位解纏方法，以對干涉對進行相位解纏并去除地形、大氣延遲等誤差。
第一代衛(wèi)星影像的形變相位容易與軌道誤差、大氣延遲以及殘余地形相位等混淆；長波段衛(wèi)星的軌道誤差能夠達到30cm，在影像上呈現(xiàn)為二階模型的變化趨勢，同時軌道誤差能引起地形估計的偏差，大氣的湍流效應也與地形線性相關(guān)。
針對第一代sar數(shù)據(jù)進行軌道、地形相關(guān)誤差去除的模型如下：
其中，x和y為距離向和方位向的坐標，h為高程；ε為隨機相位誤差；ai為待估計參數(shù)。
2.7、求解平均年形變速率及形變歷史變化值；
對差分干涉對的相位模型用最小二乘法求解可得到平均年形變速率。采用奇異值去分解(singularvaluedecompositionsvd)的方法解算去過各種誤差的解纏相位可獲得形變的歷史變化值。假設(shè)第i幅干涉圖，主輔影像獲取時間分別為ta和tb(tb＞ta)，像元x處的差分干涉相位可表示為：
式中φ(ta,x)和φ(ta,x)為在ta和tb時刻影像上的形變相位值，d(ta,x)和d(tb,x)是相對于參考時間t0的雷達視線向(los)形變量，參考時間形變量d(t0,x)≡0，λ為雷達波長。用φ表示n個時刻sar圖像上的形變相位矩陣，δφ表示m個差分干涉圖上的相位矩陣，則矩陣形式可表示為：
δφ＝aφ
其中，系數(shù)矩陣a[m×n]的每行對應一幅干涉圖，每列對應一個時間上的sar影像。用svd算法求解可得每個時間sar影像上形變對應的相位，然后可轉(zhuǎn)換為每個時間上的形變值。
上述sar數(shù)據(jù)處理融合了sbas和squeesar算法的優(yōu)勢，一方面聯(lián)合提取與分析ds點和ps點，顯著提高監(jiān)測點的數(shù)量，擴大監(jiān)測的范圍，改善相位解纏與大氣延遲去除的效果；另一方面采用小基線干涉對集合進行形變估計，在保留原始分辨率的同時降低對sar影像數(shù)量的需求。
三、對干涉結(jié)果進行精度評定、解譯以及分析遺產(chǎn)地的不穩(wěn)定機制
3.1、野外調(diào)查測量與干涉結(jié)果精度評定；
野外調(diào)查測量是必不可少的工作，一方面是要確定監(jiān)測點對應的地物，如一些樹林中的廢墟、遺址等，在普通光學影像上無法看見但能夠被高分辨率sar監(jiān)測到的地物。另一方面是對形變結(jié)果進行實地驗證，如發(fā)現(xiàn)古建筑物出現(xiàn)裂縫，地表明顯塌陷等現(xiàn)象。條件允許，可獲取實地測量數(shù)據(jù)(如水準，gps等)對干涉形變結(jié)果進行絕對精度評定。
3.2、分析遺產(chǎn)地的形變特征及引發(fā)不穩(wěn)定的因素與機制。
分析該遺產(chǎn)地的形變特征及不穩(wěn)定的因素，需要結(jié)合遺產(chǎn)地的背景環(huán)境資料，如地質(zhì)災害、地下水開發(fā)或城市建設(shè)等進行，最終得出該遺產(chǎn)地的形變特征和不穩(wěn)定因素。
以上所述，僅是本發(fā)明的較佳實施例而已，并非對本發(fā)明作任何形式上的限制，本領(lǐng)域技術(shù)人員利用上述揭示的技術(shù)內(nèi)容做出些許簡單修改、等同變化或修飾，均落在本發(fā)明的保護范圍內(nèi)。