椭球方程的矩阵形式

为了简便起见,这里的椭球指的是椭球面,如果还要包含椭球内部,则将等于号改为小于等于即可。

标准椭球方程的几何意义

考虑一个二维平面\mathbb{R}^2上的椭球方程(此时就是一个椭圆方程)是\frac{x_1^2}{r_1^2}+\frac{x_2^2}{r_2^2}=1其在坐标轴x_1,x_2上的半轴分别是r_1,r_2(r_1,r_2>0)。事实上,椭圆(椭球)可以被看作是一个单位圆(单位球)在坐标轴上进行伸缩后得到的几何图形:假设我们对坐标单位进行拉伸,x_1,x_2坐标轴上的单位长度从1分别改成r_1,r_2,由此形成一个新的矩形坐标系,新坐标系的坐标(m,n)在原坐标系的坐标是(r_1m,r_2n),于是,上面的椭圆方程实际上就是新坐标系中的一个单位圆
对于\mathbb{R}^3空间中的一个标准直角坐标系的椭球\sum_{i=1}^3\frac{x_i^2}{r_i^2}=1同样可以视作是一个新坐标系下的单位球,这个坐标系是将原标准直角坐标系的三个坐标轴的单位长度从1分别改成r_1,r_2,r_3而得到的。

现在我们给出一个\mathbb{R}^n空间中在标准直角坐标系的一个标准椭球的方程,其在x_i坐标轴方向的半轴长是r_i(r_i>0)\sum_{i=1}^n\frac{x_i^2}{r_i^2}=1若记对角矩阵\Lambda和坐标向量x分别是\Lambda=\begin{bmatrix} r_1^2\\ &r_2^2\\ &&\ddots\\ &&&r_n^2 \end{bmatrix},\quad x=\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix}则上面的标准椭球方程是x^T\Lambda^{-1} x=1但是该椭球方程仅仅表示中心点在原点、伸缩方向沿标准直角坐标系的坐标轴的椭球,要表示一个任意的椭球,需要使用坐标系的变换 。

正交矩阵和坐标系旋转

考虑\mathbb{R}^n的标准单位正交基E=[e_1,e_2,\cdots,e_n],现在将这个坐标系按原点进行任意旋转,得到一个新的坐标系V。设原来的单位正交向量e_i经旋转变为v_i,则可以得到新坐标系的单位标准正交基V=[v_1,v_2,\cdots,v_n]。我们可以用e_iv_i表示出来\begin{cases} v_1=\sum_k a_{k1}e_k=a_{11}e_1+a_{21}e_2+\cdots+a_{n1}e_n\\ v_2=\sum_k a_{k2}e_k=a_{12}e_1+a_{22}e_2+\cdots+a_{n2}e_n\\ \quad\vdots\\ v_n=\sum_k a_{kn}e_k=a_{1n}e_1+a_{2n}e_2+\cdots+a_{nn}e_n\\ \end{cases}或者使用矩阵形式写出[v_1,\dotsc,v_n]=[e_1,\dotsc,e_n]\begin{bmatrix} a_{11}&a_{12}&\cdots&a_{1n}\\ a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n1}&a_{n2}&\cdots&a_{nn} \end{bmatrix}也即V=EA设空间中任意一个点P,它在EV这两个坐标系中的坐标分别是x=[x_1,x_2,\dotsc,x_n]^Th=[h_1,h_2,\dotsc,h_n]^T,因为无论如何选取坐标系,点在空间的位置是不会发生变化的,于是就应该有Vh=Ex因为Vh=(EA)h=E(Ah)=ExE是可逆矩阵,因此有x=Ah如果将矩阵A写成列向量形式A=[a_1,a_2,\dotsc,a_n],则有v_i=Ea_i。因为v_i彼此单位正交,这就意味着v_i^Tv_j=(Ea_i)^T(Ea_j)=a_i^T(E^TE)a_j=a^T_ia_j=\begin{cases}1&i=j\\0&i\neq j\end{cases}\begin{split} A^TA&=\begin{bmatrix}a^T_1\\a^T_2\\\vdots\\a^T_n\end{bmatrix}[a_1,a_2,\dotsc,a_n]\\ &=\begin{bmatrix}a_1^Ta_1&a_1^Ta_2&\cdots&a_1^Ta_n\\ a_2^Ta_1&a_2^Ta_2&\cdots&a_2^Ta_n\\ \vdots&\vdots&\ddots&\vdots\\ a_n^Ta_1&a_n^Ta_2&\cdots&a_n^Ta_n \end{bmatrix}\\ &=\begin{bmatrix}1&0&\cdots&0\\0&1&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&1\end{bmatrix}=E_n \end{split}A是可逆矩阵,且A^{-1}=A^T,可见A是正交矩阵(同样A^T也是正交矩阵)。于是我们有h=A^{-1}x=A^Tx这个等式的意义是,将一个标准直角坐标系绕原点任意旋转后,坐标从x变为h,两个坐标之间的关系。如果将标准直角坐标系任意旋转,得到两个坐标系,对应的正交矩阵是A_1,A_2,则在标准直角坐标系下坐标为x的点在两个坐标系的下的坐标分别是h_1,h_2,应当满足关系h_1=A_1^Tx,\;h_2=A_2^Tx消去x,就有h_2=A_2^T(A_1h_1)=(A_2^TA_1)h_1易证明,两个正交矩阵的乘积仍然是正交矩阵,于是我们得出结论:空间中任意单位正交坐标系绕原点旋转后形成新的坐标系,同一个点在两个坐标系下的坐标可以使用一个正交矩阵来联系

