CN103400137B - 一种sar图像的建筑物几何参数提取方法 - Google Patents

Abstract

本发明是有关于一种SAR图像的建筑物几何参数提取方法，所述的提取方法的步骤如下：对SAR图像用滤波算法进行处理，生成滤波后图像，并生成入射角矩阵；利用线状特征检测算子从图像中检测出边界，得到边界二值图像，对边界二值图像进行数学形态学运算，得到二次散射亮线，生成真实亮度曲线，并确定方位角和入射角；确定建筑物高度、长度和权重参数的解空间，生成模拟亮度曲线集；从模拟亮度曲线集中寻找与真实亮度曲线最佳匹配的曲线，最佳匹配曲线对应的高度和长度即为该建筑物的反演参数。本发明解决了从SAR图像二次散射线提取建筑物长度偏大和常规提取建筑物高度需要大量先验知识的问题，测量所得误差小。

Description

一种SAR图像的建筑物几何参数提取方法

技术领域

本发明涉及SAR图像提取领域，特别是涉及一种SAR图像的建筑物几何参数提取方法。

背景技术

建筑物目标监测是合成孔径雷达(SyntheticApertureRadar，SAR)遥感应用的一项重要内容，在城市建设与规划、灾害损失评估、军事监视等领域具有广泛应用。SAR具有全天时、全天候、可穿透云雾快速成像的能力，在灾害应急和多云多雨地区的建筑物目标监测中具有不可替代的作用。

2007年以来，随着TerraSAR-X、Cosmo-SkyMed、Radarsat-2和TanDEM-X等卫星系统的发射成功，星载SAR图像空间分辨率达到米级，空间分辨率的不断提高已成为当前雷达遥感发展的重要方向之一。2012年11月19日，HJ-1C卫星发射成功；国家重大科技专项“高分辨率对地观测系统”和“国家自然灾害空间信息基础设施”也规划发射若干颗高分辨率SAR卫星。这些SAR卫星将为城区雷达遥感应用提供丰富的数据源，但另一方面，也给城区SAR图像应用提出了更高的要求。

由于SAR是基于测距原理通过主动发射、接收微波电磁波成像，对于地物目标的几何形态非常敏感。在SAR图像中，建筑物目标的几何形态往往体现在阴影、叠掩和多次散射特征中。现有的建筑物几何参数提取方法主要围绕这些散射特征展开。虽然高分辨率SAR卫星为城市目标三维重建和动态变化监测提供了宝贵数据源，但城市目标的高分辨率SAR图像特征非常复杂，给城区高分辨率SAR图像的应用造成很大困难。目前，基于高分辨率SAR图像的建筑物目标提取与三维重建方法仍处于非常初级的阶段，人工辅助提取水平远远优于自动算法。因此，开展城市目标高分辨率SAR目标提取与三维重建方法研究，是推动SAR遥感技术在城区监测中应用的重要方面。

在米级高分辨率SAR图像中，建筑物目标的二次散射线是一类重要特征，它是由建筑物前墙面与周边下垫面组成的二面角反射器形成一条沿墙角的亮线，对于建筑物目标的信息反演具有重要的指示意义，可以用于确定建筑物边界和反演高度。然而，在高分辨率SAR图像中，沿着二次散射亮线方向，灰度分布不均匀，甚至会出现亮线长度大于建筑物墙角线的情况，利用常规方法来提取建筑物的参数存在较大的误差。

2002年，Franceschetti等对典型建筑物的二次散射机理进行了全面深入的分析提出了采用几何光学和Kirchhoff或物理光学近似求解的理论，推导出二次散射线的后向散射截面与建筑物高度之间的定量公式如下：

$h={\mathrm{\σ}}^0\frac{{8\mathrm{\π}}^2{\cos }^2\mathrm{\θ}{\mathrm{\σ}}^2(2/l_c^2)\·\exp \left[\frac{{\tan }^2\mathrm{\θ}{\sin }^2\mathrm{\φ}}{2{\mathrm{\σ}}^2(2/l_c^2)}\right]}{{\left|S_{\mathrm{pq}}\right|}^{2}l\tan \mathrm{\θ}\cos \mathrm{\φ}(1+{\tan }^2\mathrm{\θ}{\sin }^2\mathrm{\φ})}\backslash *\mathrm{MERGEFORMAT}---\left(1\right)$

