一种用于参考条带模式InSAR的前斜视测高方法

文档序号:27383879发布日期:2021-11-15 20:51阅读:146来源:国知局
一种用于参考条带模式InSAR的前斜视测高方法
一种用于参考条带模式insar的前斜视测高方法
技术领域
1.本发明涉及干涉合成孔径雷达(insar)测高技术领域,具体涉及一种用于参考条带模式insar的前斜视测高方法。


背景技术:

2.干涉合成孔径雷达(interference synthetic aperture radar,insar)测高是指对同一区域的两幅sar图像进行干涉处理,从干涉相位中提取地面高程信息。目前,insar技术已广泛应用于机载、星载sar系统,其应用领域涵盖地标形变探测、动目标检测、海洋测绘、森林制图、洪涝监测、交通检测和冰川研究等。用于测高的传统insar系统主要工作在正侧视模式下,用于干涉的两幅sar图像中各像素点具有零多普勒中心频率或较小的多普勒中心频率;在此情况下,亚像素级的图像配准精度足以保证干涉相位精度和图像相干性。当insar系统工作在大斜视模式下时,受到多普勒中心频率影响,亚像素级的配准精度不再能够保证干涉相位精度。
3.参考条带模式sar是指通过连续地沿俯仰向调整天线波束指向,以获取不与sar平台运动轨迹平行的观测条带区域的一种新型sar成像模式。该模式与传统的条带模式相比,可以更加灵活的设置观测区域,一定程度上摆脱平台运动轨迹对观测区域的限制。为减小数据量,便于设计时变的脉冲重复间隔以应对大距离徙动的影响,参考条带模式sar需要保持波束指向垂直于条带走向,因而一般工作在大斜视模式下。
4.参考条带模式insar是指将干涉测高技术引入参考条带模式sar中,通过单航过或重复航过获取同一区域的两幅sar图像以完成干涉测高处理;为使参考条带模式insar具有干涉测高的能力,需要研究并设计一种可以用于大斜视模式的insar测高方法。


技术实现要素:

5.为解决参考条带模式insar在大斜视模式下的测高问题,本发明提供一种用于参考条带模式insar的前斜视测高方法。
6.本发明公开了一种用于参考条带模式insar的前斜视测高方法,包括:
7.步骤1、建立三维直角坐标系,读取参考条带模式insar双通道回波数据、雷达系统参数及运动参数、地面控制点参数;
8.步骤2、基于所述雷达系统参数及运动参数,采用改进的bp算法在所述三维直角坐标系中的地平面对回波数据进行成像处理;
9.步骤3、将成像处理得到的两幅sar图像进行干涉处理,并使用非线性相位滤波方法对干涉相位进行滤波;
10.步骤4、采用最小费用流方法对干涉相位进行相位解缠绕,并基于地面控制点确定绝对相位;
11.步骤5、采用前斜视模式sar的高程反演公式计算各像素点高程;
12.步骤6、采用散点插值方法,对采样网格点上的数据进行插值重采样,得到几何校
正后地面数字高程模型。
13.作为本发明的进一步改进,所述步骤1,具体包括:
14.步骤101、沿平行和垂直于飞行速度的方向及天向建立三维直角坐标系;其中,所述三维直角坐标系由满足右手定则的三个坐标轴x、y、z组成,x轴指向飞行方向,y轴垂直于地平面指向上方,z轴由右手定则确定;
15.步骤102、读取参考条带模式insar双通道回波数据、雷达系统参数及运动参数、地面控制点坐标;
16.所述参考条带模式insar双通道回波数据分别表示为e
m
和e
s
,分别对应insar系统主天线自发自收的回波信号矩阵和主天线发射辅天线接收的回波信号矩阵;e
m
和e
s
的都为n
a
×
n
r
的复矩阵,其中n
a
为方位向回波采样条数,n
r
为距离向单条回波采样的点数;
17.所述雷达系统参数及运动参数包括雷达的工作波长λ、雷达的采样率f
s
、发射的线性调频信号的基带数字波形s(n)、光速c、各方位采样时刻的回波窗开启距离r
min
(n
a
)、主天线相位中心在回波录取过程中的高度h及其空间位置坐标辅天线相位中心在回波录取过程中的空间位置坐标从主天线相位中心指向辅天线相位中心的基线矢量雷达波束矢量在场景系下的水平方位角θ
y
、雷达波束矢量在场景系下的下视角θ
p
(n
a
)、雷达天线方向图方位向3db波束宽度θ
a3db
、雷达天线方向图俯仰向3db波束宽度θ
r3db
、期望成像区域的零高度地面网格点坐标集合g(x,z)={(x
j
,0,z
j
)
t
},其中n∈{1,2,

,l
lfm
}为基带数字波形的索引,数字波形长度设为l
lfm
,n
a
∈{1,2,

,n
a
}为方位采样索引;
18.地面控制点坐标表示为
19.作为本发明的进一步改进,所述步骤2,具体包括:
20.步骤201、脉冲压缩;
21.通过在矩阵右端补零将回波信号矩阵e
m
扩充为n
a
×
(n
r
+l
lfm
)的矩阵e
mp
,通过基带数字波形s(n)的右端补零将其扩充为1
×
(n
r
+l
lfm
)的向量s
p
(i),其中i∈{1,2,

,n
r
+l
lfm
};对s
p
(i)进行快速傅里叶变换并取共轭得到匹配滤波器h
lfm
;依据如下公式逐行对e
mp
进行处理得到脉压后n
a
×
(n
r
+l
lfm
)的回波矩阵e
mcp