任意椭球的矩阵方程

我们讨论的椭球方程是以标准直角坐标系为参考系的。n维空间中的任意椭球可以通过如下步骤得到:

  1. 在标准直角坐标系中构造一个中心点在原点的n维单位超球
  2. 将该超球在各个坐标轴方向进行伸缩,得到一个正规的椭球
  3. 将该椭球绕原点进行旋转,使得它与给定的椭球方向一致
  4. 将该椭球平移至给定椭球的位置

正规椭球的矩阵方程已知是x^T\Lambda^{-1}x=1现在将该椭球绕原点进行旋转,但是这等价于将坐标系绕原点旋转,然后在新坐标系中构造出一个正规椭球,设新坐标系中椭球的坐标为h,于是该新坐标系中的正规椭球方程是h^T\Lambda^{-1}h=1现在,根据上一节的结论,我们知道如果该旋转后的椭球在原坐标系下的坐标为x,那么存在一个正交矩阵A使得h=A^Tx,于是旋转后的椭球在原坐标系下的方程是(A^Tx)^T\Lambda^{-1}(A^Tx)=x^T(A\Lambda^{-1}A^T)x=x^TPx=1现在假设所求椭球的中心坐标是x_c,我们只需要将椭球平移至该中心点即可,于是我们得到了任意椭球的方程f(x)=(x-x_c)^TP(x-x_c)=1我们现在来关注核心的矩阵P=A\Lambda^{-1}A^T其中A是正交矩阵,根据上一节和本节的推导,我们知道A的列向量代表了旋转后椭圆的各个伸缩方向的单位矢量,而\Lambda^{-1}=\mathrm{diag}(1/r_1^2,1/r_2^2,\cdots,1/r_n^2),显然\lambda_i=r_i^2\Lambda的特征值,从而\sqrt{\lambda_i}就是对应方向的伸缩系数。

