CN103698749B - A kind of method utilizing small data set SAR image sequential extraction procedures Permanent scatterers - Google Patents

Abstract

Translated fromChinese

本发明提供一种利用小数据集SAR图像序列提取永久散射体的方法。技术方案是利用包括8至15幅SAR图像的序列，首先进行SAR图像序列预处理，然后进行空域信杂比估计和时域信杂比估计，最后通过PS判定，得到PS提取结果。本发明的有益效果是：实现简单、准确度高，适用范围广。

The invention provides a method for extracting permanent scatterers by using small data set SAR image sequences. The technical solution is to use a sequence including 8 to 15 SAR images, first to perform SAR image sequence preprocessing, then to estimate the signal-to-clutter ratio in the space domain and the signal-to-clutter ratio in the time domain, and finally to obtain the PS extraction result through PS judgment. The invention has the beneficial effects of simple implementation, high accuracy and wide application range.

Description

Translated fromChinese

一种利用小数据集SAR图像序列提取永久散射体的方法A method for extracting permanent scatterers from SAR image sequences with small data sets

技术领域technical field

本发明涉及微波遥感和信号处理的交叉技术领域，特别涉及一种利用小数据集SAR（Synthetic Aperture Radar，合成孔径雷达）图像提取永久散射体的方法。The invention relates to the cross technical field of microwave remote sensing and signal processing, in particular to a method for extracting permanent scatterers by using a small data set SAR (Synthetic Aperture Radar, Synthetic Aperture Radar) image.

背景技术Background technique

PS（Permanent Scatterers，永久散射体）技术是一种典型的长时间序列差分干涉技术，通过处理星载SAR系统在不同时刻对同一区域获取的若干幅SAR图像，可获得毫米级的地表形变信息，因此被广泛用于自然灾害监测、地表形变提取、城市沉降防治等，具有重要的应用价值。PS (Permanent Scatterers, permanent scatterers) technology is a typical long-term sequence differential interferometry technology, which can obtain millimeter-level surface deformation information by processing several SAR images of the same area acquired by the spaceborne SAR system at different times. Therefore, it is widely used in natural disaster monitoring, surface deformation extraction, urban settlement prevention, etc., and has important application value.

PS技术处理流程包括PS提取、PS网络构建、PS网络参数求解和优化、地表形变场提取等步骤，其中，PS提取的准确度影响后续各步骤的处理精度和稳健度，最终影响地表形变场的提取精度。The PS technology processing flow includes PS extraction, PS network construction, PS network parameter solution and optimization, and surface deformation field extraction. Among them, the accuracy of PS extraction affects the processing accuracy and robustness of subsequent steps, and finally affects the surface deformation field. Extraction accuracy.

根据意大利学者A.Ferretti对PS的定义，PS是指在较长的时间范围内，散射相位保持稳定的点目标。现有的PS提取方法有振幅离差法和相干系数法等。振幅离差法的基本思想是，由于在特定条件下，像元的散射幅度与散射相位具有相近的稳定程度，可通过像元的振幅稳定性提取PS。这种方法通常需要处理几十幅以上的SAR图像才能获得有效的提取结果。相干系数法利用像元在整个图像序列上的平均相干系数作为其散射稳定性的评估指标，但是该算法的缺陷是容易出现的PS块状聚集现象。受SAR系统数据获取能力的限制，现实应用中往往需要利用小数据集（图像数目为8到15幅）的SAR图像序列提取区域形变场。由于现有的PS提取方法应用于小数据集SAR图像序列难以获得令人满意的提取结果，这时我们就要采用新的、有效的PS提取方法。According to the definition of PS by Italian scholar A.Ferretti, PS refers to a point target whose scattering phase remains stable in a long time range. Existing PS extraction methods include amplitude dispersion method and coherence coefficient method. The basic idea of the amplitude dispersion method is that under certain conditions, the scattering amplitude and scattering phase of the pixel have a similar degree of stability, and the PS can be extracted through the amplitude stability of the pixel. This method usually needs to process more than dozens of SAR images to obtain effective extraction results. The coherence coefficient method uses the average coherence coefficient of the pixel in the entire image sequence as the evaluation index of its scattering stability, but the defect of this algorithm is the PS block aggregation phenomenon that is easy to appear. Limited by the data acquisition capability of the SAR system, it is often necessary to extract the regional deformation field using SAR image sequences with small data sets (the number of images is 8 to 15) in practical applications. Since the existing PS extraction methods are difficult to obtain satisfactory extraction results when applied to small data set SAR image sequences, we need to adopt a new and effective PS extraction method.