其中：σ⁰为二次散射线对应的后向散射截面，S_pq为散射矩阵，p,q分别为水平、垂直极化，与墙壁和地面的介电常数有关，l和h分别为墙面的长度和高度，θ为入射角，φ为墙面方位角，σ和l_c分别为下垫面均方根高度和相关长度。

参照图1所示，该方法反演建筑物的高度的主要步骤包括：

（1）SAR图像辐射定标：利用SAR数据中提供的定标系数对图像进行辐射定标。

（2）二次散射亮线提取及后向散射截面σ⁰计算：从图像中找到建筑物对应的二次散射线，提取建筑物的长度，取整条线及其周边像素的灰度平均值作为该二次散射线的σ⁰值。

（3）参数标定：当σ⁰已知时，无法直接利用公式计算建筑物高度，需要利用已知建筑物高度、墙面等参数对公式中其他参数进行标定，该过程需要选择2-3个有已知高度的建筑物样本。

（4）高度反演：当σ⁰和其他参数都确定后，利用公式反演建筑物高度。

上述技术方案在实际使用中尚存在以下不足：

（1）该技术方案中反演高度时需要确定建筑物及下垫面的电学特性和粗糙度参数。通常情况下，这些参数难以精确测得，从而影响高度反演精度。

（2）该技术方案通过假定所有建筑物目标及下垫面的电学特性和粗糙度参数都相同来避免了大量未知参数的输入。通常情况下，这种假设不成立，因此该假设下得出的反演算法实用性较差。

（3）该技术方案中尚未考虑实际SAR图像中建筑物二次散射线大于建筑物真实长度，因此该技术方案提取建筑物长度存在偏大的情况，并且误差大小与建筑物朝向有关。

（4）该技术方案尚未考虑SAR图像中建筑物二次散射线的后向散射截面σ⁰分布不均匀性，通过取平均值计算方法会给反演结果带来误差。

由此可见，上述现有的一种SAR图像的建筑物几何参数提取方法在使用上，显然仍存在有不便与缺陷，而亟待加以进一步改进。如何能创设一种计算误差小、测量方便的新的SAR图像的建筑物几何参数提取方法,实属当前重要研发课题之一。

发明内容

本发明要解决的技术问题是提供一种SAR图像的建筑物几何参数提取方法，使其反演误差小、可快速高效的计算建筑物的几何参数，从而克服现有的提取方法误差大、测量不准确的不足。

为解决上述技术问题，本发明提供一种SAR图像的建筑物几何参数提取方法，其特征在于所述的提取方法的步骤如下：对SAR图像用滤波算法进行处理，生成滤波后图像，并生成入射角矩阵M；利用线状特征检测算子从图像中检测出边界，得到边界二值图像，对边界二值图像进行数学形态学运算，得到二次散射亮线，生成真实亮度曲线，并确定方位角和入射角；确定建筑物高度、长度和权重参数的解空间，生成模拟亮度曲线集；从模拟亮度曲线集中寻找与真实亮度曲线最佳匹配的曲线，最佳匹配曲线对应的高度和长度即为该建筑物的反演参数。

作为本发明的一种改进，所述的滤波算法采用增强Lee滤波算法

所述的线状特征检测算子采用自适应性算数均值比算子。

所述的真实亮度曲线为由测得的亮度曲线的亮度包络曲线归一化得到的归一化曲线。

所述的SAR图像经过二次散射过程，所述的二次散射过程分为散射过程A和B；其中，散射过程A为电磁波照射到墙面后，镜面反射到达到下垫面，然后经地面漫反射回传感器；散射过程B为电磁波照射到地面后镜面反射到达墙面，然后经墙面漫反射回传感器。

所述的模拟亮度曲线集由公式0<χ<L₀+w计算得出，其中I_s(χ)为二次散射亮度曲线，I₁(χ)为散射过程A的亮度曲线，I₂(χ)散射过程B的亮度曲线，ω为I₂(χ)对应的权重，