22.e
mcp
(i,:)=ifft[fft[e
mp
(i,:)]*h
lfm
]
ꢀꢀꢀ
(1)
[0023]
其中fft[
·
]与ifft[
·
]分别表示快速傅里叶变换和快速傅里叶反变换,*表示矩阵的hadamard积,(i,:)表示矩阵的第i行;
[0024]
步骤202、初始化成像网格;
[0025]
构造零初值的图像网格i(x,z)={p(x
j
,z
j
)}与地面网格点坐标集合g(x,z)一一对应;
[0026]
步骤203、对第i行回波进行逆投影处理;
[0027]
根据如下公式计算发射天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0028]
[0029]
根据如下公式计算天线系下发射天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0030][0031][0032][0033]
根据如下公式计算在天线系下的方位角θ
a
(i,j)和俯仰角θ
r
(i,j);
[0034][0035][0036]
其中表示的第k个元素;
[0037]
不同时满足以下条件的地面网格点的投影值v(i,j)=0,对其余网格点进行后续步骤的处理;
[0038][0039][0040]
若对主天线回波进行成像,则根据如下公式计算接收天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0041][0042]
若对辅天线回波进行成像,则根据如下公式计算接收天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0043][0044]
根据如下的公式计算地面网格点在回波中的索引idx(i,j);
[0045][0046]
其中||
·
||表示向量的二范数;
[0047]
不满足0≤idx(i,j)<n
r
的网格点的投影值v(i,j)=0,对其余网格点进行后续步骤的处理;
[0048]
根据如下公式计算第i个脉冲在成像网格点(x
j
,0,z
j
)
t
的投影值v(i,j);
[0049][0050]
其中表示向下取整,m为插值核点数,取正偶数;
[0051]
步骤204、对全部方位采样时刻的回波进行步骤203;
[0052]
步骤205、求和得到单视复图像;
[0053]
根据如下公式求和计算图像网格i(x,z)={p(x
j
,z
j
)};
[0054]
p(x
j
,z
j
)=∑
i
v(i,j)
ꢀꢀꢀ
(14)
[0055]
若对主天线成像则将成像结果表示为i
m
(x,z)={p
m
(x
j
,z
j
)},若对辅天线成像,则将成像结果表示为i
s
(x,z)={p
s
(x
j
,z
j
)};
[0056]
步骤206、对辅天线回波矩阵e
m
重复步骤201到步骤204的操作。
[0057]
作为本发明的进一步改进,所述步骤3,具体包括:
[0058]
步骤301、干涉处理;
[0059]
根据如下公式得到原始干涉相位;
[0060]
φ(x,z)=∠(i
m
(x,z)*conj(i
s
(x,z)))
ꢀꢀꢀ
(15)
[0061]
其中∠(
·
)表示取相位操作,conj(
·
)表示取共轭操作;
[0062]
步骤302、非线性相位滤波;
[0063]
根据王青松在ieee geoscience and remote sensing letters期刊2011年11月第8卷第6期中的论文“an efficient and adaptive approach for noise filtering of sar interferometric phase images”中提出的非线性相位滤波方法完成相位滤波操作;
[0064]
记滤波后相位为φ
filt
(x,z)。
[0065]
作为本发明的进一步改进,所述步骤4,具体包括:
[0066]
步骤401、最小费用流解相位缠绕;
[0067]
根据maori costantini在transactions on geoscience and remote sensing期刊1998年5月第36卷第3期中的论文“a novel phase unwrapping method based on network programming”中提出的最小费用最大流相位解缠绕方法完成相位解缠绕操作;
[0068]
记解缠绕后相位为φ
unwrap
(x,z);
[0069]
步骤402、计算地面控制点在主图像中聚焦位置;
[0070]
根据如下公式计算地面控制点在主图像中的聚焦位置(x
ctrl
,z
ctrl
);
[0071]
x
ctrl
=p
ctrl_x
ꢀꢀꢀ
(16)
[0072][0073]
其中,q∈{0,1},当雷达右视时q=0,当雷达左视时q=1;
[0074]
步骤403、计算地面控制点的绝对干涉相位;
[0075]
根据如下公式计算控制点的合成孔径中心时刻主天线位置根据如下公式计算控制点的合成孔径中心时刻主天线位置
[0076][0077]
根据如下公式计算控制点绝对干涉相位;
[0078][0079]
步骤404、计算解缠绕后的绝对干涉相位;
[0080]
记中离(x
ctrl
,z
ctrl
)最近的点为
[0081]
根据如下公式计算解缠绕后的绝对干涉相位
[0082][0083][0084]
其中round(
·
)表示四舍五入。
[0085]
作为本发明的进一步改进,所述步骤5,具体包括:
[0086]
步骤501、计算各地面像素网格对应的孔径中心时刻主天线相位中心位置;
[0087]
根据如下公式计算各地面像素网格g(x,z)对应的孔径中心时刻主天线相位中心位置
[0088][0089]
步骤502、计算各像素网格高程值;
[0090]
根据如下公式计算各像素网格坐标(x
j
,0,z
j
)
t
对应的高程值h(x
j
,z
j
)
[0091][0092]
作为本发明的进一步改进,所述步骤6,具体包括:
[0093]
步骤601、计算各像素网格内散射点的空间坐标;
[0094]
根据如下公式计算各像素网格点的空间坐标(x
3j
,y
3j
,z
3j
)
t