首先,易得P是一个对称矩阵P^T=(A\Lambda^{-1}A^T)^T=A\Lambda^{-1}A^T=P
此外,若a_i是正交矩阵A的第i个列向量,\xi_i是对角阵\Lambda^{-1}的第ii列元素(因此也是第i个特征值),我们断言:\xi_iP关于特征向量a_i的特征值。为了证明这一点,我们有\begin{split} &(P-\xi_iI)a_i\\ =&(A\Lambda^{-1}A^T-\xi_iI)a_i\\ =&(A\Lambda^{-1}A^T-\xi_iAA^T)a_i\\ =&A(\Lambda^{-1}-\xi_iI)A^Ta_i\\ \end{split}
因为A^Ta_i=\begin{bmatrix}a_1^T\\a_2^T\\\vdots\\a_n^T\end{bmatrix}a_i=\begin{bmatrix}a_1^Ta_i\\a_2^Ta_i\\\vdots\\a_n^Ta_i\end{bmatrix}=\begin{bmatrix}\delta_{1i}\\\delta_{2i}\\\vdots\\\delta_{ni}\end{bmatrix}=e_i其中\delta_{ij}=\begin{cases}1&i=j\\0&i\neq j\end{cases}是Kronecker符号。因此(\Lambda^{-1}-\xi_iI)A^Ta_i=(\Lambda^{-1}-\xi_iI)e_i=\eta_i此处\eta_i表示矩阵(\Lambda^{-1}-\xi_iI的第i个列向量,但显然根据定义(\Lambda^{-1}-\xi_iI的第i个列向量为0,从而(\Lambda^{-1}-\xi_iI)A^Ta_i=0,由此可知(P-\xi_iI)a_i=0从而命题可证。

现在,我们已经知道对角矩阵\Lambda^{-1}的对角线元素就是P的特征值,但是显然特征值都是正值,因此P是一个正定矩阵,又因为P还是对称是,因此P是一个对称正定矩阵。如果我们规定P是一个实对称正定矩阵,那么在线性代数中有如下定理:

实对称正定矩阵一定可以相似对角化。

换句话说,给定一个n阶实对称正定矩阵P,方程x^TPx=1就对应n维欧几里得空间中的一个椭球。我们对P一定可以相似对角化,于是得到n个单位正交特征向量,向量代表了椭球伸缩的方向,而对应的特征值的平方根倒数就是该方向的伸缩系数。而给定一个椭球,我们也可以按照上面的构造方法构造出一个实对称正定矩阵P。因此我们得到了椭球方程的矩阵形式。

P=I时,此时椭球退化为一个球。

考虑一个例子,我们希望求得一个椭圆,其中心位于(x_c,y_c)点,半长轴为a,半短轴为b,长轴沿\theta角的方向,从而短轴沿着\theta+\pi/2的方向,现在想求得该椭圆的方程。我们选取长轴和短轴方向的单位向量u_1=(\cos\theta,\sin\theta)以及u_2=(-\sin\theta,\cos\theta),显然它们彼此正交;在这两个方向上的缩放系数分别是a,b,于是可以构造矩阵A=\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix},\Lambda^{-1}=\begin{bmatrix} \frac{1}{a^2}\\&\frac{1}{b^2}\end{bmatrix}于是就有P=A\Lambda^{-1}A^T=\begin{bmatrix}\frac{\cos^2\theta}{a^2}+\frac{\sin^2\theta}{b^2}&\frac{\sin\theta\cos\theta}{a^2}-\frac{\sin\theta\cos\theta}{b^2}\\\frac{\sin\theta\cos\theta}{a^2}-\frac{\sin\theta\cos\theta}{b^2}&\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}\end{bmatrix}所以该椭圆的方程是\begin{bmatrix}x-x_c\\y-y_c\end{bmatrix}^T\begin{bmatrix}\frac{\cos^2\theta}{a^2}+\frac{\sin^2\theta}{b^2}&\frac{\sin\theta\cos\theta}{a^2}-\frac{\sin\theta\cos\theta}{b^2}\\\frac{\sin\theta\cos\theta}{a^2}-\frac{\sin\theta\cos\theta}{b^2}&\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}\end{bmatrix}\begin{bmatrix}x-x_c\\y-y_c\end{bmatrix}=1拆解开就是\begin{split}&\left(\frac{\cos^2\theta}{a^2}+\frac{\sin^2\theta}{b^2}\right)(x-x_c)^2+\left(\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}\right)(y-y_c)^2+\\\cdots&\;\;2\left(\frac{\sin\theta\cos\theta}{a^2}-\frac{\sin\theta\cos\theta}{b^2}\right)(x-x_c)(y-y_c)=1\end{split}

椭球是一个凸集