所述的I₁(χ)和对应二次散射线的长度L₁为

$I_1\left(\mathrm{\&chi;}\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&chi;}}{\min (L_0,w)}& 0<\mathrm{\&chi;}<\min (L_0,w)\\ 1& \min (L_0,w)<\mathrm{\&chi;}<\max (L_0,w)\\ 1+\frac{\max (L_0,w)}{\min (L_0,w)}-\frac{\mathrm{\&chi;}}{\min (L_0,w)}& \max (L_0,w)<\mathrm{\&chi;}<L_0+w\end{array}\right.$

L₀=h·tanθsinφ

L₁=w+h·tanθsinφ

亮度曲线I₂(χ)和对应二次散射线的长度L₂为

I₂(χ)=1,0<χ<w

L₂=w

其中，w、h和φ分别表示建筑物墙面的长度、高度和对应方位角，方位角为墙面与方位向夹角，θ为成像入射角。

所述的反演参数计算公式为 $(w_{\mathrm{est}},h_{\mathrm{est}})=\arg \underset{\mathrm{\&Omega;}}{\max }\mathrm{\&rho;}(I\left(\mathrm{\&chi;}\right),I_s\left(\mathrm{\&chi;}\right)\&Element;S_{\mathrm{\&Omega;}}),$ I(χ)为真实亮度曲线，S_Ω为模拟亮度曲线集，ρ表示两条曲线的相关性， $\mathrm{\&rho;}(X,Y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(X_i-\stackrel{\&OverBar;}{X})(Y_i-\stackrel{\&OverBar;}{Y})}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(X_i-\stackrel{\&OverBar;}{X})}^2}\sqrt{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(Y_i-\stackrel{\&OverBar;}{Y})}^2}},$

其中，X、Y分别表示两条曲线，N为曲线的点数，X_i和Y_i分别为曲线上第i点，和分别为曲线对应的均值。

采用这样的设计后，本发明至少具有以下优点：

1、本发明提出了一种利用二次散射亮度曲线与模拟亮度曲线匹配来提取建筑物长度的快速高效方法、解决了从SAR图像二次散射线提取建筑物长度偏大和常规提取建筑物高度需要大量先验知识的问题，测量所得误差小。

2、本发明通过分析建筑物二次散射形成机制，发明了一种无需先验知识直接通过建筑物几何参数来模拟二次散射亮度曲线的方法，与常规方法相比，不需要输入参数，能够自适应快速提取每一个建筑物目标几何参数。

3、本发明不同于常规方法中串行提取建筑物高度和长度参数，而是同时解算出高度和长度，充分利用了参数之间的约束关系，提高了反演精度。

4、本发明中不要求辐射定标，对图像数据要求较低，实用性强。

附图说明

上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。

图1是反演建筑物高度的现有方法流程图。

图2是本发明的反演建筑物的几何参数的流程图。

图3是本发明生成的亮度曲线示意图。

图4中为待测建筑物图.

图5为图4建筑物对应的TerraSAR-X高分辨率聚束模式地距图像。

图6为归一化的真实亮度曲线和最佳匹配曲线。

具体实施方式

参照图2、图3所示，采用本发明提供的一种基于SAR图像的建筑物几何参数提取方法对建筑物的高度和长度进行计算，其步骤如下：

（1）对SAR图像用滤波算法进行处理，生成滤波后图像S，并生成入射角矩阵M；

（2）利用线状特征检测算子从图像S中检测出边界，得到边界二值图像S_bw，对边界二值图像S_bw进行数学形态学运算，得到二次散射亮线L_db，生成真实亮度曲线，并确定方位角φ_db和入射角θ_db；

（3）确定建筑物高度、长度和权重参数的解空间，生成模拟亮度曲线集S_Ω；

（4）从模拟亮度曲线集S_Ω中寻找与真实亮度曲线最佳匹配的曲线曲线对应的高度和长度即为该建筑物的反演参数。

具体地，对各步骤进行进一步解释。

（1）对SAR图像用滤波算法进行处理，生成滤波后图像S，并生成入射角矩阵M。

本实施例中，滤波算法优选为增强Lee滤波，表达式为：