[0095]
x
3j
=x
j
ꢀꢀꢀ
(22)
[0096]
y
3j
=h(x
j
,z
j
)
ꢀꢀꢀ
(23)
[0097][0098]
其中,q∈{0,1},当雷达右视时q=0,当雷达左视时q=1;
[0099]
步骤602、散点插值得到几何校正后成像区域数字高程模型;
[0100]
根据amidror在journal of electronic imaging期刊2002年4月第11卷第2期中的论文“scattered data interpolation methods for electronic imaging systems:a survey”中提出的方法,在上一步获取的空间点云{(x
3j
,y
3j
,z
3j
)
t
}上对成像网格坐标集合g(x,z)={(x
j
,0,z
j
)
t
}中空间坐标的第二维进行散点插值;得到期望成像区域的几何校正后的数字高程模型dem={(x
j
,y
j
,z
j
)
t
}。
[0101]
与现有技术相比,本发明的有益效果为:
[0102]
本发明通过对两幅sar图像进行干涉处理、相位滤波、相位解缠、高程反演、几何校正,得到目标区域的地面数字高程模型,实现大斜视条件下的insar测高处理;其避免了传统insar测高处理方法在处理斜视数据时对主辅图像配准精度要求极高、难以保证干涉相位精度的问题,避免了传统insar测高处理方法在斜视模式下需要求解多普勒方程的问题,能够较为快速、准确的解决参考条带模式insar的数字高程模型获取问题。
附图说明
[0103]
图1为本发明一种实施例公开的用于参考条带模式insar的前斜视测高方法的流程图;
[0104]
图2为本发明一种实施例公开的雷达发射信号数字波形图;
[0105]
图3为本发明一种实施例公开的雷达回波窗开启距离变化图;
[0106]
图4为本发明一种实施例公开的雷达波束指向下视角变化图;
[0107]
图5为本发明一种实施例公开的雷达主天线bp成像结果幅度图;
[0108]
图6为本发明一种实施例公开的雷达辅天线bp成像结果幅度图;
[0109]
图7为本发明一种实施例公开的原始干涉相位图;
[0110]
图8为本发明一种实施例公开的滤波后干涉相位图;
[0111]
图9为本发明一种实施例公开的解缠绕后绝对干涉相位图;
[0112]
图10为本发明一种实施例公开的几何校正前高程图;
[0113]
图11为本发明一种实施例公开的数字高程模型图;
[0114]
图12为本发明一种实施例公开的数字高程模型三维展示图。
具体实施方式
[0115]
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
[0116]
下面结合附图对本发明做进一步的详细描述:
[0117]
如图1所示,本发明提供一种用于参考条带模式insar的前斜视测高方法,包括:
[0118]
步骤1、沿平行和垂直于飞行速度的方向及天向建立三维直角坐标系,读取参考条带模式insar的双通道回波数据,读取雷达系统参数及运动轨迹,读取地面控制点参数;
[0119]
具体包括:
[0120]
步骤101、沿平行和垂直于飞行速度的方向及天向建立三维直角坐标系;
[0121]
三维直角坐标系由满足右手定则的三个坐标轴x、y、z组成,其中x轴指向飞行方向,y轴垂直于地平面指向上方,z轴由右手定则确定,涉及到三维空间坐标的数据全部定义在该坐标系下;
[0122]
步骤102、读取参考条带模式insar双通道回波数据、读取雷达系统参数及运动参数、读取地面控制点三维坐标;
[0123]
数据、参数、坐标在本实施例中具体为:参考条带模式insar双通道回波数据分别
表示为e
m
和e
s
,分别对应insar系统主天线自发自收的回波信号矩阵和主天线发射辅天线接收的回波信号矩阵。e
m
和e
s
的都为n
a
×
n
r
的复矩阵,其中n
a
=15334为方位向回波采样条数,n
r
=3201为距离向单条回波采样的点数。雷达系统参数包括雷达的工作波长λ=0.02m、雷达的采样率f
s
=160mhz、发射的线性调频信号的基带数字波形s(n)如图2所示、光速c=3e8m/s、各方位采样时刻的回波窗开启距离r
min
(n
a
)如图3所示、主天线相位中心在回波录取过程中的高度h=10km及其空间位置坐标辅天线相位中心在回波录取过程中的空间位置坐标从主天线相位中心指向辅天线相位中心的基线矢量雷达波束矢量在场景系下的水平方位角θ
y
=45deg、雷达波束矢量在场景系下的下视角θ
p
(n
a
)如图4所示、雷达天线方向图方位向3db波束宽度θ
a3db
=1deg、雷达天线方向图俯仰向3db波束宽度θ
r3db
=2deg、期望成像区域的零高度地面网格点坐标集合g(x,z)={(x
j
,0,z
j
)
t
}为以坐标(10606.601718,0,10606.601718)
t
为中心的1000
×
1000的0.1m间隔网格,其中n∈{1,2,

,l
lfm
}为基带数字波形的索引,数字波形长度设为l
lfm
=1601,n
a
∈{1,2,

,n
a
}为方位采样索引;地面控制点坐标表示为
[0124]
步骤2、基于雷达系统参数和运动轨迹,采用改进的bp算法在所建三维直角坐标系中的地平面完成成像处理;
[0125]
具体包括:
[0126]
步骤201、脉冲压缩;
[0127]
通过在矩阵右端补零将回波信号矩阵e
m
扩充为n
a
×
(n
r
+l
lfm
)的矩阵e
mp
,通过基带数字波形s(n)的右端补零将其扩充为1
×
(n
r
+l
lfm
)的向量s
p
(i),其中i∈{1,2,

,n
r
+l
lfm
};对s
p
(i)进行快速傅里叶变换并取共轭得到匹配滤波器h
lfm
;依据如下公式逐行对e
mp
进行处理得到脉压后n
a
×
(n
r
+l
lfm
)的回波矩阵e
mcp

