复旦大学
硕士学位论文
合成孔径雷达图像反演地面高程与多时相图像变化信息自动提
取的研究
姓名：罗霖
申请学位级别：硕士
专业：电磁场与微波技术
指导教师：金亚秋
20030506
合成孔径雷达图像反演地面高程与多时相图像变化信息自动提取的研究
摘要
本研究主要包括两部分的内容。
第一部分是用合成孔径雷达（ＧＡＲ）单次飞行全极化图像数据进行地面数字高程（ＤＥＭ）反演的研究。在合成孔径雷达对地观测的全极化散射回波中，倾斜地表面的同极化或交叉极化回波的极大值所对应入射波极化的取向角沙≠０，用这种∥≠０的迁移来推算地表面坡度，并用二次相近飞行的ＳＡＲ或干涉ＳＡＲ（ＩＮＳＡＲ）图像进行地面数字高程的反演已有重要的研究。本文用全极化散射Ｍｕｅｌｌｅｒ矩阵解，将矿的这种迁移表示为ＳＡＲ图像数据中三个散射Ｓｔｏｋｅｓ参数，。，，。ｕ，的函数。通过ｇｕｌｅｒ角变换，把ｙ的迁移与倾斜地表面的水平方位角、射程角以及ＳＡＲ观测的几何结构联：系起来。在只能获得单次飞行ＳＡＲ数据的情况下，用倾斜地表水平方位排列产生的图像纹理作形态学细化算法来确定水平方位角。这样用一次飞行的ＳＡＲ图像数据，确定水平方位角和射程角及其坡度分布，用完整多重网格算法数值求解地面ｉ每程的Ｐｏｉｓｓｏｎ方程，得到地形ＤＥＭ的反演。
在第二部分中，本文提出了双阂值ＥＭ（ＥｘｐｅｃｔａｔｉｏｎＭａｘｉｍｕｍ期望极大化）与ＭＲＦ（ＭａｒｋｏｖＲａｎｄｏｍＦｉｅｌｄ马尔科夫随机场）算法，对于多时相ＳＡＲ遥感图像差值自动检测与分析城市地表的变化。不周时间获取的ＳＡＰ，图像各像素的差值表现了该像素点散射的强弱变化，用ＥＭ算法将差值图像分成：增强区、减弱区和不变区，获取两个闽值。将ＥＭ的结果作为初始值，利用ＭＲＦ的ＩＣＭ（ＩｔｅｒａｔｅｄＣｏｎｄｉｔｉｏｎａｌＭｏｄｅｓ迭代条件模式）算法迭代计算，得到关于城市地表的空间相关结构变化的分类结果。用ＥＲＳ．２ＳＡＲ对中国上海１９９６和２００２年观测图像进行了ＥＭ—ＭＲＦ自动检测与分析的实例研究，上海城市地区的实际变化证实了这种方法的有效性。
关键词合成孔径雷达ＳＡＲ，全极化Ｓｔｏｋｅｓ参数，ｙ迁移，方位角与射程角
地形ＤＥＭ，多时相ＳＡＲ，双阁值ＥＭ，ＭＲＦ，上海市
分类号ＴＰ７５１．１：０４４１．４
ＲｅｓｅａｒｃｈｏｎＤＥＭＩｎｖｅｒｓｉｏｎｆｒｏｍＳＡＩｌＩｍａｇｅｓａｎｄＡｕｔｏｍａｔｉｃＥｘｔｒａｃｔｉｏｎｏｆＭｕｌｔｉ－ｔｅｍｐｏｒａｌＩｍａｇｅｓｆｏｒＣｈａｎｇｅＤｅｔｅｃｔｉｏｎ
ＡＢＳＴＲＡＣＴ
ＴｈｉｓｄｉｓｓｅｒｔａｔｉｏｎｉｓｄｉｖｉｄｅｄＩｎｔｏｔｗｏｐａｒｔｓ．
Ｉｎｔｈｅｆｉｒｓｔｐａｒｔ，ｔｈｅｒｅｓｅａｒｃｈａｏｏｕｔｔｅｒｒａｉｎｔｏｐｏｇｒａｐｈｉｃｉｎｖｅｒｓｉｏｎｕｓｉｎｇｓｉｎｇｌｅ－ｐａｓｓｐｏｌａｒｉｍｅｔｒｉｃＳＡＲｉｍａｇｅｄａｔａｉｓｄｉｓｃｕｓｓｅｄ．ＩＩｈａｓｂｅｅｎｓｔｕｄｉｅｄｔｈａｔｔｈｅｓｈｉｆｔｏｆｐｏｌａｒｉｚａｔｉｏｎｏｒｉｅｎｔａｔｉｏｎａｎｇｌｅ∥ａｔｔｈｅｍａｘｉｍｕｍｏｆｃｏ—ｐｏｌａｒｉｚｅｄｏｒｃｒｏｓｓ—ｐｏｌａｒｉｚｅｄｂａｃｋ—ｓｃａｔｔｅｒｉｎｇｓｉｇｎａｔｕｒｅｃａｎｂｅｃｏｎｖｅｒｔｅｄｔｏｔｈｅｓｕｒｆａｃｅｓｌｏｐｅｓ．Ｉｔｉｓｔｈｅｎｕｔｉｌｉｚｅｄｔｏｇｅｎｅｒａｔｅｔｈｅｄｉｇｉｔａｌｅｌｅｖａｔｉｏｎｍａｐｐｉｎｇ（ＤＥＭ）ａｎｄｔｅｒｒａｉｎｔｏｐｏｇｒａｐｈｙｕｓｉｎｇｔｗｏ—ｐａｓｓｆｕｌｌｙｐｏｌａｒｉｍｅｔｒｉｃＳＡＰ，ｏｒｉｎｔｅｒｆｅｒｏｍｅｔｔ＇ｉｃＳＡＲ（ＩＮＳＡＲ）ｉｍａｇｅｄａｔａ．Ｔｈｉｓｐａｐｅｒ，ｕｓｉｎｇｔｈｅＭｕｅｌｌｅｒｍａｔｒｉｘｓｏｌｕｔｉｏｎｎｅｗｌｙｄｅｒｉｖｅｓｔｈｅｗｓｈｉｆｔａｓａｆｕｎｃｔｉｏｎｏｆｔｈｒｅｅＳｔｏｋｅｓｐａｒａｍｅｔｅｒｓ，１＂ｉ”Ｕ；，ｗｈｉｃｈａｒｅｍｅａｓｕｒｅｄｂｙｐｏｌａｒｉｍｅｔｒｉｃＳＡＲｉｍａｇｅｄａｔａＵｓｉｎｇｔｈｅＥｕｌｅｒａｎｇｌｅｓｔｒａｎｓｆｏｒｍａｔｉｏｎ，ｔｈｅｏｒｉｅｎｔａｔｉｏｎａｎｇｌｅ∥ｉｓｒｅｌａｔｅｄｔｏｂｏｔｈｔｈｅ
ｒａｎｇｅａｎｄａｚｉｍｕｔｈａｎｇｌｅｓｏｆｔｈｅｔｉｌｔｅｄｓｕｒｆａｃｅａｎｄｒａｄａｒｖｉｅｗｉｎｇｇｅｏｍｅｔｒｙ．Ｗｈｅｎ
ｏｎｌｙ
ａｓｉｎｇｌｅ—ｐａｓｓＳＡＲｄａｔａｉｓａｖａｉｌａｂｌｅ，ｔｈｅａｄａｐｔｉｖｅｔｈｒｅｓｈｏｌｄｉｎｇｍｅｔｈｏｄａｎｄｉｍａｇｅｍｏｒｐｈｏｌｏｇｉｃａｌｔｈｉｎｎｉｎｇａｌｇｏｒｉｔｈｍｆｏｒｌｉｎｅａｒｔｅｘｔｕｒｅｓａｒｅｐｒｏｐｏｓｅｄｔｏｆｉｒｓｔｄｅｔｅｒｍｉｎｅｔｈｅａｚｉｍｕｔｈａｎｇｌｅ．Ｔｈｅｎ，ｍａｋｉｎｇｕｓｅｏｆｆｕｌｌｍｕｌｔｉ—ｇｒｉｄａｌｇｏｒｉｔｈｍ．ｂｏｔｈｔｈｅｒａｎｇｅａｎｄａｚｉｍｕｔｈａｎｇｌｅｓａｒｅｕｔｉｌｉｚｅｄｔｏｓｏｌｖｅｔｈｅＰｏｉｓｓｏｎｅｑｕａｔｉｏｎｏｆＤＥＭｔｏ
ｐｒｏｄｕｃｅｔｈｅｔｅｒｒａｉｎｔｏｐｏｇｒａｐｈｙ．
ｉｎｔｈｅｓｅｃｏｎｄｐａｒｔ，ａｄｏｕｂｌｅ?ｔｈｒｅｓｈｏｌｄＥＭ（ＥｘｐｅｃｔａｔｉｏｎＭａｘｉｍｕｍ）ａｎｄＭＲＦ（ＭａｒｋｏｖＲａｎｄｏｍＦｉｅｌｄ）ａｌｇｏｒｉｔｔｍｌｉｓｄｅｖｅｌｏｐｅｄｔｏａｕｔｏｍａｔｉｃａｌｌｙｄｅｔｅｃｔａｎｄａｎａｌｙｓｉｓｔｈｅｃｈａｎｇｅｓｏｆｕｒｂａｎａｒｅａｓｆｏｒｍｕｌｔｉｔｅｍｐｏｒａｌＳＡＲｉｍａｇｅｓ．ＴｈｅｄｉｆｆｇｒｅｎｃｅｖａｌｕｅｏｆｅａｃｈｐｉｘｅｌｉｎｓｓｅｑｕｅｎｃｅｏｆＳＡＲｈｎａｇｅｓｅｘｈｉｂｉｔｓｔｈｅｃｈａｎｇｅｓｉｎｔｈｅｍａｇｎｉｔｕｄｅｏｆｔｈｅｓｃａｔｔｅｒｉｎｇｉｎｔｅｎｓｉｔｉｅｓ．ＴｈｅｄｉｆｆｅｒｅｎｃｅｉｍａｇｅｃａｎｂｅｄｉｖｉｄｅｄｉｎｔｏｔｈｒｅｅｒｅｇｉｏｎｓｂｙＥＭａｌｇｏｒｉｔｈｍ：ｅｎｈａｎｃｅｄｒｅｇｉｏｎ，ｗｅａｋｅｎｅｄｒｅｇｉｏｎａｎｄｃｈａｎｇｅｌｅｓｓｒｅｇｉｏｎＴａｋｉｎｇｔｈｅｏｕｔｃｏｍｅｏｆＥＭａｌｇｏｒｉｔｈｍａｓｔｈｅｉａｉｔｉａｌｖａｌｕｅ，ｔｈｅＩＣＭ（ＩｔｅｒａｔｅｄＣｏｎｄｉｔｉｏｎａｌＭｏｄｅｓ）ａｌｇｏｒｉｔｈｍｂａｓｅｏｎｔｈｅＭＲＦｔｈｅｏｒｙｉｓｐｒｏｐｏｓｅｄｔｏｏｂｔａｉｎｔｈｅｃｌａｓｓｉｆｉｃａｔｉｏｎｍａｐｏｆｔｈｅｓｐａｔｉａｌ—ｃｏｎｔｅｘｔｕａｌｃｈａｎｇｅｓｏｆｕｒｂａｒ｛ａｒｅａｓ．Ａｎｅｘａｍｐｌｅｉｓｇｉｖｅｎｔｏａｕｔｏｍａｔｉｃａｌｌｙｄｅｔｅｃｔａｎｄａｎａｌｙｓｉｓｔｈｅｃｈａｎｇｅｓｆｏｒｔｈｅ１９９６ａｎｄ２００２ＥＲＳ－２ＳＡＲｉｍａｇｅｓｉｎＳｈａｎｇｈａｉ，
Ｔ１
Ｃｈｉｎａ，ｗｈｉｃｈｃｏｎｆｉｒｍｓｔｈｅｅｆｆｅｃｔｉ＇ｖｅｎｅｓｓｏｆｏａｒＥＭ—ＭＲＦａｌｇｏｒｉｔｈｍ．
Ｋｅｙｗｏｒｄ：ｆｕｌｌｙｐｏｌａｒｉｍｅｔｒｉｃＳｔｏｋｅｓｐａｒａｍｅｔｅｒ，矿ｓｈｉｆｔ，ａｚｉｍｕｔｈａｎｇｌｅａｎｄｒａｎｇｅａｎｇｌｅ，ｔｅｒｒａｉｎｔｏｐｏｇｒａｐｈ）ｒＤＥＭ、ｍｕｌｔｉｔｅｍｐｏｒａｌＳＡＲ，ｄｏｕｂｌｅ－ｔｈｒｅｓｈｏｌｄＥＭ，ＭＲＦ，Ｓｈａｎｇｈａｉ
Ｉｎｄｅｘ：ＴＰ７５１１；０４４１．４
ＴＴＴ
第一章绪论
近年来用全极化ＳＡＲ图像对地表分类已经得到了广泛的研究。从全极化ＳｉＲ图像数据，可以得到２×２维的复散射振幅函数，。（ｐ，ｑ＝ｖ，ｈ）和４×４维的实Ｍｕｅｌ】ｅｒ矩阵Ｍ．．（ｉ，，＝１，．．４）“１。同极化和交叉极化后向散射强度是入射波极化椭圆角ｚ与取向角ｖ的函数。当地面平坦时，后向散射信号一般在入射极化取向角∥＝０处达到极大值。但是，当地面有倾斜坡度时，同极化或交叉极化后向散射信号的极大值会从∥＝０处发生迁移。这种ｖ的迁移可用来推算地表面坡度。文献［２］用实同极化和零交叉极化的散射振幅函数的假设条件推导了用散射振幅函数表示的∥的迁移。由于要同时确定水平方位角和射程角，所以需要两次或多次相近飞行的ＳＡＲ图像数据来确定表面坡度。近年来，用两次相近飞行的ＳＡＲ和ＩＮＳＡＲ数据来反演地面数字高程（ＤＥＭ）的方法已有很好的实例”。５。
但是，散射信号是随机散射介质回波功率的系综平均。测量到的散射Ｓｔｏｋｅｓ参数应该是直接与妒迁移相关的。本文角全极化散射Ｍｕｅｌｌｅｒ矩阵解，推导∥的迁移，将其表示为三个散射Ｓｔｏｋｅｓ参数』。，∽ｕ。的函数，而』，，．』。Ｌｆ，是由ＳＡＲ测量得到的。利用主坐标系和像索元的局部坐标系之间的Ｅｕｌｅｒ角变换，把入射波极化取向角矿与射程角∥、水平方位角ｙ以及ＳＡＲ观测的几何构造直接联系起来。文献［２，３］的ｖ迁移是本文推导的一种近似。在很多时候，往往只能有单次飞行的ＳＡＲ数据。本文提出倾斜地面的水平方位排列在ＳＡＲ图像中表现出来的线性纹理，用图像形态学细化算法确定图像中各像素元的地形水平方位角Ｙ，再用，。，』。，ｕ，计算∥迁移来确定射程角卢。最后用∥和ｙ来确定各像素元地形水平方位坡度与射程方位坡度。用完整多重网格算法“３数值求解地面高程ＰｏｉＳＳＯＩＩ方程，完成ＳＡＲ观测区域地表ＤＥＭ的反演。
城市与城市群体的发展是本世纪人类文明发展的重大内容，城市建筑、绿地、道路的发展、水体、植物生态、气候的变化、人Ｉ：３的迁移等等，产生的城市地域与内涵的变化都需要科学的评估。多时相的空间遥感观测为监测与评估这一发展提供了便捷可行的技术手段，它的一个直接应用就是地表变化信息的识别与分类１７－］２１。全天候高分辨率的合成孔径雷达（ＳＡＲ）近二十年的发展与空间遥感计划（如ＳＩＲ—ＣＳＡＲ，ＥＲＳ一１、２，ＲＡＤＡＲＳＡＴ－１、２）的执行能提供多年观测的数据图像．为地表变化信息的获取与评估提供了条件。
比较不同时相的海量的遥感图像数据，从直观的图像灰度值变化作人工定性的评估是十分粗糙的，难以准确地确定区域的变化，判断的闽值不可靠，极为容
易失去重要的信息。因此，如何自动获取与分析地表变化信息一直是遥感研究的重要课题。
地表变化信息的检测与分析有两类方法：有监督的方法和无监督的方法。前者需要多时相的大量实况数据作为分类过程的训练样本，后者则不需要附加信息，直接对多时相的遥感图像进行分析。本文采用无监督方法，从而并不需求大量的实况数据调查。
ＲＪｇｎｏｔ等人｛７】曾用二幅不同时间的ＳＡＲ图像之间的差值作过变化分析。ｓｉｎｇｈ【８】用红外多通道数据进行娈化矢量分析，由各对应像素计算得到的变化矢量幅值的统计来检测变化。这些：亨法相对简单易用，但不能对差值图像进行自动或非启发式的分析，对差值图像设定的阈值通常是基于经验或人为的，而阈值选取适当与否在很大程度上决定了检测结果的可靠与准确。
Ｂｒｕｚｚｏｎｅ和Ｐｒｉｅｔｏ一１近年里提出了基于Ｂａｙｅｓ理论和Ｍａｒｋｏｖ模型的自动检测区域变化的方法，通过期望极大化（ＥＭ）算法来自动获取一个阂值。他们将差值图像分成两类：变化区和无变化区。但是，他们并没有用真实ＳＡＲ图像数据来说明这一算法的可靠性与可｜亍性。同时，即使是变化区，也可以有完全不同的变化。比如楼群变成草地，后同散射可能变弱；草地上有树或变成双向反射增强的楼群，后向散射可能增强。后向散射同一ｄｂ的增强或减弱完全可以有不同的地表变化的信息内涵。随着多通道、多极化、多类遥感器多时相ＳＡＲ图像的获取，自动检测、识别、分类与评估地表变化更多的信息～定会成为重要的研究方向。
本文提出ＳＡＲ差值图像增强、减弱与不变三类区域的双闽值ＥＭ和ＭＲＹ算法，在多通道、多极化、多时相ＳＡＲ数据获取后，可很容易地推广到更多分类区域的分类检测。本文用１９９ｄ和２００２年上海市的ＥＲＳ一２ＳＡＲ图像作了实例研究，有实地情况佐证。本文的方法同样可应用于我国长江三角洲、珠江流域等城市链发展的自动识别变化信息评估。
第二章用合成孔径雷达（ＳＡＲ）单次飞行全极化图像数据进行地面数字高程（ＤＥＭ）反演
§２．１散射Ｓｔｏｋｅｓ参数表达取向角妒迁移
考虑～极化电磁波以角度印～臼，，纪）入射到地面上，其极化可用椭圆角ｚ和取向角∥来表征…。散射场可写：勾
一ｎ＊ｒ—
Ｅｔ（臼ｓ，织）。杀Ｆ（臼坤ｓ；ｌｒ一日；，纪）‘￡ｒ（万一只，Ｐ，；办７／）（１）其中２ｘ２维的复散射振幅函数Ｆ。（Ｐ，ｇ＝ｋ是）可由全极化散射测量技术得到，它包含了散射层的体散射和面散射。由式（１），可推算得到Ｓｔｏｋｅｓ矢量的实Ｍｕｅｌｌｅｒ矩阵解“１，
，ｓ（臼，，妒。）＝Ｍ（Ｏ；，妒。；丌一Ｏｉ，妒。）?，。（厅一０ｉ，妒，；ｚ，７／）（２）其中归一化的入射Ｓｔｏｋｅｓ矢量可写成
７，（ｚ，妒）＝【，。，』ｍ，Ｕ，，Ｉ，】７
＝［ｏ．ｓｏ－－ＣＯＳ２Ｚｃｏｓ２７／），０．５Ｃ＋ｃｏｓ２２＂ｃｏｓ２Ｉｆ，），～ｃｏｓ２ｚｓｉｎ２Ｐ，ｓｉｎ２ｚｒｈ”Ｍｕｅ］ｌｅｒ矩阵元素ＭＦ（ｆ，／＝１，…，４）由＜Ｌ《＞（ｐ，ｑ，ｓ，ｆ－ｖ，＾）的实部或虚部的系综平均表示，表达式可见文献［１］。应注意到，同极化项＜焉＞总是比交叉极化项＜ＬＥ＞０≠ｒ）要大得多。
同极化的后向散射系数仃。旧，丌＋妒，；口一目．，仍）可表示为…，
盯。。‘４ｚｒｃｏｓＯ，只（４）其中
ｅ＝ｏ?５ｐ。（１一ｃｏｓ２Ｚｃｏｓ２９）＋Ｉｈ（１＋ｃｏｓ２２＂ｃｏｓ２妒）＋Ｕ。ｃｏｓ２Ｚｓｉｎ２７／＋Ｖ，ｓｉｎ２ｚ］
（５）当地面平坦时，同极化后向散射仃，（ｚ，７／）的极大值一般总在妒：０处…。然而，多有研究实例表明“１３‘”，当地面颂斜时，口。的极大值会从Ｖ／：０处迁移。
在识／ａｙ＝０处，以取极大值。假设ｚ＝０的对称情形，由（ｊ）式可得到Ｑ。一‰～Ｉ。ｉ）ｓｉｎ２７／＋Ｕ。ｃｏｓ２Ｖ
．
＋０．ｓＵ把＋』。）’＋Ｏ．ｓｕ：ｓｉｎ２７／＋０．５（ｉ觚～ｗ）＇ｃｏｓ２∥№。
其中上标’表示ａ／ａ∥。由（２）式的Ｍ，比较（６）式右端各项，有
０．５（１。＋，。）－～（Ｍ１３＋Ｍ２３）ｃｏｓ２ｐ，～（Ｒｅ＜只．咒＞＋Ｒｅ＜＾，瓦＞）ｃｏｓ２９（７ａ）Ｏ．ＳＵ；～Ｍ３３ｃｏｓ２∥～Ｒｅ＜‘．．矗＋Ｌ碌＞ｃｏｓ２妒（７ｂ）（，。一，。，）～０．５（Ｍ，，十Ｍ：，）ｃｏｓ２Ｐ～ｏ．５（＜Ｉ‘。ｌ２＞＋＜Ｉ圪。Ｉ２＞）ｃｏｓ２∥
（７ｃ）ｏ．ｓ（ｚｍ一，。）１～０．Ｓ（Ｍ１３一Ｍ２３）ｃｏｓ２ｙ～０．５（Ｒｅ＜‘。，二＞－Ｒｅ＜＾。瓦＞）ｃｏｓ２∥
（８ａ）Ｕ，～０．Ｓ（Ｍ３】＋ＭＩ；２）～Ｒｅ＜‘，聪＞＋Ｒｅ＜民吒＞
（８ｂ）可看到，（７ａ）与（７ｂ）均远小于（７ｃ），而（８ａ）亦远小于（８ｂ）。因此（６）式
右端的最后三项可略去。于是，得到盯。的极大值的ｙ迁移为，
洲ｚ∥＝矗㈣
ｊ
。
ＩⅥ一Ｉ。：
可以看出，第３个Ｓｔｏｋｅｓ参数Ｕ。≠０引起了ｙ的迁移。
顺便说，如果Ｕ，和ｋ一』。都趋向于０，比如，当散射来自于均匀取向散射
体或各向同性散射介质，如厚的植被顶层时，式（９）就不能很好地确定ｖ的迁移。
进一步地，从式（２）可得到
Ｕ；～０．Ｓ（Ｍ３ｌ＋Ｍ３２）＝Ｒｅ（＜Ｅ。瓦＞＋＜厶瓦＞）（１０ａ）
，。一，。～０．Ｓ（Ｍ：２一Ｍ１，）＝０．５（＜｜，ｈ１２＞一＜ＩＥ，｜２＞）（ａｏｂ）有趣的是，如果令
Ｒｅ（＜凡瓦＞＋＜厶瓦＞）＊２■‘，ｏ．５（＜Ｉ■ｆ２＞一＜ｌＥ、．｝２＞）＝‘（氏一只），这里假定，ｖ。＝Ｅ。并定义为‘，式（９）可简化为文献［３］的结果：
恤ｚｙ２告㈣，用交叉极化的后向散射系数盯，可得到相同的结果。
§２．２用ＥｕＩｅｒ角变换得到射程角和水平方位角入射波的极化矢量‘，，；．在主坐标系中定义为…，
‘，＝－；兰ｊ呈＿＝一ｓｉｎｅ）？＋ｃ。ｓ≯，萝，ｉ，：‘，ב
‘
Ｉ￡×ｋ，ｊ
…１＋
其中入射波矢量为
ｉｊ＝ｓｉｎ８ｊｃｏｓｔｉ＋ｓｉｎｏｊｓｉｎ≯ｊｉ—ｃｏｓＯｆ２（１２）（１３）
当该像素地表面有如图１所示的局部坡度时，由局部法向量；．定义的局部坐标下的极化矢量为瓦＝警，氏＝反×￡。（】４）ｌ２６ｘＫｉｂ｜
通过两个坐标系（ｉ，萝，；）和（毛，多。．ｚｂ）之间的Ｅｕｌｅｒ角（∥，，）变换…，有
ｉ＝ｃｏｓ７ｃｏｓ犀６＋ｓｈｌ痧６一ｃｏｓ７ｓｉｎ犀６（１５ａ）
多＝一ｓｉｎｙｃｏｓ＾＋ｃｏｓ痧６＋ｓｉｎｙｓｉｎ］荛６（１５ｂ）
￡＝ｓｉｎ摩＾＋ＣＯＳ犀６（１５ｃ）
图ｉ问题的几何结构
将以上三式代入式（１２），（１３）中，可得
ｅ
￡．＝（一ｃｏｓｙｃｏｓｆｌｓｉｎ矽，一ｓｉｎ７ｃｏｓｆｌｃｏｓ妒，）ｉ６＋（一ｓｉｎｙｓｉｎ妒，＋ｃｏｓｒｃｏｓ≯，）萝６（１６）＋（ｃｏｓｙｓｉｎｆｌｓｉｎ谚＋ｓｉｎｙｓｉｎｆｌｃｏｓ妒，）ｊ６
￡ｍ＝（ｃｏｓｙｃｏｓｆｌｓｉｎ０，ｃｏｓ≯，－ｓｉｎｒｃｏｓｆｌｓｉｎＯ，ｓｉｎｅ．一ｓｉｎｌ３ｃｏｓ０。）毛
＋（ｓｉｎｙｓｉｎ０，ＣＯＳ≯，＋ＣＯＳｙｓｉｎ曰ｆｓｉｎ痧），６
＋（一ｃｏｓｙｓｉｎ∥ｓｉｎ曰，ＣＯＳ∥＋ｓｉｎｙｓｉｎｆｌｓｉｎ臼，ｓｉｎ≯，一ＣＯＳｆｌｃｏｓ０，）！＾将式（１７）代入式（１４）中，谚像素地表面的极化矢量ｈ。可表示为
．ｉ；：：≮戚ｂ＂４－ｂｙｂ
。
√口２＋６２
其中
（１７）（１８）
ａ：＝－（ｓｉｎ，ｙｓｉｎ卢０，ｂ（ｃｏｓ
ＣＯＳｓｉ。ｎ。０８：ｃ‘ｏ＋ｓ。！ｂ？：ｉｓ８ｉｎｉ鬻Ｃ。０８Ｓｉ譬ｓ。ｉ’ｎ０，ｓｉｎ庐．一ｓｉｎ卢ＣＯＳ０，，ｃ，。，
＝
７ｐ？ｊ—ｙｐ？争ｊ一０７）于是，由式（１６），（１８）可以彳焉到入射波极化的取向角ｖ为
ｃｏｓ脚ｎ‘’＝堂塑筹警幽塑
ｔａｎ。：塑望！虫殳±亟２
ｌ
ｓｉｎ８ｊ—ＣＯＳｅ？ｔＫｎｆｌｃｏｓ（），＋多？）
（２０ａ）
（２０ｂ）通常令入射庐，＝０，ｉ是射程方向，声是水平方位方向，式（２０ｂ）就简化成文献［３］的结果。注意在文献［２］中，水平方位角定义为，弓和多之间的夹角。于是，该像素地表面的射程方位坡度Ｓ。和水平方位坡度Ｓ。可定义为“３
ＳＲ＝ｔａｌｌｆｌｃｏｓｒ．
Ｓ』＝ｔｍＪｌｆｌｓｉｎ７（２１’因为在式（２０ｂ）中，用一个∥不能同时确定∥和，两个角度，所以通常需要两次或多次相近飞行的ＳＡＲ图像数据，用式（２０ｂ）建立两个方程来求解∥和ｙ，这在文献［２］中已有讨论。
§２．３ＳＡＲ图像中各像素的水平方位角
如果只有单次飞行的ＳＡＲ图像数据，除了用妒与，。，ｚ２；ｕ，的关系式（９）外，还要找到另一个已知条件，用来确定两个未知的∥和ｙ。在ＳＡＲ图像中，倾斜地表面水平方位上的排列表现出的排列纹理实际上是倾斜地表排列水平取向
图２ＳＩＲ—ＣＳＡＲ图像（Ｌ波段ｖｖ，ｈｈ，ｈｖ总功率）
从图２可以看出南北走向的丘陵起伏。为了突出地形的主要走向，我们略去了一些更为细微的局部起伏变化。注意，这里定义的水平方位角ｙ是坡度走向丸与方位多的夹角（如图ｌ所示）。
在得到各个像素元的水平方位角ｙ后，用（９）式确定的∥和（２０ｂ）式得到各个像素元的口角为ｔａｌｌ∥：＿ｊ型竺ｏ（２２）
ｓｉｎｙ十ｔａｎ∥ＣＯＳＵ』ＣＯＳｙ
其中目，由ＳＡＲ观测的几何结构给出。
由卢和ｙ得到式（２１）中的ｓ。相＆，在此基础上采用文献［１０一１３］的方法来反演地形ＤＥＭ。Ｓ。和晶的计算结果如图３ａ和３ｂ所示，
匿。３ａ水平方位坡度
图３ｂ射程方位坡度
图４由ＳＡＲ反演出的地形ＤＥＭ
图５反演出的地形高度伪彩色图
图６等高线圈
反演的ＤＥＭ与原ＳＡＲ图像所反映的直觉的地形结构以及十分简略的地面高程图“”是基本一致的，其主要地形的走向与ＳＡＲ图像的明暗纹理也是吻合的。但由于缺乏精确的地面高程实况数据，尚无法作实际地形与等高线的验证。
§２．５结论
本文提出倾斜地表散射回波极值所处的∥角迁移是３个散射Ｓｔｏｋｅｓ参数的函数。它可由ＳＡＲ全极化散射的Ｍｕｅｌｌｅｔ矩阵元测量来确定。借助Ｅｕｌｅｒ角的变换，建立了由３个Ｓｔｏｋｅｓ参：数确定的∥迁移与倾斜表面像素的射程角口和水平方位角ｙ之间的联系。在只有单次飞行ＳＡＲ数据的条件下，用自适应阈值的二值化和图像形态学细化算法来确定地表水平方位ｙ的走向，从而可以用单次飞行ＳＡＲ图像来确定每个像素的两个角度∥和ｙ，以及地表面射程方位坡度和水平方位坡度。由此用完整多重网格算法数值求解ＤＥＭ的ＰｏｉＳＳＯｌ｜１方程，得到ＤＥＭ反演的有效的新方法。本文用ＳＩＲ一：数据在中国广东惠州地区进行了实例研究。
第三章双阈值ＥＭ与ＭＲＦ算法自动分析多时相～
一图
像差值ｉｔ－ＥＲＳ２ＳＡＲ
§３．１酬算法估计变化分类的先验、条件概率
考虑不同时间（ｆ，，ｒ２）的两张已经对准的大小为，×Ｌ，图像ｘ１和』２。其差值图像爿：（Ⅳ，一ｚ．），设ｘ是一个随机变量，代表了差值图像中，×，个像素的值Ｘ。＝｛Ⅳ（ｆ，／）．１ｓｉｓ』，１≤Ｊ≤ｄ）。
检测与分析差值图像的目标；黾区分不变化类由。和变化类０９。。差值图像ｘ。的像素值的概率密度分布函数芦（∞包含了这两类的条件概率分布，即有
ｐ（Ｘ）＝ｐ（ｚＩｔＯ。）Ｐ（ｃｏ。）＋ｐ（Ｘ｝ｔＯ。）Ｐ（∞。）（１）
这里Ｐ（ｃｏ。），ｔ＇（ｃｏ，）是两类的先验概率密度函数，ｐ（Ｘ１０９。）（≈＝”，ｃ）是条件概率分布函数。用期望极大化（ＥＭ）算法对ｐ（ＸｌＣＯ。）、ｐ（ＸＩ。。）、ｐ（ｃｏ。）和ｐ（ｃｏ。）进行无监督估计。
ＥＭ算法是估计不完整性数据问题的最大似然估计法的一种普遍的方法［８１。…］。它包括了一个估计步骤和一个极大化步骤，两个步骤都需要迭代到收敛。估计步骤用受到观测值控制的参数估计作为未知变量来进行计算。极大化步骤提供了对参数的新的估计。假设ｐ（ＸＩ∞。）和ｐ（Ｘｌｆ－Ｏ。）服从Ｇａｕｓｓ分布，∞。对应的概率密度函数由均值风和方差Ｄ：来描述，同样，与ｑ相关的概率密度函数由均值∥。和方羞盯：描述。可以证明”Ｊ，估计０３。的各参数与先验概率分布有公式如下：
ｖ兰！！竺１２旦’！兰！！！塑ｌ竺！！
∥－（ｑ）：二壶生—擘坐止。孙蒹。。盟警删，
∥１一号是筮鐾枷簇』。Ｐ５（ｘ（ｉ，肋
，。…羔丝警阮护朋２渺≯“＝盛≮名等而嘉≠１Ｆ：！竺！！！！：：！：！！！Ｉ竺！！耶靠』。Ｐ。（Ｘ（ｉ，／））
上标，和ｒ＋１代表当前和下一次的迭代。类似地，这些方程也用来估计吐的先验概率和条件概率密度函数的均值和方差。
参数估计先从其初始值开始通过迭代上述公式到收敛来获得。可以证明嘲【９１，每迭代一次，估计的参数获得一个对数函数形式Ｌ（Ｏ）＝Ｉｎｐ（Ｘ。１０）的增加量，其中目＝｛Ｐ（ｃａ。），尸（脚。），以，／ａ。，盯：．口；），收敛时将获得对数函数形式的本地最大值。
参数估计的初始值可以通远分析差值图像直方图ｈ（ｘ）来获得。在直方图上确定两个阈值Ｌ和疋来获得一个比较接近ｑ，的像素子集ｓ。和另一个比较接近婢的像素子集Ｓ。。我们将Ｔ＊ｎｚ写成Ｔ＝ＭＤ（１一ｙ）和Ｔ＝ＭＤ（１＋ｙ），其中ＭＤ是＾（ｚ）的中值。（即ＭＤ＝［ｍａｘ｛ＸＤ）＋ｍｉｎ｛ＸＤ｝］／２），ｙ∈（Ｏ，１）是一个用来定义那些不容易识别的像素在Ｍ。周围变化范围的初始化参数。那么集合Ｓ。＝｛ｘ（ｉ，Ｊ）Ｉｘ（ｉ，／）＜Ｌ｝和Ｓ。：：｛ｘ（ｉ，川ｘ（ｉ，，）＞瓦）就可以用来计算两类∞。和０３．、各自的初始统计参数。
在各像素值相互独立的假设下，由Ｂａｙｅｓ最小误差原理，要使差值图像ｚ。中的像素Ｘ（ｉ，Ｊ）属于类别ｑ则须使得其后验条件概率最大，即写为
ｑ＝ａｒｇｍａ）（｛Ｐ（ｑＩｘ（ｉ，朋｝
＝ａｒｇｍａｘ｛．ｒ’（ｃｏ，）ｐ（Ｘ（ｉ，Ｊ）ｌ∞，）｝（５）
这等价于在差值图像上对两类∞。和ｇＯ。求解最大似然边界ｔ。最佳阈值瓦由如下方程给出
旦盟：旦蚓蚴
ｒ６、
Ｐ（％）ｐ（Ｘｌ∞。）…’在Ｇａｕｓｓ概率密度分布的假设下，由式（６）可得
（盯：一Ｊ；）巧＋２（∥一ｚ一∥ｃ盯：）瓦＋∥；盯：一∥：口；一２盯：盯；ｌｎｌ［矿ｃｒ。ＣＰ尸（（ｃ。。ｏ。））Ｊ］＝。（７）§３．２双阈值的ＥＭ算法
ＳＡＲ图像各像素值表示雷达回波（后向散射）的强弱，两幅不同时间的ＳＡＰ，图像的差值可表明各像素点上雁向散射增强、不变和减弱三种情况。增强与减弱表现了地表变化的不同信息。我们将差值图像各点分为三类：增强出。、不变∞。和减弱ｑ：。在ｚ。中即分别为正的点集、０周围的点集和负的点集。
于是勘的概率密度函数ｐ０Ｄ为
ｐ（Ｘ）＝ｐ（Ｘｌ珊ｄ）Ｐ（ｃｏ＿）＋ｐ（Ｘｉ∞。）２（６００）＋ｐ（Ｘｌ甜。２）ｐ（ｃｏ。２）（８）
采用ＥＭ算法，首先对增强类吼．和非增强类峨，（０９。和０３，：）得到￡，，再对
减弱类甜。２和非减弱类蛾２（国。和ＣＯ。１）得到Ｌ２，即∞。１＝卯。ＬＪｃｏ∥％２＝∞。Ｕ峨１。
将这两个结果合并得至Ｕ－－－种分类。
Ｊｌｎ¨Ｊ
不采样的象素
不采样的象素／／＼
。。—
Ｖ＼＼一／—＼
厂＼／
疋？Ｌｊ０Ｔｑ｝Ｔｃ｝蜊
氍
＊
匿
型
拭
（ａ）
差值图像量化后的灰度值
（ｂ）
图１差值图像直方图
（ａ）在差值图像直方图上的阈值划分示意图，（ｂ）在实际差值图像的归一化直方图上的三类分布示意图
如图１（ａ），直方图ｈ（ｘ）中堆＋大区（Ｘ≥０）内用两个阈值Ｌ：和ｔ，来获得一个比较接近峨，的像素子集ｓ。和另一个比较接近∞。。像素子集ｓ，将Ｌ．和ｔ，写成Ｌ．＝Ｍ。，（１一ｙ）和Ｌ．＝Ｍ。．（１＋ｙ），其中Ｍ目是ｈ（ｘ）中增强区的中值（即Ｍ。＝［ｍａｘ｛ＸＤ）／２】）。同样，在ｈ（ｘ）中减弱区（Ｚ＜０）内找到比较接近ＯＪ。，和６０。２的像素子集ｓ。２和ｓ∥其中一汀ｄ＝【ｍｉｎ｛ＸＤ｝／２］。本文中，取Ｏ．５。
图１（ｂ）为延中绿地区域实际数据后向散射差值图像的灰度直方图。为了便于比较分类，在图中用曲线大致标出了三类的概率密度分布。直方图也进行归一化，使得直方图可以与概率密度分布相一致。
在增强区内迭代计算时，不管减弱区的像素：在减小区内迭代时，不管增强区的像素。那么得到集合Ｓ。，＝｛ｚ（ｉ，Ｊ）｝０＜Ｘ（ｉ√）＜Ｌ，），Ｓ。?＝｛ｘ（ｉ，Ｊ）ＩＸ（ｉ，Ｊ）＞Ｌ｝和Ｓ。：＝｛ｘ（ｉ，Ｊ）ＩＬ：＜ｘ（ｉ，Ｊ）＜０），Ｓ。２＝（ｘ（ｉ，Ｊ）ｌｘ（ｉ，Ｊ）＜ｔ２），用来计算两类０９∥０９。ｌ和国。２，０９。：各自相关的初始统计参数。
确定了初始参数后，用式（２“）进行上述四类的迭代，收敛后得到四类的先验概率密度函数及其统计参数［Ｐ（ｅｏ。，），∥。，盯ｊ］，［Ｐ（ｃｏ。），段，，盯：］，［Ｐ（ｏａ。。），以：，仃。２：］和［户（吐：），∥。盯己］。再先后求解式（７），分别得到阈值Ｌ，和闽值ｒ，。
§３．３Ｂａｒｋｏｖ随机场（ＭＲＦ）对空间结构相关的变化分类
Ｂｒｕｚｚｏｎｅ和Ｐｒｉｅｔ００１使用了两种方法来检测多时相遥感图象上的区域变化，第一种假设差值图像各像素之间是统计独立的，第二种则认为各像素是空间结构相关的。事实上，一个像素属于６）ｋ类，那么它相邻的像素也很有可能属于功。类。使用第二种方法能改善分类的性能，得到更可靠的结果【３】。
用集合Ｃ＝｛ｃ，，１≤ｆ＜三｝表示差值图像上所有分类结果的集合，ｃｆ２｛Ｃ∥，，），１≤沤』，１≤ＪｓＪｊ，Ｃ，（２，Ｊ）∈｛０９ｃ１，ｑ：，０９。ｊ，￡：３“，，．，分别为图像的高度和宽度。
根据Ｂａｙｅｓ最小误差理论，分类结果应使后验概率为最大，即
Ｃ＾＝ａｒｇｍａｘ｛Ｐ（Ｃ，ｆＸＤ）ｊ
ＬＪｅＣ
＝ａｒｇｍａｘ｛Ｐ（Ｃ，）ｐ（石Ｄｌｃ『））（９）
Ｃ“
其中Ｚ。为差值图像，Ｐ（Ｃ，）是类别标注的先验概率，ｐ（ＸｆＣ，）］是…Ｚ４＇／ｘ……Ｉ…－。：注条件下的差值图像的像素值的联合概率密度函数。式（９）的求解需要同时估计Ｊｐ（ｃ，）和ｐ（Ｘ。ｆＣ，），计算上会过于复杂。采用Ｍａｒｋｏｖ随机场（ＭＲＦ）来建立差值