$\left\{\begin{array}{cc}{g}_{\mathrm{ij}}^{\&prime;}={\stackrel{\&OverBar;}{g}}_{\mathrm{ij}}& {\mathrm{\&sigma;}}_{\mathrm{ij}}\&le;c_u\\ g_{\mathrm{ij}}^{\&prime;}={\stackrel{\&OverBar;}{g}}_{\mathrm{ij}}\mathrm{\&beta;}+g(1-\mathrm{\&beta;})& c_u<{\mathrm{\&sigma;}}_{\mathrm{ij}}<c_{\max }\\ g_{\mathrm{ij}}^{\&prime;}=g_{\mathrm{ij}}& {\mathrm{\&sigma;}}_{\mathrm{ij}}\&GreaterEqual;c_{\max }\end{array},\backslash *\mathrm{MERGEFORMAT}---\left(2\right)\right.$

其中是Lee滤波的权函数；L为成像视数；c_u,c_max分别为均质区域和异质区域的阈值；σ_I为局部标准差。

入射角矩阵M用SAR图像近距和远距入射角通过线性内插获得，方法如下：

${\mathrm{\&theta;}}_{i,j}=\frac{j}{\mathrm{cols}}{\mathrm{\&theta;}}_{\mathrm{far}}+(1-\frac{j}{\mathrm{cols}}){\mathrm{\&theta;}}_{\mathrm{near}}\backslash *\mathrm{MERGEFORMAT}---\left(3\right)$

其中，θ_i,j为第i行第j列对应的入射角，cols为SAR图像总列数，θ_near和θ_far分别为图像近距和远距入射角。

（2）利用线状特征检测算子和数学形态学运算提取二次散射亮线L_db，生成真实亮度曲线并确定方位角φ_ab和入射角θ_db。

首先利用线状特征检测算子从图像S中检测出边界，得到边界二值图像S_bw。线状特征检测算子优选为自适应性算数均值比ROA(RatioofAverage)算子，该算子将SAR图像的区域灰度均值作为检测特征量，假定一个以给定像素为中心的固定大小的检测窗口被经过其中心且方向为a的直线分割成两个大小相等的独立部分，设u_1a和u_2a是两个区域的像素灰度均值，计算不同方向的均值比，并取其最小值r。

r取值在0和1之间。当r趋于1时。两区域均值越接近，则属于同一块区域；若r趋于0，则两区域差别越大，待检测点处于两区域边界上。

对边界二值图像S_bw进行数学形态学运算，得到二次散射亮线L_db。由于SAR图像中通常会出现平行于距离向的非二次散射亮线，因此还需要去掉平行于距离向的亮线。

参照图3所示，加粗黑色实线为提取二次散射亮线L_db，从SAR图像S中可提取出真实亮度曲线I₀(χ)，真实亮度曲线为沿着χ方向的亮度剖面曲线，方位角φ_db为L_db与方位向夹角。利用公式（3），从入射角矩阵M中得到L_db对应的入射角θ_db。由于真实亮度曲线I₀(χ)波动较大，进一步地，用其亮度包络曲线I_con(χ)代替真实亮度曲线I₀(χ)，利用公式对亮度包络曲线I_con(χ)进行归一化，得到归一化的真实亮度曲线I(χ)。

$I\left(\mathrm{\&chi;}\right)=\frac{I_{\mathrm{con}}\left(\mathrm{\&chi;}\right)-\min I_{\mathrm{con}}\left(\mathrm{\&chi;}\right)}{\max I_{\mathrm{con}}\left(\mathrm{\&chi;}\right)-\min I_{\mathrm{con}}\left(\mathrm{\&chi;}\right)}$

\*MERGEFORMAT(5)

（3）根据建筑物长度和高度的取值空间来生成模拟亮度曲线集S_Ω。

确定建筑物高度h、长度w和权重参数的解空间Ω，如下式

\*MERGEFORMAT(6)

其中，R_w、R_h、和分别为变量w、h和取值范围，下标min和max分别为最小值和最大值。

对于Ω任意元素，在成像条件(φ_db,θ_db)下，利用公式即可生成模拟亮度曲线集S_Ω。

模拟亮度曲线集S_Ω的生成过程如下：

二次散射过程分解为散射过程A和B。其中散射过程A为电磁波照射到墙面后镜面反射到达到下垫面（过程表示为S_w），然后经地面漫反射回传感器（过程表示为D_g），整个过程表示为散射过程B为电磁波照射到地面后镜面反射到达墙面（过程表示为S_g），然后经墙面漫反射回传感器（过程表示为D_w），过程表示为因此，二次散射亮度曲线I_s(χ)为散射过程A的亮度曲线I₁(χ)和散射过程B的亮度曲线I₂(χ)的加权和，假设I₂(χ)对应的权重为I₁(χ)对应权重则二次散射亮度曲线I_s(χ)表示为公式。

其中，亮度曲线I₁(χ)和对应二次散射线的长度L₁为

$I_1\left(\mathrm{\&chi;}\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&chi;}}{\min (L_0,w)}& 0<\mathrm{\&chi;}<\min (L_0,w)\\ 1& \min (L_0,w)<\mathrm{\&chi;}<\max (L_0,w)\\ 1+\frac{\max (L_0,w)}{\min (L_0,w)}-\frac{\mathrm{\&chi;}}{\min (L_0,w)}& \max (L_0,w)<\mathrm{\&chi;}<L_0+w\end{array}\right.$