发明内容Contents of the invention

本发明要解决的技术问题是，提供一种适用于小数据集的、有效的PS提取的方法。本方法利用了SAR图像数据的时空域信息，从而减少了对SAR图像数目的要求，在利用小数据集SAR图像提取PS实现突破。The technical problem to be solved by the present invention is to provide an effective PS extraction method suitable for small data sets. This method utilizes the time and space domain information of SAR image data, thereby reducing the requirement for the number of SAR images, and achieves a breakthrough in extracting PS from SAR images with small data sets.

本发明技术方案的基本思路是，由于像元信杂比和散射相位稳定度存在正相关关系，首先利用SAR图像数据的时空域信息准确估计SAR图像中各像元的信杂比，然后依据信杂比的大小判定各像元是否为PS。The basic idea of the technical solution of the present invention is that since there is a positive correlation between the signal-to-clutter ratio of the pixel and the scattering phase stability, firstly, the SNR of each pixel in the SAR image is accurately estimated by using the time-space domain information of the SAR image data, and then the signal-to-clutter ratio of each pixel in the SAR image is accurately estimated, and then The size of the clutter ratio determines whether each pixel is PS.

本发明的技术方案是，基于星载SAR系统不同时刻对同一区域获取的M(8≤M≤15)幅SAR图像序列S_i(1≤i≤M)，进行以下处理：The technical solution of the present invention is to perform the following processing on the M (8≤M≤15) SAR image sequence S_i (1≤i≤M) acquired in the same area based on the spaceborne SAR system at different times:

步骤一：SAR图像序列预处理。Step 1: SAR image sequence preprocessing.

SAR图像序列预处理包括SAR图像序列配准、重采样和辐射校正三个步骤，具体方法参照博士学位论文《基于永久散射体雷达差分干涉探测区域地表形变的研究》（西南交通大学，陈强，2006.8，第40页至第66页）进行处理。SAR image sequence preprocessing includes three steps of SAR image sequence registration, resampling and radiation correction. For the specific method, refer to the doctoral dissertation "Research on Regional Surface Deformation Based on Permanent Scatterer Radar Differential Interference Detection" (Southwest Jiaotong University, Chen Qiang, 2006.8 , pages 40 to 66) for processing.

步骤二：信杂比估计。该步骤包括空域信杂比估计和时域信杂比估计。Step 2: Signal-to-clutter ratio estimation. This step includes estimating the SNR in the space domain and estimating the SNR in the time domain.

其中，空域信杂比估计包括以下三个步骤：Among them, the airspace signal-to-clutter ratio estimation includes the following three steps:

第（一）步，SAR图像序列时域平均滤波。In the first step, the temporal average filtering of the SAR image sequence is performed.

对预处理后的M幅SAR图像序列S_i，利用下式进行时域平均滤波，得到平均幅度图For the preprocessed M SAR image sequence S_i , use the following formula to perform time-domain average filtering to obtain the average amplitude map

$\stackrel{<span><span>\&OverBar;</span>\&OverBar;</span>}{\mathrm{<span><span>S</span>S</span>}}<span><span>=</span>=</span>\frac{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>M</span>m</span>}}\underset{\mathrm{<span><span>i</span>i</span>}<span><span>=</span>=</span>\mathrm{<span><span>1</span>1</span>}}{\overset{\mathrm{<span><span>M</span>m</span>}}{\mathrm{<span><span>\&Sigma;</span>\&Sigma;</span>}}}{\mathrm{<span><span>S</span>S</span>}}_{\mathrm{<span><span>i</span>i</span>}}$

