-
文中研究的成像光谱仪,通过视场分离形成VNIR和SWIR两通道,沿轨方向存在的视场分离角,造成同一时刻两通道成像地物不同。SWIR通道受探测器限制,采用双狭缝成像,每个狭缝对应2个成像探测器,4个探测器交错排列,对应4个子视场:SWIR1-1、SWIR2-1、SWIR1-2和SWIR2-2,子视场间相互重叠拼接形成宽幅视场,如图1所示。SWIR1-1和SWIR1-2通过狭缝S1成像,相对VNIR的视场分离角相同,SWIR2-1和SWIR2-2通过狭缝S2成像,视场分离角相同。经过一段时间推扫,得到的1个VNIR和4个SWIR图像,在沿轨方向上存在失配。
运动补偿原理如图2所示,卫星从t1时刻开始运动补偿,补偿镜使观测光线相对星下点A夹角为θn,指向A1;控制补偿镜转动,t2时观测光线相对星下点B夹角-θn,指向B1;t1到t2完成一次运动补偿,t3时刻开始下一次运动补偿,补偿倍率为n时,|AB|=n*|A1B1|。推扫成像时补偿角不断变化,仪器的运动补偿镜匀角速度转动,起始时刻运动补偿角为:
$${\theta _n} = \arcsin \left( {\frac{{R\sin \alpha }}{{\sqrt {{{\left( {R + H} \right)}^2} + {R^2} - 2R\left( {R + H} \right)\cos \alpha } }}} \right)$$ (1) 式中:α=(n-1)|A1B1|/2R;R为地球半径;H为轨道高度。运动补偿时,图像在连续变化的观测角下获取,但积分时间内角度变化很小,可将运动补偿角等效为积分时间中点的倾斜观测角。第i推扫行运动补偿角为:
$$ {\theta }_{i}={\theta }_{n}-\left(2i-1\right) \cdot {\theta }_{n} \cdot \tau /\Delta T$$ (2) 式中:ΔT为成像时长;τ为积分时间。
推扫成像时各行的成像关系相互独立,以VNIR为基准建立第i行像空间坐标系Fi,原点为投影中心Si,z轴由像面中心指向Si,y轴平行于推扫行,x轴满足右手法则。像元p(i,j)的观测光线与VNIR中心光线夹角:
$$\left\{ \begin{array}{l} {\theta _{xi}} = \Delta x + {\theta _i}\\ {\theta _{yi}} = \arctan \left( {\left( {j - \Delta y - \dfrac{{{N_V} - 1}}{2}} \right) \cdot \dfrac{{\Delta d}}{f}} \right) \end{array} \right.$$ (3) 式中:Δx为视场分离角;Δy为SWIR子视场中心相对VNIR中心y方向的偏移量;NV为VNIR探测器像元个数;Δd为探测器尺寸;f为系统焦距。运动补偿下,p(i, j)的像空间观测矢量为:
$$\overrightarrow {{u_{Fi}}} = {\left[ {\begin{array}{*{20}{c}} {\tan {\theta _{xi}}}&{\tan {\theta _{yi}}}&{ - 1} \end{array}} \right]^{\rm{T}}}$$ (4) ${\vec u_{Fi}}$ 通过坐标转换得到地面坐标系观测矢量${\vec u_{Ei}}$ ,与国家标准的数字高程模型(DEM)产品求交点,得到对应地面点经纬度坐标P (Bij, Lij, Hij)[9]。推扫成像时间内,依次计算各通道的每一行所有像元的地面点坐标,建立了各通道各子视场的严格成像几何模型。 -
推扫成像时各行图像是在不同时刻得到的,运动补偿时各行的成像观测角度不同,从而导致成像过程中不同行对应的地面分辨率不同。
令轨道高度708 km,倾角98.217°,VNIR、SWIR探测器像元数:NV=2 000,NS=512,瞬时视场角42.5 μrad,视场分离角与运动补偿角如表1所示,根据VNIR和SWIR各子视场成像几何模型,确定像元p(i, j)与地面点P(Bi, j, Li, j, Hi, j)对应关系。利用相邻像元在沿轨方向x、交轨方向y上对应的地面点距离,表征地面分辨率:
$$\Delta {r_x}\left( {i,j} \right) = \frac{{{L_{i + 1,j}} - {L_{i,j}}}}{{\cos \left( {{\alpha _{i,j}} - \dfrac{\pi }{2}} \right)}},\Delta {r_y}\left( {i,j} \right) = \dfrac{{{B_{i,j + 1}} - {B_{i,j}}}}{{\cos \left( {{\beta _{i,j}} - \dfrac{\pi }{2}} \right)}}$$ (5) 式中:αi, j为沿轨方向与经度方向夹角;βi, j为交轨方向与纬度方向夹角。
表 1 视场分离角与运动补偿角
Table 1. Field separation and motion compensation angles
Field separation angle Motion compensation angle θn SWIR1 SWIR2 n=2 n=4 n=6 0.5° 1.0° 2.426° 7.238° 11.932° 运动补偿主要影响沿轨方向的地面分辨率(Ground Sampling Distance,GSD),不同补偿倍率下,VNIR和SWIR图像各行中心像元在x方向地面分辨率如图3所示。运动补偿下,VNIR和SWIR图像各行像元的地面分辨率发生非线性变化,补偿倍率越大地面分辨率变化越大。相同运动补偿倍率下,由于存在视场分离角,VNIR和SWIR的同一行图像对应地面分辨率不同。
图 3 不同行中心像元地面分辨率曲线(n为运动补偿倍率)
Figure 3. Ground sampling distance of middle pixel on different lines (n is motion compensation ratio)
VNIR通道中心行运动补偿角为0,即垂直观测;SWIR通道由于视场分离角的影响,垂直观测的推扫行与视场分离角有关。以垂直观测的图像行作为参考行m,运动补偿下沿轨方向像元到参考行像元的地面点坐标距离,与无畸变时地面距离的差值作为畸变量:
$$ {\delta }_{x}\left(i,j\right)=\frac{{L}_{i,j}-\left[{L}_{m,j}+\left(i-m\right) \cdot \Delta {r}_{L}\left(m,j\right)\right]}{\Delta {r}_{L}\left(m,j\right)}$$ (7) 式中:
$ \Delta {r_L}\left( {m,j} \right){\rm{ = }}{L_{m + 1,\;j}} - {L_{m,\;j}} $ 。VNIR和SWIR子视场图像各行中心像元沿x方向的畸变量曲线如图4所示。运动补偿下,VNIR和SWIR图像各行在沿轨方向存在非线性畸变,且运动补偿倍率越大,推扫成像行数越多,畸变越大。6倍运动补偿,成像2 000行时,VNIR最大畸变18像元,SWIR2最大畸变达28像元。同一运动补偿倍率下,对比同一时刻的VNIR和SWIR同一行图像,两者畸变大小不同,即推扫得到的VNIR和SWIR图像畸变不一致。
-
利用SWIR像元相对VNIR同一像元的地面点经纬度坐标距离,表示SWIR通道该像元的失配量。VNIR像元p(i, j)对应地面点P(Bi, j, Li, j, Hi, j),SWIR像元p′(i, j)对应地面点P′(B′i, j, L′i, j, H′i, j),则SWIR像元相对VNIR像元的失配量为:
$$\begin{split} &\left[ {\begin{array}{*{20}{c}} {\Delta {p_x}\left( {i,j} \right)} \\ {\Delta {p_y}\left( {i,j} \right)} \end{array}} \right] =\\ &\frac{1}{{\cos \left( {{\beta _{i,j}} - {\alpha _{i,j}}} \right)}}\left[ {\begin{array}{*{20}{c}} {\cos {\beta _{i,j}}}&{\sin {\beta _{i,j}}} \\ { - \sin {\alpha _{i,j}}}&{\cos {\alpha _{i,j}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dfrac{{{B_{i,j}}^\prime - {B_{i,j}}}}{{\Delta {r_x}\left( {i,j} \right)}}} \\ {\dfrac{{{L_{i,j}}^\prime - {L_{i,j}}}}{{\Delta {r_y}\left( {i,j} \right)}}} \end{array}} \right] \end{split}$$ (8) SWIR两个子视场图像沿x方向各行中心像元的失配量曲线如图5示。对于SWIR图像,沿x方向每行的失配量绝对值先减小后增大,运动补偿倍率越高,成像行数越多,失配量越大,各行的失配量非线性变化。
图 5 SWIR通道失配量变化曲线(n为运动补偿倍率)
Figure 5. Mismatch curve of SWIR channel(n is motion compensation ratio)
对比同一SWIR子视场不同运动补偿倍率下的失配曲线,补偿倍率n越大,失配量变化曲线的曲率越大,失配量的变化范围也越大。2倍运动补偿,成像2 000行时,SWIR2-1失配量变化在1像元内;6倍运动补偿时,相同条件下失配量变化达23像元。对比同一运动补偿倍率下,SWIR两个子视场图像的失配曲线,在视场分离角影响下,起始行和终止行失配量不对称。SWIR2-1的初始视场分离角较大,同一时刻两个子视场的同一行图像,SWIR2-1失配量较大,SWIR2-1图像的失配量变化范围也大于SWIR1-1。
综上所述,运动补偿造成的非线性畸变,使SWIR图像不同区域失配量不同,边缘失配大中心小。随运动补偿倍率、视场分离角和成像行数增大,图像的失配量变化范围也迅速增大。设计运动补偿倍率和视场分离角时,考虑信噪比的同时应结合两者对图像畸变和失配的影响,根据图像的应用需求具体分析。对于运动补偿下SWIR和VNIR图像配准,需要先消除VNIR和SWIR各视场图像的畸变再配准。
-
利用VNIR和SWIR的严格成像几何模型,分别校正图像畸变进行粗配准。对于残留失配量,采用基于图像的方法精配准。由于地物在不同波长下反射特性不同,同一地物在VNIR和SWIR通道成像时灰度特征存在差异。精配准时先预处理得到相似度较高的待配准图像,采用对灰度依赖较小的相位相关法配准。
-
根据几何定标得到的传感器参数、卫星姿轨数据和运动补偿角度等,建立VNIR通道和SWIR通道各子视场的成像几何模型,得到像元p(i, j)对应的地面点坐标P(Bi, j, Li, j, Hi, j)。经纬度坐标(Bi, j, Li, j, Hi, j)投影到地图坐标系,像元坐标p(i, j)和地图投影坐标P(Xi, j, Yi, j)一一对应。以VNIR中心像元的地面分辨率为采样间隔rx、ry,得到地图坐标系下几何校正后图像像元坐标p(I, J):
$$I = \frac{{{X_{i,j}} - {X_{\min }}}}{{{r_x}}},J = - \frac{{{Y_{i,j}} - {Y_{\max }}}}{{{r_y}}}$$ (8) 式中:rx、ry相等;Xmin和Ymax为校正后图像左上点像元的地图投影坐标。
采用双线性插值确定(Xi, j, Yi, j)到(I, J)的插值关系,重采样得到像元p(I, J)灰度值,依次对VNIR和SWIR子视场图像几何校正,四个SWIR子视场图像按地理坐标拼接。
-
(1)图像预处理
剔除吸收峰、边缘波段等低信噪比波段,计算VNIR和SWIR剩余有效波段图像间的灰度互相关值,分别选取相关性好的连续10个波长,进行最小噪声分量变换(MNF),变换后第一分量包含绝大多数特征且信噪比高,选取VNIR和SWIR第一分量作为参考图和待配准图。
(2)像素级失配量估计
VNIR图像fv(x, y)和SWIR图像fs(x-Δx, y-Δy)的离散傅里叶变换为FV和FS,则相位相关函数为[5]:
$$ C(x,y)={F}^{-1}\left(FS \cdot \overline{FV}/\left|FS \cdot \overline{FV}\right|\right)$$ (9) 式中:C(x,y)为脉冲函数,根据脉冲峰值位置可确定SWIR相对VNIR图像的失配量Δx和Δy。
(3)亚像素级失配量估计
传统相位相关法只能确定像元级失配量,在相位相关峰值邻域利用函数拟合可将精度扩展到亚像素,sinc函数是常用的精度较高的拟合函数[5]。将VNIR和SWIR划分成有重叠的图像块fv(xm, ym)和fs(xm−Δxm, ym−Δym),图像块间的相位相关函数为C(xm, ym)。用sinc函数在C(xm, ym)峰值邻域拟合,并在x、y方向上分离变量:
$$\begin{gathered} {C_\eta }({x_m}) = {\alpha _{xm}}\frac{{\sin \pi ({x_m} - \varDelta {x_m})}}{{\pi ({x_m} - \varDelta {x_m})}} + {\gamma _{xm}} \\ {C_\eta }({y_m}) = {\alpha _{ym}}\frac{{\sin \pi ({y_m} - \varDelta {y_m})}}{{\pi ({y_m} - \varDelta {y_m})}} + {\gamma _{ym}} \\ \end{gathered} $$ (10) 式中:αxm、αym为相关系数;γxm、γym为噪声项。通过最小二乘法,得到图像块中心在x、y方向失配量Δxm和Δym。将图像块中心失配量,采用三次样条插值得到每个像元的亚像素失配量。
根据SWIR的像素级失配量移动图像,再利用亚像素级失配量,对SWIR图像空间重采样。
-
利用东天山HyMap高光谱数据,仿真得到运动补偿下双通道高光谱辐亮度数据。为验证运动补偿的影响分析以及配准方案有效性,选择高倍率运动补偿图像:4倍补偿,SWIR1、SWIR2视场分离角0.5°和1°,成像1335行时得到VNIR和SWIR四个子视场图像,如图6所示。
建立VNIR、SWIR几何校正模型,校正图像畸变,将SWIR子视场按地理坐标拼接,如图7所示。选择VNIR波段:band14~band23,SWIR波段:band31~band40,通过MNF变换得到参考图和待配准图。利用相位相关法得到图像整体像素级失配量Δx=−1 pixel,Δy=1 pixle。将VNIR和SWIR划分成100×100图像块,中心间隔50,扩展相位相关法得到各像元的亚像素失配量Δxm:−0.4~0.4 pixel,Δym:−0.3~0.4 pixel。对SWIR图像空间重采样,完成双通道图像配准。
-
为验证配准的几何精度,选取10个特征地物作为控制点,如图7所示。利用控制点11×11邻域的灰度值求质心坐标,有效波长下质心坐标均值作为特征地物几何位置。传统方法直接基于SIFT特征配准图像[10],与文中方案配准结果对比如表2所示。运动补偿下VNIR和SWIR不同区域失配差别大,基于SIFT配准后仍存在平均3.9像元的误差,而文中方案的精度为0.3像元。
表 2 特征地物几何位置误差
Table 2. Geometric errors of feature ground objects
SIFT Proposed method Object coordinate VNIR Object coordinate SWIR Error/pixel Object coordinate VNIR Object coordinate SWIR Error/pixel |x| |y| Distance |x| |y| Distance 1 (896.49,1214.57) (891.30,1213.94) 5.2 0.6 5.2 (764.34,892.07) (764.52,892.44) 0.2 0.4 0.4 2 (807.72,801.63) (803.68,802.30) 4.0 0.8 4.1 (721.48,1308.15) (721.63,1308.43) 0.1 0.3 0.3 3 (683.18,1581.33) (678.90,1582.21) 4.1 0.9 4.2 (1080.13,591.93) (1080.47,592.32) 0.3 0.4 0.5 4 (415.23,98.09) (417.95,96.88) 2.7 1.2 2.9 (964.61,2115.28) (964.95,2114.86) 0.3 0.4 0.5 5 (385.33,592.64) (389.34,592) 4.0 0.6 4.0 (1125.58,1641.73) (1125.42,1641.49) 0.3 0.2 0.4 6 (413.81,1045.86) (416.40,1046.47) 2.6 0.6 2.7 (1217.77,1190.68) (1217.63,1190.51) 0.1 0.2 0.2 7 (230.92,252.77) (234.06.,251.92) 3.4 0.8 3.5 (1194.58,2019.39) (1194.84,2019.07) 0.3 0.3 0.4 8 (309.25,868.61) (312.89,869.13 3.6 0.5 3.7 (1275.79,1390.07) (1275.88,1390.28) 0.1 0.1 0.1 9 (168.55,1439.82) (173.13,1440.84) 4.6 1.0 4.7 (1569.27,872.26) (1569.46,872.17) 0.2 0.1 0.2 10 (93.60,1874.59) (97.04,1875.85) 3.4 1.3 3.6 (1747.45,470.69) (1747.35,470.54) 0.1 0.1 0.1 Average 3.8 0.8 3.9 Average 0.2 0.3 0.3 为验证配准后VNIR-SWIR光谱的一致性,将配准后的VNIR-SWIR辐亮度数据经过FLAASH大气校正,得到地物光谱反射率曲线,通过重叠波段的光谱反射率重合度误差评价:
$$ORE(p) = \left( {\frac{1}{{{N_o}}}\sum\limits_{i = 1}^{{N_o}} {\frac{{{R_S}({\lambda _{Si}},p)}}{{{R_V}({\lambda _{Si}},p)}}} - O{R_0}} \right) \times 100 {\text{%}} $$ (11) 式中:No为SWIR通道重叠波段数目;OR0=1为理想情况下的光谱重合度;RS(λSi, p)为SWIR通道波长λSi处像元p的光谱反射率;RV(λSi, p)为VNIR通道的光谱反射率。VNIR和SWIR通道光谱分辨率不同,分别为4.6 nm和9.7 nm,故VNIR通道的RV(λSi, p)由λSi附近波长插值得到。
分别采用文中方案和SIFT特征配准前后,VNIR图像中特征地物9的光谱反射率曲线如图8(a)和(b)所示,文中方案配准后能够得到特征地物的VNIR-SWIR连续光谱曲线,SIFT特征法配准后VNIR和SWIR光谱曲线不属于同一地物像元,因而不能衔接。图8(c)为两种方法配准后重叠波段的光谱反射率曲线,采用文中方案配准后重叠波段VNIR和SWIR的光谱反射率曲线基本重合,明显优于SIFT特征配准结果。计算10个特征地物的ORE,SIFT特征法重合度误差为41.5%,文中方案重合度误差为1.2%。
Image registration of the dual-channel spaceborne hyperspectral imager with motion compensation
-
摘要: 最新一代可见近红外(VNIR)和短波红外(SWIR)双通道星载高光谱成像仪,多采用视场分离器将VNIR和SWIR通道分离为多个子视场,同一时刻各子视场对地成像区域不同,在采用运动补偿技术提高图像信噪比时,各子视场对同一地物的观测角不同,导致图像间失配关系复杂,无法获取同一地物的VNIR-SWIR连续光谱。通过建立运动补偿下的严格成像几何模型,定量分析了双通道图像的畸变和失配规律,进而提出了各子视场分别几何定位再相位相关法配准的方案,并利用东天山区域运动补偿下星载双通道高光谱仿真数据进行验证。结果表明,传统的基于图像的配准方案精度为3.9像元,仍无法得到同一地物像元的VNIR-SWIR光谱曲线,文中方案配准精度提高到0.3像元,VNIR和SWIR重叠波段的反射率光谱重合度误差由41.5%降低到1.2%。Abstract: Latest generation of dual-channel spaceborne hyperspectral imager based on visible near infrared (VNIR) and short wave infrared (SWIR) uses the field slitter to separate VNIR and SWIR channels into several sub-fields, and each sub-field image has different ground area at the same time. When using motion compensation technology to improve signal to noise ratio of the instrument, the observation angles of each sub-field are different, which leads to more complicated mismatch of images and make it impossible to get the continuous VNIR-SWIR spectrum of the ground pixel. The rule of image distortion and dual-channel mismatch quantitatively was analyzed by Rigorous imaging model, and the registration scheme by using geometric orientation of each sub-field separately as well as the phase correlation method was proposed on this basis. Verification based on dual-channel spaceborne hyperspectral simulated data of Dongtianshan under motion compensation was performed. The result shows that registration accuracy of traditional scheme based on correlation of images reaches 3.9 pixel, which means the continuous VNIR-SWIR spectrum of the ground pixel is still unavailable. The registration accuracy of the scheme proposed by this paper reaches 0.3 pixel, and the reflectance spectrum overlap ratio error of the overlapping bands of VNIR and SWIR reduces from 41.5% to 1.2%.
-
Key words:
- hyperspectral imager /
- motion compensation /
- mismatch analysis /
- image registration /
- resampling
-
表 1 视场分离角与运动补偿角
Table 1. Field separation and motion compensation angles
Field separation angle Motion compensation angle θn SWIR1 SWIR2 n=2 n=4 n=6 0.5° 1.0° 2.426° 7.238° 11.932° 表 2 特征地物几何位置误差
Table 2. Geometric errors of feature ground objects
SIFT Proposed method Object coordinate VNIR Object coordinate SWIR Error/pixel Object coordinate VNIR Object coordinate SWIR Error/pixel |x| |y| Distance |x| |y| Distance 1 (896.49,1214.57) (891.30,1213.94) 5.2 0.6 5.2 (764.34,892.07) (764.52,892.44) 0.2 0.4 0.4 2 (807.72,801.63) (803.68,802.30) 4.0 0.8 4.1 (721.48,1308.15) (721.63,1308.43) 0.1 0.3 0.3 3 (683.18,1581.33) (678.90,1582.21) 4.1 0.9 4.2 (1080.13,591.93) (1080.47,592.32) 0.3 0.4 0.5 4 (415.23,98.09) (417.95,96.88) 2.7 1.2 2.9 (964.61,2115.28) (964.95,2114.86) 0.3 0.4 0.5 5 (385.33,592.64) (389.34,592) 4.0 0.6 4.0 (1125.58,1641.73) (1125.42,1641.49) 0.3 0.2 0.4 6 (413.81,1045.86) (416.40,1046.47) 2.6 0.6 2.7 (1217.77,1190.68) (1217.63,1190.51) 0.1 0.2 0.2 7 (230.92,252.77) (234.06.,251.92) 3.4 0.8 3.5 (1194.58,2019.39) (1194.84,2019.07) 0.3 0.3 0.4 8 (309.25,868.61) (312.89,869.13 3.6 0.5 3.7 (1275.79,1390.07) (1275.88,1390.28) 0.1 0.1 0.1 9 (168.55,1439.82) (173.13,1440.84) 4.6 1.0 4.7 (1569.27,872.26) (1569.46,872.17) 0.2 0.1 0.2 10 (93.60,1874.59) (97.04,1875.85) 3.4 1.3 3.6 (1747.45,470.69) (1747.35,470.54) 0.1 0.1 0.1 Average 3.8 0.8 3.9 Average 0.2 0.3 0.3 -
[1] Marshall M, Thenkabail P. Advantage of hyperspectral EO-1 Hyperion over multispectral IKONOS, GeoEye-1, WorldView-2, Landsat ETM+, and MODIS vegetation indices in crop biomass estimation [J]. ISPRS Journal of Photogrammetry & Remote Sensing, 2015, 108: 205-218. [2] Guanter L, Kaufmann H, Segl K, et al. The EnMAP spaceborne imaging spectroscopy mission for earth observation [J]. Remote Sensing, 2015, 7(7): 8830-8857. doi: 10.3390/rs70708830 [3] Zhang L. Perspectives on Chinese developments in spaceborne imaging spectroscopy: What to expect in the next 5–10 years[C]//Geoscience and Remote Sensing Symposium (IGARSS), 2016. [4] Gao Ya, Zhou Jialin, Hou Xue, et al. Registration algorithm for hyperspectral image based on Gaussian fitting [J]. Infrared and Laser Engineering, 2016, 45(S2): S223002. (in Chinese) [5] Wang Yun. Study on systematic geometric correction and image registration of hyperspectral Images[D]. Beijing: University of Chinese Academy of Sciences, 2011. (in Chinese) [6] Liu Zhiwen, Liu Dingsheng, Wei Jingbo. Rigorous geometric model of Tiangong-1 hyperspectral spectrometer [J]. Journal of Remote Sensing, 2014, 18(S1): 62-67. (in Chinese) [7] Wu Yinan, Li Guoning, Zhang Ke, et al. Registration model based on homologous points tracking of space camera assembly imaging [J]. Infrared and Laser Engineering, 2016, 45(3): 0326002. (in Chinese) [8] Wang Yueming, Jia Jianxin, He Zhiping, et al. Key technologies of advanced hyperspectral imaging system [J]. Journal of Remote Sensing, 2016, 20(5): 850-857. (in Chinese) [9] Jia G, Zhao H, Shang H. Pixel-size-varying method for simulation of remote sensing images [J]. Journal of Applied Remote Sensing, 2014, 8(1): 083551. (in Chinese) doi: 10.1117/1.JRS.8.083551 [10] Zhang Hongwei, Fan Xiang, Zhu Bin, et al. Dual-band infrared image registration with the introduction of outliers rejection mechanism [J]. Infrared and Laser Engineering, 2015, 44(S1): 23-28. (in Chinese)