\*MERGEFORMAT(8)

L₀=h·tanθsinφ\*MERGEFORMAT(9)

L₁=w+h·tanθsinφ\*MERGEFORMAT(10)

亮度曲线I₂(χ)和对应二次散射线的长度L₂为

I₂(χ)=1,0<χ<w\*MERGEFORMAT(11)

L₂=w\*MERGEFORMAT(12)

其中w、h和φ分别表示建筑物墙面的长度、高度和对应方位角，方位角为墙面与方位向夹角，θ为成像入射角。

（4）从模拟亮度曲线集S_Ω中寻找与真实亮度曲线最佳匹配的曲线曲线对应的高度和长度即为该建筑物的反演参数(w_est,h_est)。

将二次散射亮线L_db对应的真实亮度曲线与模拟亮度曲线集S_Ω中所有曲线进行匹配运算，真实亮度曲线优选为归一化的真实亮度曲线I(χ)，相关性最大的模拟亮度曲线对应的长度和高度即为反演值(w_est,h_est)。该过程用公式表示如

$(w_{\mathrm{est}},h_{\mathrm{est}})=\arg \underset{\mathrm{\&Omega;}}{\max }\mathrm{\&rho;}(I\left(\mathrm{\&chi;}\right),I_s\left(\mathrm{\&chi;}\right)\&Element;S_{\mathrm{\&Omega;}})\backslash *\mathrm{MERGEFORMAT}---\left(13\right)$

其中，ρ表示两条曲线的相关性，曲线X和曲线Y相关性公式如。

$\mathrm{\&rho;}(X,Y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(X_i-\stackrel{\&OverBar;}{X})(Y_i-\stackrel{\&OverBar;}{Y})}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(X_i-\stackrel{\&OverBar;}{X})}^2}\sqrt{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(Y_i-\stackrel{\&OverBar;}{Y})}^2}}\backslash *\mathrm{MERGEFORMAT}---\left(14\right)$

其中，N为曲线的点数，X_i和Y_i分别为曲线上第i点，和分别为曲线对应的均值。

参照图4-6所示,以一具体建筑物为例提取其几何参数。图4为一栋待测建筑物，其视向方向为电磁波入射方向，该建筑物的真实长度和高度为(w_ture,h_true)=(52.0m,171.0m)。参照图5所示，该建筑物对应的TerraSAR-X高分辨率聚束模式地距图像，像元大小为0.5m，框线部分为建筑物叠掩区，直线为墙面与地面形成的二次散射线，对应的方位角和入射角分别为9.6°和42°；通过上述方法得到该建筑物的归一化的真实亮度曲线和最佳匹配曲线,参照图6所示，其中实线为归一化的真实亮度曲线，虚线为最佳匹配曲线，经反演得到要提取的建筑物的长度和高度为(w_est,h_est)=(50.5m,172.0m)，长度和高度对应的误差为(e_w,e_h)=(-1.5m,1.0m)。由此可说明本发明提出的一种利用二次散射亮度曲线与模拟亮度曲线匹配方法能够有效提取建筑物长度和高度，误差小。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。