第（二）步，估计空域信杂比。对平均幅度图利用下式估计第j个像元的杂波强度J表示平均幅度图的像元个数。The second step is to estimate the airspace signal-to-clutter ratio. plot of average magnitude Use the following formula to estimate the clutter intensity of the jth pixel J for mean magnitude graph the number of pixels.

${\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>I</span>I</span>}}}_{\mathrm{<span><span>j</span>j</span>}}<span><span>=</span>=</span>\frac{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>K</span>K</span>}}\underset{\mathrm{<span><span>k</span>k</span>}<span><span>=</span>=</span>\mathrm{<span><span>1</span>1</span>}}{\overset{\mathrm{<span><span>K</span>K</span>}}{\mathrm{<span><span>\&Sigma;</span>\&Sigma;</span>}}}{{\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>jk</span>jk</span>}}}^{\mathrm{<span><span>2</span>2</span>}}<span><span>,</span>,</span>{\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>jk</span>jk</span>}}<span><span>\&Element;</span>\&Element;</span>{\mathrm{<span><span>\&Omega;</span>\&Omega;</span>}}_{\mathrm{<span><span>j</span>j</span>}}$

其中，Ω_j为第j个像元的杂波估计区域，使用空域加窗法选定（具体方法可参考《合成孔径雷达目标检测理论、算法及应用》第122到125页）。a_jk为平均幅度图第j个像元杂波估计区域Ω_j内的第k个像元的振幅，K为杂波估计区域Ω_j的像元个数。Among them, Ω_j is the clutter estimation area of the jth pixel, which is selected using the spatial windowing method (for specific methods, please refer to pages 122 to 125 of "Synthetic Aperture Radar Target Detection Theory, Algorithms, and Applications"). a_jk is the average amplitude map The amplitude of the kth pixel in the jth pixel clutter estimation area Ω_j , K is the number of pixels in the clutter estimation area Ω_j .

假设第j个像元在平均幅度图的振幅为P_j，则其空域信杂比为：Suppose the jth pixel is in the mean magnitude map The amplitude is P_j , then its spatial signal-to-clutter ratio for:

${\mathrm{<span><span>S</span>S</span>}\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>C</span>C</span>}}\mathrm{<span><span>R</span>R</span>}}_{{\mathrm{<span><span>1</span>1</span>}}_{\mathrm{<span><span>j</span>j</span>}}}<span><span>=</span>=</span>\left\{\begin{array}{cc}\mathrm{<span><span>10</span>10</span>}\mathrm{<span><span>log</span>log</span>}<span><span>(</span>(</span>\frac{{{\mathrm{<span><span>P</span>P</span>}}_{\mathrm{<span><span>j</span>j</span>}}}^{\mathrm{<span><span>2</span>2</span>}}<span><span>-</span>-</span>{\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>I</span>I</span>}}}_{\mathrm{<span><span>j</span>j</span>}}}{{\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>I</span>I</span>}}}_{\mathrm{<span><span>j</span>j</span>}}}<span><span>)</span>)</span>& {{\mathrm{<span><span>P</span>P</span>}}_{\mathrm{<span><span>j</span>j</span>}}}^{\mathrm{<span><span>2</span>2</span>}}<span><span>></span>></span>{\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>I</span>I</span>}}}_{\mathrm{<span><span>j</span>j</span>}}\\ \mathrm{<span><span>0</span>0</span>}& {{\mathrm{<span><span>P</span>P</span>}}_{\mathrm{<span><span>j</span>j</span>}}}^{\mathrm{<span><span>2</span>2</span>}}<span><span>\&le;</span>\&le;</span>{\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>I</span>I</span>}}}_{\mathrm{<span><span>j</span>j</span>}}\end{array}\right.<span><span>;</span>;</span>$