[0128]
e
mcp
(i,:)=ifft[fft[e
mp
(i,:)]*h
lfm
]
ꢀꢀꢀ
(25)
[0129]
其中fft[
·
]与ifft[
·
]分别表示快速傅里叶变换和快速傅里叶反变换,*表示矩阵的hadamard积,(i,:)表示矩阵的第i行。
[0130]
步骤202、初始化成像网格;
[0131]
构造零初值的图像网格i(x,z)={p(x
j
,z
j
)}与地面网格点坐标集合g(x,z)一一对应;
[0132]
步骤203、对第i行回波进行逆投影处理;
[0133]
根据如下公式计算发射天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0134][0135]
根据如下公式计算天线系下发射天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0136][0137][0138][0139]
根据如下公式计算在天线系下的方位角θ
a
(i,j)和俯仰角θ
r
(i,j);
[0140][0141][0142]
其中表示的第k个元素;
[0143]
不同时满足以下条件的地面网格点的投影值v(i,j)=0,对其余网格点进行后续步骤的处理;
[0144]

a
(i,j)|≤0.5deg
ꢀꢀꢀ
(32)
[0145]

r
(i,j)|≤1deg
ꢀꢀꢀ
(33)
[0146]
若对主天线回波进行成像,则根据如下公式计算接收天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0147][0148]
若对辅天线回波进行成像,则根据如下公式计算接收天线相位中心到成像网格点(x
j
,0,z
j
)
t
的矢量表示
[0149][0150]
根据如下的公式计算地面网格点在回波中的索引idx(i,j);
[0151][0152]
其中||
·
||表示向量的二范数;
[0153]
不满足0≤idx(i,j)<n
r
的网格点的投影值v(i,j)=0,对其余网格点进行后续步骤的处理;
[0154]
根据如下公式计算第i个脉冲在成像网格点(x
j
,0,z
j
)
t
的投影值v(i,j);
[0155][0156]
其中表示向下取整,m=16为插值核点数;
[0157]
步骤204、对全部方位采样时刻的回波进行步骤203;
[0158]
步骤205、求和得到单视复图像;
[0159]
根据如下公式求和计算图像网格i(x,z)={p(x
j
,z
j
)};
[0160]
p(x
j
,z
j
)=∑
i
v(i,j)
ꢀꢀꢀ
(38)
[0161]
若对主天线成像则将成像结果表示为i
m
(x,z)={p
m
(x
j
,z
j
)},若对辅天线成像,则将成像结果表示为i
s
(x,z)={p
s
(x
j
,z
j
)};
[0162]
步骤206、对辅天线回波矩阵e
m
重复步骤201到步骤204的操作;
[0163]
经过步骤201到步骤206处理的主辅天线复图像幅值如图5、图6所示。
[0164]
步骤3、将成像处理得到的两幅sar图像进行干涉处理,并使用非线性相位滤波方法对干涉相位进行滤波;
[0165]
具体包括:
[0166]
步骤301、干涉处理;
[0167]
根据如下公式得到原始干涉相位如图7所示;
[0168]
φ(x,z)=∠(i
m
(x,z)*conj(i
s
(x,z)))
ꢀꢀꢀ
(39)
[0169]
其中∠(
·
)表示取相位操作,conj(
·
)表示取共轭操作;
[0170]
步骤302、非线性相位滤波
[0171]
根据王青松在ieee geoscience and remote sensing letters期刊2011年11月第8卷第6期中的论文“an efficient and adaptive approach for noise filtering of sar interferometric phase images”中提出的非线性相位滤波方法完成相位滤波操作;
[0172]
记滤波后相位为φ
filt
(x,z),如图8所示;
[0173]
步骤4、采用最小费用流方法对干涉相位进行相位解缠绕,并基于地面控制点确定绝对相位;
[0174]
具体包括:
[0175]
步骤401、最小费用流解相位缠绕;
[0176]
根据maori costantini在transactions on geoscience and remote sensing期刊1998年5月第36卷第3期中的论文“a novel phase unwrapping method based on network programming”中提出的最小费用最大流相位解缠绕方法完成相位解缠绕操作;
[0177]
记解缠绕后相位为φ
unwrap
(x,z),如图9所示;
[0178]
步骤402、计算地面控制点在主图像中聚焦位置;
[0179]
根据如下公式计算地面控制点在主图像中的聚焦位置(x
ctrl
,z
ctrl
);
[0180]
x
ctrl
=p
ctrl_x
=10606.5518m
ꢀꢀꢀ
(40)
[0181][0182]
其中,q∈{0,1},当雷达右视时q=0,当雷达左视时q=1;
[0183]
步骤403、计算地面控制点的绝对干涉相位;
[0184]
根据如下公式计算控制点的合成孔径中心时刻主天线位置根据如下公式计算控制点的合成孔径中心时刻主天线位置
[0185][0186]
根据如下公式计算控制点绝对干涉相位;
[0187][0188]
步骤404、计算解缠绕后的绝对干涉相位;
[0189]
记中离(x
ctrl
,z
ctrl
)最近的点为
[0190]
根据如下公式计算解缠绕后的绝对干涉相位
[0191][0192][0193]
其中round(
·
)表示四舍五入;
[0194]
步骤5、采用前斜视模式sar的高程反演公式计算各像素点高程;
[0195]
具体包括:
[0196]
步骤501、计算各地面像素网格对应的孔径中心时刻主天线相位中心位置;
[0197]
根据如下公式计算各地面像素网格g(x,z)对应的孔径中心时刻主天线相位中心位置
[0198][0199]
步骤502、计算各像素网格高程值;
[0200]
根据如下公式计算各像素网格坐标(x
j
,0,z
j
)
t
对应的高程值h(x
j
,z
j
)如图10所示;
[0201][0202]
步骤6、采用散点插值方法,对采样网格点上的数据进行插值重采样,得到几何校正后地面数字高程模型;
[0203]
具体包括:
[0204]
步骤601、计算各像素网格内散射点的空间坐标;
[0205]
根据如下公式计算各像素网格点的空间坐标(x
3j
,y
3j
,z
3j
)
t