Claims (6)

1.一种SAR图像的建筑物几何参数提取方法，其特征在于所述的提取方法的步骤如下：

对SAR图像用滤波算法进行处理，生成滤波后图像，并生成入射角矩阵M；

利用线状特征检测算子从图像中检测出边界，得到边界二值图像，对边界二值图像进行数学形态学运算，得到二次散射亮线，生成真实亮度曲线，并确定方位角和入射角；

确定建筑物高度、长度和权重参数的解空间，生成模拟亮度曲线集；

从模拟亮度曲线集中寻找与真实亮度曲线最佳匹配的曲线，最佳匹配曲线对应的高度和长度即为该建筑物的反演参数；

其中，所述的SAR图像经过二次散射过程，所述的二次散射过程分为散射过程A和B；

散射过程A为电磁波照射到墙面后，镜面反射到达到下垫面，然后经地面漫反射回传感器；

散射过程B为电磁波照射到地面后镜面反射到达墙面，然后经墙面漫反射回传感器；

所述的模拟亮度曲线集由公式0＜χ＜L₀+w计算得出，其中I_s(χ)为二次散射亮度曲线，I₁(χ)为散射过程A的亮度曲线，I₂(χ)为散射过程B的亮度曲线，为I₂(χ)对应的权重，w为建筑物墙面的长度。

2.根据权利要求1所述的一种SAR图像的建筑物几何参数提取方法，其特征在于：所述的滤波算法采用增强Lee滤波算法。

3.根据权利要求1所述的一种SAR图像的建筑物几何参数提取方法，其特征在于：所述的线状特征检测算子采用自适应性算数均值比算子。

4.根据权利要求1所述的一种SAR图像的建筑物几何参数提取方法，其特征在于：所述的真实亮度曲线为由测得的亮度曲线的亮度包络曲线归一化得到的归一化曲线。

5.根据权利要求1所述的一种SAR图像的建筑物几何参数提取方法，其特征在于：所述的I₁(χ)和对应二次散射线的长度L₁为

$I_1\left(\mathrm{\&chi;}\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&chi;}}{min(L_0,w)}& 0<\mathrm{\&chi;}<min(L_0,w)\\ 1& min(L_0,w)<\mathrm{\&chi;}<max(L_0,w)\\ 1+\frac{max(L_0,w)}{min(L_0,w)}-\frac{\mathrm{\&chi;}}{min(L_0,w)}& max(L_0,w)<\mathrm{\&chi;}<L_0+w\end{array}\right.$

L₀＝h·tanθsinφ

L₁＝w+h·tanθsinφ

亮度曲线I₂(χ)和对应二次散射线的长度L₂为

I₂(χ)＝1,0＜χ＜w

L₂＝w

其中，w、h和φ分别表示建筑物墙面的长度、高度和对应方位角，方位角为墙面与方位向夹角，θ为成像入射角。

6.根据权利要求1所述的一种SAR图像的建筑物几何参数提取方法，其特征在于：

所述的反演参数计算公式为 $(w_{est},h_{est})=\arg \underset{\mathrm{\&Omega;}}{max}\mathrm{\&rho;}\left(I\right(\mathrm{\&chi;}),I_s(\mathrm{\&chi;})\&Element;S_{\mathrm{\&Omega;}}),$ I(χ)为真实亮度曲线，S_Ω为模拟亮度曲线集，ρ表示两条曲线的相关性， $\mathrm{\&rho;}(X,Y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(X_i-\stackrel{\&OverBar;}{X})(Y_i-\stackrel{\&OverBar;}{Y})}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(X_i-\stackrel{\&OverBar;}{X})}^2}\sqrt{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(Y_i-\stackrel{\&OverBar;}{Y})}^2}},$

其中，X、Y分别表示两条曲线，N为曲线的点数，X_i和Y_i分别为曲线上第i点，和分别为曲线对应的均值。