另外，时域信杂比估计包括以下两个步骤：In addition, the time-domain signal-to-clutter ratio estimation includes the following two steps:

第（一）步，估计像元的振幅离差值。The first step is to estimate the amplitude dispersion value of the pixel.

基于预处理后得到的SAR图像序列，利用下式估计第j个像元的振幅离差值，j＝1,2,…,J：Based on the SAR image sequence obtained after preprocessing, use the following formula to estimate the amplitude deviation value of the jth pixel, j=1,2,...,J:

${\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>D</span>D.</span>}}}_{\mathrm{<span><span>j</span>j</span>}}<span><span>=</span>=</span>\frac{\sqrt{\frac{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>M</span>m</span>}}\underset{\mathrm{<span><span>i</span>i</span>}<span><span>=</span>=</span>\mathrm{<span><span>1</span>1</span>}}{\overset{\mathrm{<span><span>M</span>m</span>}}{\mathrm{<span><span>\&Sigma;</span>\&Sigma;</span>}}}{{\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>ij</span>ij</span>}}}^{\mathrm{<span><span>2</span>2</span>}}<span><span>-</span>-</span>{<span><span>(</span>(</span>\frac{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>M</span>m</span>}}\underset{\mathrm{<span><span>i</span>i</span>}<span><span>=</span>=</span>\mathrm{<span><span>1</span>1</span>}}{\overset{\mathrm{<span><span>M</span>m</span>}}{\mathrm{<span><span>\&Sigma;</span>\&Sigma;</span>}}}{\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>ij</span>ij</span>}}<span><span>)</span>)</span>}^{\mathrm{<span><span>2</span>2</span>}}}}{\frac{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>M</span>m</span>}}\underset{\mathrm{<span><span>i</span>i</span>}<span><span>=</span>=</span>\mathrm{<span><span>1</span>1</span>}}{\overset{\mathrm{<span><span>M</span>m</span>}}{\mathrm{<span><span>\&Sigma;</span>\&Sigma;</span>}}}{\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>ij</span>ij</span>}}}$

其中，a_ij为第j个像元在i时刻的SAR图像振幅。Among them, a_ij is the SAR image amplitude of the jth pixel at time i.

第（二）步，利用估计第j个像元的时域信杂比如下式所示：In the second step, use Estimate the temporal signal-to-noise ratio of the jth pixel As shown in the following formula:

${\mathrm{<span><span>S</span>S</span>}\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>C</span>C</span>}}\mathrm{<span><span>R</span>R</span>}}_{{\mathrm{<span><span>2</span>2</span>}}_{\mathrm{<span><span>j</span>j</span>}}}<span><span>=</span>=</span>\mathrm{<span><span>20</span>20</span>}\mathrm{<span><span>log</span>log</span>}<span><span>(</span>(</span>\frac{{\stackrel{<span><span>^</span>^</span>}{\mathrm{<span><span>D</span>D.</span>}}}_{\mathrm{<span><span>j</span>j</span>}}^{\mathrm{<span><span>2</span>2</span>}}}{\mathrm{<span><span>2</span>2</span>}}<span><span>)</span>)</span>$

步骤三：PS判定。Step 3: PS judgment.

设定信杂比阈值为8db（decibel，分贝），根据第j个像元的空域信杂比和时域信杂比判定其是否为PS：Set the signal-to-clutter ratio threshold to 8db (decibel, decibel), according to the spatial signal-to-clutter ratio of the jth pixel and time-domain signal-to-noise ratio Determine whether it is PS:

如和均大于阈值，则判定第j个像元为PS；否则，判定该像元为非PS。这样对所有的像元进行判定后，可得到最终的PS提取结果。like and If both are greater than the threshold, it is determined that the jth pixel is PS; otherwise, it is determined that the pixel is not PS. After all the pixels are judged in this way, the final PS extraction result can be obtained.