如果P是一个实对称正定矩阵,那么可以对角化分解为A\Lambda A^T,其中A的正交矩阵(A^TA=AA^T=I),其列向量是P的特征向量,而\Lambda则是由对应特征值组成的对角矩阵\Lambda=\mathrm{diag}(\lambda_1,\lambda_2,\dotsc,\lambda_n)。因为是正定矩阵,因此对角阵的对角线元素必然是正值。于是我们可以定义P的开方为P^{\frac{1}{2}}=A\Lambda^{\frac{1}{2}}A^T其中\Lambda^{\frac{1}{2}}=\mathrm{diag}(\sqrt{\lambda_1},\sqrt{\lambda_2},\dotsc,\sqrt{\lambda_n})。这样的定义是合理的,因为
P^{\frac{1}{2}}P^{\frac{1}{2}}=(A\Lambda^{\frac{1}{2}}A^T)(A\Lambda^{\frac{1}{2}}A^T)=A\Lambda A^T=P
此外定义P^{-\frac{1}{2}}=(P^{\frac{1}{2}})^{-1}=(A\Lambda^{\frac{1}{2}}A^T)^{-1}=A\Lambda^{-\frac{1}{2}}A^T其中\Lambda^{-\frac{1}{2}}=\mathrm{diag}(\frac{1}{\sqrt{\lambda_1}},\frac{1}{\sqrt{\lambda_2}},\dotsc,\frac{1}{\sqrt{\lambda_n}})容易验证\Lambda^{-\frac{1}{2}}\Lambda^{-\frac{1}{2}}=\Lambda^{-1}我们已经知道\mathbb{R}^n空间中的任意椭球的方程是K(x)=(x-x_c)^TP(x-x_c)现在将其进行变换\begin{split} K(x)&=(x-x_c)^T(P^{\frac{1}{2}}P^{\frac{1}{2}})(x-x_c)\\ &=(x-x_c)^T((P^{\frac{1}{2}})^TP^{\frac{1}{2}})(x-x_c)\\ &=[P^{\frac{1}{2}}(x-x_c)]^T[P^{\frac{1}{2}}(x-x_c)]\\ &=u^Tu=B(u)\leqslant 1 \end{split}其中u=P^{\frac{1}{2}}(x-x_c)进一步得到x=P^{-\frac{1}{2}}u+x_c我们现在来看,集合B=\{u|u^Tu\leqslant1\}表示一个单位球,定义映射x=f(u)=P^{-\frac{1}{2}}u+x_c我们有f(B)=\{f(u)|u\in B\}=\{x|(x-x_c)^TP(x-x_c)\leqslant 1\}=K因此,给定一个这样的映射f,它可以将一个圆心在原点的单位球唯一映射为一个椭球。因为B是一个凸集,如果f是一个仿射映射,那么根据凸优化理论,B的象K也是一个凸集。

一个映射f(x)是一个仿射映射,如果可以写为f(x)=L(x)+b,其中b是一个常量,而L(x)是一个线性函数(即L(ax+by)=aL(x)+bL(y))。显然f是一个仿射函数,因此椭球K就是一个凸集。

椭球矩阵正定性和椭圆包含关系

椭球球心在原点的一个任意椭球的方程已知是x^TPx\leqslant1,其中P\in S^n_{++},表示P是一个n阶对称正定矩阵。设椭球\mathcal{A}:x^TAx\leqslant1和椭球\mathcal{B}:x^TBx\leqslant1A,B\in S^n_{++}),我们断言:
B^{-1}-A^{-1}\in S^n_{+}\iff A-B\in S^n_{+}\iff\mathcal{A}\subseteq\mathcal{B}
其中S^n_{+}表示所有对称半正定矩阵的集合,它是一个半正定锥。

我们的证明通过如下过程
\begin{array}{rlr} \mathcal{A}\subseteq\mathcal{B}&\iff \forall x\in\mathcal{A}\Rightarrow x\in\mathcal{B}&(1)\\ &\iff x^TAx\leqslant1\Rightarrow x^TBx\leqslant1&(2)\\ &\iff x^TBx\leqslant x^TAx&(3)\\ &\iff x^T(A-B)x\geqslant0&(4)\\ &\iff (A-B)\in S^n_{+}&(5)\\ &\iff B^{-1}-A^{-1}\in S^n_+&(6) \end{array}
我们首先来看(2)\iff(3):其中(3)\Rightarrow(2)是显然成立的。为了证明(2)\Rightarrow(3),也就证明:如果x^TAx\leqslant1能推出x^TBx\leqslant1,那么就有x^TBx\leqslant x^TAx,我们使用反证法,假设存在一个x_0,满足x_0^TBx_0>x_0^TAx_0,注意到A是对称正定矩阵,即二次型是一个正实数,即t=x_0^TAx_0>0,于是我们就有\left(\frac{x_0}{\sqrt{t}}\right)^TB\left(\frac{x_0}{\sqrt{t}}\right)>1,但是,因为x^TAx\leqslant1\Rightarrow x^TBx\leqslant1的逆否命题是x^TBx>1\Rightarrow x^TAx>1,这就意味着\left(\frac{x_0}{\sqrt{t}}\right)^TA\left(\frac{x_0}{\sqrt{t}}\right)>1,进而有x_0^TAx_0>t=x_0^TAx_0,这显然是不可能的,所以这样的x_0是不存在的,命题可证。