[0206]
x
3j
=x
j
ꢀꢀꢀ
(48)
[0207]
y
3j
=h(x
j
,z
j
)
ꢀꢀꢀ
(49)
[0208][0209]
其中,q∈{0,1},当雷达右视时q=0,当雷达左视时q=1;
[0210]
步骤602、散点插值得到几何校正后成像区域数字高程模型;
[0211]
根据amidror在journal of electronic imaging期刊2002年4月第11卷第2期中的论文“scattered data interpolation methods for electronic imaging systems:a survey”中提出的方法,在上一步获取的空间点云{(x
3j
,y
3j
,z
3j
)
t
}上对成像网格坐标集合g(x,z)={(x
j
,0,z
j
)
t
}中空间坐标的第二维进行散点插值;得到期望成像区域的几何校正后的数字高程模型dem={(x
j
,y
j
,z
j
)
t
},如图11与图12所示。
[0212]
本发明的优点为:
[0213]
本发明通过对两幅sar图像进行干涉处理、相位滤波、相位解缠、高程反演、几何校正,得到目标区域的地面数字高程模型,实现大斜视条件下的insar测高处理;其避免了传统insar测高处理方法在处理斜视数据时对主辅图像配准精度要求极高、难以保证干涉相位精度的问题,避免了传统insar测高处理方法在斜视模式下需要求解多普勒方程的问题,能够较为快速、准确的解决参考条带模式insar的数字高程模型获取问题。
[0214]
以上仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1