本发明的有益效果是：实现简单、准确度高，适用范围广。由于利用了SAR图像数据的空域和时域信息，本发明在提取精度上高于现有的提取方法。本发明放宽对SAR图像数目的要求，即使利用小数据集的SAR图像序列，也能获得比较准确的提取结果，因此实现简单、使用范围广，这是现有方法难以实现的。The invention has the beneficial effects of simple implementation, high accuracy and wide application range. Because the spatial domain and time domain information of the SAR image data are utilized, the extraction precision of the present invention is higher than that of the existing extraction method. The present invention relaxes the requirement on the number of SAR images, and can obtain relatively accurate extraction results even if the SAR image sequence of a small data set is used, so the implementation is simple and the application range is wide, which is difficult to achieve by existing methods.

附图说明Description of drawings

图1是本发明的原理流程示意图；Fig. 1 is a schematic flow chart of the principle of the present invention;

图2为利用本发明实施例进行实验的结果图；Fig. 2 is the result figure that utilizes the embodiment of the present invention to carry out experiment;

图3为对图2所示实验结果进行PS准确度评估的结论。Fig. 3 is the conclusion of PS accuracy evaluation on the experimental results shown in Fig. 2.

具体实施方式Detailed ways

图1是本发明的原理流程示意图。技术方案包括对基于星载SAR图像序列进行SAR图像序列预处理，然后进行空域信杂比估计和时域信杂比估计，通过PS判定，得到PS提取结果。Fig. 1 is a schematic flow chart of the principle of the present invention. The technical solution includes preprocessing the SAR image sequence based on the spaceborne SAR image sequence, and then performing the estimation of the signal-to-clutter ratio in the space domain and the signal-to-clutter ratio in the time domain, and obtaining the PS extraction result through PS judgment.

图2为利用本发明实施例进行实验的结果图。其中，SAR图像对应区域为美国的阿纳海姆市，横坐标为方位向(km)，纵坐标为距离向(km)，提取的PS用白色圆点标明。其中，利用的SAR图像序列包括12幅图像，即M=12，使用空心加窗法选定杂波估计区域，矩形窗口总大小为12像素×12像素，窗口空心区域大小5像素×5像素，即杂波估计区域的像元个数K=119=12像素×12像素-5像素×5像素。Fig. 2 is a graph showing the results of experiments conducted using the embodiment of the present invention. Among them, the corresponding area of the SAR image is Anaheim in the United States, the abscissa is the azimuth (km), the ordinate is the distance (km), and the extracted PS is marked with a white dot. Among them, the SAR image sequence used includes 12 images, that is, M=12, and the clutter estimation area is selected using the hollow window method. The total size of the rectangular window is 12 pixels×12 pixels, and the size of the window hollow area is 5 pixels×5 pixels. That is, the number of pixels in the clutter estimation area K=119=12 pixels×12 pixels−5 pixels×5 pixels.

图3为对图2所示实验结果进行PS准确度评估的结论。Fig. 3 is the conclusion of PS accuracy evaluation on the experimental results shown in Fig. 2.

采用的PS准确度指标为PS网络优化率（具体计算方法参照国际会议论文《Testing And Evaluation Of Permanent Scatterers Candidates Selection Methods》，国防科学技术大学，陈绍劲，2013.05）。The PS accuracy index used is the PS network optimization rate (the specific calculation method refers to the international conference paper "Testing And Evaluation Of Permanent Scatterers Candidates Selection Methods", National University of Defense Technology, Chen Shaojin, 2013.05).