我们接下来说明(5)\iff(6)。我们依次证明下面的命题:

  • 如果A是对称可逆矩阵,那么A^{-1}也是对称可逆矩阵:E=AA^{-1}=(AA^{-1})^T=(A^{-1})^TA^T=(A^{-1})^TA从而A^{-1}=(A^{-1})^T
  • 如果A是对称(半)正定/(半)负定矩阵,那么A^{-1}也是对称(半)正定/(半)负定矩阵。这是因为\forall x\begin{split}x^TA^{-1}x&=x^TA^{-1}AA^{-1}x\\&=\left((A^{-1})^Tx\right)^TA(A^{-1}x)\\&=(A^{-1}x)^TA(A^{-1}x)\\&=x'^TAx'\\&\lesseqgtr 0\end{split}
  • A,B合同,如果存在可逆矩阵P使得P^TAP=B。我们说,合同变换不改变矩阵的(半)正定/(半)负定性。这是因为\forall x\begin{split}x^TBx&=x^TP^TAPx\\&=(Px)^TA(Px)\\&=x'^TAx'\\&\lesseqgtr0\end{split}
  • 两个(半)正定矩阵(或两个(半)负定矩阵)的和仍然是(半)正定矩阵(或(半)负定矩阵),其证明是显然的。

一般的,我们有\begin{split}B^{-1}-A^{-1}&=A^{-1}(A-B)B^{-1}\\&=A^{-1}(A-B)(A^{-1}+B^{-1}-A^{-1})\\&=A^{-1}(A-B)A^{-1}+A^{-1}(A-B)(B^{-1}-A^{-1})\\&=A^{-1}(A-B)A^{-1}+A^{-1}(A-B)B^{-1}(A-B)A^{-1}\end{split}由于A,B都是对称正定矩阵,从而A^{-1},B^{-1}也是对称正定矩阵,A-B是对称矩阵,于是B^{-1}-A^{-1}=\underbrace{(A^{-1})^T(A-B)(A^{-1})}_{M}+\underbrace{\left((A-B)A^{-1}\right)^TB^{-1}\left((A-B)A^{-1}\right)}_{N}如果A-B是半正定矩阵,那么M,(A-B)两个矩阵合同,故而M是半正定矩阵;此外,(A-B)A^{-1}是可逆矩阵,因此N,B^{-1}合同,从而N正定。于是可得B^{-1}-A^{-1}是半正定的。同理,如果B^{-1}-A^{-1}半正定,那么A-B也半正定。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 202,607评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,047评论 2 379
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,496评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,405评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,400评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,479评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,883评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,535评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,743评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,544评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,612评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,309评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,881评论 3 306
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,891评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,136评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,783评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,316评论 2 342

推荐阅读更多精彩内容

  • 8、行列式 8.1 什么是行列式? 首先方阵才有行列式,我们先来简单回顾一下2*2和3*3的矩阵的行列式: 那行列...
    城市中迷途小书童阅读 19,648评论 2 45
  • 如果不熟悉线性代数的概念,要去学习自然科学,现在看来就和文盲差不多。”,然而“按照现行的国际标准,线性代数是通过公...
    Drafei阅读 1,532评论 0 3
  • 线性代数在科学领域有很多应用的场景,如下: 矩阵,是线性代数中涉及的内容, 线性代数是用来描述状态和变化的,而矩阵...
    zhoulujun阅读 11,952评论 3 44
  • 童话里竟没人告诉我西安的夏天会热到甩开极限几条街,你以为阳伞就能遮出阴凉吗?还不如用它来挡夏天里莫名而来的阵雨! ...
    Lu小雨在穿越阅读 107评论 0 0
  • 我是昨天晚上到的宿舍,都八点多了,弟弟打电话问我走不走,他要开车来厂里找她女朋友,可以捎上我。收拾下东西就坐上他的...
    清香的泥土阅读 334评论 3 2