利用图2实验中的SAR图像序列，分别采用现有的振幅离差法和相干系数法（参照国防科学技术大学孙希龙博士学位论文《SAR层析与差分层析成像技术研究》第30页到第40页）以及本方法提取PS，然后计算各方法的PS网络优化率指标。计算结果表明，在处理12景SAR图像时，本方法的PS准确度指标远高于现有的振幅离差法和相干系数法，这充分说明了本方法在利用小数据集SAR图像提取PS的优越性。Using the SAR image sequence in the experiment in Fig. 2, the existing amplitude dispersion method and coherence coefficient method are respectively used (refer to pages 30 to 30 of Dr. Sun Xilong's dissertation "Research on SAR Tomography and Differential Tomography" 40) and this method extracts PS, and then calculates the PS network optimization rate index of each method. The calculation results show that when processing 12 SAR images, the PS accuracy index of this method is much higher than that of the existing amplitude dispersion method and coherence coefficient method, which fully demonstrates the effectiveness of this method in extracting PS from SAR images with small data sets. Superiority.

Claims (2)

1. utilize a method for small data set SAR image sequential extraction procedures Permanent scatterers, based on Spaceborne SAR System not in the same time to the M width SAR image sequence S that the same area obtains_i, wherein SAR represents synthetic-aperture radar, 8≤M≤15, and 1≤i≤M, is characterized in that, comprises the steps:

Step one: SAR image sequence pre-service:

The pre-service of SAR image sequence comprises SAR image sequence registration, resampling and radiant correction three steps;

Step 2: signal to noise ratio is estimated:

This step comprises spatial domain signal to noise ratio and estimates and the estimation of time domain signal to noise ratio;

Wherein, spatial domain signal to noise ratio is estimated to comprise following two steps:

(1) step, the filtering of SAR image sequence time domain average:

To pretreated M width SAR image sequence S_i, utilize following formula to carry out time domain average filtering, obtain average amplitude figure

$\stackrel{\&OverBar;}{S}=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{S}_i$

(2) step, estimate spatial domain signal to noise ratio:

To average amplitude figurefollowing formula is utilized to estimate the noise intensity of a jth pixelj represents average amplitude figurepixel number;

${\hat{I}}_j=\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{a_{\mathrm{jk}}}^2,a_{\mathrm{jk}}\&Element;{\mathrm{\&Omega;}}_j$

Wherein, Ω_jfor the clutter estimation region of a jth pixel, windowing method in spatial domain is used to select; a_jkfor average amplitude figurea jth pixel clutter estimation region Ω_jthe amplitude of an interior kth pixel, K is clutter estimation region Ω_jpixel number;

Suppose that a jth pixel is at average amplitude figureamplitude be P_j, then its spatial domain signal to noise ratiofor:

${\widehat{\mathrm{SCR}}}_{1_j}=\left\{\begin{array}{cc}10\log \left(\frac{{P_j}^2-{\hat{I}}_j}{{\hat{I}}_j}\right)& {P_j}^2>{\hat{I}}_j\\ 0& {P_j}^2\&le;{\hat{I}}_j\end{array}\right.;$

Wherein, time domain signal to noise ratio is estimated to comprise following two steps:

(1) step, estimate the amplitude deviation value of pixel:

Based on the SAR image sequence obtained after pre-service, following formula is utilized to estimate the amplitude deviation value of a jth pixel, j=1,2 ..., J:

${\hat{D}}_j=\frac{\sqrt{\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{a_{\mathrm{ij}}}^2-{\left(\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ij}}\right)}^2}}{\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_{\mathrm{ij}}}$

Wherein, a_ijfor a jth pixel is at the SAR image amplitude in i moment;

(2) step, utilizesestimate the time domain signal to noise ratio of a jth pixelbe shown below:

${\widehat{\mathrm{SCR}}}_{2_j}=20\log \left(\frac{{\hat{D}}_j^2}{2}\right);$

Step 3: PS judges;

Setting signal to noise ratio threshold value is 8 decibels, according to the spatial domain signal to noise ratio of a jth pixelwith time domain signal to noise ratiojudge it whether as PS:

Aswithall be greater than threshold value, then judge that a jth pixel is as PS; Otherwise, judge that this pixel is as non-PS; After so all pixels being judged, final PS can be obtained and extract result; Wherein PS represents Permanent scatterers.

2. the method utilizing small data set SAR image sequential extraction procedures Permanent scatterers according to claim 1, it is characterized in that, when using hollow windowing method to select clutter estimation region, the total size of rectangular window is 12 pixel × 12 pixels, window hollow area size 5 pixel × 5 pixel.